MSTest: An R Package for Testing Markov Switching Models

We present the R package MSTest, which implements hypothesis testing procedures to determine the number of regimes in Markov switching models. These models have wide-ranging applications in economics, finance, and many other fields. MSTest provides several testing frameworks, including Monte Carlo likelihood ratio tests, moment-based tests, parameter stability tests, and classical likelihood ratio procedures. In addition, the package offers tools for simulating and estimating univariate and multivariate Markov switching and hidden Markov models using either the expectation–maximization algorithm or maximum likelihood estimation. The functionality of the package is demonstrated through simulation-based examples.

1 Introduction

Markov switching models were introduced by Goldfeld and Quandt (1973) but became an active area of research in economics after Hamilton (1989) proposed modeling the first difference of U.S. GNP as a nonlinear stationary process, where the nonlinearity arises from discrete shifts in the process.

These models have since been applied in a wide range of macroeconomic and financial settings. For example, they have been used to identify business cycles and provide probabilistic statements about the state of the economy (Hamilton 1989; Diebold and Rudebusch 1996; Chauvet 1998; Kim and Nelson 1999; Chauvet et al. 2002; Chauvet and Hamilton 2006; Qin and Qu 2026), to model stock market volatility using Markov switching ARCH, GARCH, and stochastic volatility models (Hamilton 1994; Gray 1996; So et al. 1998; Klaassen 2002; Haas et al. 2004; Pelletier 2006), to model interest rate dynamics (Cai 1994; Garcia and Perron 1996), to study state-dependent impulse response functions (Sims and Zha 2006; Caggiano et al. 2017), to identify structural shocks in SVAR models (Lanne et al. 2010; Herwartz and Lütkepohl 2014; Lütkepohl et al. 2021), and more recently to improve measures of core inflation by allowing for multiple inflation regimes (Le Bihan et al. 2024; Rodriguez-Rondon 2024; Ahn and Luciani 2026). Hamilton (2016) provides a detailed survey of regime switching models in macroeconomics.

Outside of macroeconomic and financial applications, these models have also been applied in climate change research (Golosov et al. 2014; Dietz and Stern 2015), environmental and energy economics (Charfeddine 2017; Cevik et al. 2021), and industrial organization (Aguirregabiria and Mira 2007; Sweeting 2013). Additionally, there is a related class of models—the hidden Markov model—which has been widely used in computational molecular biology (Baldi et al. 1994; Krogh et al. 1994), handwriting and speech recognition (Nag et al. 1986; Rabiner and Juang 1986; Rabiner and Juang 1993; Jelinek 1997), computer vision and pattern recognition (Bunke and Caelli 2001), and other machine learning applications.

Given their empirical relevance, determining the number of regimes needed to capture these nonlinearities is a key step in specifying a Markov switching model, since it is not determined endogenously during estimation. However, the asymptotic results underlying conventional testing procedures do not apply in this setting, because the required regularity conditions are violated. This has motivated a range of alternative testing procedures in the literature.

A range of procedures have been proposed for testing the null hypothesis of a linear model against a two-regime alternative (Hansen 1992, 1996a; Garcia 1998; Cho and White 2007; Marmer 2008; Carrasco et al. 2014; Kasahara et al. 2014; Dufour and Luger 2017; Qu and Zhuo 2021). Testing the null of an \(M\)-regime model against an \(M+m\)-regime alternative for \(M \geq 1\) and \(m = 1\) has been considered by Kasahara and Shimotsu (2018), who show that the parametric bootstrap test can be asymptotically valid under restrictions on the parameter space and for univariate models with fixed or predetermined regressors. Qu and Zhuo (2021) present similar bootstrap results for \(M = m = 1\) for a broader class of univariate models, though still under constrained parameter spaces. More recently, Rodriguez-Rondon and Dufour (2026) propose Monte Carlo likelihood ratio tests that accommodate both \(M \geq 1\) and \(m \geq 1\) and extend to multivariate settings, which had not been addressed previously. Their procedures are the most general available and deal transparently with the violations of regularity conditions that arise in this setting. Importantly, they do not require parametric restrictions, normality of the errors, or stationarity of the underlying process, since they do not rely on the existence of an asymptotic distribution. As a result, they apply in a wider range of settings than the parametric bootstrap, including cases where the bootstrap’s asymptotic validity has not been established. The maximized Monte Carlo version of the test further controls size in finite samples, which is particularly relevant for macroeconomic applications with quarterly data that are often short, and is robust to the identification problems that are common in Markov switching models.

Because carrying out these tests is not always straightforward, and conducting more than one can quickly become tedious, we introduce MSTest, an R package designed to test the null hypothesis of \(M\) regimes against the alternative of \(M+m\) regimes for both univariate and multivariate models. By making these testing procedures readily available to a broad audience, the package allows users to compare different tests and determine the appropriate number of regimes for economic research and policy-relevant applications.

The package also allows users to simulate and estimate univariate and multivariate Markov switching models, as well as hidden Markov models, with estimation available through the expectation–maximization (EM) algorithm or maximum likelihood estimation (MLE). For computational efficiency, MSTest relies on Rcpp (Eddelbuettel and Balamuta 2018) and RcppArmadillo (Eddelbuettel et al. 2026) and provides parallel computing options, which are especially valuable given the computational burden of testing in the presence of nuisance parameters. The package implements the methodologies of Rodriguez-Rondon and Dufour (2026), Dufour and Luger (2017), Carrasco et al. (2014), and Hansen (1992). The parametric bootstrap discussed by Qu and Zhuo (2021) and Kasahara and Shimotsu (2018) is not implemented explicitly, but it can be carried out using specific settings in combination with the local Monte Carlo likelihood ratio test of Rodriguez-Rondon and Dufour (2026).

The paper first introduces the Markov switching and hidden Markov models, then the testing procedures the package implements and how they address the regularity and identification problems specific to this setting. It then turns to the main focus of the article: the MSTest package itself—its principal functions, a typical simulation, estimation, and testing workflow, and practical guidance that complements the CRAN documentation. A brief summary concludes.

2 Markov switching models

The MSTest package considers Markov switching models in which only the mean and variance are governed by the Markov process \(S_t\). We also consider the hidden Markov model, a special case with no autoregressive coefficients. In both cases, exogenous explanatory variables may be included.

The first-order Markov process \(S_t\) that governs changes in the parameters of the Markov switching model is unobserved and evolves according to a first-order ergodic Markov chain with an (\(M \times M\)) transition probability matrix given by \[\begin{align*} \textbf{P} & = \begin{bmatrix} p_{11} & \dots & p_{M1}\\ \vdots & \ddots & \vdots \\ p_{1M} & \dots & p_{MM} \end{bmatrix} \end{align*}\] Here, \(p_{ij} = P(S_t = j \mid S_{t-1} = i)\) is the probability that state \(i\) is followed by state \(j\), \(M\) is the number of regimes, and \(S_t \in \{1, \dots, M\}\). The columns of \(\mathbf{P}\) sum to one (i.e., \(\sum_{j=1}^M p_{ij} = 1\) for all \(i\)). We can obtain the ergodic probabilities using \[\begin{align*} \pmb{\pi} = (\mathbf{A}'\mathbf{A})^{-1}\mathbf{A}'\mathbf{e}_{M+1} \hspace{0.25cm} \& \hspace{0.25cm} \mathbf{A} = \begin{bmatrix} \mathbf{I}_M-\textbf{P}\\ \pmb{1}' \end{bmatrix} \end{align*}\] where \(\mathbf{e}_{M+1}\) is the \((M+1)\)th column of \(\mathbf{I}_{M+1}\). These ergodic probabilities tell us, on average, in the long run, the proportion of time the process \(S_t\) spends in each regime.

A Markov switching model can be expressed as \[\begin{align} y_t & = \mu_{S_t} + \sum^{p}_{k=1} \phi_{k} (y_{t-k} - \mu_{S_{t-k}}) + Z_t\beta_z + \sigma_{S_t}\epsilon_t \tag{1} \end{align}\] where, in a univariate setting, \(y_t\) is a scalar, \(Z_t\) is a \((1 \times q_{z})\) vector of exogenous variables whose coefficients do not depend on the latent Markov process \(S_t\), and \(\epsilon_t\) represents the error process, which, for example, may be distributed as a \(\mathcal{N}(0,1)\). The error term is multiplied by the standard deviation \(\sigma_{S_t}\), which may either depend on the Markov process or remain constant throughout (i.e., \(\sigma\)).

This model is labeled MSARmdl in MSTest when exogenous regressors \(Z_t\) are excluded, and MSARXmdl when included. Other error distributions, such as the Student-\(t\), may be added in future versions. Currently, MSTest simulates only normally distributed errors, so we focus on this setup.

For the Markov switching model in equation (1) with \(M\) regimes, the sample log likelihood conditional on the first \(p\) observations of \(y_{t}\) is given by \[\begin{align} L_{T} (\theta) = \mathrm{log} f(y_{1}^{T}|y^{0}_{-p+1};\theta) = \sum^{T}_{t=1} \mathrm{log} f(y_{t}|\mathcal{Y}_{t-1};\theta) \tag{2} \end{align}\] where \(\mathcal{Y}_{t-1} = \sigma\)-field\(\{\dots,Z_{t-1},y_{t-2}, Z_{t}, y_{t-1}\}\) and \(\theta = (\mu_{1}, \dots, \mu_{M}, \beta, \sigma_{1}, \dots, \sigma_{M}, vec(\textbf{P}))\) and \(vec(\cdot)\) is the vectorization operator, which stacks the columns of a matrix to form a column vector.

The conditional density \(f(y_{t}\,|\,\mathcal{Y}_{t-1};\theta)\) is obtained by summing the regime-specific densities over the possible configurations of the current and previous \(p\) regimes, weighted by their conditional probabilities. The explicit expression can be found in (Hamilton 1989; Hamilton 1994; Rodriguez-Rondon and Dufour 2026).

Krolzig (1997) generalized the univariate model to the Markov switching vector autoregressive (MS-VAR) model, which in the MSTest package can be expressed as \[\begin{align} \pmb{y}_t & = \pmb{\mu}_{S_t} + \pmb{\Phi}_1(\pmb{y}_{t-1} - \pmb{\mu}_{S_{t-1}}) + \cdots + \pmb{\Phi}_p(\pmb{y}_{t-p} - \pmb{\mu}_{S_{t-p}}) + Z_{t}\pmb{\beta} + \pmb{\Sigma}^{1/2}_{S_t}\pmb{\epsilon}_t \tag{3} \end{align}\] where \(\pmb{y}_t = [y_{1,t}, \dots, y_{q,t}]'\), \(\pmb{\mu}_{S_t} = [\mu_{1,S_t}, \dots, \mu_{q,S_t}]'\), \(\pmb{\epsilon}_t = [\epsilon_{1,t}, \dots, \epsilon_{q,t}]'\), \(\pmb{\Phi}_k\) is a (\(q \times q\)) matrix containing the autoregressive parameters at lag \(k\), \(\pmb{\beta}\) is now a (\(q_z \times q\)) matrix, and \(\pmb{\Sigma}_{S_t}=\pmb{\Sigma}^{1/2}_{S_t}(\pmb{\Sigma}^{1/2}_{S_t})'\) is the (\(q \times q\)) regime-dependent covariance matrix. As in the univariate setting, the MSTest package includes versions with (MSVARXmdl) and without (MSVARmdl) exogenous regressors. More sophisticated versions of the MS-VAR model, together with their likelihood functions, are described in Krolzig (1997).

Hidden Markov models (HMMs) are a special case of the Markov switching model defined above. Because they need not be applied to time series data, they typically exclude lags of the endogenous variable \(y_t\). Removing them from (3) yields \[\begin{align} \pmb{y}_t & = \pmb{\mu}_{S_t} + Z_{t}\pmb{\beta} + \sigma_{S_t}\epsilon_t \tag{4} \end{align}\] When \(q=1\), we recover a univariate HMM from (1). This version and its multivariate counterpart are the HMMs considered in the package MSTest and are labeled as HMmdl.

As An et al. (2013) note, dependence on past observations allows richer interactions between \(y_t\) and \(S_t\), capturing more complex relationships between economic and financial variables. Including past observations is common in economic time series, often to control for stochastic trends, which may explain why Markov switching models are more popular than basic HMMs in this literature.

Markov switching models are typically estimated by the expectation–maximization (EM) algorithm (Dempster et al. 1977), Bayesian methods, or the Kalman filter in a state-space representation. They can also be estimated by maximum likelihood (MLE), which is commonly done, but only when the number of regimes \(M\) is small, since \(S_t\) is latent and the likelihood can exhibit several modes of equal height and other irregularities that make optimization harder as \(M\) grows.

In MSTest, the estimation functions described below use either the EM algorithm (control = list(method = "EM")) or MLE (control = list(method = "MLE")). Estimates can often be improved by using the EM results as starting values for the Newton-type (MLE) optimizer, a two-step procedure common in empirical work.

We omit the details of EM and MLE, since our focus is the MSTest package. Detailed treatments are given in (Hamilton 1990; Hamilton 1994), and, for MS-VAR models, in Krolzig (1997).

3 Hypothesis testing for number of regimes

The number of regimes in a Markov switching model must be specified by the researcher, since it is not determined endogenously during estimation. Testing for it is non-standard, however, because the regularity conditions underlying conventional procedures fail.

The general hypothesis of interest is \[\begin{align*} H_0: M = M_0 \quad \text{against} \quad H_1: M = M_0 + m \end{align*}\] where, in principle, both \(M_0, m \geq 1\), though most classical procedures only handle the simplest case, \(H_0: M = 1\) against \(H_1: M = 2\). Below we describe the main procedures, focusing on those in MSTest. The Monte Carlo LR procedures of Rodriguez-Rondon and Dufour (2026) are the most general, accommodating \(M_0 \geq 1\), \(m \geq 1\), non-Gaussian errors, non-stationary data, and multivariate models, though at a higher computational cost that can make other methods attractive in the simpler cases they were designed for. For each, we also note how it fits within the literature and handles identification failures and violations of regularity conditions.

Evidence on the size and power of the individual tests is available in the references given for each procedure below. In addition, Rodriguez-Rondon and Dufour (2026) provide a direct comparison of six of the seven tests implemented in MSTest, the LMC-LRT, MMC-LRT, the local and maximized moment-based tests of Dufour and Luger (2017), and the supTS and expTS tests (all but the Hansen test), under a common set of data-generating processes and show that all tests control size. Even in the simple \(M_0 = m = 1\) case, the LMC-LRT is the most powerful in the large majority of settings. For example, it has the highest power when only the mean switches, a case in which the moment-based and CHP tests can have substantially lower power, especially when the autoregressive persistence is high. All tests tend to improve when the variance switches, but even then the LMC-LRT typically leads, motivating its use even in simple settings, despite its greater computational cost.

3.1 Likelihood ratio–type tests

Hansen (1992) was the first to propose a testing procedure for Markov switching models when \(M_0 = m = 1\), and we begin with it, as it is available in MSTest (HLRTest(); see Table 4). Hansen discusses several issues that complicate the likelihood ratio approach in this setting. First, the likelihood is usually assumed locally quadratic near the null and the global optimum, but because some parameters are unidentified under the null, it is instead flat in those directions. Such unidentified nuisance parameters have been studied by (Davies 1977, 1987; Andrews and Ploberger 1994; Dufour 2006). Second, the score, usually assumed positive, can be identically \(0\) at the restricted maximum likelihood estimator of the linear (null) model. Third, some parameters, such as the transition probabilities, may take the boundary values \(0\) or \(1\), giving rise to the parameter boundary problem (Andrews 1999, 2001). In addition, the likelihood surface may have multiple local optima, so the null need not lie in the same region as the global optimum.

To address these issues, Hansen (1992) models the likelihood function as an empirical process of the unknown parameters and applies empirical process theory to obtain a bound on the asymptotic distribution of a standardized likelihood ratio statistic. The nuisance parameters, which include the second-regime parameters \((\mu_2, \sigma_2)\) and the transition probabilities \((p_{11}, p_{22})\), are not identified under the null, so the statistic is evaluated over a grid of values for these parameters and optimized accordingly. Because the underlying empirical process may be serially correlated for some parameter values, Hansen (1996a) proposes a correction that is applied when computing the asymptotic distribution. This correction is incorporated in the implementation in MSTest.

This procedure has two main drawbacks. First, it provides only a bound for the standardized likelihood ratio statistic, so it can be conservative, and the critical values reported by MSTest are those of the process \(Q\) rather than of the statistic itself. Second, it optimizes over the nuisance parameters by grid search, which becomes computationally intensive as more parameters are allowed to switch. We nonetheless include it because it is a common benchmark in the literature.

Several likelihood ratio–based procedures have followed Hansen (1992). Garcia (1998) reduce the dimensionality of the nuisance parameter space but rely on asymptotic assumptions known to be problematic in Markov switching settings (Andrews and Ploberger 1994; Hansen 1996b), and Cho and White (2007) develop a quasi-likelihood ratio test that accounts for parameters on the boundary but may not properly handle the time dependence induced by autoregressive dynamics (Cho and White 2011; Carter and Steigerwald 2012); for these reasons neither is included in MSTest. More recently, Qu and Zhuo (2021) and Kasahara and Shimotsu (2018) establish the asymptotic validity of parametric bootstrap procedures for likelihood ratio tests, the former for a broad class of univariate models and the latter for testing \(M_0 \geq 1\) against \(M_0 + 1\) regimes under fixed or predetermined regressors. While these parametric bootstrap procedures are not implemented directly in MSTest, they can be replicated using the local Monte Carlo likelihood ratio test of Rodriguez-Rondon and Dufour (2026) described below, which applies in even more general settings.

Rodriguez-Rondon and Dufour (2026) propose the maximized and local Monte Carlo likelihood ratio tests (MMC-LRT and LMC-LRT), the most general procedures currently available (implemented as MMCLRTest() and LMCLRTest(); see Table 4). They apply when both \(M_0 \ge 1\) and \(m \ge 1\), remain valid under non-stationarity, non-Gaussian errors, and parameters on the boundary, and are the only tests available for multivariate Markov switching models, including the MS-VAR and multivariate hidden Markov models discussed above. The MMC-LRT additionally controls size in finite samples and is robust to identification problems, which is especially relevant for macroeconomic applications.

Here, we provide a brief overview of the MMC-LRT and LMC-LRT procedures. A formal treatment is given in Rodriguez-Rondon and Dufour (2026). Both procedures are based on the likelihood ratio statistic \(LR_{T}=2[\bar{L}_{T}(H_{1})-\bar{L}_{T}(H_{0})]\), where \(\bar{L}_{T}(H_{1})\) and \(\bar{L}_{T}(H_{0})\) are the maximized log-likelihoods under the alternative and the null, with the log-likelihood as defined above. Under a null of \(M=1\), the parameter vector \(\theta_{0}=(\mu,\sigma^{2},\phi_{1},\ldots,\phi_{p})'\) collects the nuisance parameters, since the null distribution of \(LR_{T}\) depends on \(\theta_{0}\). Because the model is parametric, for any \(\theta_{0}\) we can simulate \(N\) i.i.d. replications \(LR_{T}^{(1)}(\theta_{0}),\ldots,LR_{T}^{(N)}(\theta_{0})\) of the statistic, and, provided the observed statistic \(LR_{T}^{(0)}\) and the simulated statistics are exchangeable for some \(\theta_{0}\in\bar{\Omega}_{0}\), the Monte Carlo \(p\)-value is

\[\begin{equation} \hat{p}_{N}[x\,|\,\theta_{0}]=\frac{N+1-R_{LR}[LR_{T}^{(0)};\,N]}{N+1} \tag{5} \end{equation}\] where \(R_{LR}[LR_{T}^{(0)};\,N] =\sum_{i=1}^{N}I[LR_{T}^{(0)}\geq LR_{T}^{(i)}(\theta _{0})]\) and \(I(C):=1\) if condition \(C\) holds, and \(I(C)=0\) otherwise. As can be seen, \(R_{LR}[LR_{T}^{(0)};\) \(N]\) simply computes the rank of the test statistic from the observed data within the generated series \(LR(N,\)\(\theta_{0})\). As shown in Rodriguez-Rondon and Dufour (2026), a critical region for this test statistic with level \(\alpha\) is given by \(\sup_{\theta_{0}\in \bar{\Omega}_{0}}\,\hat{p}_{N}[LR_{T}^{(0)}\,|\,\theta_{0}]\leq \alpha\). More precisely, if \((N+1)\alpha\) is an integer, then \[\begin{equation*} \mathbb{P}\left[ \sup \{\hat{p}_{N}[LR_{T}^{(0)}\,|\,\theta_{0}]:\theta_{0}\in \bar{\Omega}_{0}\}\leq \alpha \right] \leq \alpha \end{equation*}\] under the null hypothesis, so it is a valid test with level \(\alpha\) for \(H_{0}\). Because this result does not depend on the sample size \(T\), the test is also valid in finite samples.

This maximized Monte Carlo procedure requires searching for the maximum Monte Carlo \(p\)-value over the nuisance space \(\bar{\Omega}_{0}\). Because this space can be large, following Dufour (2006) the maximization is carried out over a consistent set \(C_T\) built around a consistent estimate \(\hat{\theta}_{0}\). MSTest lets users define \(C_T\) from a confidence interval around \(\hat{\theta}_{0}\) (\(C_T^{CI}\)), from a fixed-radius set (\(C_T^{\epsilon}\)), or from their union.

As discussed in Dufour (2006) and Rodriguez-Rondon and Dufour (2026), the solution to this optimization problem need not be unique, meaning that the maximum Monte Carlo \(p\)-value may correspond to multiple parameter vectors. For this reason, derivative-free numerical optimization methods are recommended to locate the maximum Monte Carlo \(p\)-value over the nuisance parameter space. MSTest allows users to employ the generalized simulated annealing algorithm (GenSA), the genetic algorithm (GA), and particle swarm optimization (pso); see (Dufour 2006; Scrucca 2013; Xiang et al. 2013; Zambrano-Bigiarini et al. 2013; Dufour and Neves 2019).

Choosing \(C_T\) to be the singleton \(\{\hat{\theta}_{0}\}\) yields the local Monte Carlo test (LMC-LRT), a finite-sample analogue of the parametric bootstrap that requires only a consistent estimate \(\hat{\theta}_{0}\) and, unlike the bootstrap, does not require \(N\to\infty\). As shown in Rodriguez-Rondon and Dufour (2026), both the LMC-LRT and the MMC-LRT remain valid even when no asymptotic distribution exists, and the parametric bootstrap procedures of Qu and Zhuo (2021) and Kasahara and Shimotsu (2018) can be reproduced as special cases by constraining the nuisance space accordingly.

3.2 Moment-based tests

Dufour and Luger (2017) propose an alternative approach to testing Markov switching models that avoids the statistical issues associated with likelihood ratio–type tests, is computationally less demanding, and perfectly controls the size of the test through the Monte Carlo testing methods developed in Dufour (2006) (implemented as DLMCTest() and DLMMCTest(); see Table 4). However, their method is restricted to the case where \(M_0 = m = 1\).

The moment-based test of Dufour and Luger (2017) is based on computing moments of the least squares residuals from autoregressive models estimated under the null hypothesis of \(M_0=1\). More specifically, they focus on the mean, variance, skewness, and excess kurtosis of the least squares residuals. These moments are calculated as \[\begin{equation*} M(\hat{\epsilon}) = \frac{|m_2 - m_1|}{\sqrt{s^2_1+s^2_2}}, \quad V(\hat{\epsilon}) = \frac{\vartheta_2(\hat{\epsilon})}{\vartheta_1(\hat{\epsilon})}, \quad S(\hat{\epsilon})= | \frac{\Sigma^T_{t=1} \hat{\epsilon}_{t}^3}{T(\hat{\sigma}^2)^{3/2}} |, \quad K(\hat{\epsilon})= | \frac{\Sigma^T_{t=1} \hat{\epsilon}_{t}^4}{T(\hat{\sigma}^2)^{2}} -3| \end{equation*}\] where the truncated means \(m_1, m_2\), the truncated variances \(s_1^2, s_2^2\), the variance components \(\vartheta_1, \vartheta_2\), and \(\hat{\sigma}^2 = T^{-1}\sum_{t=1}^{T}\hat{\epsilon}_t^2\) are computed from the negative and positive parts of the least squares residuals, as defined in Dufour and Luger (2017).

The testing procedure computes a test statistic for each moment, obtains the corresponding individual \(p\)-values, and combines them using one of two standard methods for aggregating independent tests: a minimum-based combination (Tippett 1931; Wilkinson 1951), or a product-based combination (Fisher 1932; Pearson 1933) (see Dufour et al. 2004; Dufour et al. 2015 for further discussion). The Monte Carlo \(p\)-value of the combined test statistic is then obtained from its rank among the values simulated under \(\eta \sim N(0, I_T)\), in the same way as for the likelihood ratio statistic above.

Because the transition probabilities, mean, and variance need not be treated as nuisance parameters—only the parameters on explanatory variables may be unidentified under the null—the test avoids much of the identification difficulty that arises in Hansen (1992), Garcia (1998), and Carrasco et al. (2014), and it extends easily to a maximized Monte Carlo framework. Although it is limited to comparing linear models against two-regime alternatives only, it is the least computationally intensive procedure available and can be computed in seconds, even in its maximized Monte Carlo version.

3.3 Optimal test for regime switching

Carrasco et al. (2014) propose a test for parameter constancy in random coefficient and Markov switching models that extends the information matrix test of White (1982). Like the moment-based test of Dufour and Luger (2017), it applies only to the case where \(M_0 = m = 1\) and only requires estimation under the null, which is attractive in Markov switching settings where estimating the alternative is computationally demanding because of nonlinearity and multiple local optima. In contrast, likelihood ratio–based procedures such as those of Hansen (1992), Garcia (1998), Cho and White (2007), Qu and Zhuo (2021), Kasahara and Shimotsu (2018), and Rodriguez-Rondon and Dufour (2026) require estimation under both hypotheses.

Although Carrasco et al. (2014) show that their test is asymptotically locally equivalent to the likelihood ratio test, simulation evidence in (Qu and Zhuo 2021; Rodriguez-Rondon and Dufour 2026) suggests that likelihood ratio–based approaches can be more powerful in some settings, such as when only the mean switches and the autoregressive persistence is high. The procedure also involves bootstrapping and optimization over nuisance parameters, which can make it more computationally intensive than the moment-based test of Dufour and Luger (2017) or the Monte Carlo likelihood ratio tests of Rodriguez-Rondon and Dufour (2026), particularly when the variance switches.

Carrasco et al. (2014) test the null of parameter constancy, \(\theta_t = \theta_0\), against the alternative \(\theta_t = \theta_0 + \eta_t\), where the switching variable \(\eta_t\) is unobservable, stationary, and may depend on nuisance parameters \(\beta\). Their statistic uses the second derivatives of the log-likelihood and the outer products of the scores, as in the information matrix test, together with an extra term that captures the serial dependence of the time-varying coefficients, so it depends on the latent process \(\eta_t\) only through its second-order properties. To handle the nuisance parameters, Carrasco et al. (2014) follow Davies (1987) and set \(\eta_t = chS_t\), where \(c\) is the amplitude of the change, \(h\) its direction, and \(S_t\) a bounded, mean-zero AR(1) process, collecting the nuisance parameters in \(\beta = (c^2, h', \rho')\). From the second-order properties of the score they obtain the supremum-type statistic \[\begin{equation*} \text{supTS} = \sup_{\{h, p:||h||=1, \underline{\rho} < \rho < \bar{\rho}\}} \frac{1}{2} (max(0, \frac{\Gamma^*_T}{\sqrt{\hat{\epsilon^{*'}}\hat{\epsilon^*}}}))^2 \end{equation*}\] where \(\Gamma^*_T\) and \(\hat{\epsilon}^*\) are normalized quantities that do not depend on \(c^2\), defined in Carrasco et al. (2014). The test bootstraps over the distribution of the nuisance parameters, for which MSTest uses a uniform prior. Because \(c^2\) need not be bounded above, they also suggest an exponential-type statistic following Andrews and Ploberger (1994), \[\begin{equation*} expTS = \int\limits_{\{\underline{\rho} \leq \rho \leq \bar\rho, ||h||<1 \}} \Psi(h,\rho) d\rho dh \end{equation*}\] where \(\Psi(h,\rho)\) is defined in Carrasco et al. (2014).

The tests of Carrasco et al. (2014) have been widely used in empirical applications (Hamilton 2005; Warne and Vredin 2006; Kahn and Rich 2007; Dufrénot et al. 2011; Morley and Piger 2012), in testing MS-GARCH models (Hu and Shin 2008), and as benchmarks (Dufour and Luger 2017; Qu and Zhuo 2021; Rodriguez-Rondon and Dufour 2026). Given their broad use and optimality properties, they are included in MSTest (CHPTest(); see Table 4).

4 The R package MSTest

The MSTest package is designed to conduct hypothesis tests for the number of regimes in Markov switching models. Because many of these procedures require estimating restricted and unrestricted models and simulating under the null, the package also provides tools for simulating and estimating Markov switching models. These features primarily support the testing procedures and are not the main focus of the package.

What distinguishes MSTest is its implementation of testing procedures that remain valid under violations of standard regularity conditions, together with its use of Rcpp and parallel computing, which make it feasible to handle nuisance parameters and Monte Carlo–based procedures in practice. Although other software exists for simulating and estimating Markov switching models, this combination of generality, computational efficiency, and modern hypothesis tests sets MSTest apart.

4.1 Data sets

The MSTest package includes three samples of U.S. real GNP and one of U.S. real GDP, all readily accessible once the package is loaded. These series have been used by Hansen (1992), Carrasco et al. (2014), Dufour and Luger (2017), and Rodriguez-Rondon and Dufour (2026), among others, to test for the number of regimes in Markov switching models. Table 1 reports the label, series, source, and sample span of each.

Table 1: U.S. real GNP and GDP data sets included in MSTest

Label	Series	Source	Sample span
hamilton84GNP	Real GNP	Hamilton (1989)	1951Q2–1984Q4
chp10GNP	Real GNP	Carrasco et al. (2014)	1951Q2–2010Q4
USGNP	Real GNP	Rodriguez-Rondon and Dufour (2026)	1947Q2–2024Q2
USRGDP	Real GDP	Rodriguez-Rondon and Dufour (2026)	1947Q2–2024Q2

These data sets can be loaded using the following commands once MSTest has been attached:

data("hamilton84GNP", package = "MSTest")
data("chp10GNP", package = "MSTest")
data("USGNP", package = "MSTest")
data("USRGDP", package = "MSTest")

Each data set has three columns: a Date column of class Date (constructed with as.Date()), the level of real GNP or GDP (GNP or RGDP), and the corresponding growth rate (GNP_gr or RGDP_gr).

4.2 Simulation

Simulation plays a central role in MSTest, as many of the testing procedures rely on Monte Carlo methods to obtain valid critical values in the presence of nuisance parameters and identification failures. The simulation functions are used primarily by the hypothesis testing procedures that build sample or asymptotic null distributions, such as LMCLRTest, MMCLRTest, DLMCTest, DLMMCTest, and CHPTest. They may also be useful to users developing new estimation or testing procedures for Markov switching models who wish to assess performance through controlled experiments.

Table 2 lists the simulation functions provided by MSTest and the classes of processes they generate. Each function takes as input a List specifying the data-generating process (DGP). We present illustrative examples below, and a complete description of each function’s arguments is available in the package documentation on CRAN.

Table 2: Simulation functions available in MSTest

Function	Description
simuNorm	Normally distributed process, optionally with exogenous regressors.
simuAR	Autoregressive process, AR(\(p\)).
simuARX	AR(\(p\)) process with exogenous regressors.
simuVAR	Vector autoregressive process, VAR(\(p\)).
simuVARX	VAR(\(p\)) process with exogenous regressors.
simuMSAR	Markov switching AR(\(p\)) process.
simuMSARX	Markov switching AR(\(p\)) process with exogenous regressors.
simuMSVAR	Markov switching VAR(\(p\)) process.
simuMSVARX	Markov switching VAR(\(p\)) process with exogenous regressors.
simuHMM	Hidden Markov model (HMM), optionally with exogenous regressors.

As an illustration, consider the simuNorm function, which generates a univariate or multivariate Gaussian process. The required input is a List containing the sample size (n), the number of series (q) (where q=1 indicates a univariate process and q>1 a multivariate one), a (\(q \times 1\)) vector of means, and a (\(q \times q\)) covariance matrix. Users may optionally specify a burnin period to discard initial observations, which defaults to zero for simuNorm but is larger for autoregressive processes to reduce sensitivity to initialization, or supply a custom matrix of innovations via the eps argument. The code below applies these functions to simulate several linear and regime-switching processes, including multivariate normal, AR, VAR, Markov switching AR and VAR, and hidden Markov models.

mdl_norm <- list(n     = 500, 
                 q     = 2,
                 mu    = c(5, -2),
                 sigma = rbind(c(5.0, 1.5),
                               c(1.5, 1.0)))
simu_norm <- simuNorm(mdl_norm)

mdl_ar <- list(n     = 500, 
               mu    = 5,
               sigma = 1,
               phi   = c(0.75))
simu_ar <- simuAR(mdl_ar)

mdl_var <- list(n     = 500, 
                p     = 1,
                q     = 2,
                mu    = c(5, -2),
                sigma = rbind(c(5.0, 1.5),
                              c(1.5, 1.0)),
                phi   = rbind(c(0.50, 0.30),
                              c(0.20, 0.70)))
simu_var <- simuVAR(mdl_var)

mdl_hmm <- list(n     = 500, 
                q     = 2,
                mu    = rbind(c(5, -2),
                              c(10, 2)),
                sigma = list(rbind(c(5.0, 1.5),
                                   c(1.5, 1.0)),
                             rbind(c(7.0, 3.0),
                                   c(3.0, 2.0))),
                k     = 2,
                P     = rbind(c(0.95, 0.10),
                              c(0.05, 0.90)))
simu_hmm <- simuHMM(mdl_hmm)

mdl_ms <- list(n     = 500, 
               mu    = c(5,10),
               sigma = c(1,1),
               phi   = c(0.75),
               k     = 2,
               P     = rbind(c(0.95, 0.10),
                             c(0.05, 0.90)))
simu_msar <- simuMSAR(mdl_ms)

mdl_msvar <- list(n     = 500, 
                  p     = 1,
                  q     = 2,
                  mu    = rbind(c(5, -2),
                                c(10, 2)),
                  sigma = list(rbind(c(5.0, 1.5),
                                     c(1.5, 1.0)),
                               rbind(c(7.0, 3.0),
                                     c(3.0, 2.0))),
                  phi   = rbind(c(0.50, 0.30),
                                c(0.20, 0.70)),
                  k     = 2,
                  P     = rbind(c(0.95, 0.10),
                                c(0.05, 0.90)))
simu_msvar <- simuMSVAR(mdl_msvar)

Figure 1: Simulated linear (left column) and Markov switching (right column) processes. Abrupt shifts between regimes are visible in the Markov switching processes.

Figure 1 shows the simulated processes. Regime switching is clearly visible in the Markov switching specifications, most notably when comparing the autoregressive process (middle left) with its Markov switching counterpart (middle right), and similar differences arise between the VAR and MS-VAR processes.

4.3 Model estimation

MSTest can estimate the ten model classes listed in Table 3, together with their MSTest labels and specifications. Although Nmdl and HMmdl use multivariate notation, both automatically detect a univariate setting when given a \((T \times 1)\) vector. Labels containing an “X” allow exogenous regressors, though these can also be included in Nmdl and HMmdl despite their names.

Table 3: Models and their specifications available in MSTest.

Model	Label	Equation
N\((\pmb{\mu} + \pmb{x}_{t} \pmb{\beta}, \pmb{\Sigma})\)	Nmdl	\(\pmb{y}_t = \pmb{\mu} + \pmb{x}_{t} \pmb{\beta} + \pmb{\Sigma}^{1/2} \pmb{\epsilon}_t\)
AR(\(p\))	ARmdl	\(y_t = \mu + \sum^p_{k=1} \phi_k (y_{t-k} - \mu) + \sigma \epsilon_t\)
ARX(\(p\))	ARXmdl	\(y_t = \mu + \sum^p_{k=1} \phi_k (y_{t-k} - \mu) + \pmb{x}_{t} \pmb{\beta} + \sigma \epsilon_t\)
VAR(\(p\))	VARmdl	\(\pmb{y}_t = \pmb{\mu} + \sum^p_{k=1}(\pmb{y}_{t-k} - \pmb{\mu})\pmb{\Phi}_k + \pmb{\Sigma}^{1/2}\pmb{\epsilon}_t\)
VARX(\(p\))	VARXmdl	\(\pmb{y}_t = \pmb{\mu} + \sum^p_{k=1}(\pmb{y}_{t-k} - \pmb{\mu})\pmb{\Phi}_k + \pmb{x}_t \pmb{\beta} + \pmb{\Sigma}^{1/2}\pmb{\epsilon}_t\)
MS-AR(\(p\))	MSARmdl	\(y_t = \mu_{S_t} + \sum^p_{k=1} \phi_k (y_{t-k} - \mu_{S_{t-k}}) + \sigma_{S_t} \epsilon_t\)
MS-ARX(\(p\))	MSARXmdl	\(y_t = \mu_{S_t} + \sum^p_{k=1} \phi_k (y_{t-k} - \mu_{S_{t-k}}) + \pmb{x}_{t} \pmb{\beta} + \sigma_{S_t} \epsilon_t\)
MS-VAR(\(p\))	MSVARmdl	\(\pmb{y}_t = \pmb{\mu}_{S_t} + \sum^p_{k=1}\pmb{\Phi}_k(\pmb{y}_{t-k} - \pmb{\mu}_{S_{t-k}}) + \pmb{\Sigma}^{1/2}_{S_t}\pmb{\epsilon}_t\)
MS-VARX(\(p\))	MSVARXmdl	\(\pmb{y}_t = \pmb{\mu}_{S_t} + \sum^p_{k=1}\pmb{\Phi}_k(\pmb{y}_{t-k} - \pmb{\mu}_{S_{t-k}}) + \pmb{x}_{t} \pmb{\beta} + \pmb{\Sigma}^{1/2}_{S_t}\pmb{\epsilon}_t\)
HMM	HMmdl	\(\pmb{y}_t = \pmb{\mu}_{S_t} + \pmb{x}_{t} \pmb{\beta} + \pmb{\Sigma}^{1/2}_{S_t}\pmb{\epsilon}_t\)

The simulation functions from the previous subsection generated a multivariate hidden Markov process (\(q = 2\) series), a Markov switching autoregressive process (\(p = 1\) lag), and a Markov switching vector autoregressive process (\(p = 1\) lag, \(q = 2\) series). Each returns a list with the simulated data, the true latent state variable \(S_t\), and other DGP elements. These series serve as inputs to the corresponding estimation functions below, for the hidden Markov model and the two Markov switching models.

control <- list(msmu   = TRUE, 
                msvar  = TRUE,
                method = "EM",
                use_diff_init = 30)
mdl_est_hmm <- HMmdl(simu_hmm[["y"]], k = 2, control = control)
summary(mdl_est_hmm)

Hidden Markov Model
              coef     s.e.
mu_1,1    5.151200 0.119030
mu_2,1   -2.007400 0.053269
mu_1,2    9.628800 0.236840
mu_2,2    1.821400 0.117870
sig_11,1  4.973200 0.375510
sig_12,1  1.388800 0.140140
sig_22,1  0.962750 0.075585
sig_11,2  7.384600 0.909040
sig_12,2  3.006800 0.422390
sig_22,2  1.738000 0.227480
p_11      0.961300 0.010523
p_12      0.038702 0.010523
p_21      0.099480 0.026777
p_22      0.900520 0.026777

log-likelihood =  -1796.469
AIC =  3620.939
BIC =  3679.943

Residuals:
       Min       1Q     Median      3Q    Max
Y1 -6.5756 -1.68620 -0.0424050 1.72610 6.8022
Y2 -2.9919 -0.75449  0.0072404 0.75641 3.2839

control <- list(msmu   = TRUE, 
                msvar  = FALSE, 
                method = "EM",
                use_diff_init = 30)
mdl_est_msar <- MSARmdl(simu_msar[["y"]], p = 1, k = 2, control = control)
summary(mdl_est_msar)

Markov Switching Autoregressive Model
           coef     s.e.
mu_1   4.994000 0.160730
mu_2  10.041000 0.186710
phi_1  0.715960 0.031774
sig    0.895540 0.057329
p_11   0.953170 0.011205
p_12   0.046832 0.011205
p_21   0.092460 0.022379
p_22   0.907540 0.022379

log-likelihood =  -797.9881
AIC =  1611.976
BIC =  1645.677

Residuals:
      Min       1Q    Median      3Q    Max
Y1 -3.723 -0.62857 -0.021467 0.61808 2.8027

control <- list(msmu   = TRUE, 
                msvar  = TRUE,
                method = "EM",
                use_diff_init = 30)
mdl_est_msvar <- MSVARmdl(simu_msvar[["y"]], p = 1, k = 2, control = control)
summary(mdl_est_msvar)

Markov Switching Vector Autoregressive Model
              coef     s.e.
mu_1,1    4.897400 0.543750
mu_2,1   -2.099200 0.509290
mu_1,2   10.290000 0.581890
mu_2,2    2.141000 0.547620
phi_1,11  0.508240 0.063077
phi_1,12  0.287540 0.075495
phi_1,21  0.204040 0.031217
phi_1,22  0.687860 0.037718
sig_11,1  5.561300 0.460000
sig_12,1  1.911700 0.193760
sig_22,1  1.204800 0.103680
sig_11,2  7.440600 0.830910
sig_12,2  3.343700 0.412780
sig_22,2  2.115700 0.239960
p_11      0.927340 0.015909
p_12      0.072658 0.015909
p_21      0.136240 0.030757
p_22      0.863760 0.030757

log-likelihood =  -1868.071
AIC =  3772.142
BIC =  3847.969

Residuals:
       Min       1Q     Median      3Q    Max
Y1 -9.8948 -1.63520  0.0415660 1.59780 7.6673
Y2 -3.7363 -0.81008 -0.0036109 0.83766 3.2878

In all three examples, method = "EM" in the control list estimates the models by the expectation–maximization (EM) algorithm, while method = "MLE" instead performs direct maximum likelihood estimation. The arguments msmu and msvar control whether the mean and variance switch across regimes, and setting either to FALSE holds the corresponding parameter constant across regimes. For the Markov switching AR example we set msvar = FALSE because its data-generating process holds the variance constant across regimes.

The option use_diff_init = 30 runs the optimizer thirty times from different initial values, retaining the highest-log-likelihood solution as the primary output while storing all runs in the trace element. All estimation functions return S3 objects with dedicated print() and summary() methods. The summary() method reports parameter estimates, log-likelihood values, information criteria (AIC and BIC), and diagnostics such as residual quantiles, which in this controlled setting can be compared directly with the true DGP values.

The summary() method also reports asymptotic standard errors for the estimated parameters, computed from the inverse of the observed information matrix (the numerically approximated Hessian of the log-likelihood) evaluated at the estimate. Large within-regime standard errors can serve as an informal indication that a model with fewer regimes may be adequate, which in turn motivates the formal testing procedures that are the focus of MSTest. These standard errors are valid under the usual regularity conditions for a given, correctly specified \(M\)-regime model. They should be interpreted with caution near the boundary of the parameter space (for example, transition probabilities close to \(0\) or \(1\)) or under weak regime separation, where the normal approximation may be poor, and MSTest returns NA in such cases. They concern inference within a fixed \(M\)-regime model and are distinct from the non-standard problem of testing the number of regimes, which is what the test procedures described below are designed to address.

Figure 2: Simulated processes (black and blue), true regimes (red), and smoothed probabilities (green) from the estimated models. The smoothed probabilities closely track the true regime paths.

Figure 2 displays the simulated series (black and blue), the true latent regime states \(S_t\) (red, dashed), and the estimated smoothed regime probabilities from the fitted models (green, solid). The smoothed probabilities track the true regime changes closely across all three models, and particularly well for the hidden Markov and Markov switching autoregressive models. Regime classification is somewhat less precise over certain intervals in the Markov switching VAR case, reflecting the greater complexity of multivariate models, which can typically be mitigated by a larger sample size or more initial values.

4.4 Hypothesis testing

This section describes the hypothesis testing functions available in MSTest. Designed for ease of use, the package in most cases requires only the series \(y_t\) (and \(Z_t\) when applicable) and the number of regimes to be tested.

Table 4 summarizes the hypothesis testing functions currently implemented in MSTest and reports, for each, the number of regimes it can test, its applicability to multivariate models, its identification robustness, its validity under non-stationarity, and its relative computational cost, as discussed above for each procedure. The DLMCTest and DLMMCTest procedures were not originally designed for non-stationary data, but Rodriguez-Rondon and Dufour (2026) provide simulation evidence that they control size and retain power under a random walk, which is why they carry an asterisk (\(^{*}\)) in Table 4. Each returns the test statistic, its \(p\)-value, and, when applicable, critical values and parameter estimates under the null and alternative hypotheses.

Table 4: Hypothesis tests available in MSTest. Multiv.: multivariate models; Ident. robust: valid under identification failure without relying on an asymptotic distribution; Non-stat.: valid under non-stationarity.

Function	Description	\(H_0\): \(M_0\) vs. \(H_1\): \(M_0{+}m\)	Multiv.	Ident. robust	Non-stat.	Cost
LMCLRTest	Local Monte Carlo LRT	\(M_0\) & \(m \ge 1\)	\(\checkmark\)	\(\checkmark\)	\(\checkmark\)	High
MMCLRTest	Maximized Monte Carlo LRT	\(M_0\) & \(m \ge 1\)	\(\checkmark\)	\(\checkmark\)	\(\checkmark\)	High
DLMCTest	Local MC moment-based test	\(M_0 = m = 1\)	–	\(\checkmark\)	\(\checkmark^{*}\)	Low
DLMMCTest	Maximized MC moment-based test	\(M_0 = m = 1\)	–	\(\checkmark\)	\(\checkmark^{*}\)	Low
CHPTest	Optimal parameter-constancy test	\(M_0 = m = 1\)	–	–	–	Med.
HLRTest	LRT empirical bound	\(M_0 = m = 1\)	–	–	–	High

The reported critical values warrant clarification. For HLRTest(), they are those of the process \(Q\) introduced in Hansen (1992), which bounds the likelihood ratio statistic from above rather than providing critical values for the LR statistic itself. For LMCLRTest(), MMCLRTest(), DLMCTest(), and DLMMCTest, they are empirical quantiles of the simulated null distribution. Throughout the package, “critical values” therefore denotes, in a generic sense, rejection thresholds derived from the relevant reference distribution.

Monte Carlo likelihood ratio test

The local Monte Carlo likelihood ratio test (LMC-LRT) is implemented via LMCLRTest() (Table 4). Because it simulates the null distribution by estimating both the restricted and unrestricted models, the estimation options accepted by the model functions can be passed to the test through mdl_h0_control and mdl_h1_control. The same structure applies to the maximized Monte Carlo likelihood ratio test discussed below. These control lists specify whether the mean and variance switch across regimes, the estimation method, the number of initial values used in optimization, and related options. A separate number of initial values for the null distribution simulation can be set via use_diff_init_sim, which defaults to the value used on the observed data (in mdl_h0_control and mdl_h1_control) and is generally recommended. Beyond these estimation controls, LMCLRTest() also requires the number of lags p, the number of regimes under the null k0, and under the alternative k1.

lmc_control = list(N = 19,
                   mdl_h0_control = list(const  = TRUE,
                                         getSE  = FALSE),
                   mdl_h1_control = list(msmu   = TRUE,
                                         msvar  = TRUE,
                                         getSE  = FALSE,
                                         method = "EM",
                                         use_diff_init = 1))
lmclrt <- LMCLRTest(simu_msvar[["y"]], p = 1, k0 = 1 , k1 = 2, control = lmc_control)
summary(lmclrt)

Restricted Model
             coef
mu_1      6.77680
mu_2     -0.62311
phi_1,11  0.55850
phi_1,12  0.23027
phi_1,21  0.18764
phi_1,22  0.63895
sig_11    8.22430
sig_12    3.98700
sig_22    2.91830

log-likelihood =  -1938.175
AIC =  3894.349
BIC =  3932.263

Unrestricted Model
             coef
mu_1,1    6.66590
mu_2,1   -0.55236
mu_1,2    7.63580
mu_2,2   -0.35457
phi_1,11  0.60883
phi_1,12  0.22478
phi_1,21  0.22281
phi_1,22  0.64515
sig_11,1  6.25930
sig_12,1  2.18750
sig_22,1  1.33840
sig_11,2 14.00300
sig_12,2  9.91160
sig_22,2  8.27010
p_11      0.75303
p_12      0.24697
p_21      0.83510
p_22      0.16490

log-likelihood =  -1901.817
AIC =  3839.633
BIC =  3915.46

Rodriguez-Rondon & Dufour (2026) Local Monte Carlo Likelihood Ratio Test
         LRT_0  0.90%  0.95%  0.99% p-value
LMC_LRT 72.716 22.728 30.894 47.549    0.05

We illustrate it on the Markov switching vector autoregressive model with \(p=1\) lag simulated in the previous subsection, testing \(H_0: M = 1\) against \(H_1: M = 2\). Since the data were generated from a two-regime model, rejection is expected. The output reports the likelihood ratio statistic, Monte Carlo critical values, and the Monte Carlo \(p\)-value, and here the null of a single regime is rejected. Because the \(p\)-value is computed from a finite number of simulations, set here to N = 19, it is discrete and can only take values on the grid \(\{1/(N+1), 2/(N+1), \ldots, N/(N+1)\}\). With \(N = 19\) the smallest attainable value is \(1/(19+1) = 0.05\), so the reported \(p\)-value of \(0.05\) is the strongest possible evidence against the null at this simulation size. As discussed in (Dufour 2006; Rodriguez-Rondon and Dufour 2026), larger N can improve power but need not be large, since these methods do not approximate asymptotic critical values. Consistent with this, Dufour et al. (2004) show that using more than \(100\) replications (including the observed statistic) has little effect on power, so \(N = 99\) is typically recommended. We use N = 19 purely for illustration.

As noted earlier, LMCLRTest() can also replicate the parametric bootstrap procedures (Kasahara and Shimotsu 2018; Qu and Zhuo 2021). Setting mdl_h1_control = list(method = "MLE") enables direct maximum likelihood estimation and parameter constraints via mle_theta_low and mle_theta_upp, which allow constraining the parameters away from the boundary of the parameter space, as required for asymptotic validity (Kasahara and Shimotsu (2018) additionally constrain the variances). Because the bootstrap approximates asymptotic critical values, a larger N is typical, though lower values remain reasonable (Qu and Zhuo (2021) use \(N = 199\) and Kasahara and Shimotsu (2018) use \(N = 299\)) given the cost of repeatedly estimating Markov switching models.

The maximized Monte Carlo likelihood ratio test (MMC-LRT) is implemented via MMCLRTest(). It shares the LMC-LRT estimation controls for the restricted and unrestricted models and likewise requires p, k0, and k1, adding several options specific to the maximization step over the nuisance parameter space. The constant eps defines the consistent set over which the maximization is performed. Setting CI_union = TRUE constructs \(C_T^\ast = C_T^{CI} \cup C_T^{\epsilon}\), while eps = 0 with CI_union = TRUE restricts the search to \(C_T^{CI}\). The optimization algorithm is selected via type. By default the search stops when a Monte Carlo \(p\)-value of \(1\) is reached (the maximum), or, by setting threshold_stop, once the test fails to reject, that is, when the \(p\)-value exceeds the significance level \(\alpha\).

mmc_control = list(N = 19,
                   eps = 0.3,
                   threshold_stop = 0.05 + 1e-6,
                   type = "GenSA",
                   CI_union = FALSE,
                   silence = TRUE,
                   mdl_h0_control = list(const  = TRUE,
                                         getSE  = FALSE),
                   mdl_h1_control = list(msmu   = TRUE,
                                         msvar  = TRUE,
                                         getSE  = FALSE,
                                         method = "EM"),
                   maxit  = 100)
mmclrt <- MMCLRTest(simu_norm[["y"]], p = 0, k0 = 1 , k1 = 2, control = mmc_control)
summary(mmclrt)

Restricted Model
          coef
mu_1    5.0774
mu_2   -1.9382
sig_11  4.9946
sig_12  1.5268
sig_22  1.0552

log-likelihood =  -1688.442
AIC =  3386.884
BIC =  3407.957

Unrestricted Model
                coef
mu_1,1    5.33240000
mu_2,1   -1.77730000
mu_1,2    4.80550000
mu_2,2   -2.10980000
sig_11,1  5.85090000
sig_12,1  1.79770000
sig_22,1  1.03890000
sig_11,2  3.93790000
sig_12,2  1.14750000
sig_22,2  1.01540000
p_11      0.06294900
p_12      0.93705000
p_21      0.99983000
p_22      0.00016602

log-likelihood =  -1684.86
AIC =  3397.72
BIC =  3456.724

Rodriguez-Rondon & Dufour (2026) Maximized Monte Carlo Likelihood Ratio Test
         LRT_0 p-value
MMC_LRT 7.1641     0.7

We illustrate the MMC-LRT on the linear bivariate model with \(M = 1\) regime simulated from a normal distribution in the previous subsection, testing \(H_0: M = 1\) against \(H_1: M = 2\) with p = 0 since the series has no autoregressive dynamics. Because the data are single-regime, we expect no rejection. We set eps = 0.3, CI_union = FALSE, and use simulated annealing via type = "GenSA", with threshold_stop = 0.05 + 1e-6 so the search stops once the null fails to be rejected. Parallel computation can be enabled by setting workers to the number of cores to use (for example, workers = 8), distributing the null distribution simulation across them. Here the MMC-LRT fails to reject the null, as indicated by the reported Monte Carlo \(p\)-value, whose resolution again depends on the number of simulations N. Even so, the MMC-LRT remains the most general likelihood ratio–based test currently available for Markov switching models, maximizing over the nuisance parameters while retaining finite-sample validity even when identification fails.

Moment-based tests

The Monte Carlo moment-based test proposed by Dufour and Luger (2017) is implemented in MSTest via the functions DLMCTest() (local version) and DLMMCTest() (maximized version). As with the other testing procedures, the user supplies the series \(y_t\) and the number of autoregressive lags \(p\), and N sets the number of Monte Carlo replications used to compute the test (default \(N = 99\)). The procedure additionally approximates the distribution of the \(p\)-value of each moment-based statistic through an auxiliary simulation, with the number of replications set by simdist_N (default \(10{,}000\)).

Below we illustrate the moment-based local Monte Carlo test on the previously simulated Markov switching autoregressive process. The test is computationally inexpensive and completes almost instantaneously. Because we set getSE = TRUE, the summary() method reports the estimated restricted model, here an AR(1) with autoregressive coefficient \(\phi_1\), together with standard errors for its coefficients, followed by the test statistics for each of the four moments (mean, variance, skewness, and excess kurtosis), the combined statistic \(F(\varepsilon)\), the corresponding critical values (quantiles of the simulated distribution), and the Monte Carlo \(p\)-values. The value of \(\phi_1\) reported in the test block equals the restricted-model estimate above, because the local test simulates the null distribution at this consistent point estimate. As expected, the null hypothesis of linearity is clearly rejected, since the data were generated from a Markov switching process with \(M = 2\) regimes.

lmc_control = list(N = 99,
                   simdist_N = 10000,
                   getSE = TRUE)
lmcmoment <- DLMCTest(simu_msar[["y"]], p = 1, control = lmc_control)

summary(lmcmoment)

Restricted Model
         coef     s.e.
mu    6.76700 0.450970
phi_1 0.84667 0.023987
sig   2.38540 0.151020

log-likelihood =  -924.9604
AIC =  1855.921
BIC =  1868.559

Dufour & Luger (2017) Moment-Based Local Monte Carlo Test
           phi_1 M(eps) V(eps)  S(eps) K(eps) F(eps)   0.90%   0.95%
LMC_min  0.84667 1.3431 16.942 0.37429  3.473      1 0.96333 0.97859
LMC_prod 0.84667 1.3431 16.942 0.37429  3.473      1 0.99822 0.99932
           0.99% p-value
LMC_min  0.99455    0.01
LMC_prod 0.99977    0.01

This computational efficiency carries over to the maximized moment-based test, implemented via DLMMCTest(), which maximizes the Monte Carlo \(p\)-value over the nuisance parameter space while remaining fast in practice. The optimization algorithm is selected via optim_type (the analogue of type in MMCLRTest()), and, as with MMCLRTest(), eps and CI_union define the consistent set over which the maximization is performed. Many of the other optimization options available there also apply here.

mmc_control <- list(N = 99,
                    getSE = TRUE,
                    eps = 0,
                    CI_union = TRUE,
                    optim_type = "GenSA",
                    silence = TRUE,
                    threshold_stop = 0.05 + 1e-6,
                    maxit = 100)
mmcmoment <- DLMMCTest(simu_msar[["y"]], p = 1, control = mmc_control)

summary(mmcmoment)

Restricted Model
         coef     s.e.
mu    6.76700 0.450970
phi_1 0.84667 0.023987
sig   2.38540 0.151020

log-likelihood =  -924.9604
AIC =  1855.921
BIC =  1868.559

Dufour & Luger (2017) Moment-Based Maximized Monte Carlo Test
           phi_1 M(eps) V(eps)  S(eps) K(eps) F(eps) p-value
MMC_min  0.84667 1.3431 16.942 0.37429  3.473      1    0.01
MMC_prod 0.84667 1.3431 16.942 0.37429  3.473      1    0.01

In the example above, we again apply the test to the simulated Markov switching autoregressive process. We set the stopping threshold to \(0.05 + 1\text{e-}6\), so that the optimization terminates once the null hypothesis fails to be rejected, although this does not occur in this case. We also set eps = 0, restricting the search to the confidence interval–based consistent set, as in Dufour and Luger (2017). The test block also reports the value of \(\phi_1\) that maximizes the Monte Carlo \(p\)-value over this set, which here coincides with the point estimate from the restricted model, although in general the two need not be equal. Once again, the null hypothesis of a linear model is rejected, in line with the results obtained from the other testing procedures.

Parameter stability test

The parameter stability test proposed by Carrasco et al. (2014) is implemented in MSTest via the function CHPTest(). As with the other testing procedures, the user must provide the series \(y_t\) and the number of autoregressive lags \(p\). Here the parameter \(N\) specifies the number of bootstrap replications used to obtain critical values and is set to \(3000\) by default, following Carrasco et al. (2014). The parameter \(\rho_b\) bounds the autocorrelation nuisance parameter \(\rho\) and is set to \(0.7\) by default, again as in their study, so that the grid search is performed over \(\rho \in [-\rho_b, \rho_b] = [-0.7, 0.7]\).

chp_control = list(N = 1000,
                   rho_b = 0.7,
                   msvar = FALSE)
pstabilitytest <- CHPTest(simu_ar[["y"]], p = 1, control = chp_control)
summary(pstabilitytest)

Restricted Model
         coef     s.e.
mu    4.91830 0.182520
phi_1 0.74989 0.029756
sig   1.03960 0.065817

log-likelihood =  -717.7421
AIC =  1441.484
BIC =  1454.122

Carrasco, Hu, & Ploberger (2014) Parameter Stability Test 

- Switch in Mean only
      test-stat  0.90%  0.95%  0.99% p-value
supTS   0.37868 1.8429 2.5673 4.1247   0.574
expTS   0.77573 1.3717 1.8293 3.4564   0.541

Above, we apply the parameter stability test to the linear autoregressive process simulated earlier. We set N = 1000 to keep the example fast and msvar = FALSE, so that the variance is held constant and the test focuses on potential instability in the mean. The summary() method reports the parameter estimates of the restricted model, which is the only model that needs to be estimated for this test, along with the results for both the supTS and expTS versions of the test statistic. As expected, the test fails to reject the null hypothesis of parameter stability, consistent with the fact that the data were generated from a linear autoregressive model without regime switching.

Stochastic likelihood ratio test

The stochastic likelihood ratio test proposed by Hansen (1992) is implemented in MSTest via the function HLRTest(). As with DLMCTest(), DLMMCTest(), and CHPTest(), the user need only provide the series \(y_t\) and the number of autoregressive lags \(p\), since this test likewise assesses only the null hypothesis of linearity, that is, a single regime, in an autoregressive model.

This procedure performs a grid search over the nuisance parameters of the unrestricted model, so most options control the construction of that grid. The argument gridsize sets the number of grid points, and mugrid_from and mugrid_by give the starting value and step size for the mean grid, with analogous arguments for the variance grid. For each point on this grid, the test estimates the restricted model over a bounded region defined by theta_null_low and theta_null_upp. Setting msvar = TRUE additionally allows regime-dependent variances, but this enlarges the nuisance parameter space and can substantially raise the computational cost.

hlrt_control  <- list(msvar          = FALSE,
                      gridsize       = 20,
                      mugrid_from    = 0,
                      mugrid_by      = 1,
                      theta_null_low = c(0,-0.99,0.01),
                      theta_null_upp = c(20,0.99,20))
hlrt <- HLRTest(simu_msar[["y"]], p = 1, control = hlrt_control)
summary(hlrt)

Restricted Model
         coef     s.e.
mu    6.76700 0.450970
phi_1 0.84667 0.023987
sig   2.38540 0.151020

log-likelihood =  -924.9604
AIC =  1855.921
BIC =  1868.559

Hansen (1992) Likelihood Ratio Bound Test -  Switch in Mean only
      test-stat 0.90 % 0.95 % 0.99 % p-value
M = 0     6.972 2.6900 2.9574 3.4479       0
M = 1     6.972 2.7176 2.9612 3.4493       0
M = 2     6.972 2.7856 3.0960 3.5397       0
M = 3     6.972 2.8679 3.1516 3.6352       0
M = 4     6.972 2.8962 3.2477 3.8034       0

Above, we illustrate the stochastic likelihood ratio test using the previously simulated Markov switching autoregressive process. The summary() method reports parameter estimates for the restricted model, which is the only model that needs to be estimated for this test, together with the standardized likelihood ratio statistic and its Monte Carlo critical values and \(p\)-values for each bandwidth used in the heteroskedasticity-robust procedure of Hansen (1996a). Across all bandwidths, the null hypothesis of linearity is strongly rejected, with Monte Carlo \(p\)-values effectively equal to zero, consistent with the data having been generated from a Markov switching process with two regimes.

5 Conclusion

The importance of testing the number of regimes in Markov switching models has motivated a substantial body of research, reflecting the statistical and computational challenges inherent in this problem. Seminal contributions such as Hansen (1992), Carrasco et al. (2014), and Dufour and Luger (2017) focus on testing the null hypothesis of a single regime (i.e., a linear model) against the alternative of two regimes. More recently, Rodriguez-Rondon and Dufour (2026) introduced a class of Monte Carlo likelihood ratio tests that allow one to test a null of \(M_0\) regimes against an alternative of \(M_0 + m\) regimes, for any \(M_0 \geq 1\) and \(m \geq 1\) and for a broad class of models, including multivariate and non-stationary ones, substantially expanding the scope of feasible hypothesis testing.

The MSTest package makes these procedures readily available by implementing the methods of these four studies within a unified, user-friendly framework, and it additionally provides tools for simulating and estimating a wide range of Markov switching and hidden Markov models to support the testing procedures. This paper has reviewed the underlying methodology and illustrated how MSTest is used in practice to conduct inference on the number of regimes. By lowering the computational and implementation barriers to regime testing and offering multiple complementary procedures with differing assumptions and robustness properties, the package lets users select methods well suited to their empirical setting and make informed decisions about model specification.

Computational details

All results in this paper were obtained using R 4.4.0 (R Core Team 2024) and version 0.1.9 of the package MSTest, available from the Comprehensive R Archive Network (CRAN). For computational efficiency, MSTest depends on Rcpp 1.1.1 (Eddelbuettel and Balamuta 2018) and RcppArmadillo 15.2.4.1 (Eddelbuettel et al. 2026), and it uses GenSA 1.1.14.1 (Xiang et al. 2013), GA 3.2.4 (Scrucca 2013), and pso 1.0.4 (Bendtsen 2022) for numerical optimization.

All tables and figures can be reproduced with the file article.R, distributed with the MSTest package (inst/examples/article.R) and also available from the package’s public GitHub repository. Because the Monte Carlo procedures are simulation-based, exact reproduction requires the random seed set in article.R, though results can vary slightly across hardware, BLAS/LAPACK implementations, and operating systems. On Darwin 25.6.0 (arm64), the full script runs in approximately three minutes of elapsed time, which varies with hardware and parallelization settings.

5.1 Supplementary materials

Supplementary materials are available in addition to this article. It can be downloaded at RJ-2026-044.zip

5.2 CRAN packages used

MSTest, Rcpp, RcppArmadillo, GenSA, GA, pso

5.3 CRAN Task Views implied by cited packages

HighPerformanceComputing, NumericalMathematics, Optimization