Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
108,734 characters · 24 sections · 237 citation commands
mbMSTest: An R-Package for Testing Markov Switching Models
\thispagestyle{empty}
\setcounter{page}{1}
Markov switching models were first introduced by goldfeld1973markov, but they were later popularized and became an active area of research in economics after hamilton89 proposed modeling the first difference of U.S. GNP as a nonlinear stationary process rather than a linear stationary process, as was typically done. The nonlinearity here arises from discrete shifts in the process.
These models have now been considered in various macroeconomic and financial applications. For example, they have been used in the identification of business cycles to provide probabilistic statements about the state of the economy (see chauvet1998econometric; diebold1996measuring; kimnel1999; chauvet2006dating; quqin2021), in modeling stock market volatility with Markov switching ARCH, GARCH, and Stochastic Volatility models (see hamilton1994; gray1996; klaassen2002improving; haas2004new; pelletier2006regime; so1998stochastic), in modeling interest rate dynamics (see cai1994markov; garper1996), in considering state-dependent impulse response functions (see sims2006were; caggiano2017), in the identification of structural shocks in SVAR models (see lanne2010structural; herwartz2014structural; lutkepohl2021testing), and more recently in improving measures of core inflation with multiple inflation regimes (see rodron_inf_2024). hamilton2016 provides a detailed survey of regime switching in macroeconomics.
Outside of macroeconomic and financial applications, these models have also been applied in climate change research (see golosov2014optimal; dietz2015endogenous), environmental and energy economics (see cevik2021renewable; charfeddine2017impact), industrial organization (see aguirregabiria2007sequential; sweeting2013dynamic), and health economics (see hernandez2016switching; anser2021impact). Additionally, there is a related model—the Hidden Markov model—which has various applications in computational molecular biology (see krogh1994hidden; baldi1994hidden), handwriting and speech recognition (see rabiner1986introduction; nag1986script; RabJuanFundamentals; jelinek1997statistical), computer vision and pattern recognition (see bunke2001hidden), and other machine learning applications.
Given their empirical relevance, it is important to determine the number of regimes needed to properly capture the nonlinearities present in the data, as this is not determined endogenously when estimating Markov switching models. However, the asymptotic results of conventional hypothesis testing procedures do not apply in this setting because the regularity conditions required for such results are violated. Consequently, alternative hypothesis testing procedures have been proposed in the literature.
Notable contributions to testing the null hypothesis of a linear model against a model with two regimes include hansen92, hansen96, garcia98, chowhite07, marmer2008, chp14, kso2014modqlr, dufourluger17, and quzhuo2021likelihood. Testing the null of an $M$-regime model against the alternative of an $M+m$-regime model for $M \geq 1$ and $m = 1$ has been considered by kasshi2018, where the authors show that the parametric bootstrap test can be asymptotically valid when imposing certain restrictions on the parameter space and specifically considering univariate models with fixed or predetermined regressors. quzhuo2021likelihood present similar results regarding the parametric bootstrap for $M = m = 1$ but for a broader set of univariate models, however, still containing the parameter space. More recently, rodrondufour_mcmstest propose Monte Carlo likelihood ratio tests that handle cases where both $M \geq 1$ and $m \geq 1$ and also consider multivariate settings, which had not been addressed previously. Their test procedures are the most general procedure available and deal transparently with issues related to violations of regularity conditions. Importantly, they do not require parametric restrictions, normality of errors, or stationary processes, as the existence of an asymptotic distribution is unnecessary. This makes the tests applicable in more cases than the parametric bootstrap and allows for use in settings where the asymptotic validity of the parametric bootstrap has not been established. The maximized Monte Carlo version of their test even controls the test size in finite samples, which is particularly relevant for many macroeconomic applications using quarterly data. Additionally, this version of their test is robust to identification problems, which are common when working with Markov switching models.
Testing the number of regimes that a Markov switching model should include is an important step when deciding on the model's specification. However, performing these tests is not necessarily trivial, and conducting more than one test can quickly become tedious. As a result, we introduce {\fontseries{m}\fontseries{b}\selectfont MSTest}, an \proglang{R} package that can be used to test the null hypothesis of $M$ regimes against the alternative hypothesis of $M+m$ regimes for both univariate and multivariate models. The purpose of this \proglang{R} package is to enable econometricians to determine the number of regimes that a model should include for a given process by making hypothesis testing procedures readily available to a general audience. It aims to facilitate the comparison of different test procedures and the determination of the number of regimes in a model for economic research. The package also allows users to simulate and estimate univariate and multivariate Markov switching models, as well as hidden Markov models. Estimation is provided through the use of the expectation maximization (EM) algorithm or maximum likelihood estimation (MLE). {\fontseries{m}\fontseries{b}\selectfont MSTest} utilizes Rcpp (eddetal2018) and RcppArmadillo (eddsan2024) for computational efficiency. This is especially important given the computational burden of testing in the presence of nuisance parameters, as is the case here. The {\fontseries{m}\fontseries{b}\selectfont MSTest} package includes the methodologies presented in rodrondufour_mcmstest, dufourluger17, chp14, and hansen92. The parametric bootstrap discussed by quzhuo2021likelihood and kasshi2018 is not explicitly provided but can be performed using specific settings while employing the local Monte Carlo likelihood ratio test of rodrondufour_mcmstest. These testing procedures, along with other notable contributions previously mentioned, are discussed in more detail in Section ((ref)) of this paper.
In Section (ref), we describe Markov switching and hidden Markov models for which the hypothesis test are implemented. In Section (ref), we describe some theoretical aspects of the hypothesis testing procedures the {\fontseries{m}\fontseries{b}\selectfont MSTest} package offers. That is, we discuss the identification difficulties in more detail and the methodologies developed in each framework to test the hypothesis of interest. Section (ref) describes the package in more detail and the functions available to test for the number of regimes. This section is meant to compliment the {\fontseries{m}\fontseries{b}\selectfont MSTest} documentation available through CRAN, by giving a general overview of available functions and can also serve as a short manual for the {\fontseries{m}\fontseries{b}\selectfont MSTest} package. Section (ref) includes an empirical example where we apply the test procedures included in {\fontseries{m}\fontseries{b}\selectfont MSTest} to the U.S. GNP data of hamilton89, the extended data used in chp14 and dufourluger17 and a further extension of the data ranging from 1951Q2 to 2024Q2. We compare the different models and provide tables with: test statistics, critical values and p-values. Finally we provide brief concluding remarks in section (ref).
The {\fontseries{m}\fontseries{b}\selectfont MSTest} package examines Markov switching models where only the mean and variance are governed by the Markov process $S_t$ and so here, we describe such models. We also consider a specific case of the Markov switching model—the Hidden Markov model—in which no autoregressive coefficients are included as explanatory variables. In both cases, other exogenous explanatory variables may be included.
We begin by describing the first-order Markov process $S_t$ that governs the changes in the parameters of the Markov switching model. We assume that the process $S_t$ is unobserved and evolves according to a first-order ergodic Markov chain with a ($M \times M$) transition probability matrix given by
where for example $p_{ij} = P(S_t = j | S_{t-1} = i)$ is the probability of state $i$ being followed by state $j$ and $M$ is the total number of regimes. Specifically, if we consider $M$ regimes, the process takes integer values $S_t=\{1, \dots, M\}$. Additionally, the columns of the transition matrix $\textbf{P}$ must sum to one in order to have a well defined transition matrix (i.e., $\sum^M_{j=1} p_{ij} = 1$).
Considering the example in hamilton1994 where the Markov process has only two regimes, we only need a $(2 \times 2)$ transition matrix to summarize the transition probabilities $\textbf{P}$ as follows:
We can also obtain the ergodic probabilities, $\pi = (\pi_1, \pi_2)'$, which are given by
in a setting with two regimes. More generally, for any number of $M$ regimes we could use
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.
The Markov switching model considered in the {\fontseries{m}\fontseries{b}\selectfont MSTest} package can be expressed as
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 Markov switching autoregressive model is labeled as “\code{MSARmdl}” in the {\fontseries{m}\fontseries{b}\selectfont MSTest} package when $Z_t$ is excluded, or as “\code{MSARXmdl}” when exogenous regressors, $Z_t$, are included. These are the versions most commonly used in various economic and financial applications, as well as other time series-related applications. Other error distributions, such as a Student-t distribution, may also be considered in future versions of the package. Currently, the {\fontseries{m}\fontseries{b}\selectfont MSTest} package only considers Markov switching autoregressive models with normally distributed errors when simulating the processes, so we focus on this setup in this paper.
Continuing with the example where a Markov switching model given by equation ((ref)) has $M=2$ regimes, such that $S_t = \{1, 2\}$, the sample log likelihood conditional on the first $p$ observations of $y_{t}$ is given by
where $\mathscr{Y}_{t-1} = \sigma$-field$\{\dots,Z_{t-1},y_{t-2}, Z_{t}, y_{t-1}\}$ and $\theta = (\mu_{1}, \mu_{2}, \beta, \sigma_{1}, \sigma_{2}, vec(\textbf{P}))$ and $vec(\cdot)$ is the vectorization operator which stacks the columns of a matrix to form a column vector. Here,
and more specifically
where we set
and $\text{Pr}(S^{*}_{t}=s^{*}_{t}|\mathscr{Y}_{t-1};\theta)$ is the probability that this occurs. Note that as in, rodrondufour_mcmstest, here we denote the latent variable that determines the regimes at time $t$ as $S_{t}$ and let $s_{t}$ denote the (observed) realization of $S_{t}$.
krolzig1997markov generalized the the univariate Markov switching autoregressive model to the multivariate setting and hence introduced the Markov switching Vector Autoregressive (MS-VAR) model. The MS-VAR model considered in the {\fontseries{m}\fontseries{b}\selectfont MSTest} package can be expressed as
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 with the univariate setting, the {\fontseries{m}\fontseries{b}\selectfont MSTest} package also includes a version without exogenous regressors $Z_{t}$, labeled “\code{MSVARmdl}”, and a version which allows the inclusion of exogenous regressors, labeled “\code{MSVARXmdl}”. More sophisticated versions of the MS-VAR model, as well as their likelihood functions, are described in krolzig1997markov and we direct the interested reader to consider this reference to learn more about these models and their components.
Hidden Markov models (HMM) can be shown to be a special case of the more general Markov switching model we defined above. Specifically, Hidden Markov models don't necessarily have to be applied to time series data and for this reason they typically do not include lags of the endogenous variable, $y_{t}$.
For example, we can recover a hidden Markov from ((ref)) by simply excluding lags of $y_t$ as explanatory variables giving
When $q=1$, we recover a univariate Hodden Markov model from ((ref)). This version and its multivariate counterpart are the HMMs considered in the package {\fontseries{m}\fontseries{b}\selectfont MSTest} and are labeled as “\code{HMmdl}”.
As described by an2013identifiability, the dependence on past observations allows for more general interactions between $y_t$ and $S_t$ which can be used to model more complicated causal links between economic or financial variables of interest. Including past observations is a very common practice in economic time series applications as a way to control for stochastic trends, which may explain why Markov-switching models are more popular than the basic HMM in this literature.
Typically, Markov switching models are estimated using the Expectation-Maximization (EM) algorithm (see dempster_maximum_1977), Bayesian methods, or with the Kalman filter when employing the state-space representation of the model. In very simple cases, Markov switching models can be estimated using Maximum Likelihood Estimation (MLE). However, since the Markov process $S_{t}$ is latent and, more importantly, because the likelihood function may exhibit several modes of equal height along with other unusual features that complicate MLE estimation, this approach is less commonly used.
The \proglang{R} package {\fontseries{m}\fontseries{b}\selectfont MSTest} enables estimation of the models described above via the EM algorithm by setting “\code{control = list(method = `EM')}” or via MLE by setting “\code{control = list(method = `MLE')}” in the estimation functions, which are detailed below in Section (ref). In practice, empirical estimates can sometimes be improved by using the EM algorithm results as initial values in a Newton-type optimization algorithm. This two-step estimation procedure is used to obtain results presented in the empirical section of rodrondufour_mcmstest, and in other related works.
We omit a detailed explanation of the EM algorithm and MLE, as our focus is on describing the {\fontseries{m}\fontseries{b}\selectfont MSTest} package. For the interested reader, the estimation of a Markov switching model via the EM algorithm and MLE is described in detail in hamilton1990 and hamilton1994, and for Markov-switching VAR models, in krolzig1997markov.
When estimating a Markov switching model, it is essential to determine the number of regimes to be estimated, as this is not decided endogenously during the estimation process. However, it is well understood in the literature that when considering testing for the number of regimes in a Markov switching model, conventional hypothesis testing procedures are no longer valid as the typical regularity conditions needed for asymptotic validity are not met. To address these challenges, various studies have proposed alternative methods that yield valid testing procedures. In this section, we describe some key procedures, focusing on those included in the \proglang{R} package {\fontseries{m}\fontseries{b}\selectfont MSTest}. As mentioned in the introduction, the purpose of {\fontseries{m}\fontseries{b}\selectfont MSTest} is to make the most useful of these procedures accessible to a general audience. Here, we explain how these procedures fit within the existing literature on hypothesis testing for Markov switching models and briefly review how these procedures circumvent these violations of regularity conditions.
In general, the hypothesis test of interest is
where both \(M_0, m \geq 1\). However, most available test procedures can only address the case when $M_0=m=1$ and so
In this case, a linear model (one regime) is being considered under the null hypothesis and is being compared against a Markov switching model with two regimes under the alternative hypothesis. Currently, the most general procedure able to deal with settings where \(M_0, m \geq 1\) are the Monte Carlo likelihood ratio tests described in rodrondufour_mcmstest but this flexibility comes at the cost of, relatively speaking, being computationally intensive. For this reason {\fontseries{m}\fontseries{b}\selectfont MSTest} also includes other procedures such as the moment-based tests of dufourluger17 and the parameter stability test of chp14, which are computationally efficient and useful when only considering the more simple case of $M_0=m=1$. Moreover, the package includes the standardized likelihood ratio test of hansen92, which is a significant contribution and has frequently served as a benchmark for evaluating other test procedures.
hansen92 was the first to propose a testing procedure for Markov switching models when $M_0=m=1$, so we begin with a review of this procedure, which is available in {\fontseries{m}\fontseries{b}\selectfont MSTest}. The author provides a thorough description of the problems that plague the likelihood ratio approach for testing the number of regimes in a Markov switching model. First, it is typically assumed that the likelihood function is locally quadratic in the region where the null hypothesis and the globally optimal estimated parameters are found. However, as the author notes, since some parameters are not identified under the null, this region is likely flat with respect to those unidentified parameters rather than quadratic. Unidentified nuisance parameters under the null hypothesis are issues that have been considered in davies1977, davies87, andrewsploberger94, and dufour2006. Second, it is commonly assumed that the score is positive; however, as described, it can be identically $0$ under the restricted maximum likelihood estimator (MLE) of a linear model (the null hypothesis). Third, some parameters, such as the transition probabilities, can take values of $0$ and $1$, which leads to the parameter boundary problem discussed in andrews1999 and andrews2001. Additionally, the likelihood surface can have multiple local optima, meaning the null hypothesis may not lie in the same region as the global optimum.
hansen92 introduces a new approach that does not require the conventional assumptions associated with likelihood ratio tests. Instead, the authors model the likelihood function as an empirical process of the unknown parameters. They utilize empirical process theory to establish a bound for the asymptotic distribution of the standardized likelihood ratio test. hansen92 formulates the hypothesis as follows:
where $\alpha_0$ represents the parameter values under the null and $\alpha$ the parameter values under the alternative. They begin by decomposing the likelihood ratio as such:
where $R_n(\alpha) = E[LR_n(\alpha)]$ is the expectation of the likelihood ratio function and $Q_n(\alpha) = \Sigma_{i=1}^n q_i(\alpha) = [l_i(\alpha) - l_i(\alpha_0)] - E[l_i(\alpha) - l_i(\alpha_0)]$ is the deviations from the mean. Fluctuations in Q play an important role in identifying an optimum as:
and by using the fact that \(R_n(\alpha) \leq 0 \) for all \(\alpha\) when the null hypothesis is true, we can see that \(\frac{1}{\sqrt{n}} LR_n (\alpha) \leq \frac{1}{\sqrt{n}} Q_n (\alpha)\) and so it follows that,
Thus, the distribution of the empirical process Q can provide a bound for the asymptotic distribution of the LR statistic. The test statistic is further standardized:
where
As suggested by equation ((ref)), they resolve the issue of nuisance parameters by evaluating the standardized LR statistic for different parameter values $\alpha$. Specifically, they evaluate the test statistic over a grid of different parameter values and optimize with respect to those nuisance parameter values. To be clear, hansen92 sets $\alpha = (\mu_2, \sigma_2, p_{11},p_{22})$ as the vector of nuisance parameters, which includes the second state parameters and transition parameters. The first regime parameters $\theta=(\mu_1,\sigma_1,\phi_1,\dots,\phi_p)$ are fully identified. They further split $\alpha$ into $\beta =(\mu_2, \sigma_2)$ and $\gamma = (p_{11},p_{22})$, where $\beta$ is treated as a parameter of interest and set to be identical to $\mu_1$ and $\sigma_1$ under the null. However, the process $Q$ may also be serially correlated for some values of $\alpha$, and so hansen96 adds a correction that should be used for calculating the asymptotic distribution of the test statistic, which is also included in the implementation of this test procedure in {\fontseries{m}\fontseries{b}\selectfont MSTest}.
There are two main drawbacks to consider when using this likelihood ratio procedure for testing Markov switching models. The first drawback is that this test only provides a bound for the standardized LR statistic, which can be conservative. Therefore, it is important to remember that the critical values provided by this test in the {\fontseries{m}\fontseries{b}\selectfont MSTest} package are not necessarily those of the standardized LR statistic but rather of the process $Q$, which provides a bound for the standardized LR statistic. The second drawback is that it involves optimizing the value of the nuisance parameters through a grid search. Although this process is manageable when analyzing only a few variables, as it requires switching between regimes in addition to the transition parameters $p_{11}$ and $p_{22}$, it can quickly become computationally intensive for models that consider more parameters as switching between regimes. Despite some of these drawbacks, we include this test in {\fontseries{m}\fontseries{b}\selectfont MSTest} because its properties are well understood and it has often been used as a benchmark for comparison in the literature.
garcia98 also reviewed the problem of testing the number of regimes in a Markov switching models using a likelihood ratio approach. garcia98 builds on hansen92's approach but differs in that they only treat $\gamma$ as nuisance parameters over which we must optimize. This change simplifies some of the computational burden. The $\beta$ parameter remains a parameter of interest but is incorporated into $\theta$, the identified parameters. Although this is a significant contribution and the authors provide valuable insights, this test is not included in this version of {\fontseries{m}\fontseries{b}\selectfont MSTest}. One reason is that the author assumes that the LR test can be expressed as the supremum of a chi-square functional asymptotically under the null hypothesis, a claim that hansen96inference and andrewsploberger94 suggest cannot be made for Markov switching models. However, it may still be included in future versions of {\fontseries{m}\fontseries{b}\selectfont MSTest} for completeness.
chowhite07 also considers hypothesis testing when the null hypothesis is a linear model and the alternative is a Markov switching model with two regimes. They address a difficulty not covered by hansen92 or garcia98, specifically the case where parameters lie on the boundary of the parameter space. The authors classify the null hypothesis into two mutually exclusive subsets: one where $p \in (0,1)$ and another where $p = 0$ or $p = 1$, where $p$ is a transition probability (e.g., $p_{11}$). The second case presents the boundary parameter problem. Building on the work of andrews1999 and andrews2001, the authors develop a QLR test statistic that accounts for this boundary issue. Additionally, chowhite07 employs a normal mixture model framework to formulate the QLR test within this context. They argue that this approach allows them to disregard certain time series dependence properties implied by the Markov regime-switching process. The authors describe the QLR test as sensitive to the mixture aspect of the regime-switching process, claiming it delivers a test with strong power under the alternative hypothesis. A significant finding of their work is that the critical values of the asymptotic distribution of the test statistic are influenced by the consideration of the boundary problem. In other words, they demonstrate the importance of addressing the boundary problem when examining the asymptotic distribution of the test statistic. However, carsteig2012 highlight a critical flaw in this test, arguing that it may overlook the time dependency of the Markov chain but fails to account for the time dependencies present in an autoregressive model when parameters change across regimes. The authors acknowledge this limitation in chowhi2011. For this reason, the test proposed by chowhite07 is not currently included in {\fontseries{m}\fontseries{b}\selectfont MSTest}, as autoregressive models are particularly relevant in many economic and financial applications.
quzhuo2021likelihood presents a novel characterization of the conditional regime probabilities for a family of likelihood ratio-based tests and establishes the asymptotic distribution for the test statistics. Like hansen92, garcia98, and chowhite07, they derive an approximation of the likelihood ratio as an empirical process, where $\{(p_{11},p_{22}): \epsilon \leq p_{11},p_{22} \leq 1-\epsilon \ \& \ p_{11}+p_{22}\geq 1+\epsilon\}$ and $\epsilon$ is a small constant. They also provide a finite sample refinement to correct some of the over-rejections that can occur in specific cases. As a result, they are able to study the asymptotic null distribution and find that, even though the null hypothesis has only one regime and thus parameters do not switch, the nuisance parameters can affect the limiting distribution, which will depend on which parameters are allowed to switch. Furthermore, quzhuo2021likelihood describe how some of these results can explain why specific bootstrap procedures may be inconsistent (e.g., when including weakly exogenous regressors) and why standard information criteria such as the BIC can be sensitive to the hypothesis and the model structure. Although the \proglang{R} package {\fontseries{m}\fontseries{b}\selectfont MSTest} currently does not include the test procedure proposed by quzhuo2021likelihood, future versions will likely incorporate this procedure, as it represents a significant contribution and provides the best approximation to the asymptotic distribution of the likelihood ratio test statistic, making it useful for users interested in these asymptotic results.
Another likelihood ratio-based test is that of kasshi2018. In kasshi2015, the authors introduce a re-parameterization and higher-order expansion of the likelihood ratio function. In kasshi2018, they apply this technique along with the Difference in Quadratic Mean (DQM) approximation introduced by liushao2003 to address some of the issues discussed above while estimating and studying the asymptotic distribution of the likelihood ratio test statistic for Markov switching models. Additionally, they do this in a more general setting where $M_0 \geq 1$ but $m=1$ still, allowing them to consider a null hypothesis with more than one regime. This makes their approach more general than the likelihood ratio test procedures discussed so far in that regard. In doing so, like quzhuo2021likelihood, they also show that the parametric bootstrap procedure is an asymptotically valid test procedure, but in this case, they only consider the scenario where fixed or predetermined regressors are included in the model while examining Markov switching autoregressive models. The parametric bootstrap procedure is not directly available in {\fontseries{m}\fontseries{b}\selectfont MSTest}, but it can be implemented using specific options within the Local Monte Carlo likelihood ratio test (LMC-LRT) proposed by rodrondufour_mcmstest, which is included. This process is further detailed in section (ref), where we discuss the test procedures available within the {\fontseries{m}\fontseries{b}\selectfont MSTest} package. Notably, the Monte Carlo procedures of rodrondufour_mcmstest are applicable in even more settings (e.g., $m>1$ and more) and so they are discussed next.
In rodrondufour_mcmstest, the authors propose the Maximized Monte Carlo likelihood ratio test (MMC-LRT) and the Local Monte Carlo likelihood ratio test (LMC-LRT), which can be used to compare very general Markov switching models. These procedures represent the most general type of testing methods, as they can be applied to hypothesis testing in its broadest form when both $M_{0}$ and $m$ are greater than or equal to 1. Furthermore, as described in rodrondufour_mcmstest, these Monte Carlo likelihood ratio tests can be utilized when the process $y_{t}$ is non-stationary, when the model exhibits non-Gaussian errors, when parameters take values at the boundary, and in multivariate settings, which include the previously discussed Markov switching VAR and multivariate Hidden Markov models. To be more precise, the MMC-LRT and LMC-LRT are the only test procedures available in the literature and in the \proglang{R} package {\fontseries{m}\fontseries{b}\selectfont MSTest} that can be employed to test multivariate Markov switching models. Additionally, the MMC-LRT procedure is valid in finite samples and is robust to identification issues, features that are empirically relevant when considering macroeconomic applications of Markov switching models, as discussed in rodrondufour_mcmstest.
Here, we summarize the MMC-LRT and LMC-LRT procedures but readers interested in further details are referred to the more formal description provided in rodrondufour_mcmstest and dufour2006 for even further details on the Monte Carlo techniques used in these procedures. For simplicity of exposition, we use an example where we are interested in a null hypothesis of a linear model (i.e., only $M_{0}=1$ regime) and an alternative hypothesis of $M_{0}+m=2$ regimes and consider an autoregressive model where only the mean and variance are subject to change. Since this is a likelihood ratio-based approach, the log-likelihood values under the null and alternative hypothesis are required. The log-likelihood for the model under the alternative (and under the null hypothesis if $M_{0}>1$) is given by ((ref)) - ((ref)):
where
Here, the subscript of $1$ underscores the fact that $\theta _{1}$ is the parameter vector under the alternative hypothesis. The set $\Omega$ satisfies any theoretical restrictions we may wish to impose on $\theta _{1}$ [such as $\sigma _{1}>0$ and $\sigma _{2}>0$]. On the other hand, the log-likelihood under the null hypothesis ($M_{0}=1$) is given by
where
Note that $\bar{\Omega}_{0}$ has lower dimension than $\Omega $. The null and alternative hypotheses can be written as:
where $\delta _{1}=(\mu _{1}$, $\sigma _{1})$ and\ $\delta _{2}=(\mu _{2}$, $\sigma _{2})$. Clearly, $H_{0}$ is a restricted version of $H_{1}$: for each $\theta _{0}\in \bar{\Omega}_{0}$, we can find $\theta _{1}$ such that
where $\Omega _{0}$ is the subset of vectors $\theta _{1}\in \Omega $ such that $\theta _{1}$ satisfies $H_{0}$. Under $H_{0}$, the vector $\theta_{0}\in \bar{\Omega}_{0}$ is a nuisance parameter: the null distribution of any test statistic for $H_{0}$ depends on $\theta _{0}\in \bar{\Omega}_{0}$. In this problem, the null distribution of the test statistic, $LR_{T}$, is in fact completely determined by $\theta _{0}\in \bar{\Omega}_{0}$. As in garcia98 and the parametric bootstrap procedure describe in quzhuo2021likelihood and kasshi2018, it is assumed that the null hypothesis depends only on the mean, variance, and autoregressive coefficients. The likelihood ratio statistic for testing $H_{0}$ against $H_{1}$ can then written as
where
Since the model is parametric, we can generate a vector $N\;$i.i.d replications of $LR_{T}$ for any given value of $\theta _{0}\in \bar{\Omega}_{0}$:
As discussed in rodrondufour_mcmstest, the main assumptions required are that the random variables $LR_{T}^{(0)}, \,LR_{T}^{(1)}(\theta _{0}),\,\ldots \,,LR_{T}^{(N)}(\theta _{0})$ are exchangeable for some $\theta _{0}\in \bar{\Omega}_{0}$ each with distribution function $F[x\,|\,\theta _{0}]$ (i.e., they are $i.i.d$). From here, we can compute the Monte Carlo $p$-value which is given by
where
and $I(C):=1$ if condition $C$ holds, and $I(C)=0$ otherwise. As can be seen from ((ref)), $R_{LR}[LR_{T}^{(0)};$ \thinspace $N]$ simply computes the rank of the test statistic from the observed data within the generated series $LR(N,$\thinspace $\theta _{0})$. Then, as shown in rodrondufour_mcmstest, a critical region for this test statistic with level $\alpha $ is then given by
More precisely, if $(N+1)\alpha $ is an integer, then
under the null hypothesis and so it is a valid test with level $\alpha $ for $H_{0}$ and this result does not depend on the sample size $T$ and so it is also valid in finite samples.
This is the Maximized Monte Carlo likelihood ratio test and it requires searching for the maximum Monte Carlo p-value over the nuisance parameter space $\bar{\Omega}_{0}$. Since this space can be very large and grows as the number of autoregressive components and the number of regimes increases, the authors propose a more efficient alternative which involves searching over a consistent set $C_T$ [as originally proposed in dufour2006]. A consistent set can be defined using the consistent point estimate. For example, let $\hat{\theta}_{0}$ be the consistent point estimate of $\theta_{0}$. Then, we can define
where $c$ is a fixed positive constant that does not depend on $T$ and $\left\Vert \cdot \right\Vert $ is the Euclidean norm in $\mathbb{R}^{k}$. A consistent set of interest may also be $C_{T}^{\ast }=C_{T}^{CI}\cup C_{T}^{\epsilon }$ where
Hence, $C_{T}^{CI}$ is defined by a confidence interval based on consistent point estimates, while $C_{T}^{\epsilon }$ is determined using a fixed constant $\epsilon$ that is independent of $T$. The union of these two sets allows for values that may lie outside the confidence interval for some parameters and within it for others, depending on the choice of $\epsilon$. The {\fontseries{m}\fontseries{b}\selectfont MSTest} \proglang{R} package enables users to define the consistent set $C_T$ by specifying only a fixed positive constant $\epsilon$, only the confidence interval, or the union of both.
As discussed in dufour2006 and rodrondufour_mcmstest, the solution to this optimization problem may not be unique, meaning that the maximum $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 within the nuisance parameter space. The package {\fontseries{m}\fontseries{b}\selectfont MSTest} allows users to use the Generalized Simulated Annealing algorithm, Genetic Algorithm, and Particle Swarm algorithm [see xiaetal2013, zametal2013, scrucca2013, dufour2006, and dufour2019finite].
Finally, as described in rodrondufour_mcmstest, we can define $C_{T}$ as the singleton set $C_{T}={\hat{\theta}_{0}}$, resulting in the Local Monte Carlo Likelihood Ratio Test (LMC-LRT). Here, the consistent set includes only the consistent point estimate $\hat{\theta}_{0}$, so the Monte Carlo $p$-value depends solely on $\hat{\theta}_{0}$. This LMC version of the test can be interpreted as the finite-sample analogue of the parametric bootstrap. In this context, asymptotic validity pertains to $\hat{\theta}_{0}$ converging to the true parameter $\theta _{0}$ as the sample size increases, rather than to the asymptotic validity of critical values emphasized in studies like hansen92, garcia98, chowhite07, quzhuo2021likelihood, and kasshi2018. Specifically, akin to the parametric bootstrap, the LMC procedure is valid only as $T\rightarrow \infty$. However, unlike the parametric bootstrap, a large number of simulations (i.e., $N\rightarrow \infty$) is not necessary, as the procedure does not approximate asymptotic critical values or assume asymptotic convergence of the test statistic distribution but instead uses critical values derived from the sample distribution. This design allows for computational efficiency since it eliminates the need for extensive simulations to obtain asymptotically valid critical values. Further, as discussed in rodrondufour_mcmstest, The procedure is valid even in cases where an asymptotic distribution does not exist. This is another feature which makes both the MMC-LRT and the LMC-LRT procedures more general than the parametric bootstrap. Still, it is worth noting that the parametric bootstrap procedure discussed in quzhuo2021likelihood and kasshi2018 can be implemented by constraining the parameter space of the transition probabilities based on the assumptions in those studies. Using a larger number of simulations in these contexts can approximate asymptotic critical values where prior research has demonstrated the parametric bootstrap’s validity. Nevertheless, as outlined in rodrondufour_mcmstest, the LMC-LRT and MMC-LRT procedures remain the most general likelihood ratio test procedures currently available.
dufourluger17 propose a different way to test Markov switching models that also avoids the statistical issues described above for likelihood ratio type tests. Additionally, their test is less costly computationally in comparison to all test mentioned above, the parameter stability test discussed next, and allows the econometrician to perfectly control the level of the test through the use of the Monte Carlo test methods described in dufour2006 and used in rodrondufour_mcmstest. However, their proposed method can only deal with the case where $M_0=m=1$. The moment-based test of dufourluger17 involves computing moments of the least-square residuals of autoregressive models under the null hypothesis. More Specifically, they focus on the mean, variance, skewness and excess kurtosis of the least-square residuals. These moments are calculate as
where, \(m_1 = \frac{\Sigma^T_{t=1} \hat{\epsilon_t} \mathbbm{1}[\hat{\epsilon_t} < 0]}{\Sigma^T_{t=1}\mathbbm{1}[\hat{\epsilon_t} < 0]}\), \(m_2 = \frac{\Sigma^T_{t=1} \hat{\epsilon_t} \mathbbm{1}[\hat{\epsilon_t} > 0]}{\Sigma^T_{t=1}\mathbbm{1}[\hat{\epsilon_t} > 0]}\), \(s^2_1 = \frac{\Sigma^T_{t=1} (\hat{\epsilon_t}-m_1)^2 \mathbbm{1}[\hat{\epsilon_t} < 0]}{\Sigma^T_{t=1}\mathbbm{1}[\hat{\epsilon_t} < 0]}\) and \(s^2_2 = \frac{\Sigma^T_{t=1} (\hat{\epsilon_t}-m_2)^2 \mathbbm{1}[\hat{\epsilon_t} > 0]}{\Sigma^T_{t=1}\mathbbm{1}[\hat{\epsilon_t} > 0]}\)
where, \(\vartheta_1 = \frac{\Sigma^T_{t=1} \hat{\epsilon_t}^2 \mathbbm{1}[\hat{\epsilon_t}^2 < \hat{\sigma}^2]}{\Sigma^T_{t=1} \mathbbm{1}[\hat{\epsilon_t}^2 < \hat{\sigma}^2]} \), \(\vartheta_2 = \frac{\Sigma^T_{t=1} \hat{\epsilon_t}^2 \mathbbm{1}[\hat{\epsilon_t}^2 > \hat{\sigma}^2]}{\Sigma^T_{t=1} \mathbbm{1}[\hat{\epsilon_t}^2 > \hat{\sigma}^2]} \) and \(\hat{\sigma}^2=T^{-1}\Sigma^T_{t=1} \hat{\epsilon_t}^2\)
and
The testing procedure involves calculating the test statistic for each moment, obtaining the individual p-values and using two different methods of combining independent test statistics. The first method is based on the min of the p-values and was suggested by tippett1931 and wilkinson1951. Here, the test statistic becomes,
where for example, \(\hat{G_M}[M(\hat{\epsilon})] = 1 - \hat{F_M}[M(\hat{\epsilon})]\) is the Monte Carlo p-value of \(M(\hat{\epsilon})\). The second method of combining the test statistics involves taking the product of them. This method of combining test statistics was suggested by fisher1932 and pearson1933. In this case the the test statistic becomes,
Interested readers should see dufour2004 and dufour2014 which provide further discussion of these methods of combining test statistics. Finally, the Monte Carlo p-value of the combined statistics is given by
and
where \(R_{F_{min}}\) and \(R_{F_{\times}}\) are the ranks of \(F_{min}(\hat{\epsilon})\) and \(F_{\times}(\hat{\epsilon})\) in \(F_{min}(\hat{\eta}_1)\), ..., \(F_{min}(\hat{\eta}_{N-1})\) and \(F_{\times}(\hat{\eta}_1),\) ...,\( F_{\times}(\hat{\eta}_{N-1})\) respectively, when ordered. Also, \(\hat{\eta} = \eta - \bar{\eta}\) and \(\eta \sim N(0,I_T)\).
The computational efficiency of this test makes it easily extendable to the use of Maximized Monte Carlo when nuisance parameters are present. Furthermore, this test is not subject to the same level of statistical difficulties such as unidentified parameters under the null as in hansen92, garcia98 and chp14. This is because, transition probability parameters, the mean and the variance do not need to be treated as nuisance parameters. Parameters of explanatory variables are the only ones which may be unidentified under the null and so only these are treated as nuisance parameters. garcia98 also reduced the nuisance parameter space by treating only the transition probabilities $p_{11}$ and $q_{22}$ as nuisance parameters, however, this is an even further reduction of the nuisance parameter space and makes the treatment of autoregressive models with more lags more tractable. Although the moment-based test can only be used to compare linear models against Markov switching models with two regimes, it is the least computationally intensive procedure available and takes only seconds to compute, even when considering the MMC version of the test, and so it is included in the \proglang{R} package {\fontseries{m}\fontseries{b}\selectfont MSTest}.
chp14 proposes a test that can be described as an optimal test for the consistency of parameters in random coefficient and Markov switching models. This test can be understood as an extension of white1982's information matrix test. As suggested by the authors, it shares several advantages, such as the need to estimate the model only under the null hypothesis, which, as we saw, is also the case for the moment-based test of dufourluger17. This feature of needing to estimate only the restricted model is particularly advantageous. In contrast, the likelihood ratio tests proposed by hansen92, garcia98, chowhite07, quzhuo2021likelihood, kasshi2018, and rodrondufour_mcmstest all require estimating the model under both the null and alternative hypotheses. The presence of non-linearity in estimating Markov switching models introduces multiple local optima, necessitating numerical procedures and making likelihood ratio test procedures relatively more computationally intensive. Moreover, the authors invoke the Neyman-Pearson lemma to prove the optimality of their test, demonstrating that it is asymptotically locally equivalent to the likelihood ratio test. However, simulation results presented in rodrondufour_mcmstest and discussion in quzhuo2021likelihood suggest their LRT approach have better power in certain settings, such as when only the mean is subject to change and persistence is high. Further, this method also involves a bootstrap procedure while searching over the nuisance parameter space, which can make obtaining asymptotic critical values more computationally intensive than the moment-based approach of dufourluger17, for example.
The authors formulate the hypothesis in the following way:
where the switching variable \(\eta_t\) is unobservable, stationary and may depend on nuisance parameters \(\beta\). Their test makes use of the second derivatives of the log-likelihood and the outer products of the scores as in the information matrix test with the addition of an extra term, which captures the serial dependence of the time-varying coefficients. This means that the form of the test depends on the latent process \(\eta_t\) only through its second-order properties. Additionally, the distribution of \( \eta_t \) is assumed to exist even under the null, but does not play a role with regards to the distribution of the data \((y_T, y_{T-1},y_{T-2},...y_1)\) under the null. That is, under the null, they are mutually exclusive.
The authors first propose a Sup-type test as in davies87 to combat the presence of nuisance parameters. They set \(\eta_t = chS_t \), where $c$ is a scalar specifying the amplitude of the change, $h$ a vector specifying the direction of the alternative and \(S_t\) is a Markov-chain, which follows an autoregressive process such as \(S_t = \rho S_{t-1} + e_t\), where \(e_t\) is i.i.d. U[-1,1] and \(-1 <\rho < 1\) so that \(S_t\) is bounded by support \(( -1/(1-|\rho|), 1/(1-|\rho|)\) and has zero mean. Letting \(\beta =(c^2, h',\rho ') \) be the vector of nuisance parameters, we can write
which allows us to get the expression
as in chp14, where \(\mu^*_{2,t}(\beta,\theta) = \mu_{2,t}(\beta,\theta)/c^2\), \(\Gamma^*_T = \Gamma^*_T(\beta,\theta) = \sum\limits_{t} \mu^*_{2,t}(\beta,\theta)/\sqrt{T}\) and \(\hat{\epsilon^*}\) are the residuals from regressing \(\mu^*_{2,t}(\beta,\theta)\) on \( l_t^{(1)}(\hat{\theta})\) so that \(\Gamma^*_T \) and \(\hat{\epsilon^*}\) are not dependent on \(c^2\). As previously mentioned, this methodology involves bootstrapping over the distributions of the nuisance parameters. As a result, one must choose a prior distribution for the nuisance parameters. The most commonly used distribution in this case is the uniform distribution, and this is what is implemented in the {\fontseries{m}\fontseries{b}\selectfont MSTest} package. However, since the parameter \(c^2\) is not necessarily bounded from above, a uniform distribution may not always be appropriate. As a result, chp14 also suggest using an Exponential-type test as in andrewsploberger94. They propose the following statistic:
where
These tests proposed by chp14 have been used in the empirical applications of hamilton05, warnevredin06, kahnrich07, morleypiger12, and DMPF11, in testing MS-GARCH models by hushin08, and in dufourluger17, quzhuo2021likelihood, and rodrondufour_mcmstest as a benchmark to compare their test procedures. Due to their wide use and optimality results, the {\fontseries{m}\fontseries{b}\selectfont MSTest} package also includes these test procedures.
The MSTest package is designed for conducting hypothesis testing for the number of regimes in Markov switching models within the \proglang{R} programming environment. Since many of the test procedures included here require estimating the restricted and unrestricted models and simulating the restricted model, it also makes available simulation and estimation of Markov switching models. However, it should be noted that these features are simply a by-product and not the focus of this package. Specifically, they were designed with the purpose of being compatible with the hypothesis testing functions.
Several other platforms also offer functionalities for Markov switching models. for example, in \proglang{MATLAB}, the Econometrics Toolbox provides functions for estimating state-space models, including Kalman filters, though it may lack hypothesis testing for these models. In \proglang{Python}, the {\fontseries{m}\fontseries{b}\selectfont statsmodels} library also offers support for estimating state-space models, including Markov switching models, but again does not include dedicated functions for hypothesis testing of Markov switching model. Similarly, in \proglang{STATA}, while commands like \code{mswitch} exist for estimating Markov switching models, they do not cover the extensive range of testing procedures that {\fontseries{m}\fontseries{b}\selectfont MSTest} provides.
For these reasons, {\fontseries{m}\fontseries{b}\selectfont MSTest} stands out in its implementation of novel testing procedures that are robust to the violation of regularity condition. Additionally, the integration of Rcpp enhances computational efficiency, allowing users to handle nuisance parameters effectively. Hence, while there are alternative tools available across platforms for simulating and estimating Markov switching models, {\fontseries{m}\fontseries{b}\selectfont MSTest} offers a unique combination of flexibility, computational efficiency, and availability of hypothesis testing procedures that sets it apart.
The R package {\fontseries{m}\fontseries{b}\selectfont MSTest} includes three samples of U.S. real GNP and one sample of U.S. real GDP, all of which can easily be accessed once the package is loaded. Specifically, it provides the original sample used in hamilton89, the sample ending in 2010 first considered in chp14, and a more complete sample ranging from the second quarter of 1947 to the second quarter of 2024. The U.S. real GDP data also covers this extended period. These samples have been used by hansen92, chp14, dufourluger17, and rodrondufour_mcmstest, among others, to test for the number of regimes in a Markov switching model and to showcase the performance of their proposed test procedures. Table (ref) provides the label used to identify each sample in the \proglang{R} package {\fontseries{m}\fontseries{b}\selectfont MSTest} and describes the span of each sample.
These data sets can be loaded using the following commands in \proglang{R} once {\fontseries{m}\fontseries{b}\selectfont MSTest} has been loaded:
All three data sets include three columns: (1) \code{Date}, (2) \code{GNP} or \code{RGDP}, and (3) \code{GNP_gr} or \code{RGDP_gr}. The first column is of \code{Date} type, defined using the \code{as.Date()} function. The second column contains the levels of U.S. real GNP or GDP, and the third column contains the growth rate of U.S. real GNP or GDP.
This section describes a set of functions available in {\fontseries{m}\fontseries{b}\selectfont MSTest} that allow users to simulate different types of processes. These functions are utilized by the hypothesis testing procedures, specifically those that involve using simulation to obtain the (sample or asymptotic) null distribution of the test statistic, such as \code{LMCLRTest}, \code{MMCLRTest}, \code{DLMCTest}, \code{DLMMCTest}, and \code{CHPTest}. Users interested in developing new estimation or testing procedures for Markov switching models may find these functions useful for testing the performance of their procedures and comparing them to those available in {\fontseries{m}\fontseries{b}\selectfont MSTest}.
Table (ref) provides the labels used to identify each simulation function in the {\fontseries{m}\fontseries{b}\selectfont MSTest} package and describes the processes they can simulate. Each function requires a \code{List} object as input, containing values for the data-generating process. Below, we present examples for some of the included simulation functions, which should give users an idea of the elements the input \code{List} must include. As with all functions in {\fontseries{m}\fontseries{b}\selectfont MSTest}, a more exhaustive description of each function's inputs is provided in the {\fontseries{m}\fontseries{b}\selectfont MSTest} documentation available through CRAN.
For example, as described in Table (ref), simulating a normally distributed process can be done using the \code{simuNorm} function. In this case, the function requires the user to input a \code{List} containing: the sample size (\code{n}), the number of series to be simulated (\code{q})—where \code{q=1} indicates a univariate process and \code{q>1} indicates a multivariate process—a ($q \times 1$) vector of means for each series, and a ($q \times q$) covariance matrix. Users may also specify the number of additional observations to simulate and discard at the beginning using the \code{burnin} parameter, which for \code{simuNorm} defaults to $0$. For autoregressive type processes however, the default \code{burnin} value is higher to avoid any dependence on the initialization used in simulating the process. Users also have the option to provide a ($(n+\text{burnin}) \times q$) matrix of errors if they prefer not to use normally distributed errors. This is done by defining the element \code{eps} in the input \code{list} with that matrix. As an example, we can simulate a multivariate normal process, an autoregressive process, a vector autoregressive process, a Markov switching autoregressive process, a hidden Markov process, and a Markov switching vector autoregressive process using the following code:
These simulated processes are shown in Figure (ref). From here, we can already observe the regime switching nature of some of the processes. It is especially more obvious when comparing the two middle charts which plot an autoregressive process (left) and Markov switching autoregressive process (right).
Here, we briefly describe the set of functions available in {\fontseries{m}\fontseries{b}\selectfont MSTest} that allow users to estimate different types of Markov switching models. In total, there are ten distinct types of models that can be estimated. These models are listed in Table (ref), along with their labels in {\fontseries{m}\fontseries{b}\selectfont MSTest} and the corresponding equations for each model. Although \code{Nmdl} and \code{HMmdl} use multivariate notation, if a ($T \times 1$) vector is provided as input, the function will detect that this is a univariate setting. Also, as may be apparent from their equations, the labels that include an `X' in their names are those that allow for the inclusion of exogenous regressors. Although the \code{Nmdl} and \code{HMmdl} functions do not follow this nomenclature, including exogenous regressors within these functions is always possible.
Above, we used the simulation functions to simulate a multivariate Hidden Markov process with $q=2$ series, a Markov switching autoregressive process with $p=1$ lag, and a Markov switching vector autoregressive process with $p=1$ lag and $q=2$ series. The output of these functions is a list that includes the simulated process, the true state variable $S_{t}$, and other elements of the data-generating process (DGP) that were provided as input. Below, we provide an example where the simulated processes are used as input for estimating the respective model. This is demonstrated for the Hidden Markov model and the two Markov switching models.
In each case, we set \code{method="EM"} in the \code{control} \code{List} object, which is used to specify options when estimating these models, to utilize the expectation maximization algorithm for optimization. As previously mentioned and described in the package documentation, the user can also set \code{method="MLE"} if they wish to employ maximum likelihood estimation. Additionally, we can specify whether we want the mean and variance to change according to the regime by setting \code{msmu=TRUE} and \code{msvar=TRUE}, respectively. Setting either of these to false would result in a model where that parameter is constant across regimes. The option \code{use_diff_init=30} is used to estimate the model thirty times with different initial values each time. The model with the highest log-likelihood value is kept as the main output, but the output \code{List} of these functions provides results for all iterations under the output \code{List} \code{trace}. The output \code{List} also contains various elements, all of which are described in the package documentation. Importantly, the outputs are \code{S3} objects, and the \code{print()} and \code{summary()} methods are provided for each, as can be seen in the example code above. Specifically, we see that the \code{summary()} method can be used to display the parameter estimates, which in this controlled setting can be compared to the true values, along with other characteristics such as the log-likelihood, AIC, BIC, and quantiles of the residuals.
The Figure (ref) plots the simulated processes in black and blue, the true regime states $S_{t}$, which were also an output of the simulation functions, in red (dashed), and the model-estimated smoothed probabilities in green (dashed). From this, we can see that the detection of changes in regime is captured quite well in the estimation process. For some periods, there is, relatively speaking, a bit more difficulty with the Markov switching VAR model, but increasing sample sizes or considering more initial values could likely help in this regard.
In this section, we describe the hypothesis testing functions available in the {\fontseries{m}\fontseries{b}\selectfont MSTest} package. The {\fontseries{m}\fontseries{b}\selectfont MSTest} package has been designed with ease of use in mind and hence, in most cases, only requires the user to provide the variable $y_t$ (and $Z_t$ when relevant) and specify the number of regimes to test for.
Table (ref) shows all the functions currently available in the {\fontseries{m}\fontseries{b}\selectfont MSTest} package and briefly describes them. All functions return the test statistic, critical values of the test statistic distribution when available, the p-value, and parameter estimates under the null and alternative hypothesis depending on which were used for testing. It is important to note that, even though they are called critical values, the values returned by \code{HLRTest()} are the critical values of the process $Q$ discussed in hansen92, which provide a bound for the LR but are not the critical values of the LR test for Markov switching. Likewise, the values returned by \code{LMCLRTest()}, \code{MMCLRTest()}, \code{DLMCTest()}, and \code{DLMMCTest()}, described as critical values, are in fact the percentiles of the simulated null distribution and hence we are using the term critical values loosely here.
We begin by discussing the Local Monte Carlo Likelihood Ratio test. As indicated in Table (ref), this test can be implemented in the {\fontseries{m}\fontseries{b}\selectfont MSTest} package using the function \code{LMCLRTest()}. Since this test requires estimating both the restricted and unrestricted models to simulate the null distribution, the options we would typically pass into the estimation functions using a \code{List} can be passed into the options of this testing procedure using \code{mdl_h0_control} and \code{mdl_h1_control}. The same applies to the Maximized Monte Carlo Likelihood Ratio test, which we describe next. Specifically, we can set options such as whether the mean or variance switches according to the regime, the number of initial values to use when estimating the models, and other related settings. However, the \code{LMCLRTest()} function allows users to specify a different number of initial values for estimating models during null distribution simulation through the \code{use_diff_init_sim} option. By default, this is set to match the value used to estimate the model with observed data, as specified in the \code{mdl_h0_control} and \code{mdl_h1_control} \code{List}s, which is generally recommended.
Apart from the options that can be set, the \code{LMCLRTest()} function also requires specifying the number of lags, \code{p}, the number of regimes under the null hypothesis, \code{k0}, and the number of regimes under the alternative hypothesis, \code{k1}, to be compared. Below, we provide an example where we apply the LMC-LRT procedure to the Markov switching autoregressive model with $p=1$ lag, that we simulated above. In this example, we test the null hypothesis $H_0: M=1$ against the alternative $H_1: M=2$. Given that the simulated data is a Markov switching model with $M=2$ regimes, we would expect the null hypothesis to be rejected, which is indeed the result obtained.
As previously mentioned, the Local Monte Carlo Likelihood Ratio test function, \code{LMCLRTest()}, can be used to replicate the parametric bootstrap approach discussed in quzhuo2021likelihood and kasshi2018. To do this, we start by setting \code{mdl_h1_control = list(method = "MLE")}. In those studies, MLE is used, and using MLE in {\fontseries{m}\fontseries{b}\selectfont MSTest} also allows us to constrain the parameter space. To set parameter constraints, we define the vectors \code{mle_theta_low} and \code{mle_theta_upp} in the \code{mdl_h1_control} \code{List}. In addition to constraining the parameter space for the transition probabilities, kasshi2018 also impose restrictions on the variance, which can be implemented here. Additionally, we may wish to increase the value of \code{lmc_control = list(N)} to improve the approximation of the asymptotic critical values. For example, quzhuo2021likelihood set $N=199$, while kasshi2018 set $N=299$ for similar purposes. Ideally, for a bootstrap procedure, this value should be higher. However, given the computational demands of estimating Markov switching models, using lower values can be reasonable.
The Maximized Monte Carlo Likelihood Ratio test can be conducted using the \code{MMCLRTest()} function. This function shares the same options as the LLMC-LRT procedure for estimating both the restricted and unrestricted models. Similarly, users must specify the number of lags, \code{p}, the number of regimes under the null hypothesis, \code{k0}, and the number of regimes under the alternative hypothesis, \code{k1}, for comparison. However, \code{MMCLRTest()} also offers additional settings related to the testing procedure, as this method is more intricate. For instance, users can set the fixed constant $c$, used to define the consistent set over which to search, by setting \code{eps}. To use the set $C_{T}^{\ast} = C_{T}^{CI} \cup C_{T}^{\epsilon}$, as described earlier, \code{CI_union=TRUE} should be specified. To use only $C_{T}^{CI}$, one can set \code{eps=0} and \code{CI_union=TRUE}. Additionally, users can select the numerical optimization algorithm for this search via \code{type}. By default, the algorithm stops when a p-value of $1$ is reached, as this is the highest possible p-value. Alternatively, users can choose to stop when the test fails to reject, that is, upon reaching any value above the test level $\alpha$, by using \code{threshold_stop}.
Below, we provide an example where the MMC-LRT procedure is used to test a linear autoregressive model with only $M=1$ regime. Here, we set \code{eps=0.3} and \code{CI_union=TRUE} to use the consistent set $C_{T}^{\ast }$ and set \code{type="pso"} to employ the particle swarm algorithm (see psopack). Other available optimization algorithms include: Simulated Annealing using the “GenSA” package introduced by xiaetal2013, Particle Swarm algorithm using the “PSO” package introduced byzametal2013 and Genetic algorithm using the “GA” introduced by scrucca2013. Importantly, we also set \code{workers=8}, which allows us to use a parallel version of the test to improve computational efficiency. This option is also available for the \code{LMCLRTest()} function. Specifically, the null distribution is simulated using $8$ different workers in this case. Note that, as shown in the example, the user must register a parallel pool and close the cluster in order to make use of this functionality.
As expected, we fail to reject the null hypothesis of no Markov switching (i.e., $M=1$) in this case since the true DGP is one of a linear autoregressive process. Given that we set \code{threshold_stop = 0.05 + 1e-6}, the algorithm converged quickly. We could continue searching for the parameter values under the null that give the maximum p-value, but for the purpose of exposition, we set the \code{threshold} parameter so that we stop searching once the test fails to reject the null hypothesis.
The Monte Carlo moment-based test proposed by dufourluger17 is invoked using \code{DLMCTest()} (local version) and \code{DLMMCTest()} (maximized version). Like all other functions, these require $y_{t}$ and $p$, the number of autoregressive lags, to be specified by the user. The parameter $N$ sets the number of Monte Carlo replications. The default value for $N$ is $99$, as dufour2004 shows that having more than $100$ replications (i.e., including the value obtained from the observed series) has a minimal effect on power. This test also involves approximating the distribution of the p-value for each statistic through a separate round of simulation. The number of replications used in this approximation is determined by $N2$, which is set to $10,000$ by default. For the moment-based local Monte Carlo test, these are the only required and available options.
Below, we provide an example of this test using the previously simulated Markov switching autoregressive data. Results are obtained very quickly. The output of the \code{summary()} function presents the nuisance parameter value (in this case, $\phi_1$), the test statistic for each of the four moments, the combined test statistic under $F(e)$, the critical values, and the Monte Carlo p-values. Here, we observe a clear rejection of the null hypothesis, which is expected since the process used is indeed a Markov switching process with $M=2$ regimes.
The computational efficiency of this procedure extends well to the Moment-based maximized Monte Carlo test proposed by dufourluger17, which again can be performed very quickly using the {\fontseries{m}\fontseries{b}\selectfont MSTest} package. Here, more options are available regarding the optimization over the nuisance parameter space. Here, \code{optim_type} is used to determine the numerical optimization algorithm to be used. As with the MMC-LRT procedure, we can use \code{eps} and \code{CI_union} to define the consistent set over which to maximize the p-value.
The computational efficiency of this procedure extends well to the moment-based maximized Monte Carlo test proposed by dufourluger17, which can also be performed very quickly using the {\fontseries{m}\fontseries{b}\selectfont MSTest} package, despite the maximization over the nuisance parameter space. In this case, like with the MMC-LRT procedure, more options are available regarding the optimization over this space. Specifically, \code{optim_type} is used to determine the numerical optimization algorithm to be employed. Similar to the MMC-LRT procedure, we can utilize \code{eps} and \code{CI_union} to define the consistent set over which to maximize the p-value. In fact, many of the same options available for the \code{MMCLRTest()} and available here also.
Above, we provide an example of this function, where we again utilize the Markov switching autoregressive process that was previously simulated. We set the \code{threshold} parameter to be $0.05 + 1e-6$ so that we stop searching if the test fails to reject the null hypothesis, which should not be the case. This time, we also set \code{eps=0} so that the search is performed over the confidence interval, as in dufourluger17. Here, we reject the null hypothesis of a linear model once more, a result that is consistent with the tests used so far.
The test proposed by chp14 can be employed by using CHPTest(), where again $y_{t}$ and $p$ must be provided by the user. In this case, $N$ is the number of bootstraps and is set to $3000$ by default, as in their work, and $\rho$ is used to set the bound for one of the nuisance parameters. This value is set to $0.7$ by default, as in their work also. By setting $\rho = 0.7$, the grid search occurs over the parameter space $\rho \in [-0.7, 0.7]$.
Above, we provide an example of this test using the linear autoregressive process simulated previously again. In this case, we also set \code{msvar=TRUE} to consider a test where changes in the variance are also considered. The \code{summary()} function provides estimates of the restricted model, as this is the only model needed to perform the test, and displays results from the supTS and expTS versions of the test. As expected, we also find that the test fails to reject the null hypothesis of a linear model.
Finally, in the {\fontseries{m}\fontseries{b}\selectfont MSTest} package, the test proposed by hansen92 can be invoked using the command HLRTest(). As before, this function requires the user to provide $y_{t}$, the variable of interest. This test is designed to assess the null hypothesis of linearity in autoregressive models, so a value for the number of autoregressive components $p$ must again be provided. The test performs a grid search over the nuisance parameters, meaning most of the available options are related to configuring this grid. Specifically, users can set the \code{gridsize}, as well as the starting point and step size for both the mean and variance grids. These values pertain to the unrestricted model. While keeping these values fixed, the test procedure optimizes the values of the restricted model; thus, the user can utilize \code{theta_null_low} and \code{theta_null_upp} to define the optimization space. It is important to note that while setting \code{msvar = TRUE} is possible, it will require more computation time as the function must optimize over a larger parameter space.
Above, we provide an example of this test using the previously simulated Markov switching autoregressive process. In this case, we also set \code{msvar=TRUE} to consider a test where changes in variance are enabled. The \code{summary()} function provides estimates of the restricted model, as this is the only model needed to perform the test in this instance. This output also displays results from the different types of standardizations mentioned in hansen96. In all cases, we find that the test rejects the null hypothesis of a linear model, yielding a p-value that is very small, approximately $0$.
In this section, we apply the full suite of tests available in {\fontseries{m}\fontseries{b}\selectfont MSTest} to three empirical data sets provided in the package, specifically to three samples of U.S. GNP growth rates covering distinct periods of interest. According to kimnel1999, a model with two regimes should capture the structural decline in business cycle volatility that began in the mid-1980s, a phenomenon known as the Great Moderation. Additionally, as discussed by hamilton89, the presence of business cycle fluctuations support a two-regime model to to reflect shifts in the growth rate of U.S. GNP. These characteristics make U.S. GNP growth an ideal empirical example for the tests provided by {\fontseries{m}\fontseries{b}\selectfont MSTest}. If the model accurately captures these structural shifts, by allowing both the mean and variance to vary by regime we should expect the tests to at least reject the null hypothesis of one regime (i.e., a linear model) for the second and third sample, which covers both the high and low-volatility periods and various recessionary and expansionary periods.
We conduct the tests using two model specifications. In the first, only the mean changes by regime, as initially suggested by hamilton89. IN the second, and an alternative specification where both the mean and variance vary by regime is considered. Table (ref) presents the test results for all three samples across these two specifications. The first panel shows results when only the mean changes by regime. Here, we observe that the tests fail to reject the null hypothesis of a linear model for the first sample but reject it for the third, extended sample. For the second sample, covering 1951Q2 to 2010Q4, the supTS, expTS, and H-LRT tests fail to reject the null hypothesis of a linear model, while the LMC-LRT and MMC-LRT procedures reject it. Simulation evidence from rodrondufour_mcmstest suggests that in settings where only the mean changes, the LMC-LRT and MMC-LRT tests generally have higher power and may provide more reliable results. The second panel shows results when both the mean and variance vary by regime. Here, all tests are consistent. That is, they suggest that there is insufficient evidence to reject the null hypothesis of a linear model in the first sample, but find sufficient evidence against this null hypothesis in the two larger samples.
While it would be valuable to further investigate the potential presence of a third regime in the second and third samples, this analysis is thoroughly conducted in rodrondufour_mcmstest. In their work, the authors also consider controlling for the Great Moderation and COVID periods by treating these as known structural breaks in the mean, offering a more comprehensive analysis. As they employ the {\fontseries{m}\fontseries{b}\selectfont MSTest} package for these tests, we direct interested readers to their paper for these results.
The importance of testing the number of regimes in Markov switching models has led to numerous contributions, each addressing the statistical and computational challenges inherent in this problem. Notable works include hansen92, chp14, and dufourluger17, which focus on testing the null hypothesis of a single regime (i.e., a linear model) versus the alternative of two regimes. More recently, rodrondufour_mcmstest introduced a set of Monte Carlo test procedures for testing a null hypothesis of $M_0$ regimes against an alternative of $M_0+m$ regimes, applicable when both $M_0 \geq 1$ and $m \geq 1$. The \proglang{R} package {\fontseries{m}\fontseries{b}\selectfont MSTest} makes these test procedures available, implementing the methods from these four studies. This paper reviews these procedures and explains how {\fontseries{m}\fontseries{b}\selectfont MSTest} can be used to perform these tests, as well as other package features such as simulation and model estimation. The goal of {\fontseries{m}\fontseries{b}\selectfont MSTest} is to provide researchers with a user-friendly tool for conducting these tests, facilitating inference, and determining the appropriate model specifications for their data.
The results in this paper were obtained using \proglang{R} 4.4.0 r2024 with the packages {\fontseries{m}\fontseries{b}\selectfont MSTest} 0.1.3 MSTest2024, {\fontseries{m}\fontseries{b}\selectfont Rcpp} 1.0.13 eddetal2018, {\fontseries{m}\fontseries{b}\selectfont RcppArmadillo} 14.0.2.1 eddsan2024, {\fontseries{m}\fontseries{b}\selectfont GenSA} 1.1.14.1 xiaetal2013, {\fontseries{m}\fontseries{b}\selectfont pso} 1.0.4 psopack, {\fontseries{m}\fontseries{b}\selectfont foreach} 1.5.2 foreachpack, and {\fontseries{m}\fontseries{b}\selectfont doParallel} 1.0.17 doparpack. Computations were performed on Apple macOS Sonoma Version 14.2.1 aarch64-apple-darwin20. Code for the computations is available in the R script article.R, available in the GitHub repository at \url{https://github.com/roga11/MSTest/tree/main/inst/examples}. \proglang{R} itself and all packages used are available from the Comprehensive \proglang{R} Archive Network (CRAN) at \url{https://CRAN.R-project.org/}.
This work was supported by the Fonds de recherche sur la société et la culture Doctoral Research Scholarships (B2Z).