hdMTD: High-Dimensional MTD Modeling in R

• hdMTD is an R package for fitting and simulating high-dimensional Mixture Transition Distribution models that capture long-range dependencies in categorical time series.
• It implements various lag-selection methods—including BIC-based, CUT, FS, and FSC approaches—to estimate sparse lag sets and reduce the exponential parameter growth of high-order Markov chains.
• The package also provides EM-based parameter estimation and exact stationary simulation, with empirical illustrations demonstrating improved performance over traditional models.

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 (Gripp et al., 1 Sep 2025).

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$0-dimensional transition with

$d$1

where $d$2 satisfy

$d$3

Here $d$4 is a distribution on $d$5 independent of the past, and each $d$6 is a distribution on $d$7 for lag $d$8 and past symbol $d$9. In sampling terms, one first draws a lag index $\{X_t\}_{t\in\mathbb{Z}}$0 with probability $\{X_t\}_{t\in\mathbb{Z}}$1; if $\{X_t\}_{t\in\mathbb{Z}}$2, one samples from $\{X_t\}_{t\in\mathbb{Z}}$3, and otherwise from $\{X_t\}_{t\in\mathbb{Z}}$4.

This representation is parsimonious because it requires $\{X_t\}_{t\in\mathbb{Z}}$5 lag-specific matrices and mixture weights rather than a full $\{X_t\}_{t\in\mathbb{Z}}$6-lag joint transition law. In the multimatrix case, the parameter count is linear in $\{X_t\}_{t\in\mathbb{Z}}$7, up to $\{X_t\}_{t\in\mathbb{Z}}$8. The package is centered on the multimatrix MTD, in which each lag has its own matrix $\{X_t\}_{t\in\mathbb{Z}}$9, 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 $A$0 can be large, possibly of order $A$1, 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 $A$2, because estimating high-dimensional joint distributions becomes infeasible when most length-$A$3 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 $A$4,

$A$5

where $A$6 is total variation distance and $A$7 differ only at lag $A$8. For an MTD,

$A$9

If $d$0 or all rows of $d$1 are identical, then $d$2, so lag $d$3 is irrelevant. The relevant lag set is therefore

$d$4

Once $d$5 is known, the transition law depends only on $d$6, not on the full $d$7-lag past.

2. Mathematical formulation and likelihood structure

hdMTD assumes a stationary Markov chain of order $d$8 on a finite state space $d$9. Each lag kernel $$P(X_t = a \mid X_{t-1} = x_{-1}, \dots, X_{t-d} = x_{-d}) = P(a \mid x_{-1}, \dots, x_{-d}),$$0 can be represented by a $$P(X_t = a \mid X_{t-1} = x_{-1}, \dots, X_{t-d} = x_{-d}) = P(a \mid x_{-1}, \dots, x_{-d}),$$1 stochastic matrix,

$$P(X_t = a \mid X_{t-1} = x_{-1}, \dots, X_{t-d} = x_{-d}) = P(a \mid x_{-1}, \dots, x_{-d}),$$2

Inference is built from empirical joint and conditional probabilities computed from a sample $$P(X_t = a \mid X_{t-1} = x_{-1}, \dots, X_{t-d} = x_{-d}) = P(a \mid x_{-1}, \dots, x_{-d}),$$3.

For $$P(X_t = a \mid X_{t-1} = x_{-1}, \dots, X_{t-d} = x_{-d}) = P(a \mid x_{-1}, \dots, x_{-d}),$$4 and $$P(X_t = a \mid X_{t-1} = x_{-1}, \dots, X_{t-d} = x_{-d}) = P(a \mid x_{-1}, \dots, x_{-d}),$$5, the joint counts are

$$P(X_t = a \mid X_{t-1} = x_{-1}, \dots, X_{t-d} = x_{-d}) = P(a \mid x_{-1}, \dots, x_{-d}),$$6

with empirical joint probabilities

$$P(X_t = a \mid X_{t-1} = x_{-1}, \dots, X_{t-d} = x_{-d}) = P(a \mid x_{-1}, \dots, x_{-d}),$$7

If $$P(X_t = a \mid X_{t-1} = x_{-1}, \dots, X_{t-d} = x_{-d}) = P(a \mid x_{-1}, \dots, x_{-d}),$$8, the empirical conditional probabilities are

$$P(X_t = a \mid X_{t-1} = x_{-1}, \dots, X_{t-d} = x_{-d}) = P(a \mid x_{-1}, \dots, x_{-d}),$$9

For a lag set $d$0, hdMTD also uses

$d$1

and

$d$2

Given a subset $d$3, the nonparametric log-likelihood is

$d$4

and the BIC objective used in the package is

$d$5

where $d$6 is the number of free parameters under the chosen MTD structure and $d$7 is a user-specified constant.

This formulation is significant because it separates two statistical tasks. The first is structural recovery of $d$8, 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 (Gripp et al., 1 Sep 2025).

The BIC-based algorithm, available through hdMTD_BIC() or hdMTD(..., method="BIC"), selects a lag set $d$9 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 $|A|^d(|A|-1)$0 and determines whether each lag $|A|^d(|A|-1)$1 is relevant by comparing empirical conditional distributions over pairs of $|A|^d(|A|-1)$2-compatible pasts. The rule is based on

$|A|^d(|A|-1)$3

with threshold

$|A|^d(|A|-1)$4

where

$|A|^d(|A|-1)$5

for tunings $|A|^d(|A|-1)$6, $|A|^d(|A|-1)$7, and $|A|^d(|A|-1)$8 with $|A|^d(|A|-1)$9. If any compatible pair exceeds the threshold, lag $d$0 is kept; otherwise it is cut. Its complexity is $d$1, and consistency is proved when $d$2 and $d$3, even if $d$4. 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 $d$5. At step $d$6, it adds the lag maximizing

$d$7

This quantity measures how strongly the distribution of $d$8 changes with $d$9, conditional on the currently selected lag set. The complexity is

$d$0

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 $d$1, the probability that $d$2 is high for large $d$3, 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 $d$4; 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 $d$5 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 $d$6, 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 $d$7. The E-step computes expected counts for component usage and lag-specific transitions, and the M-step updates $d$8, $d$9, and $d$00 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 $d$01. 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 $d$02, 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() (Gripp et al., 1 Sep 2025). Given an MTD object with $d$03, it uses a backward randomization algorithm inspired by Comets et al. (2002). For each time $d$04, one samples $d$05 with $d$06. To simulate $d$07, one samples a lag index $d$08. If $d$09, $d$10 is sampled from $d$11 and the recursion stops. Otherwise, if $d$12 has not yet been generated, the procedure is called recursively at time $d$13, after which $d$14 is sampled from $d$15. Because $d$16, 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 $d$17, relevant lags $d$18, and weights $d$19, $d$20, $d$21, and $d$22. 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 $d$23. As the sample size $d$24 increases, FS’s mean total variation error $d$25 approaches Oracle’s $d$26, FS recovers $d$27 exactly in almost all replicates for $d$28, and the Naive estimator using lags $d$29 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, $d$30 °C, about $d$31; and Category 2, $d$32 °C, about $d$33. After reversing time order, hdMTD_FS(Temp12, d = 400, l = 3) returns 1 364 6. The interpretation given is lag $d$34 for yesterday’s regime, lag $d$35 for a weekly pattern, and lag $d$36 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 $d$37 for the independent model, $d$38 for the second-order Markov chain, and $d$39 for the FS-based MTD; sensitivity approximately $d$40, $d$41, and $d$42; precision approximately $d$43, $d$44, and $d$45; and F1 score approximately $d$46, $d$47, and $d$48. 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 $d$49 small and, in the cited framework, growing at most like $d$50. This is the regime in which lag selection remains statistically and computationally manageable even when the maximum order $d$51 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$52. 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 $d$53 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 (Assländer et al., 2022). 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.