Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
73,238 characters · 14 sections · 70 citation commands
\thispagestyle{empty}
Abstract: { Maximum likelihood estimation of large Markov-switching vector autoregressions (MS-VARs) can be challenging or infeasible due to parameter proliferation. To accommodate situations where dimensionality may be of comparable order to or exceeds the sample size, we adopt a sparse framework and propose two penalized maximum likelihood estimators with either the Lasso or the smoothly clipped absolute deviation (SCAD) penalty. We show that both estimators are estimation consistent, while the SCAD estimator also selects relevant parameters with probability approaching one. A modified EM-algorithm is developed for the case of Gaussian errors and simulations show that the algorithm exhibits desirable finite sample performance. In an application to short-horizon return predictability in the US, we estimate a 15 variable 2-state MS-VAR(1) and obtain the often reported counter-cyclicality in predictability. The variable selection property of our estimators helps to identify predictors that contribute strongly to predictability during economic contractions but are otherwise irrelevant in expansions. Furthermore, out-of-sample analyses indicate that large MS-VARs can significantly outperform "hard-to-beat" predictors like the historical average. }
JEL Classifications: C13, C32, C55, G12
Key words: high-dimensional time series, Markov regime-switching, oracle property, SCAD, stock return predictability.
\setcounter{footnote}{0}
We consider the asymptotic properties of the maximum likelihood estimator (MLE) for high-dimensional Markov-switching (MS) vector autoregressive (VAR) models in a double asymptotic framework where both the sample size and number of parameters are allowed to diverge. Specifically, we study the following:
where $y_t \in \mathbb{R}^{d_T}$ is a vector of endogenous time series, $x_t \in \mathbb{R}^{d^*_T}$ allows for the inclusion of exogenous variables, and $\varepsilon_t$ is independently and identically distributed with mean zero and positive definite variance-covariance matrix $\Sigma_T(S_t)$. The parameters $A_{j,T}(S_t)$, $B_{j,T}(S_t)$, and $\Sigma_T(S_t)$ are respectively $d_T \times d_T$, $d_T \times d_T^*$, and $d_T \times d_T$ matrices that depend on an underlying state $\{S_t\}$, which is modeled as a latent first-order Markov chain on a finite and discrete state space taking values from $\{1 , \ldots, M\}$ with $M$ known and fixed. The transition probabilities of the chain are given by $P(S_t = j| S_{t-1} = i) = p_{i\shortrightarrow j}$, and initial distributions are given by $P(S_1 = j) = p_{j}$. We define the precision matrix to be the inverse of the variance-covariance matrix, $Q_T(S_t) \equiv \Sigma_T(S_t)^{-1}$.
By including the transition matrix and initial probability distributions as estimable parameters, we have a total of $K_T \equiv M(p_T d_T^2 + q_T d_T d_T^* + d_T(d_T+1)/2) + M^2$ parameters. Furthermore, the dimensions of the model, namely $d_T$, $d_T^*$, $p_T$, and $q_T$, are allowed to increase with the sample size subject to regularity conditions. Nonetheless, it is easy to see that, especially with many regimes, such a set-up could easily lead to a high-dimensional problem where the number of parameters exceed the sample size.
To address this, we assume that the model in (ref) is sparse. This means that only a small subset of parameters, relative to the sample size, are truly non-zero, while the rest are regarded as irrelevant. Effectively, this means that the 'true' model with only the relevant parameters is a low-dimension system. A key question is thus whether we are able to asymptotically recover the 'true' model given the full data without a priori knowledge of the sparsity pattern and the latent states.
Low-dimensional MS models, where the number of parameters are small and fixed, are ubiquitous in macroeconometrics and empirical finance. Univariate and multivariate MS models have been used in investigations of the business cycle hamilton1989new, monetary policy sims2006were, exchange rates and currency crisis cerra2005did, ichiue2011regime, forecast combinations elliott2005optimal, equity returns henkel2011time, and asset allocation ang2002international, guidolin2007asset, guidolin2008international among other applications\footnote{See the surveys by hamilton2016macroeconomic and ang2012regime for a more extensive list of theoretical and empirical applications in macroeconomics and finance respectively}.
Many empirical investigations in macroeconomics and finance are fundamentally questions about high-dimensional endogenous systems. For example, bianchi2019modeling study the stock returns of 83 S&P100 firms, within the framework of Markov-switching Seemingly Unrelated Regressions, to address questions on network centrality and systemic risk. If we consider larger models or when the number of variables are allowed to diverge, the modeling of such systems may lead to parameter proliferation which causes estimation to be computationally infeasible or unstable. On the other hand, forecasting or policy analysis that restricts the system to only a few variables may suffer from significant omitted variable bias chan2020reducing, koop2013large, banbura2010large.
There are at least four traditional responses to this problem: (i) using aggregated data to conserve on degrees of freedom (e.g. stock indices instead of individual prices); (ii) limiting the number of endogenous variables either in an ad-hoc fashion or by invoking economic theory; (iii) factor-augmented approaches introduced by bernanke2005measuring and related work on factor analysis by stock2002forecasting and bai2008large, and (iv) Bayesian estimation when the number of parameters is relatively large but still smaller than the sample size to obtain more stable parameter estimates. Note that (i) and (ii) are compromises and might entail a change in research question, while (iii) does indeed work with high-dimensional data and can be extended to the regime-switching setting as in liu2016regime.
More recently, kock2015oracle introduced an alternative to factor models in estimating stationary high-dimensional VAR(p) models by means of shrinkage. The authors showed that penalization via the (adaptive) least absolute shrinkage and selection operator (Lasso) allows one to asymptotically recover the sparsity patterns in the coefficient matrices (here, $A_{j,T}$). han2015direct show similar theoretical results under milder conditions but with the number of lags fixed at 1, while zhu2020nonconcave provide similar results for the smoothly clipped absolute deviation (SCAD) penalty.
The model in (ref) can be viewed as a generalization of kock2015oracle by introducing dependency of the multivariate process and its sparsity pattern on a latent state, while allowing for exogenous covariates. We propose two shrinkage type estimators for high-dimensional MS-VARs using either the Lasso or the SCAD penalty, which is a folded concave penalty\footnote{See fan2014strong for a characterization.}. We include SCAD penalization in our investigation because it requires only mild conditions for model selection consistency as opposed to the Lasso\footnote{The Lasso requires a strong irrepresentability condition zhao2006model, which may not hold in some empirical settings.}. More practically, folded concave penalties have been documented to yield better finite sample performance in empirical applications.
The extension for MS is notably non-trivial as it classically requires some form of MLE, which is not only computationally more intensive compared to the equation-by-equation Lasso framework in kock2015oracle, but also theoretically more challenging. Furthermore, the inclusion of folded concave penalties to a optimization problem generally introduces multiple local optima. Existing asymptotic results on high-dimensional penalized MLE, for example kwon2012large for the SCAD penalty, establish oracle properties for a theoretic local optimum\footnote{The oracle property is satisfied when the penalized MLE is asymptotically equivalent to the estimator obtained from maximizing the likelihood with irrelevant parameters and penalty terms excluded.}. However, there is no guarantee that a computed solution to our penalized likelihood function is indeed the desired local optimum. For that to be true unequivocally, we require at least an assumption of strict concavity which is only possible when the number of parameters is smaller than the sample size (see theorem 3 of kim2008smoothly or theorem 2 of kwon2012large), and hence is not directly applicable here.
To deal with this without deviating from our high-dimensional setting, we adopt a similar estimation strategy as in fan2014strong, which relies on the local linear approximation (LLA) algorithm for the folded concave penalties introduced by zou2008one. This approach does not attempt to find the desired local optimum, but instead seeks a lower probability bound to the claim that the computed solution is indeed the desired optimum. As the SCAD penalty requires an initial estimator, theorem 1 of fan2014strong claims that if we can find an initial estimate that is asymptotically close to the true parameter vector, in a sense to be made clear in \thref{fsconst}, then the LLA algorithm delivers the oracle estimator in a single step with probability approaching one. Borrowing from their terminology, we will call such an initial estimator "localizable". However, given that the likelihood function of MS models tend to be general and potentially multimodal, we discipline our investigation by requiring the likelihood surface at the true parameter to be locally concave and the oracle solution to be unique. Such a restriction is nonetheless considerably weaker than the strict concavity condition required in kim2008smoothly or kwon2012large. Taken together, this suggests a two-stage approach for the SCAD estimator, similar to that in li2015model and chen2020time, with the search of a localizable initial estimator in the first stage, and the SCAD-penalized MLE subsequently.
Our results open up the possibility of studying Markov regime-switching dynamics in high-dimensional macroeconomic and financial systems. To illustrate, we extend the investigation on short-horizon stock return predictability considered in henkel2011time to incorporate 14 aggregate return predictors from welch2008comprehensive, up from the original 4, which mitigates potential omitted variable bias. This results in a relatively high-dimensional system with 711 estimable parameters. Penalized maximum likelihood estimation with the SCAD penalty yields the often reported counter-cyclicality in return predictability. Furthermore, the variable selection property of our estimators helps to identify predictors that contribute strongly to predictability during economic contractions but are otherwise irrelevant in expansions.
The rest of this paper is organized as follows. Section 2 describes the penalized maximum likelihood problems while Section 3 provides an EM algorithm for the special case of Gaussian errors. Section 4 contains the key asymptotic results. Section 5 consider Monte Carlo experiments to assess the finite sample properties of the proposed estimators, while section 6 applies it to the problem of short-horizon return predictability. Finally, section 7 concludes. All proofs are collected in the appendix.
In this section, we describe the maximum likelihood problem in a manner amenable to the derivation of theoretical results. This construction is however, not convenient for computation, which will be accomplished in section (ref). There, we will rely on the Expectation-Maximization (EM) algorithm.
Before proceeding, we define some notation. For this section and the next, we suppress the dependence of the parameters on the sample size for notational convenience.
To begin, set $\phi \equiv (\theta_A^\top, \theta_B^\top, Q^{ND\top}, Q^{D\top}, \pi^\top)^\top$ to include all parameters across all states. Let $A(s) = [A_1(s), \ldots, A_{p_T}(s)]$, and $A = [A(1), \ldots, A(M)]$, then $\theta_A = vech(A)$, which is a $M p_T d_T^2 \times 1$ vector. $\theta_B$ is a $Mq_Td_Td_T^* \times 1$ vector that is defined analogously. Recall that the precision matrix is given as $Q(s) = \Sigma(s)^{-1}$. Let $Q^{ND}$ be the vector containing unique off-diagonal elements of the precision matrix across all states i.e. $Q^{ND} = (q_{21}(1), q_{31}(1), \ldots, q_{d_T1}(1), q_{23}(1), \ldots, q_{d_T(d_T-1)}(1), q_{21}(2), \ldots, q_{d_T(d_T-1)}(M))^\top$, and $Q^{D}$ includes the diagonal elements: $Q^{D} =(q_{11}(1),\ldots,q_{d_Td_T}(1),\ldots, q_{d_Td_T}(M))^\top$. $Q = (Q^{ND\top}, Q^{D\top})^\top$ is $M d_T (d_T +1)/2 \times 1$. Finally, $\pi = (p_{1 \shortrightarrow 1}, p_{1 \shortrightarrow 2}, \ldots, p_{1 \shortrightarrow M},\ldots, p_{M \shortrightarrow M}, p_1, \ldots, p_M)^\top$ is a vector of transition probabilities and initial distributions.
Written this way, it is clear that our characterization of sparsity is equivalent to saying that the vectors $\theta_A$ , $\theta_B$, and $Q^{ND}$ are sparse. Formally, we assume that only $h^A$, $h^B$, and $h^{Qnd}$ elements of $\theta_A$, $\theta_B$, and $Q^{ND}$ respectively, are non-zero and that $h^A + h^B + h^{Qnd} \ll T$. Sparsity cannot be imposed on $Q^{D}$ to maintain positive definiteness.
To begin, let $\mathcal{I}^{t-1}_{t-v} = (y_{t-v}, \ldots, y_{t-1}, x_{t-v}, \ldots, x_{t-1} )$ be the information set from time $t-v$ to $t-1$ for some integer $v$, $\Phi_{s} = (A(s), B(s), Q(s))$ be the state-specific VAR parameters, and define $\mathcal{Y}_T = (y_1, \ldots, y_T)$ and analogously for $\mathcal{X}_T$. Without loss of generality, assume that $p_T = \max\{p_T, q_T\}$ is the largest lag, and thus we can denote the conditional density of $y_t$ from (ref) as $g(y_t| \mathcal{I}^{t-1}_{t-p_T}; \Phi_{S_t})$.
Then, it can be shown that the (non-penalized) conditional likelihood is given as
Here we emphasize that $p_{i \shortrightarrow j}$ and the initial distributions $p_j$ are functions of the parameter vector $\phi$. As prefaced earlier, direct optimization of (ref) is challenging because the transition probabilities $p_{i \shortrightarrow j}$ are highly non-linear functions of the parameters. To deal with this issue, we will rely on a modified EM algorithm in section (ref).
We propose solving the following optimization problem with Lasso penalty for VAR coefficients and graphical Lasso (gLasso) penalty for the precision matrix friedman2008sparse to induce sparsity,
where $a_{mn}(s)$ refers to the element in the $m^{th}$ row and $n^{th}$ column of $A(s)$, and analogously for $b_{mn}(s)$ and $q_{mn}(s)$ for $B(s)$ and $Q(s)$ respectively. $\lambda^{Lasso}$ and $\lambda^{gLasso}$ are penalty terms. Note that only the non-diagonal elements are penalized in the sparse precision matrix estimation. It is possible to set $\lambda^{Lasso}$ to be different for parameters in $A(s)$, and $B(s)$ as long as the penalties are proportional to one another, however, doing so may introduce significant computational costs when tuning the penalty terms. Nonetheless, we call the solution to (ref) the Lasso estimator.
To use the LLA zou2008one for the SCAD problem, we require a localizable initial estimate, $\tilde{\phi}$. As we show later in section (ref), the Lasso estimate above satisfies this property. Hence, we re-optimize the log-likelihood with the Lasso estimate as $\tilde{\phi}$ with the LLA for the SCAD penalty. This utilizes the first derivative of the penalty which is defined as
where $(m)_+ = m$ if $m > 0$ and $0$ otherwise, $\lambda$ is a penalty parameter, and $a>2$ is a constant. Here, we set $a=3.7$ as suggested in fan2001variable. The optimization problem is given by
where $\lambda$ and $\lambda^*$ are penalty terms, and parameters with tilde are from the initial estimate $\tilde{\phi}$.
A key theoretical advantage of using the SCAD penalty instead of simply stopping once we have obtained the Lasso estimate, is that we can, under relatively mild conditions, attain selection consistency, or in other words, asymptotically recover the true sparse support of the parameters. On the other hand, penalization with Lasso requires strong irrepresentability conditions zhao2006model. However, if prediction is the main goal of estimation, then both Lasso and SCAD are applicable.
In this section, we propose an EM algorithm to solve (ref) and (ref) similar to that of monbet2017sparse. For concreteness and parsimony, we focus only on the scenario where $\varepsilon_t \sim^{i.i.d} N(0, \Sigma(S_t))$ in this section, which implies that the conditional density $g(y|\cdot)$ is Gaussian. Not only is this assumption standard in the VAR literature, it also induces a closed form problem which simplifies the estimation of the VAR parameters.
The EM algorithm was proposed to estimate models with incomplete or hidden data. baum1970maximization applied the algorithm to estimate hidden Markov models, which is a general class that encompasses many MS models in econometrics. Intuitively, the algorithm works not by directly optimizing the likelihood function in (ref), which as mentioned earlier is a highly non-linear function of the parameters, but instead optimizes a constructed auxiliary function (label it $\Omega(\cdot)$) to derive a monotonically increasing lower bound on the value of the original likelihood.
Formally, the algorithm works iteratively with two steps per iteration: expectation (E) and maximization (M). With each iteration applied to $\Omega(\cdot)$, we obtain updates on the parameters. The corresponding likelihood $\mathcal{L}(\mathcal{Y}_T|\mathcal{X}_T;\phi)$ is guaranteed to be non-decreasing with each update and eventually finds a stationary point\footnote{This can mean a local or global maximum, or a saddle point.} in the likelihood surface dempster1977maximum. Although the statistical guarantees on the EM algorithm have been developed for the low-dimensional context, it is not difficult to show that Theorem 1 of dempster1977maximum (monotonicity of the algorithm) will still hold with penalization.
Since both estimators in section (ref) entail a maximization problem, we can apply the EM algorithm to optimize either (ref) or (ref). The E-step for both problems will be similar, while the M-step is different because of differences in penalization.
To begin, we define the auxiliary function in the E-step.
E-step. We consider the complete-data log-likelihood given as $\ell(\mathcal{S}_T,\mathcal{Y}_T;\mathcal{X}_T, \phi)$ where $\mathcal{S}_T = (S_1, \ldots, S_T)$. Note that the incomplete-data likelihood in (ref) can be written as
and thus
The log-likelihood is called 'complete' because it treats the latent state as observable data. However, since we do not actually observe it, we consider its conditional expectation, $E[\ell(\mathcal{S}_T,\mathcal{Y}_T;\mathcal{X}_T, \phi) \\| \mathcal{I}_T, \phi^{(j-1)}]$, where $\phi^{(j-1)}$ is the estimate from the previous $(j-1)^{th}$ EM iteration and is treated as pre-determined during the current iteration. Since we have additive penalties, we can construct the auxiliary function as
where pen(N) for N $\in \{\text{Lasso, SCAD}\}$ refers to either group of penalties. Ignoring penalties for now, we have
By construction, pen(N) would be relevant to the optimization of VAR parameters in $Z_1(s)$. On the other hand, $Z_2(s,s^{'})$ and $Z_3(s)$ are only a function of $\pi$ and thus can be maximized by standard optimization procedures as this is a low-dimensional problem.
Since we have Gaussianity of $g(y_t|\cdot)$, we can show that
where
and
Note that $Z_2(s,s^{'}; \pi)$ contains the smoothed probabilities $P(S_t=s^{'}|S_{t-1}=s, \mathcal{I}_T, \phi^{(j-1)})$, which can be computed via an iterative backward-forward recursion (see for e.g. hamilton1990analysis).
M-step. Here, we maximize the auxiliary function in (ref) with respect to $\phi$. Suppose we are executing the EM algorithm for the N estimator (i.e. either (ref) or (ref)). Then, for all states $s$, we solve
where pen(N,$s$) refers to the $s^{th}$ state penalties for the N estimator. For computational efficiency, it might be convenient to further separate (ref) into two parts. To illustrate, consider a M-step for the Lasso. Firstly, fix the coefficient values $A(s)$ and $B(s)$, and estimate
Next, given an estimate of $Q(s)$, we optimize
This sequential partitioning helps with computation because we can now individually apply fast (block) coordinate descent algorithms friedman2007pathwise, friedman2008sparse to each part of the problem. Note that these algorithms will also work with the SCAD penalty as defined in (ref).
The M-step provides an update of the parameters $\phi^{(j)}$. It can be shown that repeated updating of $\phi^{(j)}$ will lead the respective original likelihoods in (ref) or (ref) to either increase or remain constant in value, but not decrease. Hence, the maximizer can be found subject to some termination condition on the implied increments of the original penalized likelihood.
To implement the algorithm, we require an initial vector of parameters $\phi^{(0)}$. We recommend using reasonable randomly generated values for the Lasso problem, and subsequently the computed Lasso optimum as $\phi^{(0)}$ for the SCAD problem.
Similar to kock2015oracle, we propose selecting $\lambda^{Lasso}$ and $\lambda^{gLasso}$ for the Lasso estimator, and $\lambda$ and $\lambda^*$ for the SCAD estimator, using a modified version of the Bayes Information Criterion (BIC). Specifically, for either Lasso or SCAD, choose $\lambda_1$ and $\lambda_2$ to minimize
where
and parameters with 'hats' indicate estimates. Following wang2009shrinkagetuning, $C_T$ is set at $\log K_T$, where we recall that $K_T$ is the total number of parameters in the system, to help obtain consistency of the BIC in high-dimensional regressions, while $l$ is the number of estimates that are identified as non-zero in the system. Simulation results in section (ref) indicate that tuning parameter selection with BIC can yield consistent estimates.
We establish theoretical properties for each estimator in a double asymptotic framework, where we allow the sample size $T \rightarrow \infty$ and the dimensions of the candidate models $p_T, q_T, d_T,$ and $d_T^*$ to diverge at appropriate rates. By extension, it is natural to allow the dimensions of the true non-zero parameters $h^A_T$, $h^B_T$, and $h^{Qnd}_T$ to grow albeit subject to stricter restrictions. Consequently, the dimension of the true parameter vector, $\phi_T$, may extend to infinity.
Traditionally, theoretical results on consistency in MS models with finite state spaces rely on limit theorems from random matrix theory to attain some generalized form of the Kullback-Leibler divergence leroux1992maximum, francq1998ergodicity, whereby consistency follows from an identification condition. This strategy is challenging in the high-dimensional context because of the dependence of $\phi_T$ on the sample size. Instead, we adopt a modified approach to fan2004nonconcave to establish consistency and oracle properties under a diverging parameter framework.
Recall that $K_T$ is the total number of VAR and transition matrix parameters from (ref). Let $\phi^*_T \in \Theta_T$ be a vector of true (sparse) parameters where $\Theta_T \subset \mathbb{R}^{K_T}$ is an open subset. In addition, let $K^{Sp}_T = h_T^A + h_T^B + h_T^{Qnd}$ be the number of non-zero parameters out of those that are subject to penalization (i.e. the VAR parameters), and let $K_T^* = K^{Sp}_T + Md_T + M^2$ be the total number of non-zero parameters in the true parameter vector where $M^2$ and $Md_T$ are from the transition matrix and the diagonals of variance-covariance matrices respectively. Without loss of generality, we assume that $\phi^*_T$ can be re-arranged in the following form
We impose the following regularity conditions for deriving our asymptotic results. Let $\mathcal{L}(\mathcal{Y}_T| \mathcal{X}_T; \phi_T) = \mathcal{L}_T(\phi_T)$. In addition, $\nabla^{k}$ represents the $k^{th}$ derivative with respect to $\phi_T$, and $\nabla^{k}_{j_1, \ldots, j_k}$ represent the $k^{th}$ derivative with respect to the $j_1, \ldots, j_k^{th}$ element in $\phi_T$.
Remarks. (ref)-(ref) are standard assumptions in the MS literature on asymptotic normality bickel1998asymptotic. The assumption of local concavity in (ref) is much weaker than that of strict concavity commonly assumed in high-dimensional regularized MLE problems. We need the likelihood to be (locally) concave between the localizable initial estimator and the true parameter to apply theorem 1 of fan2014strong so that the SCAD estimation initialized with Lasso estimates delivers the oracle result. The assumption of a locally concave neighborhood is not uncommon in the theoretical literature on high-dimensional problems. For example, in the context of a two-step estimation procedure for high-dimensional sparse principal components analysis, jankova2018biased assumes that a rough initial estimate can be found in a locally convex neighborhood of the population parameter. (ref)(i) is commonly assumed in the SCAD literature and helps us to establish selection consistency, while (ref)(ii) is essentially identical to assumption 2(ii) in li2015model. (ref)(ii) is a modification of the restricted eigenvalue condition in bickel2009simultaneous, and this particular formulation follows assumption 2(iii) in li2015model. (ref) restricts the number of true parameters $K_T^*$ to only increase at a rate slower than $T^{1/4}$. Furthermore, $\lambda^{lasso} \rightarrow 0$ but converges slower than $1/\sqrt{T}$, and $\lambda \rightarrow 0$ but slower than $\sqrt{K_T^*}\lambda^{lasso}$. (ref) is a standard assumption on the minimum signal strength of relevant parameters. This condition is required for proving the oracle property of the SCAD estimator and mirrors the characterization in zhu2020nonconcave.
Our first result states that the Lasso estimator is consistent. Let $\|\cdot\|$ denote the $\ell_2$ norm.
We obtain the final equality because of (A8). This means that the Lasso estimator exhibits estimation consistency, although we do not guarantee that it is able to asymptotically identify the relevant parameters.
Next, define $\hat{S}=\{j: |\hat{\phi}_{T,j}| > 0\}$ where $\hat{\phi}_{T,j}$ is the $j^{th}$ element in $\hat{\phi}_{T}$. Hence, $\hat{S}$ is the index set of estimated non-zero parameters from the SCAD procedure. Define $S_0$ for $\phi_T^*$ in an analogous manner (i.e. the index set of truly relevant population parameters). The following result shows that the SCAD estimator, initialized by the Lasso estimates, achieves not just estimation consistency, but also selection consistency.
Part (1) of \thref{consistency} states that the SCAD estimator is consistent since $\sqrt{K_T^*/T} \rightarrow 0$ by (A8). The second part claims that the sparsity property holds. In other words, we are able to exactly distinguish the parameters that are truly non-zero from those that are irrelevant in the theoretical limit. Note that these results hold even if the dimensions of the parameter vectors diverge.
Next, we discuss the asymptotic normality of the SCAD estimator. To do so, define $\hat{\phi_T}^{(\neg 0)}$ to be the estimate $\hat{\phi}_T$ with the zeroes removed, and similarly for $\phi_T^{*(\neg0)}$.
We remark that the pre-multiplication of $G_T$ helps with the exposition since $\phi_T^*$ can be diverging in dimension. \thref{norm} implies that it is asymptotically justifiable to apply the same statistical inference for the estimate obtained from maximizing the problem with a priori knowledge on the sparsity pattern, to the SCAD-penalized solution $\hat{\phi}_T$. Here, we note that this convergence in distribution is a pointwise result, and is not guaranteed to hold uniformly with respect to the parameter vector leeb2005model. Uniform inference for model selection in multivariate time series is however still a nascent area of research masini2020machine and is beyond the scope of this paper.
This section studies the finite sample properties of the proposed EM algorithm in estimating large MS-VARs using both Lasso and SCAD penalization schemes\footnote{For SCAD, we initialize the algorithm with Lasso estimates as described earlier.}. The number of endogenous variables considered are $d = 10$ and $16$. We consider three experiments\footnote{The numerical experiments here are similar to those considered in kock2015oracle.} with 2 states ($M=2$) throughout:
In all experiments, we exclude an intercept and set the transition probabilities to be $p_{1\rightarrow 1} = p_{2\rightarrow 2} = 0.8$.
500 iterations are generated for sample sizes $T=100,200,$ and $300$. To evaluate our procedures, we consider the following metrics. The first three metrics are concerned with selection consistency. True model included looks at the share of iterations in which the estimate includes the true model for both states (i.e. truly non-zero coefficients are estimated as non-zero). Selected variables is the number of non-zero parameters estimated by the system. As it may be the case that the estimation sets a truly non-zero coefficient to zero, it is informative to study the share of truly non-zero parameters that are identified correctly as non-zero by the algorithm. Subsequently, we consider estimation consistency as measured by the root mean squared error (RMSE) of the parameters. This is given by $\sqrt{\frac{1}{200} \sum_{i=1}^{200} \| \hat{\phi}_T(i) - \phi_T^* \|}$ where $\hat{\phi}_T(i)$ is the estimated parameter vector for iteration $i$ containing all parameters in the system. $RMSE_{VAR}, RMSE_{COV},$ and $RMSE_{p}$ are similarly defined RMSE measures for estimated VAR coefficients, variance-covariance parameters, and transition probabilities respectively.
The results of the experiments are presented in Table (ref). Looking at the metrics for estimation consistency, we see that all $RMSE$ measures are declining as the sample size increases for both Lasso and SCAD. It is interesting to note that SCAD appears to perform better for larger sample sizes in terms of estimation accuracy as measured by $RMSE$. Furthermore, we see that, in most cases, both estimators get better at including the true model while the number of selected parameters fall, which provides evidence of model selection consistency. We also note that SCAD has a slight advantage over the Lasso in selection accuracy for larger sample sizes as seen in experiments 2 and 3.
It is well established that short-horizon stock return predictability (both in- and out-of-sample) exhibits significant time-variation chen2012testing, rapach2013forecasting. Specifically, henkel2011time (henceforth HMN) and dangl2012predictive show that return predictability is correlated with business cycles in a distinctively counter-cyclical fashion. In particular, HMN estimate a MS-VAR(1) with 2 states, one state corresponding to an 'expansion' regime and the other to a 'recession', for the period of 1953 to 2008 with excess returns, dividend yield, short rate, term spread, and default spread as endogenous variables for the US. With the estimated system, they calculate the adjusted $R^2$ ($\overline{R}^2$) from a predictive regression of lagged predictors on one-month ahead excess returns (the first equation from the VAR\footnote{Excess returns are ordered first in the system.}) for both model-implied periods of recession and expansion. The authors find that $\overline{R}^2$ is close to 0 during expansions while it averages around 0.175 for recessions. This led them to argue that return predictability is negligible during expansions and is present exclusively during recessions.
A major drawback of this strategy is that the predictors were selected in a relatively ad-hoc manner, which may yield significant omitted variable bias. The framework precludes the possibility that changes in return predictability are due to predictors beyond that of the chosen variables. This is a significant problem because if, as the authors argued, aggregate predictor variables are jointly determined by the "micro-level objectives of firms and central banks" which are in turn driven by business cycles, we should expect this counter-cyclical relationship between predictor and excess returns to potentially manifest in any predictor that relate to the "micromotives" of economic agents, which the literature on return predictability is not short of.
We approach this problem by incorporating 14 predictors as considered by welch2008comprehensive in their study on return predictability: dividend price ratio (d/p), dividend yield (d/y), earnings price ratio (e/p), dividend payout ratio (d/e), stock variance (svar), book-to-market ratio (b/m), net equity expansion (ntis), treasury bill rate (tbl), long term rate of returns (ltr), long term yield (lty), term spread (tms), default yield spread (dfy), default return spread (dfr), and inflation (infl). Note that d/p, tbl, tms, and dfy overlap with the predictors in HMN. This results in a 2 state MS-VAR(1) system with 15 endogenous variables including excess returns (r) ordered first, with 711 estimable parameters including an intercept\footnote{To maintain the assumption of stationarity, we detrend all time series by subtracting a moving average of the past 12 months following ang2007stock. Furthermore, as is customary in the literature on Lasso, we rescale all variables to have zero mean and a standard deviation of 1.}. Additionally, we update the investigation period to cover April 1953 to December 2018 (sample size of 787).
Given the large number of parameters involved, the system is estimated with the SCAD penalty initialized with Lasso estimates. As a diagnostic check on the EM algorithm, we first verify that the estimated system yields regimes consistent with the interpretation of expansions and recessions. Figure (ref) compares the smoothed probabilities for state 1 with that of an NBER-based recession indicator from the St. Louis FRED database (USREC). We see that when the probability of state 1 is close to unity, the NBER-based index often indicates a recession (value of 1), which provides evidence that state 1 corresponds to that of a recession. More formally, if we classify recessions as having a state 1 probability of greater than 0.5 and subsequently compare the resultant series with the NBER-based indicator, we get an agreement rate of 70%, which is comparable to the 77% obtained in HMN.
Several other features of our system are similar. The estimated transition probability matrix yields $p_{1\rightarrow 1} = 0.79$ (recession to recession) and $p_{2 \rightarrow 2} = 0.88$ (expansion to expansion), which are close to the respective estimates of $0.80$ and $0.91$ in HMN. We also replicate the finding of increased volatility for all predictors during recession relative to periods of expansion as reported in table (ref).
Panel A of Table (ref) reports the $\overline{R}^2$ derived from using the first equation of the MS-VAR (since $r$ is ordered first) in predicting one-month ahead excess returns conditional on state. These results present strong evidence in favor of the counter-cyclicality of return predictability given the larger magnitude of $\overline{R}^2$ during recessions (0.126) compared to expansion (0.112), which is in line with HMN. However, our estimate of $\overline{R}^2$ during expansions is at least 4 times larger than theirs, which is consistent with the findings of dangl2012predictive that return predictability might still be present during booms, but are indeed stronger during busts.
A key merit of our regularized MS-VAR in this context is the ability to identify predictors that contribute the most to this counter-cyclical behavior of return predictability. Looking at Panel B of Table (ref), we see that 2 predictors are relevant to both states: d/p, and svar. More importantly, we see that tbl drops out during times of crisis, which is reasonable given that the policy rate has recently been set constantly close to zero during protracted periods of economic slowdown. The inclusion of 6 new relevant predictors during recessions (e/p, ltr, tms, dfr, infl and lty) indicate that they contribute to return predictability during bad times only, which is a finding that is challenging to replicate in small MS-VARs with few endogenous variables. Interestingly, of all the selected variables, only d/p, tbl and tms were included in the MS-VAR of HMN, which suggests a potential omitted variable bias problem. Capturing these dynamics matter because they may help to consolidate the narrative on return predictability. For example, the non-relevance of e/p during expansions and its relevance during recessions is consistent with the idea of "accounting conservatism", as argued by HMN, which meant that balance sheets may be quicker in reporting bad news while reacting slower during good times. Slower disclosure and propagation of reported earnings information during expansions may hinder the predictive power of earnings ratios during such a period, while a faster response with bad news during crisis may lead to more accurate reporting of earnings, and contribute to predictive power.
Next, we study the out-of-sample (OOS) forecasting properties of the proposed large MS-VAR with an OOS period of January 2010 to December 2018. The model is estimated with an expanding window and one-step ahead forecasts are constructed as the weighted average of predictions from both states. To compare our results, we consider also the original 5 variable MS-VAR estimated in HMN. In addition, we form forecasts of $r$ using an expanding window historical average and include predictions from an ARMA(1,1) model. Specifically, the historical average is often hailed as notoriously difficult to beat welch2008comprehensive. To formally compare the results, we employ two forecast comparison tests: Diebold and Mariano test diebold1995comparing and the 'Reality check' test white2000reality.
Panel C of Table (ref) reports the results. Notably, the proposed MS-VAR achieves the lowest mean squared forecast errors (MSFE) followed by the MS-VAR of HMN. Although the difference between the proposed MS-VAR and that of HMN is not statistically significant, we note that the large MS-VAR performs significantly better than both the historical average and ARMA(1,1), as indicated by both tests.
Additionally, we consider an out-of-sample period encompassing the 2008 financial crisis in Panel D. Here, we expect regime-switching methods to perform better given that the possibility of having a recession state is explicitly accounted for. Indeed, we see that, relative to the historical average benchmark, this is the case for the SCAD-penalized MS-VAR but not for the original MS-VAR in HMN. One conjecture for the difference is that, as suggested in Panel B, many of the predictors that are relevant for predictability during recessions are included in the large MS-VAR but not the HMN model. This provides some evidence that prediction performance can be improved by considering more predictors. This is particularly the case because the large MS-VAR continues to perform the best in terms of achieving the lowest MSFE although this difference in predictive accuracy does not appear to be statistically significant.
In this paper, we have proposed two new shrinkage type estimators to handle parameter proliferation in sparse high-dimensional MS-VARs. Theoretically, we have shown that both the Lasso and SCAD estimators are estimation consistent, while the latter has the added benefit of selecting relevant variables with high probability. Consequently, the SCAD estimator exhibits the oracle property in that it is asymptotically equivalent to an estimator that assumes a priori knowledge of the sparsity pattern in the system. Results from numerical experiments show that the proposed EM algorithm is able to handle large MS-VARs well and the finite sample performance of the estimators provides support for our theoretical results. The empirical investigation on the counter-cyclicality of return predictability highlights the flexibility of the proposed estimation in incorporating many endogenous predictors and the merit of allowing for variable selection in regime-switching applications. Furthermore, the significantly better OOS performance of our model suggests that sizeable improvements to stock return prediction can be attained with a larger pool of predictors. Notably, this suggests that our framework can be generalized to other applications where regime-switching is of interest but where high-dimensionality may be a limiting factor, such as those frequently encountered in monetary policy, asset allocation, and other macroeconomic or financial systems.
\processdelayedfloats