Package 'MSTest'

Description

This package implements hypothesis testing procedures that can be used to identify the number of regimes in a Markov-Switching model.

Author(s)

Gabriel Rodriguez-Rondon, [email protected] (Maintainer)

Jean-Marie Dufour, [email protected]

References

Carrasco, Marine, Liang Hu, and Werner Ploberger. 2014. “Optimal test for Markov switching parameters.” Econometrica 82 (2): 765–784.

Dempster, A. P., N. M. Laird, and D. B. Rubin. 1977. “Maximum Likelihood from Incomplete Data via the EM Algorithm.” Journal of the Royal Statistical Society. Series B 39 (1): 1–38..

Dufour, Jean-Marie, and Richard Luger. 2017. “Identification-robust moment-based tests for Markov switching in autoregressive models.” Econometric Reviews 36 (6-9): 713–727.

Kasahara, Hiroyuk, and Katsum Shimotsu. 2018. “Testing the number of regimes in Markov regime switching models.” arXiv preprint arXiv:1801.06862.

Krolzig, Hans-Martin. 1997. “The markov-switching vector autoregressive model.”. Springer.

Hamilton, James D. 1989. “A new approach to the economic analysis of nonstationary time series and the business cycle.” Econometrica 57 (2): 357–384.

Hamilton, James D. 1994. "Time series analysis". Princeton university press.

Hansen, Bruce E. 1992. “The likelihood ratio test under nonstandard conditions: testing the Markov switching model of GNP.” Journal of applied Econometrics 7 (S1): S61–S82.

Rodriguez-Rondon, G., & Dufour, J.-M. 2026a. "Monte Carlo Likelihood-Ratio Tests for Markov Switching Models." Bank of Canada Staff Working Paper, No. 2026-23. doi: 10.34989/swp-2026-23.

Rodriguez-Rondon, G., & Dufour, J.-M. 2026b. "MSTest: An R-Package for Testing Markov Switching Models." Bank of Canada Staff Working Paper, No. 2026-7. doi: 10.34989/swp-2026-7.

Qu, Zhongjun, and Fan Zhuo. 2021. “Likelihood Ratio-Based Tests for Markov Regime Switching.” The Review of Economic Studies 88 (2): 937–968.

Description

Runs the Local Monte Carlo likelihood ratio test of Rodriguez-Rondon & Dufour (2026a) twice on the same observed series and the same pre-drawn simulation draws, for a k0 >= 2 null in the boundary (possibly absorbing regime) setting: once at the fitted null-model parameters, and once at their boundary projection – the absorbing regime's self-transition probability set to exactly 1 (its column's remaining entries to 0), everything else at its estimate. The regimes are first relabeled canonically (canonicalize_regimes: transient regimes first, absorbing regime last), and every simulated sample starts with the chain initialized in the transient regime with no burn-in (the change-point-like family of Chib, 1998) – which is why burnin = 0 is required. The test rejects only if BOTH component p-values are below the chosen level (reported as pval = max(pval_lmc, pval_fix)). At p = 0 the first observation's regime is a draw from the transient regime's transition distribution (the hidden Markov simulator applies a transition before the first emission). Multivariate models are supported: the boundary projection edits only the absorbing regime's transition-matrix column, so the fixed simulation parameter keeps the fitted regime covariances – estimates of an admissible model, positive definite by construction.

Usage

aBfixLMCTest(Y, p, k0, k1, Z = NULL, control = list())

Arguments

Value

List of class aBfixLMCTest with the same structure as aURfixLMCTest (component p-values pval_lmc, pval_fix, composition pval, the two null distributions, and theta_fix – here the boundary projection). mdl_h0 is returned in canonical regime order.

References

Rodriguez-Rondon, G., & Dufour, J.-M. 2026a. "Monte Carlo Likelihood-Ratio Tests for Markov Switching Models." Bank of Canada Staff Working Paper, No. 2026-23. doi: 10.34989/swp-2026-23.

Chib, S. 1998. "Estimation and Comparison of Multiple Change-Point Models." Journal of Econometrics 86(2): 221-241.

Description

This function estimates an autoregresive model with p lags. This can be used for the null hypothesis of a linear model against an alternative hypothesis of a Markov switching autoregressive model with k regimes.

Usage

ARmdl(Y, p, control = list())

Arguments

Value

List of class ARmdl (S3 object) with model attributes including:

• y: a (T-p x 1) matrix of observations.

• X: a (T-p x p + const) matrix of lagged observations with a leading column of 1s if const=TRUE or not if const=FALSE.

• x: a (T-p x p) matrix of lagged observations.

• fitted: a (T-p x 1) matrix of fitted values.

• resid: a (T-p x 1) matrix of residuals.

• mu: estimated mean of the process.

• beta: a ((1 + p) x 1) matrix of estimated coefficients.

• intercept: estimate of intercept.

• phi: estimates of autoregressive coefficients.

• stdev: estimated standard deviation of the process.

• sigma: estimated variance of the process.

• theta: vector containing: mu, sigma, and phi.

• theta_mu_ind: vector indicating location of mean with 1 and 0 otherwise.

• theta_sig_ind: vector indicating location of variance with 1 and 0 otherwise.

• theta_var_ind: vector indicating location of variance with 1 and 0 otherwise. This is the same as theta_sig_ind in ARmdl.

• theta_phi_ind: vector indicating location of autoregressive coefficients with 1 and 0 otherwise.

• stationary: Boolean indicating if process is stationary if TRUE or non-stationary if FALSE.

• n: number of observations after lag transformation (i.e., n = T-p).

• p: number of autoregressive lags.

• q: number of series. This is always 1 in ARmdl.

• k: number of regimes. This is always 1 in ARmdl.

• control: List with model options used.

• logLike: log-likelihood.

• AIC: Akaike information criterion.

• BIC: Bayesian (Schwarz) information criterion.

• Hess: Hessian matrix. Approximated using hessian and only returned if getSE=TRUE.

• info_mat: Information matrix. Computed as the inverse of -Hess. If matrix is not PD then nearest PD matrix is obtained using nearest_spd. Only returned if getSE=TRUE.

• nearPD_used: Boolean determining whether nearPD function was used on info_mat if TRUE or not if FALSE. Only returned if getSE=TRUE.

• theta_se: standard errors of parameters in theta. Only returned if getSE=TRUE.

References

Basawa, I. V., Mallik, A. K., McCormick, W. P., Reeves, J. H., & Taylor, R. L. 1991. "Bootstrapping Unstable First-Order Autoregressive Processes." The Annals of Statistics 19(2): 1098-1101.

See Also

MSARmdl

Examples

set.seed(1234)
# Define DGP of AR process
mdl_ar <- list(n     = 500, 
               mu    = 5,
               sigma = 2,
               phi   = c(0.5,0.2))

# Simulate process using simuAR() function
y_simu <- simuAR(mdl_ar)

# Set options for model estimation
control <- list(const  = TRUE, 
                getSE  = TRUE)

# Estimate model
y_ar_mdl <- ARmdl(y_simu$y, p = 2, control)
y_ar_mdl

Description

This function estimates an ARX model with p lags. This can be used for the null hypothesis of a linear model against an alternative hypothesis of a Markov switching autoregressive model with k regimes.

Usage

ARXmdl(Y, p, Z, control = list())

Arguments

Value

List of class ARmdl (S3 object) with model attributes including:

• y: a (T-p x 1) matrix of observations.

• X: a (T-p x p + const) matrix of lagged observations with a leading column of 1s if const=TRUE or not if const=FALSE.

• x: a (T-p x p) matrix of lagged observations.

• fitted: a (T-p x 1) matrix of fitted values.

• resid: a (T-p x 1) matrix of residuals.

• mu: estimated mean of the process.

• beta: a ((1 + p + qz) x 1) matrix of estimated coefficients.

• betaZ: a (qz x q) matrix of estimated exogenous regressor coefficients.

• intercept: estimate of intercept.

• phi: estimates of autoregressive coefficients.

• stdev: estimated standard deviation of the process.

• sigma: estimated variance of the process.

• theta: vector containing: mu, sigma, and phi.

• theta_mu_ind: vector indicating location of mean with 1 and 0 otherwise.

• theta_sig_ind: vector indicating location of variance with 1 and 0 otherwise.

• theta_var_ind: vector indicating location of variance with 1 and 0 otherwise. This is the same as theta_sig_ind in ARmdl.

• theta_phi_ind: vector indicating location of autoregressive coefficients with 1 and 0 otherwise.

• stationary: Boolean indicating if process is stationary if TRUE or non-stationary if FALSE.

• n: number of observations after lag transformation (i.e., n = T-p).

• p: number of autoregressive lags.

• q: number of series. This is always 1 in ARmdl.

• k: number of regimes. This is always 1 in ARmdl.

• control: List with model options used.

• logLike: log-likelihood.

• AIC: Akaike information criterion.

• BIC: Bayesian (Schwarz) information criterion.

• Hess: Hessian matrix. Approximated using hessian and only returned if getSE=TRUE.

• info_mat: Information matrix. Computed as the inverse of -Hess. If matrix is not PD then nearest PD matrix is obtained using nearest_spd. Only returned if getSE=TRUE.

• nearPD_used: Boolean determining whether nearPD function was used on info_mat if TRUE or not if FALSE. Only returned if getSE=TRUE.

• theta_se: standard errors of parameters in theta. Only returned if getSE=TRUE.

References

Basawa, I. V., Mallik, A. K., McCormick, W. P., Reeves, J. H., & Taylor, R. L. 1991. "Bootstrapping Unstable First-Order Autoregressive Processes." The Annals of Statistics 19(2): 1098-1101.

See Also

MSARmdl

Examples

set.seed(1234)
# Define DGP of AR process
mdl_ar <- list(n     = 500, 
               mu    = 5,
               sigma = 2,
               phi   = c(0.5,0.2))

# Simulate process using simuAR() function
y_simu <- simuAR(mdl_ar)

# Set options for model estimation
control <- list(const  = TRUE, 
                getSE  = TRUE)

# Estimate model
y_ar_mdl <- ARmdl(y_simu$y, p = 2, control)
y_ar_mdl

Description

Runs the Local Monte Carlo likelihood ratio test of Rodriguez-Rondon & Dufour (2026a) twice: on the observed series in levels with p autoregressive lags, and on its first difference with p - 1 lags (at p = 1 the differenced null has no lagged term). The test rejects only if BOTH component p-values are below the chosen level (an intersection-union composition, reported as pval = max(pval_levels, pval_diff)). Differencing imposes the null unit root without estimating it: at a one-regime null with p = 1 the differenced component's statistic does not depend on any estimated autoregressive parameter, and its unit-root guarantee is exact in finite samples – including with exogenous regressors present. For p >= 2 the differenced model keeps lagged differences, and the component is asymptotically valid at the unit root rather than exact. Exogenous regressors are kept IN LEVELS for the differenced component (rows 2 to T of Z) – they enter the observation equation additively, outside the mean adjustment, so differencing the series does not difference them. Both components are full LMCLRTest runs sharing the same control (including mc_seed, so their simulation streams start from the same seed position); each holds its own pre-drawn innovations fixed internally. Multivariate series run the same way – the difference is taken per series and both components are ordinary multivariate tests – but the unit-root validity argument is a univariate construction: a multivariate non-stationary null spans cointegration ranks 0 to q - 1, and first-differencing all series imposes the no-cointegration case, so at q > 1 the composition's unit-root guarantee is conditional on that maintained assumption.

Usage

aURdiffLMCTest(Y, p, k0, k1, Z = NULL, control = list())

Arguments

Value

List of class aURdiffLMCTest with pval_levels, pval_diff (the component p-values), pval (their maximum; reject only if it is below the level), lmc_levels and lmc_diff (the two full LMCLRTest objects), and control.

References

Rodriguez-Rondon, G., & Dufour, J.-M. 2026a. "Monte Carlo Likelihood-Ratio Tests for Markov Switching Models." Bank of Canada Staff Working Paper, No. 2026-23. doi: 10.34989/swp-2026-23.

Description

Runs the Maximized Monte Carlo likelihood ratio test of Rodriguez-Rondon & Dufour (2026a) twice: on the observed series in levels with p autoregressive lags (typically with the unit-root search box of MMCLRTest's box_include control, so the search set contains the unit root), and on its first difference with p - 1 lags, each with its own consistent-set search box around its own fitted null model. The test rejects only if BOTH component p-values are below the chosen level (reported as pval = max(pval_levels, pval_diff)). Multivariate series are handled as in aURdiffLMCTest, with the same scope note: first-differencing all series imposes the no-cointegration case, so at q > 1 the unit-root guarantee is conditional on that maintained assumption.

Usage

aURdiffMMCTest(Y, p, k0, k1, Z = NULL, control = list(), control_diff = NULL)

Arguments

Value

List of class aURdiffMMCTest with pval_levels, pval_diff, pval (their maximum), mmc_levels and mmc_diff (the two full MMCLRTest objects), control, and control_diff.

References

Rodriguez-Rondon, G., & Dufour, J.-M. 2026a. "Monte Carlo Likelihood-Ratio Tests for Markov Switching Models." Bank of Canada Staff Working Paper, No. 2026-23. doi: 10.34989/swp-2026-23.

Description

Runs the Local Monte Carlo likelihood ratio test of Rodriguez-Rondon & Dufour (2026a) twice on the same observed series and the same pre-drawn simulation draws: once with the null distribution simulated at the fitted null-model parameters (the ordinary LMC-LRT), and once with the null-simulation parameter's long-run autoregressive coefficient moved to exactly 1 (everything else unchanged), so the second component's simulated samples have an exact unit root – the Monte Carlo counterpart of imposing the unit root in the resampling scheme (Basawa, Mallik, McCormick, Reeves & Taylor, 1991). The test rejects only if BOTH component p-values are below the chosen level (an intersection-union composition, reported as pval = max(pval_lmc, pval_fix)), which protects the level when the data may carry a unit root. For p > 1 the long-run coefficient sum(phi) is moved to 1 while the short-run (augmented Dickey-Fuller) coefficients are held fixed – linear arithmetic on the fitted vector, nothing re-estimated; the resulting vector is reported as theta_fix. Exogenous regressors, when present, keep their estimated coefficients in the unit-root-fixed simulation parameter. Estimation on the simulated samples stays free (unrestricted) in both components, identically to the observed side. At a one-regime null with p = 1 and no exogenous regressors the imposed-root component's guarantee at the unit root is exact in finite samples: the fixed simulation parameter carries no estimated autoregressive coefficient, and the remaining location-scale nuisance is absorbed by invariance (the level cancels at the unit root, the scale by equivariance). For p >= 2 the short-run coefficients enter the fixed parameter at their estimates and the component is asymptotically valid at the unit root, as it is with regressors present; it is the differencing construction that retains finite-sample exactness at p = 1 with regressors – see aURdiffLMCTest when exactness with Z is needed. Univariate only: the augmented Dickey-Fuller map has no single multivariate counterpart (a VAR's unit-root null spans cointegration ranks 0 to q - 1, so a fixed point would impose an unstated rank choice); for a multivariate series use aURdiffLMCTest, whose differencing construction needs no fixed point.

Usage

aURfixLMCTest(Y, p, k0, k1, Z = NULL, control = list())

Arguments

Value

List of class aURfixLMCTest with mdl_h0, mdl_h1, LRT_0 (the shared observed statistic), pval_lmc and pval_fix (the component p-values), pval (their maximum; reject only if it is below the level), LRN and LRN_fix (the two simulated null distributions), LRN_cv and LRN_fix_cv, theta_fix (the unit-root-fixed parameter vector), and control.

References

Rodriguez-Rondon, G., & Dufour, J.-M. 2026a. "Monte Carlo Likelihood-Ratio Tests for Markov Switching Models." Bank of Canada Staff Working Paper, No. 2026-23. doi: 10.34989/swp-2026-23.

Basawa, I. V., Mallik, A. K., McCormick, W. P., Reeves, J. H., & Taylor, R. L. 1991. "Bootstrapping Unstable First-Order Autoregressive Processes." The Annals of Statistics 19(2): 1098-1101.

Basawa, I. V., Mallik, A. K., McCormick, W. P., & Taylor, R. L. 1991. "Bootstrap Test of Significance and Sequential Bootstrap Estimation for Unstable First Order Autoregressive Processes." Communications in Statistics - Theory and Methods 20(3): 1015-1026.

Description

Relabels the regimes of a fitted Markov-switching model in ascending order of fitted persistence (P[j,j]), so that state 1 is the least persistent (transient) regime and the last state is the most persistent. Estimation is invariant to regime relabeling, so a fitted model's labels are otherwise arbitrary; this makes them deterministic. Never called automatically by any estimation or testing function.

Usage

canonicalize_regimes(mdl, by = "persistence")

Arguments

Value

An object of the same class as mdl, with every field on the regime axis permuted to the canonical order (theta, mu, sigma, stdev, P, pinf, St, intercept, beta, standard errors when present, and the constraint-binding diagnostics). logLike, AIC, BIC, resid and fitted are unchanged (a relabeling does not change the likelihood). Fields with no regime axis (phi, betaZ, y, Z, n, p, q, k, control, the theta_*_ind position maps, theta_0, init_used, trace) are unchanged.

References

Rodriguez-Rondon, G., & Dufour, J.-M. 2026a. "Monte Carlo Likelihood-Ratio Tests for Markov Switching Models." Bank of Canada Staff Working Paper, No. 2026-23. doi: 10.34989/swp-2026-23.

Description

This data is the extension of the GNP series used in CHP (2014), Econometrica. This series ranges from 1951Q2 to 2010Q4.

Usage

chp10GNP

Format

Date

Vector of dates

GNP

US GNP series

GNP_gr

log difference of US GNP series

Source

https://www.econometricsociety.org/content/supplement-optimal-test-markov-switching-parameters

References

Carrasco, Marine, Liang Hu, and Werner Ploberger. 2014. “Optimal test for Markov switching parameters.” Econometrica 82 (2): 765–784.

Description

Parametric bootstrap described on p. 771 of Carrasco, Hu & Ploberger (2014): simulate under the null at the fitted parameters, re-estimate, and recompute the test statistic, entirely in C++ (no calls back into R per replication).

Usage

CHPbootCV(mdl, rho_b, N, msvar, burnin = 100L)

Arguments

Value

A (N x 2) matrix with the bootstrapped supTS statistic in the first column and expTS in the second.

References

Carrasco, Marine, Liang Hu, and Werner Ploberger. 2014. “Optimal test for Markov switching parameters.” Econometrica 82 (2): 765–784.

Description

This function performs the CHP (2014) parameter stability test as outline in Carrasco, M., Hu, L. and Ploberger, W. (2014). Original source code can be found here.

Usage

CHPTest(Y, p, control = list())

Arguments

Value

List of class CHPTest (S3 object) with model attributes including:

• mdl_h0: List with restricted model attributes. This will be of class ARmdl (S3 object). See ARmdl.

• supTS: supTS test statistic value.

• expTS: expTS test statistic value.

• supTS_N: A (N x 1) vector with simulated supTS test statistics under null hypothesis.

• expTS_N: A (N x 1) vector with simulated expTS test statistics under null hypothesis.

• pval_supTS: P-value for supTS version of parameter stability test.

• pval_expTS: P-value for expTS version of parameter stability test.

• supTS_cv: Vector with 90%, 95%, and 99% bootstrap critical values for supTS version of parameter stability test.

• expTS_cv: Vector with 90%, 95%, and 99% bootstrap critical values for expTS version of parameter stability test.

• control: List with test procedure options used.

References

Carrasco, Marine, Liang Hu, and Werner Ploberger. 2014. “Optimal test for Markov switching parameters.” Econometrica 82 (2): 765–784.

Examples

# load data used in Hamilton 1989 
y84 <- as.matrix(hamilton84GNP$GNP_gr)


# Set test procedure options
control = list(N = 1000, 
               rho_b = 0.7, 
               msvar = FALSE)

# perform test with switch in mean only on Hamilton 1989 data

  mdl_84_msmu <- CHPTest(y84, p = 4, control = control)
  summary(mdl_84_msmu)

Description

This function performs the Local Monte Carlo moment-based test for Markov switching autoregressive models proposed in Dufour & Luger (2017).

Usage

DLMCTest(Y, p, control = list())

Arguments

Value

List of class DLMCTest (S3 object) with attributes including:

• mdl_h0: List with restricted model attributes. This will be of class ARmdl if p>0 or Nmdl otherwise (S3 objects). See ARmdl or Nmdl.

• theta: Value of nuisance parameters. Specifically, these are the consistent estimates of nuisance parameters as discussed in Dufour & Luger (2017) LMC procedure.

• S0: A (1 x 4)) matrix containing the four moment-based test statistics defined in (11) - (14) in Dufour & Luger (2017).

• F0_min: Test statistic value for min version of Local Monte Carlo moment-based test.

• F0_prod: Test statistic value for prod version of Local Monte Carlo moment-based test.

• FN_min: A (N x 1) vector with simulated test statistics for min version of Local Monte Carlo moment-based test under null hypothesis.

• FN_prod: A (N x 1) vector with simulated test statistics for prod version of Local Monte Carlo moment-based test under null hypothesis.

• pval_min: P-value for min version of Local Monte Carlo moment-based test.

• pval_prod: P-value for prod version of Local Monte Carlo moment-based test.

• FN_min_cv: Vector with 90%, 95%, and 99% Monte Carlo critical values for min version of Local Monte Carlo moment-based test.

• FN_prod_cv: Vector with 90%, 95%, and 99% Monte Carlo critical values for prod version of Local Monte Carlo moment-based test.

• control: List with test procedure options used.

References

Dufour, J. M., & Luger, R. 2017. "Identification-robust moment-based tests for Markov switching in autoregressive models." Econometric Reviews, 36(6-9), 713-727.

Examples

set.seed(1234)
# load data used in Hamilton 1989 and extended data used in CHP 2014 
y84 <- as.matrix(hamilton84GNP$GNP_gr)
y10 <- as.matrix(chp10GNP$GNP_gr)

# Set test procedure options
lmc_control = list(N = 99,
                   simdist_N = 10000,
                   getSE = TRUE)

# perform test on Hamilton 1989 data
lmc_gnp84 <- DLMCTest(y84, p = 4, control = lmc_control)
summary(lmc_gnp84)

Description

This function performs the maximized Monte Carlo moment-based test for Markov switching autoregressive models proposed in Dufour & Luger (2017).

Usage

DLMMCTest(Y, p, control = list())

Arguments

Value

List of class DLMCTest (S3 object) with attributes including:

• mdl_h0: List with restricted model attributes. This will be of class ARmdl if p>0 or Nmdl otherwise (S3 objects). See ARmdl or Nmdl.

• theta_max_min: Value of nuisance parameters when min version of p-value is maximized as discussed in Dufour & Luger (2017) MMC procedure.

• theta_max_prod: Value of nuisance parameters when prod version of p-value is maximized as discussed in Dufour & Luger (2017) MMC procedure.

• theta_low: Lower bound on nuisance parameter values used when searching for maximum p-value.

• theta_upp: Upper bound on nuisance parameter values used when searching for maximum p-value.

• S0_min: A (1 x 4)) matrix containing the four moment-based test statistics defined in (11) - (14) in Dufour & Luger (2017) when theta_min is used.

• S0_prod: A (1 x 4)) matrix containing the four moment-based test statistics defined in (11) - (14) in Dufour & Luger (2017) when theta_prod is used.

• F0_min: Test statistic value for min version of Maximized Monte Carlo moment-based test.

• F0_prod: Test statistic value for prod version of Maximized Monte Carlo moment-based test.

• FN_min: A (N x 1) vector with simulated test statistics for min version of Maximized Monte Carlo moment-based test under null hypothesis.

• FN_prod: A (N x 1) vector with simulated test statistics for prod version of Maximized Monte Carlo moment-based test under null hypothesis.

• pval_min: Maximum p-value for min version of Maximized Monte Carlo moment-based test.

• pval_prod: Maximum p-value for prod version of Local Monte Carlo moment-based test.

• FN_min_cv: Vector with 90%, 95%, and 99% Monte Carlo critical values for min version of Local Monte Carlo moment-based test.

• FN_prod_cv: Vector with 90%, 95%, and 99% Monte Carlo critical values for prod version of Local Monte Carlo moment-based test.

• control: List with test procedure options used.

• optim_min_output: List with optimization output for min version of Maximized Monte Carlo moment-based test.

• optim_prod_output: List with optimization output for prod version of Maximized Monte Carlo moment-based test.

References

Dufour, J. M., & Luger, R. 2017. "Identification-robust moment-based tests for Markov switching in autoregressive models." Econometric Reviews, 36(6-9), 713-727.

Examples

set.seed(1234)
# load data used in Hamilton 1989 and extended data used in CHP 2014 
y84 <- as.matrix(hamilton84GNP$GNP_gr)
y10 <- as.matrix(chp10GNP$GNP_gr)

# Set test procedure options
mmc_control <- list(N = 99,
                    simdist_N = 10000,
                    getSE = TRUE,
                    eps = 0, 
                    CI_union = TRUE,
                    lambda = 100,
                    stationary_ind = TRUE,
                    optim_type = "GenSA",
                    silence = FALSE,
                    threshold_stop = 1,
                    maxit = 200)


# perform test on Hamilton 1989 data

  mmc_gnp84 <- DLMMCTest(y84, p = 4, control = mmc_control)
  summary(mmc_gnp84)

Description

This data set is used in Hansen (1992) to test the US GNP model proposed by Hamilton (1989). This series ranges from 1951Q2 to 1984Q4.

Usage

hamilton84GNP

Format

Date

Vector of dates

GNP

US GNP series

GNP_gr

US GNP log difference

Source

https://www.ssc.wisc.edu/~bhansen/progs/jae_92.html

References

Hansen, Bruce E. 1992. “The likelihood ratio test under nonstandard conditions: testing the Markov switching model of GNP.” Journal of applied Econometrics 7 (S1): S61–S82.

Hamilton, James D. 1989. “A new approach to the economic analysis of nonstationary time series and the business cycle.” Econometrica 57 (2): 357–384.

Description

This function performs Hansen's likelihood ratio test as described in Hansen (1992). Original source code can be found here.

Usage

HLRTest(Y, p, control = list())

Arguments

Value

List of class HLRTest (S3 object) with model attributes including:

• mdl_h0: List with restricted model attributes. This will be of class ARmdl (S3 object). See ARmdl.

• LR0: Likelihood ratio test statistic value.

• LRN: A (N x 1) vector with simulated LRT statistics under null hypothesis.

• pval: P-value.

• LR_cv: A (nwband x 3) matrix with 90%, 95%, and 99% critical values in each column respectively.

• coef: Vector of coefficients from restricted model and grid search that maximized standardized LRT.

• control: List with test procedure options used.

References

Hansen, Bruce E. 1992. “The likelihood ratio test under nonstandard conditions: testing the Markov switching model of GNP.” Journal of applied Econometrics 7 (S1): S61–S82.

Hansen, Bruce E. 1996. “Erratum: The likelihood ratio test under nonstandard conditions: testing the Markov switching model of GNP.” Journal of applied Econometrics 7 (S1): S61–S82.

Examples

# --------------------------- Use simulated process ----------------------------
set.seed(1234)
# Define DGP of MS AR process
mdl_ms2 <- list(n     = 200, 
                mu    = c(5,1),
                sigma = c(1,1),
                phi   = c(0.5),
                k     = 2,
                P     = rbind(c(0.90, 0.10),
                              c(0.10, 0.90)))

# Simulate process using simuMSAR() function
y_ms_simu <- simuMSAR(mdl_ms2)

hlrt_control  <- list(ix          = 1, 
                      gridsize    = 7,
                      pgrid_from  = 0.05,
                      pgrid_by    = 0.05,
                      pgrid_to    = 0.95,
                      mugrid_from = 0,
                      mugrid_by   = 1)


  hlrt <- HLRTest(y_ms_simu$y, p = 1, control = hlrt_control)
  summary(hlrt)

Description

This function estimates a Hidden Markov model with k regimes.

Usage

HMmdl(Y, k, Z = NULL, control = list())

Arguments

Value

List of class HMmdl (S3 object) with model attributes including:

• y: a (T x q) matrix of observations.

• fitted: a (T x q) matrix of fitted values.

• resid: a (T x q) matrix of residuals.

• mu: a (k x q) matrix of estimated means of each process.

• beta: if q=1, this is a ((1 + qz) x k) matrix of estimated coefficients. If q>1, this is a list containing k separate ((1 + qz) x q) matrix of estimated coefficients for each regime.

• betaZ: a (qz x q) matrix of estimated exogenous regressor coefficients.

• intercept: a (k x q) matrix of estimated intercept of each process. If Z is Null, this is the same as mu.

• stdev: If q=1, this is a (k x 1) matrix with estimated standard. If q>1, this is a List with k (q x q) matrices with estimated standard deviation on the diagonal.

• sigma: If q=1, this is a (k x 1) matrix with variances. If q>1, this is a List with k (q x q) estimated covariance matrix.

• theta: vector containing: mu and vech(sigma).

• theta_mu_ind: vector indicating location of mean with 1 and 0 otherwise.

• theta_sig_ind: vector indicating location of variance and covariances with 1 and 0 otherwise.

• theta_var_ind: vector indicating location of variances with 1 and 0 otherwise.

• theta_P_ind: vector indicating location of transition matrix elements with 1 and 0 otherwise.

• n: number of observations (same as T).

• q: number of series.

• k: number of regimes in estimated model.

• P: a (k x k) transition matrix.

• pinf: a (k x 1) vector with limiting probabilities of each regime.

• St: a (T x k) vector with smoothed probabilities of each regime at each time t.

• deltath: double with the maximum relative change in vector theta on the last iteration.

• iterations: number of EM iterations performed to achieve convergence (if less than maxit).

• converged: Boolean. TRUE if the conv convergence criterion was met before maxit iterations, FALSE otherwise.

• theta_0: vector of initial values used.

• init_used: number of different initial values used to get a finite solution. See description of input maxit_converge.

• msmu: Boolean. If TRUE model was estimated with switch in mean. If FALSE model was estimated with constant mean.

• msvar: Boolean. If TRUE model was estimated with switch in variance. If FALSE model was estimated with constant variance.

• control: List with model options used.

• logLike: log-likelihood.

• AIC: Akaike information criterion.

• BIC: Bayesian (Schwarz) information criterion.

• Hess: Hessian matrix. Approximated using hessian and only returned if getSE=TRUE.

• info_mat: Information matrix. Computed as the inverse of -Hess. If matrix is not PD then nearest PD matrix is obtained using nearest_spd. Only returned if getSE=TRUE.

• nearPD_used: Boolean determining whether nearPD function was used on info_mat if TRUE or not if FALSE. Only returned if getSE=TRUE.

• theta_se: standard errors of parameters in theta. Only returned if getSE=TRUE.

• trace: List with Lists of estimation output for each initial value used due to use_diff_init > 1.

References

Dempster, A. P., N. M. Laird, and D. B. Rubin. 1977. “Maximum Likelihood from Incomplete Data via the EM Algorithm.” Journal of the Royal Statistical Society. Series B 39 (1): 1–38..

Böhning, D., E. Dietz, R. Schaub, P. Schlattmann, and B. G. Lindsay. 1994. “The distribution of the likelihood ratio for mixtures of densities from the one-parameter exponential family.” Annals of the Institute of Statistical Mathematics 46 (2): 373–388.

McLachlan, G. J., and T. Krishnan. 2008. The EM Algorithm and Extensions. 2nd ed. Hoboken, New Jersey: John Wiley & Sons.

Hamilton, James D. 1990. “Analysis of time series subject to changes in regime.” Journal of econometrics, 45 (1-2): 39–70.

Krolzig, Hans-Martin. 1997. “The markov-switching vector autoregressive model.”. Springer.

See Also

Nmdl

Examples

set.seed(1234)
# ----- Univariate ----- #
# Define DGP 
mdl_hmm <- list(n     = 1000, 
                q     = 1,
                mu    = as.matrix(c(5,-2)),
                sigma = list(as.matrix(5.0),
                             as.matrix(7.0)),
                k     = 2,
                P     = rbind(c(0.90, 0.10),
                              c(0.10, 0.90)))

# Simulate process using simuHMM() function
y_hmm_simu <- simuHMM(mdl_hmm)

# Set options for model estimation
control <- list(msmu   = TRUE, 
                msvar  = TRUE,
                method = "EM",
                use_diff_init = 1)


# Estimate model
y_hmm_mdl <- HMmdl(y_hmm_simu$y, k = 2, control = control)
summary(y_hmm_mdl)



# ----- Multivariate normal process ----- #
# Define DGP 
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.90, 0.10),
                              c(0.10, 0.90)))

# Simulate process using simuHMM() function
y_hmm_simu <- simuHMM(mdl_hmm)

# Set options for model estimation
control <- list(msmu   = TRUE, 
                msvar  = TRUE,
                method = "EM",
                use_diff_init = 1)

# Estimate model

  y_hmm_mdl <- HMmdl(y_hmm_simu$y, k = 2, control = control)
  summary(y_hmm_mdl)

Description

This function performs the Local Monte Carlo likelihood ratio test (LMC-LRT) proposed in Rodriguez-Rondon & Dufour (2026a). As discussed in their work, this test can be applied in very general settings and can be used to compare varioous regimes under the null and under the alternative.

Usage

LMCLRTest(Y, p, k0, k1, Z = NULL, control = list())

Arguments

Value

List of class LMCLRTest (S3 object) with attributes including:

• mdl_h0: List with restricted model attributes.

• mdl_h1: List with unrestricted model attributes.

• LRT_0: Value of test statistic from observed data.

• LRN: A (N x 1) vector of test statistics from data simulated under the null hypothesis.

• pval: P-value of Local Monte Carlo Likelihood Ratio Test.

• LRN_cv: Vector with 90%, 95%, and 99% Monte Carlo simulated critical values (from vector LRN). These are not asymptotic critical values.

• control: List with test procedure options used.

Numerical policy for the LR statistic

For nested models the LR statistic is non-negative by construction (the null fit embeds into the alternative parameter space with equal likelihood), so a negative computed value can only be a numerical optimization artifact. Both the observed statistic and every simulated statistic are therefore floored at 0 (identically on both sides, which preserves the exchangeability underlying the Monte Carlo p-value); a warning reports the magnitude when it is non-negligible. For k0 > 1, the null log-likelihood is additionally floored at the one-regime (linear) fit on both sides. Ties at 0 are handled by the randomized tie-breaking rule of Dufour (2006) already implemented in MCpval; ties at the bottom of the support can only make the test conservative. The regime variances are bounded away from zero during estimation (control em_variance_constraint of the model constructors, default 0.01 of the one-regime fit variance), because the Markov-switching likelihood is otherwise unbounded and degenerate solutions produce spurious likelihood-ratio spikes; the statistic is therefore a constrained-likelihood ratio, computed by the same rule for the observed and the simulated samples, which is all the Monte Carlo argument requires. Transition probabilities are not constrained.

References

Rodriguez-Rondon, G., & Dufour, J.-M. 2026a. "Monte Carlo Likelihood-Ratio Tests for Markov Switching Models." Bank of Canada Staff Working Paper, No. 2026-23. doi: 10.34989/swp-2026-23.

Examples

set.seed(1234)
# Define DGP of MS AR process
mdl_ms2 <- list(n     = 200, 
                mu    = c(5,10),
                sigma = c(1,4),
                phi   = c(0.5),
                k     = 2,
                P     = rbind(c(0.90, 0.10),
                              c(0.10, 0.90)))

# Simulate process using simuMSAR() function
y_ms_simu <- simuMSAR(mdl_ms2)


# ------ MS-AR example ----- #
# Set test procedure options
lmc_control = list(N = 19,
                   burnin = 100,
                   converge_check = NULL,
                   mdl_h0_control = list(const  = TRUE, 
                                         getSE  = FALSE),
                   mdl_h1_control = list(msmu   = TRUE, 
                                         msvar  = TRUE,
                                         getSE  = FALSE,
                                         method = "EM",
                                         use_diff_init = 1))


# lmctest <- LMCLRTest(y_ms_simu$y, p = 1, k0 = 1 , k1 = 2, control = lmc_control)
# summary(lmctest)

Description

This function computes the Monte Carlo P-value.

Usage

MCpval(test_stat, null_vec, type = "geq")

Arguments

Value

MC p-value of test

References

Dufour, Jean-Marie 2006. "Monte Carlo tests with nuisance parameters: A general approach to finite-sample inference and nonstandard asymptotics". Journal of Econometrics, 133(2), 443-477.

Dufour, Jean-Marie, and Richard Luger. 2017. "Identification-robust moment-based tests for Markov switching in autoregressive models". Econometric Reviews, 36(6-9), 713-727.

Description

This function performs the Maximized Monte Carlo likelihood ratio test (MMC-LRT) proposed in Rodriguez-Rondon & Dufour (2026a).

Usage

MMCLRTest(Y, p, k0, k1, Z = NULL, control = list())

Arguments

Value

List of class MMCLRTest (S3 object) with attributes including:

• mdl_h0: List with restricted model attributes.

• mdl_h1: List with unrestricted model attributes.

• mdl_h0_mmc: List with restricted model attributes at the maximizing nuisance parameter values.

• mdl_h1_mmc: List with unrestricted model attributes (same as mdl_h1).

• LRT_0: Value of test statistic from observed data.

• pval: Maximized p-value of the Maximized Monte Carlo Likelihood Ratio Test.

• theta_h0: Maximizing nuisance parameter vector for the restricted model.

• theta_h1: Parameter vector of the unrestricted model.

• control: List with test procedure options used.

• mmc_optimout: Optimization output object returned by the selected optimizer (pso, GenSA, or GA). For k0 >= 2 the search itself runs over a reduced parameter vector (the transition matrix contributes only its free entries), and the optimizer's reported solution ($par for pso/GenSA, @solution for GA) is mapped back to the full theta coordinates used everywhere else in the return value; other diagnostic fields inside mmc_optimout (e.g. pso's trace.stats$x, GenSA's trace.mat) are left in the reduced search coordinates.

• pval_0: Monte Carlo p-value at the initial (estimated) nuisance parameter values theta_0 – the Local MC point of the search. Present for every estimate-centered search; the reported pval satisfies pval >= pval_0 by construction (enforced directly, not merely relied on from the optimizer's own search). With box_boundary = "face" the fitted point lies outside the searched set, so this field is replaced by pval_proj: the evaluation at the unclamped boundary projection of the fitted parameters (reported as theta_face_proj), with pval >= pval_proj enforced the same way.

• theta_low, theta_upp: The realized search box (full parameter coordinates, named) after every widening and clamp – stored so it can be verified from saved output whether the box contained a target value.

• optim_diagnostics: List with type (the optimizer, or "early_stop"), convergence (the optimizer's own exit code where it reports one), feval (exact number of objective evaluations, counted directly), optimizer_feval (the optimizer's own reported call count where one exists, as a cross-check), pd_penalty_hits (number of candidates excluded because an assembled regime covariance was not positive definite; NA when that gate is not active, i.e. for univariate models or when neither box_boundary nor box_include is in use), margin_penalty_hits (number of candidates excluded by the stationarity-margin gate of box_match_constraints; NA when no margin is in force), and evaluated_range (per-coordinate min/max of the candidates the optimizer evaluated, in the reduced search coordinates, named by the searched positions; NULL when the search was skipped by threshold_stop).

• boundary_mode, absorbing_regime, theta_face_proj, seed_repair: Present only when box_boundary is used – the mode, the absorbing regime's index after canonical labeling (always k0), and (face mode) the unclamped boundary projection whose evaluation is pval_proj together with seed_repair (whether the clamped search seed's transient covariance needed the positive-definiteness repair, and the off-diagonal shrinkage factor applied). With box_boundary set, mdl_h0 is returned in canonical regime order (see canonicalize_regimes).

Numerical policy for the LR statistic

For nested models the LR statistic is non-negative by construction (the null fit embeds into the alternative parameter space with equal likelihood), so a negative computed value can only be a numerical optimization artifact. Both the observed statistic and every simulated statistic are therefore floored at 0 (identically on both sides, which preserves the exchangeability underlying the Monte Carlo p-value); a warning reports the magnitude when it is non-negligible. For k0 > 1, the null log-likelihood is additionally floored at the one-regime (linear) fit on both sides. Ties at 0 are handled by the randomized tie-breaking rule of Dufour (2006) already implemented in MCpval; ties at the bottom of the support can only make the test conservative. The regime variances are bounded away from zero during estimation (control em_variance_constraint of the model constructors, default 0.01 of the one-regime fit variance), because the Markov-switching likelihood is otherwise unbounded and degenerate solutions produce spurious likelihood-ratio spikes; the statistic is therefore a constrained-likelihood ratio, computed by the same rule for the observed and the simulated samples, which is all the Monte Carlo argument requires. Transition probabilities are not constrained.

References

Rodriguez-Rondon, G., & Dufour, J.-M. 2026a. "Monte Carlo Likelihood-Ratio Tests for Markov Switching Models." Bank of Canada Staff Working Paper, No. 2026-23. doi: 10.34989/swp-2026-23.

Examples

set.seed(1234)
# Define DGP of MS AR process
mdl_ms2 <- list(n     = 200, 
                mu    = c(5,10),
                sigma = c(1,4),
                phi   = c(0.5),
                k     = 2,
                P     = rbind(c(0.90, 0.10),
                              c(0.10, 0.90)))

# Simulate process using simuMSAR() function
y_ms_simu <- simuMSAR(mdl_ms2)


# Set test procedure options
mmc_control = list(N = 19,
                   burnin = 100,
                   converge_check = NULL,
                   eps = 0.1,
                   CI_union = TRUE,
                   silence = FALSE,
                   threshold_stop = 0.05 + 1e-6,
                   type = "pso",
                   maxit = 50,
                   mdl_h0_control = list(const  = TRUE, 
                                         getSE  = TRUE),
                   mdl_h1_control = list(msmu   = TRUE, 
                                         msvar  = TRUE,
                                         getSE  = TRUE,
                                         method = "EM",
                                         use_diff_init = 1))
  
# MMCtest <- MMCLRTest(y_ms_simu$y, p = 1 , k0 = 1, k1 = 2, control = mmc_control)
# summary(MMCtest)

Description

This function estimates a Markov-switching autoregressive model

Usage

MSARmdl(Y, p, k, control = list())

Arguments

Value

List of class MSARmdl (S3 object) with model attributes including:

• y: a (T x 1) matrix of observations.

• X: a (T-p x p + const) matrix of lagged observations with a leading column of 1s.

• x: a (T-p x p) matrix of lagged observations.

• fitted: a (T x 1) matrix of fitted values.

• resid: a (T x 1) matrix of residuals.

• intercept: a (k x 1) vector of estimated intercepts of each process.

• mu: a (k x 1) vector of estimated means of each process.

• beta: a ((1 + p) x k) matrix of estimated coefficients.

• phi: estimates of autoregressive coefficients.

• stdev: a (k x 1) vector of estimated standard deviation of each process.

• sigma: a (k x 1) estimated covariance matrix.

• theta: vector containing: mu and vech(sigma).

• theta_mu_ind: vector indicating location of mean with 1 and 0 otherwise.

• theta_sig_ind: vector indicating location of variances with 1 and 0 otherwise.

• theta_var_ind: vector indicating location of variances with 1 and 0 otherwise. This is the same as theta_sig_ind in MSARmdl.

• theta_P_ind: vector indicating location of transition matrix elements with 1 and 0 otherwise.

• stationary: Boolean indicating if process is stationary if TRUE or non-stationary if FALSE.

• n: number of observations (same as T).

• p: number of autoregressive lags.

• q: number of series. This is always 1 in MSARmdl.

• k: number of regimes in estimated model.

• P: a (k x k) transition matrix.

• pinf: a (k x 1) vector with limiting probabilities of each regime.

• St: a (T x k) vector with smoothed probabilities of each regime at each time t.

• deltath: double with the maximum relative change in vector theta on the last iteration.

• iterations: number of EM iterations performed to achieve convergence (if less than maxit).

• converged: Boolean. TRUE if the conv convergence criterion was met before maxit iterations, FALSE otherwise.

• theta_0: vector of initial values used.

• init_used: number of different initial values used to get a finite solution. See description of input maxit_converge.

• msmu: Boolean. If TRUE model was estimated with switch in mean. If FALSE model was estimated with constant mean.

• msvar: Boolean. If TRUE model was estimated with switch in variance. If FALSE model was estimated with constant variance.

• control: List with model options used.

• logLike: log-likelihood.

• AIC: Akaike information criterion.

• BIC: Bayesian (Schwarz) information criterion.

• Hess: Hessian matrix. Approximated using hessian and only returned if getSE=TRUE.

• info_mat: Information matrix. Computed as the inverse of -Hess. If matrix is not PD then nearest PD matrix is obtained using nearest_spd. Only returned if getSE=TRUE.

• nearPD_used: Boolean determining whether nearPD function was used on info_mat if TRUE or not if FALSE. Only returned if getSE=TRUE.

• theta_se: standard errors of parameters in theta. Only returned if getSE=TRUE.

• trace: List with Lists of estimation output for each initial value used due to use_diff_init > 1.

References

Dempster, A. P., N. M. Laird, and D. B. Rubin. 1977. “Maximum Likelihood from Incomplete Data via the EM Algorithm.” Journal of the Royal Statistical Society. Series B 39 (1): 1–38..

Böhning, D., E. Dietz, R. Schaub, P. Schlattmann, and B. G. Lindsay. 1994. “The distribution of the likelihood ratio for mixtures of densities from the one-parameter exponential family.” Annals of the Institute of Statistical Mathematics 46 (2): 373–388.

McLachlan, G. J., and T. Krishnan. 2008. The EM Algorithm and Extensions. 2nd ed. Hoboken, New Jersey: John Wiley & Sons.

Hamilton, James D. 1990. “Analysis of time series subject to changes in regime.” Journal of econometrics, 45 (1-2): 39–70.

Basawa, I. V., Mallik, A. K., McCormick, W. P., Reeves, J. H., & Taylor, R. L. 1991. "Bootstrapping Unstable First-Order Autoregressive Processes." The Annals of Statistics 19(2): 1098-1101.

See Also

ARmdl

Examples

# --------------------------- Use simulated process ----------------------------
set.seed(1234)
# Define DGP of MS AR process
mdl_ms2 <- list(n     = 500, 
                mu    = c(5,10),
                sigma = c(1,4),
                phi   = c(0.5),
                k     = 2,
                P     = rbind(c(0.90, 0.10),
                              c(0.10, 0.90)))

# Simulate process using simuMSAR() function
y_ms_simu <- simuMSAR(mdl_ms2)

# Set options for model estimation
control <- list(msmu   = TRUE, 
                msvar  = TRUE, 
                method = "EM",
                use_diff_init = 1)

# Estimate model

  ms_mdl <- MSARmdl(y_ms_simu$y, p = 1, k = 2, control = control)
  summary(ms_mdl)

Description

This function estimates a Markov-switching autoregressive model with exogenous regressors.

Usage

MSARXmdl(Y, p, k, Z, control = list())

Arguments

Value

List of class MSARXmdl (S3 object) with model attributes including:

• y: a (T x 1) matrix of observations.

• X: a (T-p x p + const) matrix of lagged observations with a leading column of 1s.

• x: a (T-p x p) matrix of lagged observations.

• fitted: a (T x 1) matrix of fitted values.

• resid: a (T x 1) matrix of residuals.

• intercept: a (k x 1) vector of estimated intercepts of each process.

• mu: a (k x 1) vector of estimated means of each process.

• beta: a ((1 + p + qz) x k) matrix of estimated coefficients.

• betaZ: a (qz x q) matrix of estimated exogenous regressor coefficients.

• phi: estimates of autoregressive coefficients.

• stdev: a (k x 1) vector of estimated standard deviation of each process.

• sigma: a (k x 1) estimated covariance matrix.

• theta: vector containing: mu and vech(sigma).

• theta_mu_ind: vector indicating location of mean with 1 and 0 otherwise.

• theta_sig_ind: vector indicating location of variances with 1 and 0 otherwise.

• theta_var_ind: vector indicating location of variances with 1 and 0 otherwise. This is the same as theta_sig_ind in MSARmdl.

• theta_P_ind: vector indicating location of transition matrix elements with 1 and 0 otherwise.

• stationary: Boolean indicating if process is stationary if TRUE or non-stationary if FALSE.

• n: number of observations (same as T).

• p: number of autoregressive lags.

• q: number of series. This is always 1 in MSARmdl.

• k: number of regimes in estimated model.

• P: a (k x k) transition matrix.

• pinf: a (k x 1) vector with limiting probabilities of each regime.

• St: a (T x k) vector with smoothed probabilities of each regime at each time t.

• deltath: double with the maximum relative change in vector theta on the last iteration.

• iterations: number of EM iterations performed to achieve convergence (if less than maxit).

• converged: Boolean. TRUE if the conv convergence criterion was met before maxit iterations, FALSE otherwise.

• theta_0: vector of initial values used.

• init_used: number of different initial values used to get a finite solution. See description of input maxit_converge.

• msmu: Boolean. If TRUE model was estimated with switch in mean. If FALSE model was estimated with constant mean.

• msvar: Boolean. If TRUE model was estimated with switch in variance. If FALSE model was estimated with constant variance.

• control: List with model options used.

• logLike: log-likelihood.

• AIC: Akaike information criterion.

• BIC: Bayesian (Schwarz) information criterion.

• Hess: Hessian matrix. Approximated using hessian and only returned if getSE=TRUE.

• info_mat: Information matrix. Computed as the inverse of -Hess. If matrix is not PD then nearest PD matrix is obtained using nearest_spd. Only returned if getSE=TRUE.

• nearPD_used: Boolean determining whether nearPD function was used on info_mat if TRUE or not if FALSE. Only returned if getSE=TRUE.

• theta_se: standard errors of parameters in theta. Only returned if getSE=TRUE.

• trace: List with Lists of estimation output for each initial value used due to use_diff_init > 1.

References

Dempster, A. P., N. M. Laird, and D. B. Rubin. 1977. “Maximum Likelihood from Incomplete Data via the EM Algorithm.” Journal of the Royal Statistical Society. Series B 39 (1): 1–38..

Böhning, D., E. Dietz, R. Schaub, P. Schlattmann, and B. G. Lindsay. 1994. “The distribution of the likelihood ratio for mixtures of densities from the one-parameter exponential family.” Annals of the Institute of Statistical Mathematics 46 (2): 373–388.

McLachlan, G. J., and T. Krishnan. 2008. The EM Algorithm and Extensions. 2nd ed. Hoboken, New Jersey: John Wiley & Sons.

Hamilton, James D. 1990. “Analysis of time series subject to changes in regime.” Journal of econometrics, 45 (1-2): 39–70.

Basawa, I. V., Mallik, A. K., McCormick, W. P., Reeves, J. H., & Taylor, R. L. 1991. "Bootstrapping Unstable First-Order Autoregressive Processes." The Annals of Statistics 19(2): 1098-1101.

See Also

ARmdl

Examples

# --------------------------- Use simulated process ----------------------------
set.seed(1234)
# Define DGP of MS AR process
mdl_ms2 <- list(n     = 500, 
                mu    = c(5,10),
                sigma = c(1,4),
                phi   = c(0.5),
                k     = 2,
                P     = rbind(c(0.90, 0.10),
                              c(0.10, 0.90)))

# Simulate process using simuMSAR() function
y_ms_simu <- simuMSAR(mdl_ms2)

# Set options for model estimation
control <- list(msmu   = TRUE, 
                msvar  = TRUE, 
                method = "EM",
                use_diff_init = 1)

# Estimate model

  ms_mdl <- MSARmdl(y_ms_simu$y, p = 1, k = 2, control = control)
  summary(ms_mdl)

Description

This function estimates a Markov-switching vector autoregressive model

Usage

MSVARmdl(Y, p, k, control = list())

Arguments

Value

List of class MSVARmdl (S3 object) with model attributes including:

• y: a (T-p x q) matrix of observations.

• X: a (T-p x p*q + const) matrix of lagged observations with a leading column of 1s.

• x: a (T-p x p*q) matrix of lagged observations.

• fitted: a (T x q) matrix of fitted values.

• resid: a (T-p x q) matrix of residuals.

• intercept: a (k x q) matrix of estimated intercepts of each process.

• mu: a (k x q) matrix of estimated means of each process.

• beta: a list containing k separate ((1 + p) x q) matrix of estimated coefficients for each regime.

• phi: estimates of autoregressive coefficients.

• Fmat: Companion matrix containing autoregressive coefficients.

• stdev: List with k (q x q) matrices with estimated standard deviation on the diagonal.

• sigma: List with k (q x q) matrices with estimated covariance matrix.

• theta: vector containing: mu and vech(sigma).

• theta_mu_ind: vector indicating location of mean with 1 and 0 otherwise.

• theta_sig_ind: vector indicating location of variance and covariances with 1 and 0 otherwise.

• theta_var_ind: vector indicating location of variances with 1 and 0 otherwise.

• theta_P_ind: vector indicating location of transition matrix elements with 1 and 0 otherwise.

• stationary: Boolean indicating if process is stationary if TRUE or non-stationary if FALSE.

• n: number of observations (same as T).

• p: number of autoregressive lags.

• q: number of series.

• k: number of regimes in estimated model.

• P: a (k x k) transition matrix.

• pinf: a (k x 1) vector with limiting probabilities of each regime.

• St: a (T x k) vector with smoothed probabilities of each regime at each time t.

• deltath: double with the maximum relative change in vector theta on the last iteration.

• iterations: number of EM iterations performed to achieve convergence (if less than maxit).

• converged: Boolean. TRUE if the conv convergence criterion was met before maxit iterations, FALSE otherwise.

• theta_0: vector of initial values used.

• init_used: number of different initial values used to get a finite solution. See description of input maxit_converge.

• msmu: Boolean. If TRUE model was estimated with switch in mean. If FALSE model was estimated with constant mean.

• msvar: Boolean. If TRUE model was estimated with switch in variance. If FALSE model was estimated with constant variance.

• control: List with model options used.

• logLike: log-likelihood.

• AIC: Akaike information criterion.

• BIC: Bayesian (Schwarz) information criterion.

• Hess: Hessian matrix. Approximated using hessian and only returned if getSE=TRUE.

• info_mat: Information matrix. Computed as the inverse of -Hess. If matrix is not PD then nearest PD matrix is obtained using nearest_spd. Only returned if getSE=TRUE.

• nearPD_used: Boolean determining whether nearPD function was used on info_mat if TRUE or not if FALSE. Only returned if getSE=TRUE.

• theta_se: standard errors of parameters in theta. Only returned if getSE=TRUE.

• trace: List with Lists of estimation output for each initial value used due to use_diff_init > 1.

References

Dempster, A. P., N. M. Laird, and D. B. Rubin. 1977. “Maximum Likelihood from Incomplete Data via the EM Algorithm.” Journal of the Royal Statistical Society. Series B 39 (1): 1–38..

Böhning, D., E. Dietz, R. Schaub, P. Schlattmann, and B. G. Lindsay. 1994. “The distribution of the likelihood ratio for mixtures of densities from the one-parameter exponential family.” Annals of the Institute of Statistical Mathematics 46 (2): 373–388.

McLachlan, G. J., and T. Krishnan. 2008. The EM Algorithm and Extensions. 2nd ed. Hoboken, New Jersey: John Wiley & Sons.

Krolzig, Hans-Martin. 1997. “The markov-switching vector autoregressive model.”. Springer.

See Also

VARmdl

Examples

set.seed(123)
# Define DGP of MS VAR process
mdl_msvar2 <- list(n     = 200, 
                   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.90, 0.10),
                                 c(0.10, 0.90)))

# Simulate process using simuMSVAR() function
y_msvar_simu <- simuMSVAR(mdl_msvar2)

# Set options for model estimation
control <- list(msmu   = TRUE, 
                msvar  = TRUE,
                method = "EM",
                use_diff_init = 1)
                
# Estimate model

  y_msvar_mdl <- MSVARmdl(y_msvar_simu$y, p = 1, k = 2, control = control)
  summary(y_msvar_mdl)

Description

This function estimates a Markov-switching vector autoregressive model with exogenous regressors.

Usage

MSVARXmdl(Y, p, k, Z, control = list())

Arguments

Value

List of class MSVARXmdl (S3 object) with model attributes including:

• y: a (T-p x q) matrix of observations.

• X: a (T-p x p*q + const) matrix of lagged observations with a leading column of 1s.

• x: a (T-p x p*q) matrix of lagged observations.

• resid: a (T-p x q) matrix of residuals.

• fitted: a (T x q) matrix of fitted values.

• intercept: a (k x q) matrix of estimated intercepts of each process.

• mu: a (k x q) matrix of estimated means of each process.

• beta: a list containing k separate ((1 + p + qz) x q) matrix of estimated coefficients for each regime.

• betaZ: a (qz x q) matrix of estimated exogenous regressor coefficients.

• phi: estimates of autoregressive coefficients.

• Fmat: Companion matrix containing autoregressive coefficients.

• stdev: List with k (q x q) matrices with estimated standard deviation on the diagonal.

• sigma: List with k (q x q) matrices with estimated covariance matrix.

• theta: vector containing: mu and vech(sigma).

• theta_mu_ind: vector indicating location of mean with 1 and 0 otherwise.

• theta_sig_ind: vector indicating location of variance and covariances with 1 and 0 otherwise.

• theta_var_ind: vector indicating location of variances with 1 and 0 otherwise.

• theta_P_ind: vector indicating location of transition matrix elements with 1 and 0 otherwise.

• stationary: Boolean indicating if process is stationary if TRUE or non-stationary if FALSE.

• n: number of observations (same as T).

• p: number of autoregressive lags.

• q: number of series.

• k: number of regimes in estimated model.

• P: a (k x k) transition matrix.

• pinf: a (k x 1) vector with limiting probabilities of each regime.

• St: a (T x k) vector with smoothed probabilities of each regime at each time t.

• deltath: double with the maximum relative change in vector theta on the last iteration.

• iterations: number of EM iterations performed to achieve convergence (if less than maxit).

• converged: Boolean. TRUE if the conv convergence criterion was met before maxit iterations, FALSE otherwise.

• theta_0: vector of initial values used.

• init_used: number of different initial values used to get a finite solution. See description of input maxit_converge.

• msmu: Boolean. If TRUE model was estimated with switch in mean. If FALSE model was estimated with constant mean.

• msvar: Boolean. If TRUE model was estimated with switch in variance. If FALSE model was estimated with constant variance.

• control: List with model options used.

• logLike: log-likelihood.

• AIC: Akaike information criterion.

• BIC: Bayesian (Schwarz) information criterion.

• Hess: Hessian matrix. Approximated using hessian and only returned if getSE=TRUE.

• info_mat: Information matrix. Computed as the inverse of -Hess. If matrix is not PD then nearest PD matrix is obtained using nearest_spd. Only returned if getSE=TRUE.

• nearPD_used: Boolean determining whether nearPD function was used on info_mat if TRUE or not if FALSE. Only returned if getSE=TRUE.

• theta_se: standard errors of parameters in theta. Only returned if getSE=TRUE.

• trace: List with Lists of estimation output for each initial value used due to use_diff_init > 1.

References

Dempster, A. P., N. M. Laird, and D. B. Rubin. 1977. “Maximum Likelihood from Incomplete Data via the EM Algorithm.” Journal of the Royal Statistical Society. Series B 39 (1): 1–38..

Böhning, D., E. Dietz, R. Schaub, P. Schlattmann, and B. G. Lindsay. 1994. “The distribution of the likelihood ratio for mixtures of densities from the one-parameter exponential family.” Annals of the Institute of Statistical Mathematics 46 (2): 373–388.

McLachlan, G. J., and T. Krishnan. 2008. The EM Algorithm and Extensions. 2nd ed. Hoboken, New Jersey: John Wiley & Sons.

Krolzig, Hans-Martin. 1997. “The markov-switching vector autoregressive model.”. Springer.

See Also

VARmdl

Examples

set.seed(123)
# Define DGP of MS VAR process
mdl_msvar2 <- list(n     = 200, 
                   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.90, 0.10),
                                 c(0.10, 0.90)))

# Simulate process using simuMSVAR() function
y_msvar_simu <- simuMSVAR(mdl_msvar2)

# Set options for model estimation
control <- list(msmu   = TRUE, 
                msvar  = TRUE,
                method = "EM",
                use_diff_init = 1)
                
# Estimate model

  y_msvar_mdl <- MSVARmdl(y_msvar_simu$y, p = 1, k = 2, control = control)
  summary(y_msvar_mdl)

Description

This function estimates a univariate or multivariate normally distributed model. This can be used for the null hypothesis of a linear model against an alternative hypothesis of a HMM with k regimes.

Usage

Nmdl(Y, Z = NULL, control = list())

Arguments

Value

List of class Nmdl (S3 object) with model attributes including:

• y: a (T x q) matrix of observations.

• fitted: a (T x q) matrix of fitted values.

• resid: a (T x q) matrix of residuals.

• mu: a (1 x q) vector of estimated means of each process.

• intercept: a (1 x q) vector of estimated intercept of each process. If Z is NULL, it is the same as mu.

• beta: a ((1 + p + qz) x q) matrix of estimated coefficients.

• betaZ: a (qz x q) matrix of estimated exogenous regressor coefficients (If Z is provided).

• stdev: a (q x 1) vector of estimated standard deviation of each process.

• sigma: a (q x q) estimated covariance matrix.

• theta: vector containing: mu, betaZ (if matrix Z is provided), and vech(sigma).

• theta_mu_ind: vector indicating location of mean with 1 and 0 otherwise.

• theta_sig_ind: vector indicating location of variance and covariances with 1 and 0 otherwise.

• theta_var_ind: vector indicating location of variances with 1 and 0 otherwise.

• n: number of observations (same as T).

• q: number of series.

• k: number of regimes. This is always 1 in Nmdl.

• control: List with model options used.

• logLike: log-likelihood.

• AIC: Akaike information criterion.

• BIC: Bayesian (Schwarz) information criterion.

• Hess: Hessian matrix. Approximated using hessian and only returned if getSE=TRUE.

• info_mat: Information matrix. Computed as the inverse of -Hess. If matrix is not PD then nearest PD matrix is obtained using nearest_spd. Only returned if getSE=TRUE.

• nearPD_used: Boolean determining whether nearPD function was used on info_mat if TRUE or not if FALSE. Only returned if getSE=TRUE.

• theta_se: standard errors of parameters in theta. Only returned if getSE=TRUE.

Examples

set.seed(1234)

# ----- Univariate ----- # 
# Define DGP 
mdl_norm <- list(n     = 1000, 
                 q     = 1,
                 mu    = as.matrix(5),
                 sigma = as.matrix(5.0))

# Simulate process using simuNorm() function
y_norm_simu <- simuNorm(mdl_norm)

# estimate parameters
y_norm_mdl <- Nmdl(y_norm_simu$y)
summary(y_norm_mdl)


# ----- Multivariate ----- # 
# Define DGP 
mdl_norm <- list(n     = 1000, 
                 q     = 2,
                 mu    = c(5, -2),
                 sigma = rbind(c(5.0, 1.5),
                               c(1.5, 1.0)))

# Simulate process using simuNorm() function
y_norm_simu <- simuNorm(mdl_norm)

# estimate parameters
y_norm_mdl <- Nmdl(y_norm_simu$y)
summary(y_norm_mdl)

Description

This function simulates an autoregresive process.

Usage

simuAR(mdl_h0, burnin = 100)

Arguments

Value

List with simulated autoregressive series and its DGP parameters.

Examples

set.seed(1234)
# Define DGP of AR process
mdl_ar <- list(n     = 500, 
               mu    = 5,
               sigma = 2,
               phi   = c(0.5,0.2))

# Simulate process using simuAR() function
y_simu <- simuAR(mdl_ar)

plot(y_simu)

Description

This function simulates an autoregresive X process.

Usage

simuARX(mdl_h0, burnin = 100)

Arguments

Value

List with simulated autoregressive series and its DGP parameters.

Examples

set.seed(1234)
# Define DGP of AR process
mdl_ar <- list(n     = 500, 
               mu    = 5,
               sigma = 2,
               phi   = c(0.5,0.2))

# Simulate process using simuAR() function
y_simu <- simuAR(mdl_ar)

plot(y_simu)

Description

This function simulates a Hidden Markov Model process.

Usage

simuHMM(mdl_h0, burnin = 100)

Arguments

Value

List with simulated series and its DGP parameters.

Examples

set.seed(1234)

# ----- Univariate ----- #
# Define DGP 
mdl_hmm <- list(n     = 1000, 
                q     = 1,
                mu    = as.matrix(c(5,
                                   -2)),
                sigma = list(as.matrix(5.0),
                             as.matrix(7.0)),
                k     = 2,
                P     = rbind(c(0.90, 0.10),
                              c(0.10, 0.90)))

# Simulate process using simuHMM() function
y_hmm_simu <- simuHMM(mdl_hmm)

plot(y_hmm_simu)

# ----- Multivariate ----- #
# Define DGP 
mdl_hmm <- list(n     = 1000, 
                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.90, 0.10),
                              c(0.10, 0.90)))

# Simulate process using simuHMM() function
y_hmm_simu <- simuHMM(mdl_hmm)


plot(y_hmm_simu)

Description

This function simulates a Markov-switching autoregressive process.

Usage

simuMSAR(mdl_h0, burnin = 100)

Arguments

Value

List with simulated Markov-switching autoregressive process and its DGP properties.

Examples

set.seed(1234)
# Define DGP of MS AR process
mdl_ms2 <- list(n     = 500, 
                mu    = c(5,10),
                sigma = c(1,2),
                phi   = c(0.5, 0.2),
                k     = 2,
                P     = rbind(c(0.90, 0.10),
                              c(0.10, 0.90)))

# Simulate process using simuMSAR() function
y_ms_simu <- simuMSAR(mdl_ms2)

plot(y_ms_simu)

Description

This function simulates a Markov-switching autoregressive process.

Usage

simuMSARX(mdl_h0, burnin = 100)

Arguments

Value

List with simulated Markov-switching autoregressive process and its DGP properties.

Examples

set.seed(1234)
# Define DGP of MS AR process
mdl_ms2 <- list(n     = 500, 
                mu    = c(5,10),
                sigma = c(1,2),
                phi   = c(0.5, 0.2),
                k     = 2,
                P     = rbind(c(0.90, 0.10),
                              c(0.10, 0.90)))

# Simulate process using simuMSAR() function
y_ms_simu <- simuMSAR(mdl_ms2)

plot(y_ms_simu)

Description

This function simulates a Markov-switching vector autoregressive process.

Usage

simuMSVAR(mdl_h0, burnin = 100)

Arguments

Value

List with simulated vector autoregressive series and its DGP parameters.

Examples

set.seed(1234)
# Define DGP of MS VAR process
mdl_msvar2 <- list(n     = 1000, 
                   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.90, 0.10),
                                 c(0.10, 0.90)))

# Simulate process using simuMSVAR() function
y_msvar_simu <- simuMSVAR(mdl_msvar2)

plot(y_msvar_simu)

Description

This function simulates a Markov-switching VARX process.

Usage

simuMSVARX(mdl_h0, burnin = 100)

Arguments

Value

List with simulated vector autoregressive series and its DGP parameters.

Examples

set.seed(1234)
# Define DGP of MS VAR process
mdl_msvar2 <- list(n     = 1000, 
                   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.90, 0.10),
                                 c(0.10, 0.90)))

# Simulate process using simuMSVAR() function
y_msvar_simu <- simuMSVAR(mdl_msvar2)

plot(y_msvar_simu)

Description

This function simulates a normally distributed process.

Usage

simuNorm(mdl_h0, burnin = 0)

Arguments

Value

List with simulated series and its DGP parameters.

Examples

set.seed(1234)
# Define DGP 
mdl_norm <- list(n     = 1000, 
                 q     = 2,
                 mu    = c(5, -2),
                 sigma = rbind(c(5.0, 1.5),
                               c(1.5, 1.0)))

# Simulate process using simuNorm() function
y_norm_simu <- simuNorm(mdl_norm)

plot(y_norm_simu)

Description

This function simulates a vector autoregresive process.

Usage

simuVAR(mdl_h0, burnin = 100)

Arguments

Value

List with simulated vector autoregressive series and its DGP parameters.

Examples

set.seed(1234)
# Define DGP of VAR process
mdl_3var2 <- list(n     = 1000, 
                  p     = 2,
                  q     = 3,
                  mu    = c(5, -2, 1),
                  sigma = rbind(c(5.0, 1.5, 2.5),
                                c(1.5, 1.0, 1.5),
                                c(2.5, 1.5, 4.2)),
                  phi   = rbind(c(0.70, 0.30, 0.35,  -0.50, -0.20,   0.25),
                                c(0.20, 0.40, 0.35,  -0.30,  0.30,   0.25),
                                c(0.20, 0.30, 0.50,  -0.30, -0.20,  -0.40)))

# Simulate process using simuVAR() function
y3var2_simu <- simuVAR(mdl_3var2)

plot(y3var2_simu)

Description

This function simulates a vector autoregressive process with exogenous regressors.

Usage

simuVARX(mdl_h0, burnin = 100)

Arguments

Value

List with simulated vector autoregressive series and its DGP parameters.

Examples

set.seed(1234)
# Define DGP of VAR process
mdl_3var2 <- list(n     = 1000, 
                  p     = 2,
                  q     = 3,
                  mu    = c(5, -2, 1),
                  sigma = rbind(c(5.0, 1.5, 2.5),
                                c(1.5, 1.0, 1.5),
                                c(2.5, 1.5, 4.2)),
                  phi   = rbind(c(0.70, 0.30, 0.35,  -0.50, -0.20,   0.25),
                                c(0.20, 0.40, 0.35,  -0.30,  0.30,   0.25),
                                c(0.20, 0.30, 0.50,  -0.30, -0.20,  -0.40)))

# Simulate process using simuVAR() function
y3var2_simu <- simuVAR(mdl_3var2)

plot(y3var2_simu)

Description

This data is used in Rodriguez-Rondon & Dufour (2026a). The series ranges from 1947Q2 to 2024Q2.

Usage

USGNP

Format

Date

Vector of dates

GNP

US GNP series

GNP_gr

log difference of US GNP series

Source

https://fred.stlouisfed.org/graph/?g=UKDQ

References

U.S. Bureau of Economic Analysis, Gross National Product [GNP], retrieved from FRED, Federal Reserve Bank of St. Louis; https://fred.stlouisfed.org/series/GNP, September, 2024.

Rodriguez-Rondon, G., & Dufour, J.-M. 2026a. "Monte Carlo Likelihood-Ratio Tests for Markov Switching Models." Bank of Canada Staff Working Paper, No. 2026-23. doi: 10.34989/swp-2026-23.

Description

This data is used in Rodriguez-Rondon & Dufour (2026a). The series ranges from 1947Q2 to 2024Q2.

Usage

USRGDP

Format

Date

Vector of dates

RGDP

US Real GDP series

RGDP_gr

log difference of US Real GDP series

Source

https://fred.stlouisfed.org/series/GDPC1

References

U.S. Bureau of Economic Analysis, Real Gross Domestic Product [GDPC1], retrieved from FRED, Federal Reserve Bank of St. Louis; https://fred.stlouisfed.org/series/GDPC1, September, 2024.

Rodriguez-Rondon, G., & Dufour, J.-M. 2026a. "Monte Carlo Likelihood-Ratio Tests for Markov Switching Models." Bank of Canada Staff Working Paper, No. 2026-23. doi: 10.34989/swp-2026-23.

Description

This function estimates a vector autoregresive model with p lags. This can be used for the null hypothesis of a linear model against an alternative hypothesis of a Markov switching vector autoregressive model with k regimes.

Usage

VARmdl(Y, p, control = list())

Arguments

Value

List of class VARmdl (S3 object) with model attributes including:

• y: a (T-p x q) matrix of observations.

• X: a (T-p x p*q + const) matrix of lagged observations with a leading column of 1s if const=TRUE or not if const=FALSE.

• x: a (T-p x p*q) matrix of lagged observations.

• fitted: a (T-p x q) matrix of fitted values.

• resid: a (T-p x q) matrix of residuals.

• mu: a (1 x q) vector of estimated means of each process.

• beta: a ((1 + p) x q) matrix of estimated coefficients. .

• intercept: estimate of intercepts.

• phi: a (q x p*q) matrix of estimated autoregressive coefficients.

• Fmat: Companion matrix containing autoregressive coefficients.

• stdev: a (q x 1) vector of estimated standard deviation of each process.

• sigma: a (q x q) estimated covariance matrix.

• theta: vector containing: mu, vech(sigma), and vec(t(phi)).

• theta_mu_ind: vector indicating location of mean with 1 and 0 otherwise.

• theta_sig_ind: vector indicating location of variance and covariances with 1 and 0 otherwise.

• theta_var_ind: vector indicating location of variances with 1 and 0 otherwise.

• theta_phi_ind: vector indicating location of autoregressive coefficients with 1 and 0 otherwise.

• stationary: Boolean indicating if process is stationary if TRUE or non-stationary if FALSE.

• n: number of observations after lag transformation (i.e., n = T-p).

• p: number of autoregressive lags.

• q: number of series.

• k: number of regimes. This is always 1 in VARmdl.

• Fmat: matrix from companion form. Used to determine is process is stationary.

• control: List with model options used.

• logLike: log-likelihood.

• AIC: Akaike information criterion.

• BIC: Bayesian (Schwarz) information criterion.

• Hess: Hessian matrix. Approximated using hessian and only returned if getSE=TRUE.

• info_mat: Information matrix. Computed as the inverse of -Hess. If matrix is not PD then nearest PD matrix is obtained using nearest_spd. Only returned if getSE=TRUE.

• nearPD_used: Boolean determining whether nearPD function was used on info_mat if TRUE or not if FALSE. Only returned if getSE=TRUE.

• theta_se: standard errors of parameters in theta. Only returned if getSE=TRUE.

See Also

MSVARmdl

Examples

# ----- Bivariate VAR(1) process ----- #
set.seed(1234)
# Define DGP of VAR process
mdl_var <- list(n     = 1000, 
                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)))

# Simulate process using simuVAR() function
y_simu <- simuVAR(mdl_var)

# Set options for model estimation
control <- list(const  = TRUE, 
                getSE  = TRUE)

# Estimate model
y_var_mdl <- VARmdl(y_simu$y, p = 2, control = control)
summary(y_var_mdl)

Description

This function estimates a vector autoregresive model with p lags. This can be used for the null hypothesis of a linear model against an alternative hypothesis of a Markov switching vector autoregressive model with k regimes.

Usage

VARXmdl(Y, p, Z, control = list())

Arguments

Value

List of class VARmdl (S3 object) with model attributes including:

• y: a (T-p x q) matrix of observations.

• X: a (T-p x p*q + const) matrix of lagged observations with a leading column of 1s if const=TRUE or not if const=FALSE.

• x: a (T-p x p*q) matrix of lagged observations.

• fitted: a (T-p x q) matrix of fitted values.

• resid: a (T-p x q) matrix of residuals.

• mu: a (1 x q) vector of estimated means of each process.

• beta: a ((1 + p + qz) x q) matrix of estimated coefficients.

• betaZ: a (qz x q) matrix of estimated exogenous regressor coefficients.

• intercept: estimate of intercepts.

• phi: a (q x p*q) matrix of estimated autoregressive coefficients.

• Fmat: Companion matrix containing autoregressive coefficients.

• stdev: a (q x 1) vector of estimated standard deviation of each process.

• sigma: a (q x q) estimated covariance matrix.

• theta: vector containing: mu, vech(sigma), and vec(t(phi)).

• theta_mu_ind: vector indicating location of mean with 1 and 0 otherwise.

• theta_sig_ind: vector indicating location of variance and covariances with 1 and 0 otherwise.

• theta_var_ind: vector indicating location of variances with 1 and 0 otherwise.

• theta_phi_ind: vector indicating location of autoregressive coefficients with 1 and 0 otherwise.

• stationary: Boolean indicating if process is stationary if TRUE or non-stationary if FALSE.

• n: number of observations after lag transformation (i.e., n = T-p).

• p: number of autoregressive lags.

• q: number of series.

• k: number of regimes. This is always 1 in VARmdl.

• Fmat: matrix from companion form. Used to determine is process is stationary.

• control: List with model options used.

• logLike: log-likelihood.

• AIC: Akaike information criterion.

• BIC: Bayesian (Schwarz) information criterion.

• Hess: Hessian matrix. Approximated using hessian and only returned if getSE=TRUE.

• info_mat: Information matrix. Computed as the inverse of -Hess. If matrix is not PD then nearest PD matrix is obtained using nearest_spd. Only returned if getSE=TRUE.

• nearPD_used: Boolean determining whether nearPD function was used on info_mat if TRUE or not if FALSE. Only returned if getSE=TRUE.

• theta_se: standard errors of parameters in theta. Only returned if getSE=TRUE.

See Also

MSVARmdl

Examples

# ----- Bivariate VAR(1) process ----- #
set.seed(1234)
# Define DGP of VAR process
mdl_var <- list(n     = 1000, 
                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)))

# Simulate process using simuVAR() function
y_simu <- simuVAR(mdl_var)

# Set options for model estimation
control <- list(const  = TRUE, 
                getSE  = TRUE)

# Estimate model
y_var_mdl <- VARmdl(y_simu$y, p = 2, control = control)
summary(y_var_mdl)

Description

Builds deterministic starting values for the k1-regime (alternative) model by embedding a fitted k0-regime (null) model into the larger parameter space. Each embedding duplicates one or more regimes of the null fit (with the transition-matrix column split so that the likelihood at the exact embedding equals the null model's log-likelihood) and then applies a small deterministic antisymmetric perturbation to the duplicated regimes (means shifted by +/- c_pert regime standard deviations when the mean switches; variances scaled by (1 + c_pert)^{+/-1} when the variance switches) so that the EM algorithm can leave the symmetric fixed point. One starting vector is returned per way of distributing the k1 - k0 extra regimes among the k0 regimes of the null fit. These starts are used in addition to the random starting values when init_method = "warmstart" in LMCLRTest and MMCLRTest.

Usage

warmstart_theta(mdl_h0, k1, msmu = NULL, msvar = NULL, c_pert = 0.5)

Arguments

Value

List of numeric vectors, each a valid starting value for the k1-regime model (in the parameter ordering used by the corresponding model constructor). Starting vectors containing non-finite values are dropped; the list may be empty if the null fit is degenerate.

References

Kasahara, H. & Shimotsu, K. 2018. "Testing the number of regimes in Markov regime switching models." arXiv:1801.06862.

Dufour, J.-M. 2006. "Monte Carlo tests with nuisance parameters: A general approach to finite-sample inference and nonstandard asymptotics." Journal of Econometrics, 133(2), 443-477.

Rodriguez-Rondon, G. & Dufour, J.-M. 2026. "Monte Carlo Likelihood-Ratio Tests for Markov Switching Models." Bank of Canada Staff Working Paper 2026-23.