---
title: 'hdMTD: High-Dimensional MTD Modeling in R'
url: https://www.emergentmind.com/topics/hdmtd
type: topic
---

# hdMTD: High-Dimensional MTD Modeling in R

hdMTD is an R package for fitting and simulating high-dimensional Mixture Transition Distribution (MTD) models, a parsimonious class of high-order Markov chains designed to capture long-range dependencies in categorical time series. In an MTD, the transition probability is written as a convex mixture of lower-order conditional distributions that depend on single lags, which replaces the exponential parameter growth of a general order-\(d\) Markov chain with a representation whose parameter count is linear in \(d\). The package implements non-parametric lag selection, maximum likelihood estimation of transition probabilities, expectation-maximization estimation of MTD parameters, and exact simulation from the stationary distribution, with illustrations on simulated data and daily temperature records from Brazil [2509.01808].

## 1. Model class and conceptual foundations

Let \( \{X_t\}_{t\in\mathbb{Z}} \) be a categorical time series taking values in a finite state space \(A\). A general order-\(d\) Markov chain specifies
\[
P(X_t = a \mid X_{t-1} = x_{-1}, \dots, X_{t-d} = x_{-d}) = P(a \mid x_{-1}, \dots, x_{-d}),
\]
and requires a transition probability for every combination of \(d\) past states. The parameter count grows like \(|A|^d(|A|-1)\), which is exponential in \(d\).

An MTD model of order \(d\) replaces the full \(d\)-dimensional transition with
\[
P(a \mid x_{-1},\dots,x_{-d})
=
\lambda_0 p_0(a) + \sum_{j=-d}^{-1} \lambda_j\, p_j(a \mid x_j), \quad a\in A,
\]
where \(\lambda_j \in [0,1]\) satisfy
\[
\lambda_0 + \sum_{j=-d}^{-1} \lambda_j = 1.
\]
Here \(p_0(\cdot)\) is a distribution on \(A\) independent of the past, and each \(p_j(\cdot \mid b)\) is a distribution on \(A\) for lag \(j\) and past symbol \(b\in A\). In sampling terms, one first draws a lag index \(j\in\{0,-1,\dots,-d\}\) with probability \(\lambda_j\); if \(j=0\), one samples from \(p_0(\cdot)\), and otherwise from \(p_j(\cdot \mid X_{t+j})\).

This representation is parsimonious because it requires \(d\) lag-specific matrices and mixture weights rather than a full \(d\)-lag joint transition law. In the multimatrix case, the parameter count is linear in \(d\), up to \(d+(|A|-1)(1+d|A|)\). The package is centered on the multimatrix MTD, in which each lag has its own matrix \(p_j\), but it also supports the single-matrix case, where all lags share the same matrix.

The high-dimensional regime considered by hdMTD is the regime in which the maximum lag \(d\) can be large, possibly of order \(n\), the sample size. This is motivated by settings in which dependence may occur simultaneously at short and long scales, such as daily, weekly, and yearly lags. Traditional nonparametric Markov estimation typically forces \(d \approx O(\log n)\), because estimating high-dimensional joint distributions becomes infeasible when most length-\(d\) pasts are unobserved or rarely observed. hdMTD adopts the sparse-lag perspective: the relevant lag set is assumed to be small even when the maximal order is large.

A central quantity is the oscillation at lag \(j\),
\[
\delta_j = \max \big\{ d_{TV}\big(P(\cdot \mid x), P(\cdot \mid y)\big) :\ x,y\in A^d,\ x_k = y_k\ \forall k\neq j \big\},
\]
where \(d_{TV}\) is total variation distance and \(x,y\) differ only at lag \(j\). For an MTD,
\[
\delta_j = \lambda_j \cdot \max_{b,c\in A} d_{TV}\big(p_j(\cdot \mid b), p_j(\cdot \mid c)\big).
\]
If \(\lambda_j=0\) or all rows of \(p_j\) are identical, then \(\delta_j=0\), so lag \(j\) is irrelevant. The relevant lag set is therefore
\[
\Lambda = \{j\in\{-d,\dots,-1\} : \delta_j > 0\}.
\]
Once \(\Lambda\) is known, the transition law depends only on \(x_\Lambda\), not on the full \(d\)-lag past.

## 2. Mathematical formulation and likelihood structure

hdMTD assumes a stationary Markov chain of order \(d\) on a finite state space \(A\). Each lag kernel \(p_j\) can be represented by a \(|A|\times |A|\) stochastic matrix,
\[
p_j = \{p_j(a\mid b): a,b\in A\}.
\]
Inference is built from empirical joint and conditional probabilities computed from a sample \(S_n=(X_1,\dots,X_n)\).

For \(S\subseteq\{-d,\dots,-1\}\) and \(j\in\{-d,\dots,0\}\setminus S\), the joint counts are
\[
N_{n,j}(x_S,b,a)
=
\sum_{t=d+1}^{n}
\mathbb{I}\{X_{t+k} = x_k,\, k\in S;\ X_{t+j}=b;\ X_t = a\},
\]
with empirical joint probabilities
\[
\hat{\pi}_{n,j}(x_S,b,a) = \frac{N_{n,j}(x_S,b,a)}{n-d}.
\]
If \(\hat{\pi}_{n,j}(x_S,b)>0\), the empirical conditional probabilities are
\[
\hat{P}_{n,j}(a \mid x_S,b)
=
\frac{\hat{\pi}_{n,j}(x_S,b,a)}{\hat{\pi}_{n,j}(x_S,b)},
\qquad
\hat{\pi}_{n,j}(x_S,b)=\sum_{a\in A}\hat{\pi}_{n,j}(x_S,b,a).
\]

For a lag set \(S\), hdMTD also uses
\[
\bar{N}_n(x_S) = \sum_{a} N_n(x_S,a), \qquad \hat{\pi}_n(x_S) = \sum_{a}\hat{\pi}_n(x_S,a),
\]
and
\[
\hat{P}_n(a \mid x_S) =
\begin{cases}
\frac{\hat{\pi}_n(x_S,a)}{\hat{\pi}_n(x_S)}, & \hat{\pi}_n(x_S)>0 \\
|A|^{-1}, & \text{otherwise}.
\end{cases}
\]

Given a subset \(S\), the nonparametric log-likelihood is
\[
\ell(S) = \sum_{x_S\in A^S} \sum_{a\in A} N_n(x_S,a)\,\log \hat{P}_n(a \mid x_S),
\]
and the BIC objective used in the package is
\[
\text{BIC}(S;n) = -\,\ell(S) + \text{Pen}(S,n), \quad \text{Pen}(S,n)=\theta(S)\,\log n\,\xi,
\]
where \(\theta(S)\) is the number of free parameters under the chosen MTD structure and \(\xi\) is a user-specified constant.

This formulation is significant because it separates two statistical tasks. The first is structural recovery of \(\Lambda\), the set of relevant lags. The second is estimation of the MTD parameters once a lag set has been chosen. hdMTD provides dedicated procedures for both tasks, rather than treating the problem as a single global search over all high-order Markov transition tables.

## 3. Lag selection, transition estimation, and EM inference

The package implements four lag-selection procedures: a BIC-based algorithm, the CUT estimator, the FS estimator, and the FSC estimator [2509.01808].

The BIC-based algorithm, available through `hdMTD_BIC()` or `hdMTD(..., method="BIC")`, selects a lag set \(S \subseteq S^*\) minimizing a penalized log-likelihood over subsets whose sizes are restricted by `minl` and `maxl`. It allows multimatrix versus single-matrix structure and models with or without an independent part. Its main advantage is classical likelihood-based model selection; its main limitation is the exponential growth of the search space, so tractability depends on restricting the candidate set `S`, the size range, or both.

The CUT estimator, implemented as `hdMTD_CUT()`, starts from a candidate set \(S\) and determines whether each lag \(j\in S\) is relevant by comparing empirical conditional distributions over pairs of \(S\setminus\{j\}\)-compatible pasts. The rule is based on
\[
d_{TV}\big(\hat{P}_n(\cdot\mid x_S),\hat{P}_n(\cdot\mid y_S)\big) > t_n(x_S,y_S),
\]
with threshold
\[
t_n(x_S,y_S)=s_n(x_S)+s_n(y_S),
\]
where
\[
s_n(x_S) = \sqrt{
\frac{\alpha(1+\xi)}{2\bar{N}_n(x_S)}
\sum_{a\in A}
\sqrt{\frac{\mu}{\mu - \psi(\mu)}\Big(\hat{P}_n(a\mid x_S) + \frac{\alpha}{\bar{N}_n(x_S)}\Big)} } + \frac{\alpha |A|}{6 \bar{N}_n(x_S)},
\]
for tunings \(\alpha>0\), \(\xi>0\), and \(\mu\in(0,3)\) with \(\mu > \psi(\mu)=e^\mu-\mu-1\). If any compatible pair exceeds the threshold, lag \(j\) is kept; otherwise it is cut. Its complexity is \(O(|A|^2\,|S|\,(n-d))\), and consistency is proved when \(|S| = O(\log n)\) and \(\Lambda\subseteq S\), even if \(d=\beta n\). In practice the defaults are `alpha=0.05`, `mu=1`, and `xi=0.5`.

The forward stepwise estimator `hdMTD_FS()` greedily builds a lag set of size \(\ell\). At step \(k\), it adds the lag maximizing
\[
\hat{\nu}_{n,j,S} =
\sum_{x_S\in A^S} \sum_{b\in A} \sum_{c\in A}
\frac{\hat{\pi}_{n,j}(x_S,b)\,\hat{\pi}_{n,j}(x_S,c)}{\hat{\pi}_n(x_S)}
\, d_{TV}\big( \hat{P}_{n,j}(\cdot\mid x_S,b), \hat{P}_{n,j}(\cdot\mid x_S,c) \big).
\]
This quantity measures how strongly the distribution of \(X_t\) changes with \(X_{t+j}\), conditional on the currently selected lag set. The complexity is
\[
O(|A|^3 \,\ell\, (n-d)(d-(\ell-1)/2)).
\]
The order of the output encodes importance: the first selected lag is the most predictive among the candidates considered at the initial step. The stated theory implies that if \(\ell \ge |\Lambda|\), the probability that \(\Lambda \subseteq \hat{\Lambda}^{\text{FS}}_n\) is high for large \(n\), so FS is naturally used to build a superset of the relevant lag set.

The FSC estimator `hdMTD_FSC()` combines FS and CUT with sample splitting. One half of the sample is used to run FS and produce a candidate set \(S\); the other half is used to run CUT on that set. This is intended to avoid double-using the same data and yields a fully data-driven estimator of \(\Lambda\) in sparse high-order settings.

Once a lag set has been selected, hdMTD provides two forms of estimation. The function `probs()` returns empirical conditional probabilities \(\hat{P}_n(a\mid x_S)\), which are the maximum likelihood estimates under the nonparametric Markov assumption. The function `MTDest()` fits an MTD model by EM. In the multimatrix case, the latent variable is the component index \(L_t\in\{0\}\cup S\). The E-step computes expected counts for component usage and lag-specific transitions, and the M-step updates \(\lambda_j\), \(p_0(a)\), and \(p_j(a\mid b)\) by normalized expected counts. Iteration stops when the increase in log-likelihood is less than `M` (default `M=0.01`) or the number of iterations reaches `nIter` (default 100). The output may include updated mixture weights, lag matrices, the independent distribution, the number of iterations, log-likelihood increments, and estimated oscillations when `oscillations=TRUE`.

## 4. Software structure and workflow in R

hdMTD exposes both model-construction functions and inference functions. The package adopts a specific lag-indexing convention: outputs are given as positive integers, but these represent negative time lags. Thus `c(1,15,30)` represents lags \(-1,-15,-30\). The input sample `X` is ordered from latest to oldest so that lag indexing matches the internal convention.

| Function | Role |
|---|---|
| `MTDmodel()` | Create an MTD object from `Lambda`, `A`, `lam0`, `lamj`, `p0`, and `pj` |
| `hdMTD()` | Generic interface with `method = "FS"`, `"CUT"`, `"BIC"`, or `"FSC"` |
| `hdMTD_FS()`, `hdMTD_CUT()`, `hdMTD_BIC()`, `hdMTD_FSC()` | Method-specific lag-selection procedures |
| `probs()` | Empirical conditional probabilities for a chosen lag set |
| `oscillation()` | True or empirical oscillations for lags |
| `MTDest()` | EM estimation of MTD parameters |
| `perfectSample()` | Exact simulation from the stationary distribution |

`MTDmodel(Lambda, A, lam0, lamj, p0, pj, single_matrix, indep_part)` creates an MTD object and returns the transition matrix over all combinations of states at the relevant lags, the component matrices `pj`, and the weights `lambdas`. If `pj` is not provided, hdMTD samples the matrices uniformly and normalizes them. The options `single_matrix=TRUE` and `indep_part=FALSE` correspond to the single-matrix variant and the model without independent part.

A typical workflow begins with data preparation, including categorical encoding and reversal of time order. Lag selection then depends on prior information. If there is no prior candidate set, the package documentation recommends `hdMTD_FS(X, d, l)` with a reasonably large `d` and modest `l`, optionally followed by `hdMTD_FSC(X, d, l)` to prune irrelevant lags. If there is a candidate set `S`, one may use `hdMTD_BIC(X, d, S, minl, maxl, ...)` or `hdMTD_CUT(X, d, S, alpha, mu, xi)`. After selecting \(\hat{\Lambda}\), empirical transition probabilities can be obtained with `probs(X, S, matrixform=TRUE/FALSE)`. Full MTD parameters are then estimated through `MTDest()`, and a complete model object can be reconstructed with `MTDmodel()`.

This organization makes a clear distinction between nonparametric transition estimation and structured MTD estimation. The former uses empirical conditional probabilities directly; the latter estimates the convex-mixture representation itself. That distinction is important in practice because a chosen lag set can support either a nonparametric Markov model on the reduced lag set or an MTD model on the same set.

## 5. Exact simulation and empirical illustrations

hdMTD includes an exact simulation procedure from the stationary distribution through `perfectSample()` [2509.01808]. Given an MTD object with \(\lambda_0>0\), it uses a backward randomization algorithm inspired by Comets et al. (2002). For each time \(s\le t\), one samples \(L_s\in\{0\}\cup\Lambda\) with \(P(L_s=j)=\lambda_j\). To simulate \(X_t\), one samples a lag index \(j\). If \(j=0\), \(X_t\) is sampled from \(p_0(\cdot)\) and the recursion stops. Otherwise, if \(X_{t+j}\) has not yet been generated, the procedure is called recursively at time \(t+j\), after which \(X_t\) is sampled from \(p_j(\cdot\mid X_{t+j})\). Because \(\lambda_0>0\), the recursion terminates in finite time almost surely. The result is a sample exactly distributed according to the stationary distribution, without burn-in or convergence diagnostics.

A simulated binary example uses \(A=\{0,1\}\), relevant lags \(\Lambda=\{-30,-15,-1\}\), and weights \(\lambda_0 = 0.01\), \(\lambda_{-1} = 0.39\), \(\lambda_{-15} = 0.30\), and \(\lambda_{-30} = 0.30\). With `d=40` and `l=4`, `hdMTD_FS(X, d = 40, l = 4)` returns `30 15 1 27`, recovering the three relevant lags and adding one extra lag because four lags were requested. In the same example, `hdMTD_BIC(X, d = 40, minl = 4, maxl = 4)` returns `1 15 17 30`. The paper also reports empirical comparisons among FS, a “Naive” estimator, and an “Oracle” estimator for a model with true \(\Lambda = \{-5,-1\}\). As the sample size \(m\) increases, FS’s mean total variation error \(\bar{\Delta}_{FS}(m)\) approaches Oracle’s \(\bar{\Delta}_{Oracle}(m)\), FS recovers \(\{-5,-1\}\) exactly in almost all replicates for \(m=10\,000\), and the Naive estimator using lags \(\{-5,-4,-3,-2,-1\}\) has much higher error.

The main real-world illustration is based on hourly maximum temperatures from Brasília (INMET), 2010–2024, aggregated into daily mean of hourly maximum temperatures and discretized into two categories: Category 1, \([12.046, 20.923)\) °C, about \(26.1\%\); and Category 2, \([20.923, 29.8]\) °C, about \(73.9\%\). After reversing time order, `hdMTD_FS(Temp12, d = 400, l = 3)` returns `1 364 6`. The interpretation given is lag \(-1\) for yesterday’s regime, lag \(-6\) for a weekly pattern, and lag \(-364\) for an annual pattern. CUT applied to the FS output keeps all three lags, and BIC on that candidate set chooses the full set `c(1, 6, 364)` as best.

The resulting FS-based MTD is compared with an independent model and a classical second-order Markov chain selected by BIC among orders 1–6. In a prediction experiment on the last 366 days, averaged over 1000 simulations, the reported metrics are: accuracy approximately \(65.1\%\) for the independent model, \(82.4\%\) for the second-order Markov chain, and \(83.5\%\) for the FS-based MTD; sensitivity approximately \(26.7\%\), \(58.7\%\), and \(63.0\%\); precision approximately \(17.5\%\), \(49.8\%\), and \(52.5\%\); and F1 score approximately \(21.1\%\), \(53.9\%\), and \(57.2\%\). The stated interpretation is that the sparse high-order model captures long-range seasonality and weekly structure that are not represented by the selected second-order chain.

## 6. Assumptions, limitations, and related usages of the acronym

The package is built for stationary MTD chains on a finite state space. Continuous data must therefore be discretized. Its theoretical guarantees are tied to sparsity of the relevant lag set, with \(|\Lambda|\) small and, in the cited framework, growing at most like \(\log n\). This is the regime in which lag selection remains statistically and computationally manageable even when the maximum order \(d\) is large.

Several limitations follow directly from that setup. If the dependence structure cannot be well represented as a convex mixture of single-lag conditional distributions, hdMTD may misrepresent the dynamics; the package documentation explicitly notes interactions among multiple lags simultaneously as an example. Even though the MTD parametrization is parsimonious, empirical probabilities can still become sparse for very large state spaces or extremely high \(d\). BIC search can be computationally expensive when the candidate set is large and the admissible subset sizes are broad, which is why FS and FSC are recommended in high-dimensional settings.

Practical interpretation also requires care. Positive lag indices in outputs represent negative time lags; oscillation values \(\delta_j\) quantify the strength of each lag’s effect; and the order of lags in FS output reflects selection order and incremental importance, whereas CUT and FSC outputs do not encode such a ranking. This suggests that hdMTD is best viewed not simply as a high-order Markov package, but as a workflow for sparse lag discovery, reduced-set transition estimation, structured MTD fitting, and stationary simulation.

A potential source of confusion is nomenclature. In a distinct context, the term “hdMTD” has been used informally to denote “hybrid-state driven MT dynamics” or “hybrid-state driven MT disposition” in quantitative magnetization transfer imaging, although the acronym is not used explicitly in that paper [2207.08259]. That usage refers to hybrid-state spin physics and generalized Bloch modeling rather than to mixture transition distribution models. In the statistical context of the R package, hdMTD denotes high-dimensional Mixture Transition Distribution models and the software environment built for their inference and simulation.

Source: https://www.emergentmind.com/topics/hdmtd