| Title: | Hypothesis Testing for Markov Switching Models |
|---|---|
| Description: | Implementation of hypothesis testing procedures described in Hansen (1992) <doi:10.1002/jae.3950070506>, Carrasco, Hu, & Ploberger (2014) <doi:10.3982/ECTA8609>, Dufour & Luger (2017) <doi:10.1080/07474938.2017.1307548>, and Rodriguez Rondon & Dufour (2022) <https://grodriguezrondon.com/files/RodriguezRondon_Dufour_2024_MonteCarlo_LikelihoodRatioTest_MarkovSwitchingModels_20241015.pdf> that can be used to identify the number of regimes in Markov switching models. |
| Authors: | Gabriel Rodriguez Rondon [cre, aut], Jean-Marie Dufour [aut] |
| Maintainer: | Gabriel Rodriguez Rondon <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 0.1.3 |
| Built: | 2024-10-31 12:29:48 UTC |
| Source: | CRAN |
This package implements hypothesis testing procedures that can be used to identify the number of regimes in a Markov-Switching model.
Gabriel Rodriguez Rondon, [email protected] (Maintainer)
Jean-Marie Dufour, [email protected]
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, Gabriel and Jean-Marie Dufour. 2022. "Simulation-Based Inference for Markov Switching Models” JSM Proceedings, Business and Economic Statistics Section: American Statistical Association.
Rodriguez Rondon, Gabriel and Jean-Marie Dufour. 2022. “Monte Carlo Likelihood Ratio Tests for Markov Switching Models.” Unpublished manuscript.
Rodriguez Rondon, Gabriel and Jean-Marie Dufour. 2022. “MSTest: An R-package for Testing Markov-Switching Models.” Unpublished manuscript.
Qu, Zhongjun, and Fan Zhuo. 2021. “Likelihood Ratio-Based Tests for Markov Regime Switching.” The Review of Economic Studies 88 (2): 937–968.
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.
ARmdl(Y, p, control = list())ARmdl(Y, p, control = list())
Y |
A |
p |
Integer determining the number of autoregressive lags. |
control |
List with model options including:
|
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.
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_mdlset.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
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.
ARXmdl(Y, p, Z, control = list())ARXmdl(Y, p, Z, control = list())
Y |
A |
p |
Integer determining the number of autoregressive lags. |
Z |
A |
control |
List with model options including:
|
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.
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_mdlset.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
Carrasco, Hu, & Ploberger 2010 GNP data
chp10GNPchp10GNP
This data is the extension of the GNP series used in CHP (2014), Econometrica. This series ranges from 1951Q2 to 2010Q4.
Vector of dates
US GNP series
log difference of US GNP series
https://www.econometricsociety.org/content/supplement-optimal-test-markov-switching-parameters
Carrasco, Marine, Liang Hu, and Werner Ploberger. 2014. “Optimal test for Markov switching parameters.” Econometrica 82 (2): 765–784.
This bootstrap procedure is described on pg. 771 of CHP 2014.
CHPbootCV(mdl, rho_b, N, msvar)CHPbootCV(mdl, rho_b, N, msvar)
mdl |
List containing model attributes (see |
rho_b |
Number determining bounds for distribution of |
N |
Number of bootstrap simulations. |
msvar |
Boolean indicator. If |
Bootstrap critical values
Carrasco, Marine, Liang Hu, and Werner Ploberger. 2014. “Optimal test for Markov switching parameters.” Econometrica 82 (2): 765–784.
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.
CHPTest(Y, p, control = list())CHPTest(Y, p, control = list())
Y |
A ( |
p |
Integer determining the number of autoregressive lags. |
control |
List with test procedure options including:
|
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.
Carrasco, Marine, Liang Hu, and Werner Ploberger. 2014. “Optimal test for Markov switching parameters.” Econometrica 82 (2): 765–784.
# 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)# 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)
This function performs the Local Monte Carlo moment-based test for Markov switching autoregressive models proposed in Dufour & Luger (2017).
DLMCTest(Y, p, control = list())DLMCTest(Y, p, control = list())
Y |
Series to be tested |
p |
Number of autoregressive lags. |
control |
List with test procedure options including:
|
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.
Dufour, J. M., & Luger, R. 2017. "Identification-robust moment-based tests for Markov switching in autoregressive models." Econometric Reviews, 36(6-9), 713-727.
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)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)
This function performs the maximized Monte Carlo moment-based test for Markov switching autoregressive models proposed in Dufour & Luger (2017).
DLMMCTest(Y, p, control = list())DLMMCTest(Y, p, control = list())
Y |
Series to be tested |
p |
Number of autoregressive lags. |
control |
List with test procedure options including:
|
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.
Dufour, J. M., & Luger, R. 2017. "Identification-robust moment-based tests for Markov switching in autoregressive models." Econometric Reviews, 36(6-9), 713-727.
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)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)
Hamilton 1984 & Hansen 1992 GNP data
hamilton84GNPhamilton84GNP
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.
Vector of dates
US GNP series
US GNP log difference
https://www.ssc.wisc.edu/~bhansen/progs/jae_92.html
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.
This function performs Hansen's likelihood ratio test as described in Hansen (1992). Original source code can be found here.
HLRTest(Y, p, control = list())HLRTest(Y, p, control = list())
Y |
A ( |
p |
Integer determining the number of autoregressive lags. |
control |
List with test procedure options including:
|
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.
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.
# --------------------------- 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)# --------------------------- 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)
This function estimates a Hidden Markov model with k regimes.
HMmdl(Y, k, Z = NULL, control = list())HMmdl(Y, k, Z = NULL, control = list())
Y |
a |
k |
integer determining the number of regimes to use in estimation. Must be greater than or equal to |
Z |
an otpional |
control |
List with model options including:
|
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 maximum absolute difference in vector theta between last iteration.
iterations: number of EM iterations performed to achieve convergence (if less than maxit).
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.
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..
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.
This function performs the Local Monte Carlo likelihood ratio test (LMC-LRT) proposed in Rodriguez-Rondon & Dufour (2024). 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.
LMCLRTest(Y, p, k0, k1, Z = NULL, control = list())LMCLRTest(Y, p, k0, k1, Z = NULL, control = list())
Y |
Series to be tested. Must be a ( |
p |
Number of autoregressive lags. Must be greater than or equal to 0. |
k0 |
Number of regimes under null hypothesis. Must be greater than or equal to 1. |
k1 |
Number of regimes under alternative hypothesis. Must be greater than |
Z |
Exogenous regressors. Optional input and default is NULL. When used, it should be a ( |
control |
List with test procedure options including:
|
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.
Rodriguez-Rondon, Gabriel and Jean-Marie Dufour. 2022. "Simulation-Based Inference for Markov Switching Models” JSM Proceedings, Business and Economic Statistics Section: American Statistical Association.
Rodriguez-Rondon, Gabriel and Jean-Marie Dufour. 2024. “Monte Carlo Likelihood Ratio Tests for Markov Switching Models.” Unpublished manuscript.
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 = TRUE), mdl_h1_control = list(msmu = TRUE, msvar = TRUE, getSE = TRUE, method = "EM", use_diff_init = 1)) lmctest <- LMCLRTest(y_ms_simu$y, p = 1, k0 = 1 , k1 = 2, control = lmc_control) summary(lmctest)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 = TRUE), mdl_h1_control = list(msmu = TRUE, msvar = TRUE, getSE = TRUE, method = "EM", use_diff_init = 1)) lmctest <- LMCLRTest(y_ms_simu$y, p = 1, k0 = 1 , k1 = 2, control = lmc_control) summary(lmctest)
This function computes the Monte Carlo P-value.
MCpval(test_stat, null_vec, type = "geq")MCpval(test_stat, null_vec, type = "geq")
test_stat |
Test statistic under the alternative (e.g. |
null_vec |
A ( |
type |
String determining type of test. options are: "geq" for right-tail test, "leq" for left-tail test, "abs" for absolute value test and "two-tail" for two-tail test. |
MC p-value of test
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.
This function performs the Maximized Monte Carlo likelihood ratio test (MMC-LRT) proposed in Rodriguez-Rondon & Dufour (2024).
MMCLRTest(Y, p, k0, k1, Z = NULL, control = list())MMCLRTest(Y, p, k0, k1, Z = NULL, control = list())
Y |
Series to be tested. Must be a ( |
p |
Number of autoregressive lags. Must be greater than or equal to 0. |
k0 |
Number of regimes under null hypothesis. Must be greater than or equal to 1. |
k1 |
Number of regimes under alternative hypothesis. Must be greater than |
Z |
Exogenous regressors. Optional input and default is NULL. When used, it should be a ( |
control |
List with test procedure options including:
|
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.
Rodriguez-Rondon, Gabriel and Jean-Marie Dufour. 2022. "Simulation-Based Inference for Markov Switching Models” JSM Proceedings, Business and Economic Statistics Section: American Statistical Association.
Rodriguez-Rondon, Gabriel and Jean-Marie Dufour. 2024. “Monte Carlo Likelihood Ratio Tests for Markov Switching Models.” Unpublished manuscript.
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 = 100, 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)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 = 100, 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)
This function estimates a Markov-switching autoregressive model
MSARmdl(Y, p, k, control = list())MSARmdl(Y, p, k, control = list())
Y |
(T x 1) vector with observational data. |
p |
integer for the number of lags to use in estimation. Must be greater than or equal to |
k |
integer for the number of regimes to use in estimation. Must be greater than or equal to |
control |
List with model options including:
|
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 maximum absolute difference in vector theta between last iteration.
iterations: number of EM iterations performed to achieve convergence (if less than maxit).
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.
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..
Hamilton, James D. 1990. “Analysis of time series subject to changes in regime.” Journal of econometrics, 45 (1-2): 39–70.
# --------------------------- Use simulated process ---------------------------- 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 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)# --------------------------- Use simulated process ---------------------------- 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 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)
This function estimates a Markov-switching autoregressive model
MSARXmdl(Y, p, k, Z, control = list())MSARXmdl(Y, p, k, Z, control = list())
Y |
a |
p |
integer for the number of lags to use in estimation. Must be greater than or equal to |
k |
integer for the number of regimes to use in estimation. Must be greater than or equal to |
Z |
a |
control |
List with model options including:
|
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 + 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 maximum absolute difference in vector theta between last iteration.
iterations: number of EM iterations performed to achieve convergence (if less than maxit).
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.
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..
Hamilton, James D. 1990. “Analysis of time series subject to changes in regime.” Journal of econometrics, 45 (1-2): 39–70.
# --------------------------- Use simulated process ---------------------------- 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 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)# --------------------------- Use simulated process ---------------------------- 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 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)
This function estimates a Markov-switching vector autoregressive model
MSVARmdl(Y, p, k, control = list())MSVARmdl(Y, p, k, control = list())
Y |
( |
p |
integer for the number of lags to use in estimation. Must be greater than or equal to |
k |
integer for the number of regimes to use in estimation. Must be greater than or equal to |
control |
List with optimization options including:
|
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 maximum absolute difference in vector theta between last iteration.
iterations: number of EM iterations performed to achieve convergence (if less than maxit).
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.
List with model characteristics
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..
Krolzig, Hans-Martin. 1997. “The markov-switching vector autoregressive model.”. Springer.
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)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)
This function estimates a Markov-switching vector autoregressive model
MSVARXmdl(Y, p, k, Z, control = list())MSVARXmdl(Y, p, k, Z, control = list())
Y |
( |
p |
integer for the number of lags to use in estimation. Must be greater than or equal to |
k |
integer for the number of regimes to use in estimation. Must be greater than or equal to |
Z |
a |
control |
List with optimization options including:
|
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.
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 maximum absolute difference in vector theta between last iteration.
iterations: number of EM iterations performed to achieve convergence (if less than maxit).
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.
List with model characteristics
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..
Krolzig, Hans-Martin. 1997. “The markov-switching vector autoregressive model.”. Springer.
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)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)
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.
Nmdl(Y, Z = NULL, control = list())Nmdl(Y, Z = NULL, control = list())
Y |
a |
Z |
an otpional |
control |
List with model options including:
|
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.
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)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)
This function simulates an autoregresive process.
simuAR(mdl_h0, burnin = 100)simuAR(mdl_h0, burnin = 100)
mdl_h0 |
List containing the following DGP parameters
|
burnin |
Number of simulated observations to remove from beginning. Default is |
List with simulated autoregressive series and its DGP parameters.
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)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)
This function simulates an autoregresive X process.
simuARX(mdl_h0, burnin = 100)simuARX(mdl_h0, burnin = 100)
mdl_h0 |
List containing the following DGP parameters
|
burnin |
Number of simulated observations to remove from beginning. Default is |
List with simulated autoregressive series and its DGP parameters.
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)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)
This function simulates a Hidden Markov Model process.
simuHMM(mdl_h0, burnin = 100)simuHMM(mdl_h0, burnin = 100)
mdl_h0 |
List containing the following DGP parameters
|
burnin |
Number of simulated observations to remove from beginning. Default is |
List with simulated series and its DGP parameters.
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)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)
This function simulates a Markov-switching autoregressive process.
simuMSAR(mdl_h0, burnin = 100)simuMSAR(mdl_h0, burnin = 100)
mdl_h0 |
List containing the following DGP parameters
|
burnin |
Number of simulated observations to remove from beginning. Default is |
List with simulated Markov-switching autoregressive process and its DGP properties.
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)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)
This function simulates a Markov-switching autoregressive process.
simuMSARX(mdl_h0, burnin = 100)simuMSARX(mdl_h0, burnin = 100)
mdl_h0 |
List containing the following DGP parameters
|
burnin |
Number of simulated observations to remove from beginning. Default is |
List with simulated Markov-switching autoregressive process and its DGP properties.
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)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)
This function simulates a Markov-switching vector autoregressive process.
simuMSVAR(mdl_h0, burnin = 100)simuMSVAR(mdl_h0, burnin = 100)
mdl_h0 |
List containing the following DGP parameters
|
burnin |
Number of simulated observations to remove from beginning. Default is |
List with simulated vector autoregressive series and its DGP parameters.
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)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)
This function simulates a Markov-switching VARX process.
simuMSVARX(mdl_h0, burnin = 100)simuMSVARX(mdl_h0, burnin = 100)
mdl_h0 |
List containing the following DGP parameters
|
burnin |
Number of simulated observations to remove from beginning. Default is |
List with simulated vector autoregressive series and its DGP parameters.
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)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)
This function simulates a normally distributed process.
simuNorm(mdl_h0, burnin = 0)simuNorm(mdl_h0, burnin = 0)
mdl_h0 |
List containing the following DGP parameters
|
burnin |
Number of simulated observations to remove from beginning. Default is |
List with simulated series and its DGP parameters.
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)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)
This function simulates a vector autoregresive process.
simuVAR(mdl_h0, burnin = 100)simuVAR(mdl_h0, burnin = 100)
mdl_h0 |
List containing the following DGP parameters
|
burnin |
Number of simulated observations to remove from beginning. Default is |
List with simulated vector autoregressive series and its DGP parameters.
This function simulates a vector autoregresive process.
simuVARX(mdl_h0, burnin = 100)simuVARX(mdl_h0, burnin = 100)
mdl_h0 |
List containing the following DGP parameters
|
burnin |
Number of simulated observations to remove from beginning. Default is |
List with simulated vector autoregressive series and its DGP parameters.
US GNP data 1947Q2 - 2024Q2
USGNPUSGNP
This data is used in Rodriguez-Rondon & Dufour (2024). The series ranges from 1947Q2 to 2024Q2.
Vector of dates
US GNP series
log difference of US GNP series
https://fred.stlouisfed.org/graph/?g=UKDQ
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, Gabriel and Jean-Marie Dufour. 2024. “Monte Carlo Likelihood Ratio Tests for Markov Switching Models.” Unpublished manuscript.
US Real GDP data 1947Q2 - 2024Q2
USRGDPUSRGDP
This data is used in Rodriguez-Rondon & Dufour (2024). The series ranges from 1947Q2 to 2024Q2.
Vector of dates
US Real GDP series
log difference of US Real GDP series
https://fred.stlouisfed.org/series/GDPC1
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, Gabriel and Jean-Marie Dufour. 2024. “Monte Carlo Likelihood Ratio Tests for Markov Switching Models.” Unpublished manuscript.
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.
VARmdl(Y, p, control = list())VARmdl(Y, p, control = list())
Y |
a |
p |
integer determining the number of autoregressive lags. |
control |
List with model options including:
|
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.
# ----- 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)# ----- 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)
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.
VARXmdl(Y, p, Z, control = list())VARXmdl(Y, p, Z, control = list())
Y |
a |
p |
integer determining the number of autoregressive lags. |
Z |
a |
control |
List with model options including:
|
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.
# ----- 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)# ----- 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)