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.
126,918 characters · 19 sections · 44 citation commands
Combining Shrinkage and Sparsity in Conjugate Vector Autoregressive Models0.75cm
\thispagestyle{empty}
\setcounter{page}{1}
This paper deals with estimating vector autoregressive (VAR) models of the following form,
where $\bm y_t = (y_{1t}, \dots, y_{mt})'$ denotes an $m$-dimensional vector of time series measured in time $t=1,\dots,T$, $\bm A_j~(j=1,\dots,p)$ is an $(m \times m)$-dimensional matrix of coefficients associated with the $j^{th}$ lag of $\bm y_t$, $\bm C$ is an $m$-dimensional intercept vector, and $\bm \varepsilon_t \sim \mathcal{N}(\bm 0, \bm \Sigma)$ is a Gaussian shock vector with zero mean and an $(m \times m)$-dimensional variance-covariance matrix $\bm \Sigma$. For further convenience, let $\bm a = \text{vec}\{(\bm A_1, \dots, \bm A_p , \bm C)' \}$ denote a vector of dimension $k=m(mp+1)$ of vectorized coefficients with $a_i ~(i=1,\dots, k)$ denoting its $i^{th}$ element. This model class has been extensively used for forecasting and policy analysis in central banks \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see][]{alessi2014central} as well as a natural starting point for unveiling stylized time series facts in order to estimate theoretical models \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see, inter alia,][]{hall2012information}.
Conditional on the first $p$ observations, estimation of the model in (ref) can be carried out using ordinary least squares (OLS). In this case, however, overfitting issues arise, translating into imprecise out-of-sample forecasts. As a potential solution, the Bayesian literature uses informative priors to push the system towards a prior model. For instance, the widely used Minnesota prior assumes that the elements in $\bm y_t$ follow a random walk a priori \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{doan1984forecasting, litterman1986forecasting, giannone2015prior}. Theoretically inspired restrictions stemming from structural models can also be used to inform parameter estimates and thus improve inference \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{ingram1994supplanting, del2004priors}. The key feature of these priors is that they are conjugate, implying that the likelihood and the prior feature the same distributional from. This yields closed-form solutions for key posterior quantities and, if simulation-based techniques are necessary, greatly improves estimation speed.
In VARs conjugacy requires that each equation in the system features the same set of predictors, potentially leading to model misspecification \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see, for example,][]{george2008bayesian,koop2013forecasting}. This translates into a Kronecker structure in the likelihood, prior, and the resulting posterior distribution, implying that inversion of the posterior variance-covariance matrix of the coefficients is computationally cheap. In contrast, models based on non-conjugate priors allow for different predictors across equations by specifying the prior on the VAR coefficients independently of $\bm \Sigma$. This, however, is computationally much more demanding, since the convenient Kronecker structure is lost.\footnote{For recent solutions that allow for estimating large-scale VARs under non-conjugate priors, see carriero2019large.}
Apart from reduced flexibility in terms of covariate selection across equations, typical shrinkage priors push many VAR coefficients towards zero. Under continuous shrinkage priors, however, this implies that the probability of observing a coefficient that exactly equals zero is zero \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see, for example,][]{griffin2010inference,carvalho2010horseshoe, bhattacharya2015dirichlet, huber2017adaptive, huber2019inducing}. Spike and slab priors allow for shrinking coefficients exactly to zero. These priors rely on an additional set of auxiliary binary indicators that determine whether a coefficient is zero or non-zero. In large models with $k$ covariates (such as the VAR models we consider), estimating these indicators is cumbersome since the number of potential models is $2^k$. In such a situation, Markov chain Monte Carlo (MCMC) techniques often fail to explore this vast model space and convergence issues arise \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see][]{polson2010shrink}.
In recent contributions, \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{hahncarvalho2015dss} and \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{bhattacharya2018signal} propose a way to circumvent insufficient zero shrinkage/variable selection in light of an increasing amount of predictors (i.e. the curse of dimensionality problem). They estimate a large-scale regression model under a suitable shrinkage prior and then post-process a point estimator (the posterior mean) such that the distance between the fit of the model based on the shrinkage prior and a model based on a sparse estimator (i.e. with coefficients set equal to zero) is minimized while accounting for a penalty term that depends on the $L1$-norm of the coefficients. This approach, labeled decoupled shrinkage and selection (DSS), yields a sparse estimator and is analogous to solving a LASSO-type problem. One key disadvantage, however, is the non-automatic nature of this approach. A semi-automatic approach that is similar in nature is described in \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{bhattacharya2018signal}. In this framework, an optimization problem to efficiently set coefficients associated with irrelevant predictors to zero is solved. But instead of performing cross validation, this task depends only on a single tuning parameter that needs to be chosen by the researcher. These techniques that combine shrinkage and sparsity have been shown to work well in a wide range of applications ranging from finance \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{puelz2017variable, puelz2019portfolio} to macroeconomics \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{huber2019inducing}.
In this paper, we deal with both issues discussed above by proposing a fully conjugate VAR model coupled with the prior proposed in \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{kadiyala1997numerical} and \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{koop2013forecasting}, that allows for different covariates across equations, sparsity in terms of the VAR coefficients, and data-based zero restrictions on the covariance parameters in $\bm \Sigma$. Instead of post-processing posterior mean/median estimates, we follow huber2019inducing in sparsifying each draw from the posterior distribution. This yields an approximate posterior distribution for a sparse vector of coefficients which can be used for uncertainty quantification. The key advantage is that this significantly reduces estimation uncertainty if the data generating process is sparse. For example, if there is strong evidence that $a_{i}$ equals zero, our proposed framework is capable of selecting this restriction consistently across different draws from the posterior distribution of $a_i$. This implies that the posterior variance of $a_i$, in the limiting case that each draw of $a_i$ is set to zero, also equals zero. In terms of forecasting, the reduced estimation uncertainty could then lead to more precise predictions, especially in situations where $k$ is large.
The merits of our proposed approach are illustrated by means of two applications. In the first application, we use synthetic data obtained from a set of different data generating processes (DGPs) that differ in terms of sparsity, model size, and number of observations. Across DGPs, we find that i.) our framework successfully detects zero values in both the VAR coefficients and the error variance-covariance matrices and ii.) it outperforms other Bayesian VARs (BVARs) in terms of root mean square errors (RMSEs). In the second application, we forecast US output, inflation, and short-term interest rates using the dataset compiled in mccracken2016fred. We find that applying the additional sparsification step often improves point and density forecasts. In turbulent times (such as the period of the global financial crisis), however, our results also show that using sparsification could harm the accuracy of density forecasts by underestimating the predictive variance. Nevertheless, these accuracy losses are never substantial and forecasts are still competitive.
The remainder of the paper is structured as follows. (ref) introduces the conjugate Bayesian VAR while (ref) discusses the techniques to achieve sparsity in this model. (ref) documents model features when applied to synthetic data. (ref) summarizes the results of the forecast exercise with real data. Finally, (ref) concludes the findings of the paper and an appendix provides details on data.
Before discussing prior implementation, it is worth noting that (ref) can be rewritten as a standard regression model,
with $\bm x_t=(\bm y_{t-1}, \dots, \bm y_{t-p}, 1)'$ denoting an $n(=pm+1)$-dimensional vector of explanatory variables. In terms of full-data matrices $\bm Y$ (with $t^{th}$ row $\bm y'_t$) and $\bm X$ (with $t^{th}$ row $\bm x'_t$), the model reads
where $\bm A = (\bm A_1, \dots, \bm A_p , \bm C)'$ and $\bm E$ is a $(T \times m)$-dimensional matrix of stacked shocks with $t^{th}$ row given by $\bm \varepsilon'_t$.
The model in (ref) features $k$ parameters in $\bm a$ and $w =m (m+1)/2$ free parameters in $\bm \Sigma$. If $m$ and $p$ become large, the number of parameters sharply increases, making precise estimation almost impossible. To deal with this issue, Bayesian econometricians rely on informative priors. This implies that more weight is placed on the prior and the resulting posterior distribution of $\bm a$ and $\bm \Sigma$ will be strongly influenced by the prior model (such as the random walk).
The general form of the conjugate prior in VARs assumes dependence between $\bm a$ and $\bm \Sigma$ and is given by
where $\bm a_0$ denotes a $k$-dimensional prior mean vector and $\bm V_0(\bm \delta)$ is a prior variance-covariance matrix that depends on a lower dimensional set of $q$ hyperparameters in $\bm \delta$. In what follows, we assume that $\bm y_t$ is stationary and thus a common choice for the prior mean would be $\bm a_0 = \bm 0$.
For $\bm V_0(\bm \delta)$, we use a variant of the conjugate Minnesota prior \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{kadiyala1997numerical, koop2013forecasting} that can be implemented using a set of dummy observations \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[][]{banbura2010large} that are then concatenated to $\bm Y$ and $\bm X$:
with $\bm V_0(\bm \delta) = ( \hbox{ \vbox{ \hrule height 0.7 pt \kern -1.8ex \hbox{ \kern-0.05em \ensuremath{\bm X} \kern-0.05em } } } ' \hbox{ \vbox{ \hrule height 0.7 pt \kern -1.8ex \hbox{ \kern-0.05em \ensuremath{\bm X} \kern-0.05em } } } )^{-1}$. Here, $\boldsymbol{J}_p = (1,\dots,p)'$, $\pi$ is a hyperparameter that determines the prior variance on the intercepts, and $\phi_i~(i=1,\dots,m)$ represents the prior mean associated with the coefficient on the first, own lag of a given variable (which is consequently set equal to zero). In addition, we let $\hat{\sigma}^2_i$ denote the OLS variances obtained by estimating $m$ univariate AR($p$) models for each element in $\bm y_t$. Finally, $\theta_1$ is a hyperparameter that controls the overall tightness of the prior. Lower values of $\theta_1$ imply a stronger prior belief, effectively pushing the elements in $\bm a$ towards $\bm a_0$. For $\pi$, we simply set it equal to a large value to render the prior on the intercept weakly informative. This prior setup implies that $\bm \delta = (\theta_1, \pi)'$ is a $2$-dimensional vector.
The final ingredient is a conjugate prior on $\bm \Sigma$. Here, conjugacy implies a prior on $\bm \Sigma$ that does not depend on $\bm a$ and follows an inverted Wishart distribution:
We let $s_0$ denote the prior degrees of freedom and $\bm S_0$ a prior scaling matrix. The main shortcoming of this prior choice is that shrinking specific covariances in $\bm \Sigma$ to zero is impossible. For instance, even if there exists significant evidence that contemporaneous relations across elements in $\bm y_t$ equal zero, this prior is not capable of selecting such restrictions and the resulting posterior estimate of $\bm \Sigma$ (and of its inverse) will be non-sparse.
In the general case (i.e. with any form of the prior hyperparameters), one can show that the conditional posterior of $\bm a$ is given by
with
Here, we let $\bm A_0$ denote an $(n \times m)$-dimensional matrix reshaped such that $\bm a_0 = \text{vec}(\bm A_0)$.
Using the Minnesota dummies, the posterior moments can be obtained by applying Theil-Goldberger mixed estimation \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{theil1961pure}:
where $ \hbox{ \vbox{ \hrule height 0.7 pt \kern 0.2ex \hbox{ \kern-0.05em \ensuremath{\bm Y} \kern-0.05em } } } =(\bm Y', \hbox{ \vbox{ \hrule height 0.7 pt \kern -1.8ex \hbox{ \kern-0.05em \ensuremath{\bm Y} \kern-0.05em } } } )'$ and $ \hbox{ \vbox{ \hrule height 0.7 pt \kern 0.2ex \hbox{ \kern-0.05em \ensuremath{\bm X} \kern-0.05em } } } = (\bm X, \hbox{ \vbox{ \hrule height 0.7 pt \kern -1.8ex \hbox{ \kern-0.05em \ensuremath{\bm X} \kern-0.05em } } } )'$ denote full-data matrices augmented with dummy observations.
Under the prior in (ref), the posterior distribution also follows an inverted Wishart distribution,
The posterior degrees of freedom are denoted by $s_1= T + s_0$ and $\bm S_1$ represents the $(m \times m)$-dimensional posterior scaling matrix, obtained by using the Minnesota-specific dummy observations:
A key advantage of conjugacy is the Kronecker structure in (ref), which implies that $\mathcal{\bm V} = \bm \Sigma \otimes \hbox{ \vbox{ \hrule height 0.7 pt \kern 0.2ex \hbox{ \kern-0.05em \ensuremath{\bm V} \kern-0.05em } } } $ is a block-diagonal matrix and computing the inverse or the Cholesky factor is computationally cheap. By contrast, if $\mathcal{\bm V}$ were a full $(k \times k)$ matrix, computation would quickly become cumbersome and impossible even for moderate values of $m$ and $p$. One further advantage of the conjugate prior is that the one-step-ahead predictive density and the marginal likelihood (ML) are available in closed form \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see, for instance,][]{zellner1985bayesian}. This implies that if interest centers on one-step-ahead forecasts, no posterior simulation is required.\footnote{For higher-order forecasts or other quantities such as impulse responses, Monte Carlo simulation is necessary.}
Unfortunately, the conjugate prior has two important shortcomings. First, each equation must include the same set of covariates (see (ref)), a feature that could be unappealing if the researcher wishes to introduce theoretically motivated restrictions across equations. Second, the structure of the prior variance-covariance matrix implies that for each equation $j=1,\dots,m$, the prior variance is given by $\sigma^2_{jj} \bm V_0(\bm \delta)$, with $\sigma^2_{jj}$ denoting the $(j,j)^{th}$ element of $\bm \Sigma$. One consequence of this is that across equations, the prior variances are proportional to each other. This implies that it is not possible to discriminate between coefficients on own (which we define as lags of the $j^{th}$ endogenous variable in equation $j$) and other (defined as the lags of the $i^{th}$ endogenous variable for $i \neq j$ within equation $j$) lags. The methods we discuss in the next section allow for different treatment of own and other lags in a simple way.
From a forecaster's perspective, heavily parameterized models, such as large-scale VARs, have another important shortcoming. The continuous shrinkage prior described in (ref) implies that the probability of observing exact zeros in $\bm a$ equals zero. One could ask whether it makes a big difference to zero out different $a_i$'s as opposed to setting them close to zero. Setting them close, but not exactly to zero essentially implies that there exists a lower bound of accuracy one can achieve under the specific prior distribution \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{huber2019inducing}. For small-scale systems, this has negligible implications on predictive accuracy. However, if $k$ is large (i.e. of order $1,000$ or $10,000$), parameter uncertainty adds up and potentially dominates the predictive variance. To see this point, notice that under the conjugate prior, the one-step-ahead predictive density follows a multivariate $t$-distribution \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see][]{koop2013forecasting} with predictive variance given by
with $\text{Var}( \bullet )$ denoting the variance operator, $x_{i T+1}$ is the $i^{th}$ element of $\bm x_{T+1}$ and $v_{ij}$ referring to the $(i,j)^{th}$ element in $ \hbox{ \vbox{ \hrule height 0.7 pt \kern 0.2ex \hbox{ \kern-0.05em \ensuremath{\bm V} \kern-0.05em } } } $.
Here, it can be seen that the predictive variance depends on the variance of the reduced-form shocks in $\bm \varepsilon_{T+1}$ and parameter uncertainty arising from the term in the parenthesis. Equation ((ref)) indicates that if $n$ increases, posterior uncertainty rises even if the $v_{ij}$'s associated with irrelevant predictors are small. This point clearly highlights the difference between sparsity and shrinkage, namely the fact that under a sparse model, $v_{ij}$ would be equal to zero if and only if the relevant predictor is excluded from the model. In the next subsections, we will show how this lower bound on accuracy (determined by small but non-zero values of $v_{ij}$) can be removed.
Equation ((ref)), moreover, highlights that the uncertainty associated with coefficients potentially adds up in large-scale models and the variance implied by the reduced-form shocks further influences the predictive variance. Without additional restrictions, these two sources can act in opposite directions. In the case of a large-scale model, the variance-covariance matrix might be underestimated due to overfitting while uncertainty surrounding parameter estimates is too large. The second effect is mainly driven by the fact that in VAR models and with standard macroeconomic datasets, covariates are often highly correlated and this, in combination with insufficient shrinkage, inflates variance estimates of the regression coefficients. Since these two sources play a crucial role in forming accurate forecasts, it is imperative to treat both of those carefully.
Since obtaining a sparse representation of $\bm a$ is unfeasible in high dimensions due to the necessity to explore a model space of cardinality $2^k$, we follow a different route that combines shrinkage and sparsity. Our approach follows hahncarvalho2015dss and bhattacharya2018signal and is based on manipulating an estimator $\hat{\bm a}$ ex-post by solving the following optimization problem,
with $\bm Z = (\bm I_m \otimes \bm X)$, $\bm \alpha$ being a sparse $k$-dimensional vector and $\left\| \bm m \right\|_2$ denoting the Euclidean norm of a vector $\bm m$. Equation ((ref)) consists of two components. The first part measures the Euclidean distance between the fit of an unrestricted model, estimated using the shrinkage prior described in (ref), and a sparse model determined by $\bm \alpha$. The second part is a penalty term that penalizes non-zero values in $\bm \alpha$, with $\kappa_j$ denoting variable-specific penalties. In light of large $k$ (which is almost always the case in moderately-sized VARs), choosing the tuning parameters $\kappa_j$ by means of cross-validation becomes computational prohibitive.
To circumvent the necessity to employ cross-validation, we adopt the signal adaptive variable selection (SAVS) estimator proposed in bhattacharya2018signal. We rewrite (ref) in terms of the $j^{th}$ column of $\bm Z$, $\bm Z_j$, and solve the optimization problem in (ref) for each covariate individually, adopting the coordinate descent algorithm \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{friedman2007pathwise}. This yields the following solution to the optimization problem in (ref),\footnote{Strictly speaking, this is the solution obtained after the first iteration of the optimization algorithm, which, conditional on initializing the algorithm at the posterior mean, already indicates convergence at this stage.}
for $j=1,\dots,k$, with $\text{sign}(c)$ returning the sign of a real number $c$ and $c_+ = \text{max}\{c, 0\}$.
We set the penalty term as follows:
which depends on the non-sparse estimate $\hat{a}_j$ and two hyperparameters $\lambda>0$ and $\zeta\ge 1$. Setting $\zeta \ge 1$ implies that smaller values of $\hat{a}_j$ receive a larger penalty and are likely to be zeroed out by the SAVS algorithm.
A typical choice, proposed in bhattacharya2018signal, sets $\lambda=1$ and $\zeta=2$. We show below that, in simulations, this choice works well. The approach stipulated in hahncarvalho2015dss is obtained by setting $\zeta = 1$ while inferring $\lambda$ by visually inspecting the posterior output. More specifically, hahncarvalho2015dss suggest choosing $\lambda$ such that the variation-explained by a sparsified linear predictor (which is akin to a standard R$^2$) statistically equals that of the non-sparsified model. For carrying out structural analysis, this poses no problem since it needs to be done once. However, if the researcher is interested in assessing forecasting accuracy, this procedure has to be repeated sequentially over a hold-out period. This makes the non-automatic nature of the approach problematic.
As stated above, the natural conjugate Minnesota prior is not capable of discriminating between own and lags other variables. In this paper, we modify the SAVS estimator accordingly. In what follows, we replace $\lambda$ with a lag-wise parameter that increases the weight associated with coefficients on higher order lags of $\bm y_t$ and impose a stronger penalty on coefficients related to the other lags within a given equation. Moreover, we do not sparsify the diagonal elements of $\bm A_1$. These parameters are specified for each $\bm A_l$ such that
for $l =1, \dots, p; i=1,\dots,m; j=1,\dots,m$. Here, we assume that $\lambda$ is some lag-invariant scaling parameter and $\lambda_{l,ij}$ increases quadratically with the lag order. Note that for the first, own lag of a given equation, we set the penalty equal to zero, capturing the notion that this covariate is crucial and should never be set equal to zero \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{banbura2010large}. For coefficients on lags of other variables we increase the penalty slightly by multiplying $\lambda$ with $l^2$ instead of $(l-1)^2$.\footnote{In the empirical application, we specify the penalty on the intercept term equal to zero.}
Up to this point, we have focused attention on obtaining a sparse representation of the VAR coefficients. In large dimensions, $\bm \Sigma$ also contains $w$ free elements and, without using more sophisticated shrinkage techniques, the existing estimate would be prone to overfitting. As a potential remedy, we propose post-processing the estimates of the precision matrix $\bm \Sigma^{-1}$ (i.e. the inverse of $\bm \Sigma$). friedman2008sparse and, more recently, bashir2018post propose methods to ex-post sparsify precision matrices using the graphical lasso. We follow this literature and specify a loss function similar to (ref) that aims to strike a balance between model fit and parsimony.
Let $\bm \Omega$ be a sparse estimate of $\bm \Sigma^{-1}$ with elements given by $\omega_{ij}$. The loss function is then given by
with $\hat{\bm S}$ denoting an estimate of the variance-covariance matrix, $\rho_{ij}$ referring to a parameter-specific penalty and $\log \det(\bullet)$ being the log-determinant while $\text{tr}(\bullet)$ denotes the trace of a square matrix. The term $\text{tr}\left( \bm \Omega \hat{\bm S} \right) - \log \left(\det \left(\bm \Omega \right)\right)$ measures the (negative) expected fit whereas $\sum_{i \neq j} \rho_{ij} |{\omega}_{ij}|$ constitutes a penalty term that penalizes non-zero precision parameters in $\bm \Omega$. Similarly to (ref), (ref) aims to find a sparse precision matrix that describes the data well while being parsimonious.\footnote{Note that, if $\hat{\omega}^*_{ij}$ with $i \neq j$, the $(i,j)^{th}$ element of $\hat{\bm \Omega}^*$, is set to zero, then $y_{it}$ and $y_{jt}$ exhibit no contemporaneous relationship.}
Optimizing (ref) is challenging and suitable penalty parameters need to be defined. We follow friedman2008sparse in adopting the coordinate descent algorithm and state (ref) as a set of independent soft-threshold problems which can be solved for each off-diagonal element, respectively.
To determine the penalty parameter, we follow glasso and use:
with $|\hat{s}^*_{ij}|$ denoting the absolute size of the $(i,j)^{th}$ element of $\hat{\bm S}^{-1}$ and $\varpi$ is a scalar penalty parameter while $\kappa \ge 1$ controls the penalty on small precision parameters. (ref) nests the specification stipulated in bashir2018post if we set $\kappa = 1$, $\hat{s}_{ij}$ to an initial estimate of the $(i,j)^{th}$ element of the precision matrix, and cross-validate $\varpi$.
It is worth discussing a promising alternative approach to regularization of precision matrices. One could also regularize $\bm \Sigma^{-1}$ by stating it as a set of nodewise regressions \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{meinshausen2006high}. Exploiting the triangular decomposition of the precision matrix one can treat each node as an independent lasso problem and use the other endogenous variables as covariates. This strategy would imply that one replaces the optimization problem in (ref) by a set of $m$ independent (node-specific) problems as outlined in (ref). As noted by friedman2008sparse and banerjee2008model, however, this approach is a special case of the graphical lasso and thus closely related.
Before discussing our posterior simulation algorithm, it is worth noting that up to this point, the different sparsification techniques have been proposed such that some estimate (i.e. the posterior mean/median) is used and then ex-post sparsified. This technique provides a sparse point estimator of $\bm a$ and $\bm \Sigma$ but is not capable of controlling for posterior uncertainty conditional on zeroing out the $v_{ij}$'s.
Following huber2019inducing, we sparsify each draw from the joint posterior distribution of $\bm a$ and $\bm \Sigma$. Drawing from the joint posterior is easily achieved. Let $\bm a^{(r)}$ and $\bm \Sigma^{(r)}$ denote the $r^{th}$ draw from the posterior, then we first sample $\bm \Sigma^{(r)}$ from its marginal posterior distribution ((ref)) and, conditional on this draw, we sample $\bm a^{(r)}$ from ((ref)). Given this pair of draws, the corresponding loss functions becomes:
Equations ((ref)) and ((ref)) indicate that we search for an optimal action that minimizes the loss (instead of the expected loss) for each draw. This guarantees that the corresponding sparse estimates associated with the $r^{th}$ draw of $\bm {\hat a}^{* (r)}$ and $\hat{\bm \Omega}^{* (r)}$ are optimal. To be consistent with the definition of the (variable-specific) penalty parameters in ((ref)) and ((ref)), we replace the corresponding point estimators with the draws of $a_j$ and $s_{ij}$.
This approach is similar in nature to woody2019model who perform (approximate) uncertainty quantification around sparse estimators. As opposed to our approach, woody2019model estimate a regression model using MCMC and then project each draw into the sparse posterior for the optimal model. This is very similar to our strategy with the main exception that we base our inferences on all sparsified MCMC draws. More precisely, while woody2019model rely on a single optimal model (selected using the posterior mean) to project the non-sparse posterior draws into the sparse regression with $q$ selected covariates, our approach allows for uncertainty about this set of $q$ regressors. In the simulation exercise, we will show that the corresponding (sparse) point estimate is close to the one of the traditional approach and thus works well empirically.
Our approach can be viewed as an approximate algorithm to draw from the joint posterior of sparsified coefficients $p(\hat{\bm a}^*, \hat{\bm \Omega}^{*}| \bm Y, \bm X )$. The optimization problem in ((ref)) is solved using the SAVS estimator. This approach, however, has been developed under the assumption that the point estimate used is the posterior mean/median. bhattacharya2018signal show that, using any of these implies that the gradient descent algorithm quickly (after one iteration) converges. Using draws from the posterior of $\bm a$ instead yields similar favorable properties of the optimization algorithm, leading to convergence after one iteration.\footnote{More precise results based on averages of the loss functions of the optimization routine are available from the first author upon request.}
As opposed to the approach proposed in hahncarvalho2015dss, our method allows for uncertainty quantification and computation of non-linear functions of the parameters such as impulse responses or higher order predictive distributions. Moreover, it allows for selecting appropriate submodels (defined through inclusion/exclusion of covariates and/or covariance relations in $\bm \Sigma$). Applying sparsification to each draw implicitly yields a sparse estimator of $\bm a$ and $\bm \Sigma$ which can be viewed as a specific restricted version of the non-sparsified model. Since we average across these different sparse estimators, we effectively average across different models. Doing this is similar to Bayesian Model Averaging. In contrast, the traditional method can be viewed as approximate Bayesian Model Selection with the shortcoming that uncertainty across models is not taken into consideration.
As mentioned in the introductory section, applying the SAVS algorithm to sparsify draws from the joint posterior during MCMC potentially implies that point estimators such as the posterior mean of $\bm \hat{\bm a}^*$ and $\hat{\bm \Sigma}^*$ are non-sparse. However, this strongly depends on the information contained in the posterior distribution; if there is significant information that a given coefficient is equal to zero, the corresponding point estimator of the sparsified coefficient could also be exactly zero.
We use a set of different data generating processes (DGPs) that vary in terms of dimension ($m \in \{3, 10, 30\}$), length of the time series ($T \in \{80, 240\})$ and whether the model is sparse, moderately sparse or dense to analyze whether sparsification improves estimation accuracy. All simulated VAR models feature five lags ($p = 5$), with the coefficient matrix $\bm A_j$ (for $j = \{1, \dots, 5\}$) being drawn from $\mathcal{N}(0, (\xi/j)^2)$. In the case of $m = 3$ we set $\xi = 0.3$ and, for stability reasons, we define $\xi = 0.2$ for $m=10$ and $\xi = 0.1$ for $m=30$. Moreover, we add $0.25$ to the diagonal elements of $\bm A_1$. To capture that higher lag orders become less important, we rescale the variance of the Gaussian by $1/j^2$ (for $j=2, \dots, 5$). Similarly, the non-zero off-diagonal elements of the lower Cholesky factor of $\bm \Sigma$ are sampled from $\mathcal{N}(0, \xi^2)$, while the elements of diag ($\bm \Sigma$) are all non-zero and set to $0.25$.
We obtain DGPs that feature different levels of sparsity by randomly setting off-diagonal elements in $\bm A_j$ (for $j = \{1, \dots, 5\}$) and in the lower Cholesky factor of $\bm \Sigma$ to zero. As stated above, we consider three levels of sparsity. The dense model features around $10\%$ zeroes in the coefficients while the moderately sparse model features around $60\%$ zeroes. Finally, we also consider an extremely sparse DGP with approximately $90\%$ zeroes. The dense DGP turns out to be a challenging case for our model. This is because it features a large number of non-zero but small coefficients (especially for $\bm A_j$ with $j>1$) which might be erroneously set equal to zero.
To assess the sensitivity of the results with respect to different choices of $\lambda$ and $\varpi$, we compute a range of sparse models and benchmark it to the non-sparse competitor. This non-sparse competitor is a Minnesota-prior BVAR with hyperparameters obtained by optimizing the marginal likelihood of the model over a grid \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see, for example,][]{carriero2019large}. Moreover, we consider the stochastic search variable selection (SSVS) prior as competitor \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{george1993variable, george1997approaches}. This model assumes a mixture of Gaussians prior to introduce sparsity but has the severe drawback of being non-conjugate and thus challenging to estimate in large dimensions. For the SSVS we follow george2008bayesian and re-scale the spike and slab component with the OLS coefficient variances denoted by $\hat{v}_{ii}$. The corresponding spike variance is then given by $0.01 \times \hat{v}_{ii}$ while the slab variance is considerably larger with $100 \times \hat{v}_{ii}$.
Finally, we add two additional competing specifications. The first one (labeled SAVS-Median) is the traditional SAVS approach stipulated in bhattacharya2018signal which sparsifies the posterior median. This model allows us to assess whether sparsifying the posterior median yields similar insights (in terms of point estimates) as our approach. The second one, labeled CDA, sparsifies each draw from the joint posterior using a standard coordinate descent algorithm (i.e. without stopping after the first iteration). This specification serves to illustrate whether using more iterations yields similar insights compared to stopping after the first iteration.
Table (ref) shows (relative) mean absolute errors (MAE) between the posterior median of the coefficients for the sparsified BVAR and the true parameter values, averaged across $150$ replications per DGP. Note that all numbers in the table feature a numerical standard error which is relatively small. Nevertheless, these findings need to be interpreted with some caution and we aim to focus on results that imply substantial differences relative to the benchmark after taking into account the simulation-induced variation. All MAEs are divided by the MAEs of non-sparse competitor.
The upper panel of the table presents the results for the VAR coefficients while the lower panel of (ref) displays the MAEs associated with the covariance parameters. In order to investigate how differing values of $\lambda$ and $\varpi$ impact estimation accuracy, we also estimate the model over a grid of values for $\lambda \in \{0.01, 0.1, 0.5, 1\}$ and set $\varpi = \lambda/10$. Before proceeding, it is worth noting that we estimate all VAR models with five lags.
Considering the upper panel of (ref), a few results are worth emphasizing. First, we observe that sparsification pays off in terms of achieving lower estimation errors. This improvement strongly depends on the true level of sparsity, with strong accuracy gains if the DGP is very sparse and the data sample is short. Especially when $T$ is small relative to the number of parameters, sparsification improves against the traditional Bayesian VAR model.
Second, the sensitivity of estimation accuracy with respect to $\lambda$ varies with the level of sparsity. For instance, we find slightly more pronounced differences if the DGP is either dense or moderately dense but, as long as $\lambda$ is set greater to $0.01$ we find only small differences in MAEs. It is worth noting, however, that these small differences across different penalty terms are often insignificant. Once we take into account numerical standard errors the specific choice of $\lambda$ (as long as it is not set too small) seems to play a minor role with differences being smaller than ten percent in MAE terms (and thus often within one standard deviation) for most models considered. This also provides some evidence that the specific choice proposed in bhattacharya2018signal (i.e. $\lambda =1$) works well in most circumstances.
Third, for large models we find that the SAVS estimator yields substantial gains, improving upon the shrinkage-only estimator by large margins. These improvements even arise if the DGP is characterized by relatively few zeros in the VAR coefficients. This finding is not surprising given the fact that the absolute number of zeros increases with the dimension of the parameter space and the small but negligible posterior estimates under the Minnesota BVAR have a detrimental effect on estimation accuracy.
Fourth, we find that the conjugate VARs in combination with SAVS very often improve upon the non-conjugate and more flexible VAR coupled with the SSVS prior. For large-scale models, we even find that the SSVS prior performs rather poorly and this might be caused by mixing issues in the indicators that determine the mixture of Gaussian to be chosen. In smaller-sized models, accuracy differences decrease but still favor our proposed sparsified model.
Finally, and importantly, comparing the performance of the SAVS - Median and CDA approaches with the corresponding model based on setting $\lambda=1$ reveals no differences in estimation accuracy. This remarkable result shows that sparsifying each draw from the joint posterior yields (almost) identical sparsified point estimators and provides evidence that using only a single iteration of the coordinate descent algorithm seems to be sufficient when compared to using a stopping rule (and thus potentially many more iterations).
The lower panel of (ref) provides similar but more mixed insights. Sparsification of the variance-covariance matrix sometimes yields accuracy improvements over its non-sparsified counterpart. These improvements range from being small (or in some rare cases even negative) to very large (in the case the DGP is sparse and the model is moderately large). We conjecture that the somewhat smaller improvements in predictive accuracy arise from the Wishart-distributed prior imposed on $\bm \Sigma^{-1}$. This prior, by construction, is not capable of discriminating between relevant and irrelevant covariance parameters since it uniformly pushes the posterior estimate of $\bm \Sigma^{-1}$ towards a diagonal matrix. Again, we find no discernible differences between SAVS-median and CDA to the respective model with $\lambda=1$.
Considering the performance of the SSVS prior shows that it provides more accurate estimates of $\bm \Sigma$ but at substantially higher computational costs. In most instances where the SSVS prior improves upon our SAVS-based model, these performance gains are often small (i.e. below ten to 15 percent in MAE terms). So given that the additional costs of applying SAVS to the posterior draws of $\bm \Sigma$ are negligible, we can recommend adding this additional step to further improve estimation accuracy.
We stressed one key advantage in (ref), namely that sparsification reduces estimation uncertainty by zeroing out the coefficient under scrutiny during posterior simulation. Thus, while the discussion in the previous paragraphs highlights that using sparsification improves estimation performance in terms of point estimators, we now investigate its consequences on the posterior variance of $\bm a$. (ref)(a) and (b) are heatmaps that show the absolute distance between the posterior median and the true coefficients (left panel) as well as a corresponding heatmap that presents the corresponding posterior standard deviation of the parameter under scrutiny (right panel). These heatmaps are created for a single realization from the sparse DGP with $T= 240$ and $m = 30$.
Considering the heatmaps reveals that sparsification improves estimation accuracy and accurately detects zeroes, as evidenced by the abundance of white cells in the figure. The slight bias along the main diagonal (which also exists under the shrinkage-only model) stems from the informative prior that is centered around zero. However, note that even with a high degree of shrinkage, the corresponding estimate of $\bm a$ with the Minnesota prior is quite dense. By contrast, applying SAVS yields a very sparse coefficient matrix and, in addition, a sparsified estimate of the variance-covariance specification.
The left panel of (ref) shows that sparsity is, not surprisingly, accompanied by appreciable decreases in posterior variance. White cells imply that the posterior standard deviation is (close to) zero. This is often the case if we post-process the posterior using SAVS. As expected, we see that the coefficients associated with the first lag often feature considerable posterior uncertainty. For higher lag orders, however, SAVS effectively reduces posterior uncertainty and, in light of the strong accuracy gains for the point estimator, leads to a much more favorable bias-variance relationship. For a standard Minnesota prior without SAVS, we observe a more dense heatmap with a great deal of purple shaded areas. We would like to stress that even for the Minnesota prior these standard deviations are often small (especially compared to some weakly informative prior). But these small elements in $v_{ij}$ could add up (see (ref)) and thus be deleterious for predictive accuracy. And it is precisely this problem which we try to circumvent by applying SAVS.
On the variance-covariance matrix we, again, see more white cells under the sparse model. This shows that we reduce estimation uncertainty. Considering the right panel of (ref)(b), moreover, implies that we also successfully reduce posterior uncertainty around the free elements in $\bm \Sigma$.
To sum up, our simulation exercise shows that using SAVS almost always improves estimates of $\bm a$ and $\bm \Sigma$. These accuracy gains increase strongly with the degree of sparsity in the DGP. At a first glance, we find some differences across the competing values of $\lambda$. As a general recommendation, we advise to stick to the standard setup proposed in bhattacharya2018signal and set $\lambda = 1$. This choice works reasonably well and is always close to the optimal value of $\lambda$ (and it often coincides with it).
In simulations, we have shown that our approach yields more precise parameter estimates if the DGP is sparse. To illustrate the merits of our approach for forecasters in central banks and other policy institutions, we now carry out a US macroeconomic forecasting exercise. Our empirical work relies on the dataset discussed in mccracken2016fred. Several recent papers have analyzed this dataset using various shrinkage and sparsification techniques \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{giannone2017economic, cross2020macroeconomic} and find mixed evidence for sparsity. Nevertheless, using this application we aim to illustrate that even when there is little evidence in favor of sparsity, the proposed framework is still capable of improving upon a model that relies solely on shrinkage priors.
In this application we use the quarterly variant of the mccracken2016fred dataset that spans the period from $1959$:Q$1$ up to $2018$:Q$4$. To investigate whether combining shrinkage and sparsity improves predictive performance, we split the sample in two parts. The first part, called the training sample, runs from $1959$:Q$1$ to $1989$:Q$4$. This initial period is used to compute the $h$-step-ahead predictive distribution (for $h \in \{1, 4, 8\}$). After obtaining the predictive density for all $h$, we expand the initial estimation period by one quarter (i.e. to $1990:$Q$1$) and repeat this procedure until we reach the penultimate point in the sample (i.e. $2018$:Q$3$). The period from $1990$:Q$1$ to $2018$:Q$4$ is consequently labeled the hold-out or verification period. For each quarter in the hold-out period, we evaluate the predictive accuracy of the models using root-mean-squared forecast errors (RMSEs), log predictive likelihoods (LPLs), and normalized forecast errors.\footnote{As alternative metric, we also computed compute continuous rank probability scores (CRPS). Compared to the LPLs, the CRPS yield qualitatively similar insights. For the sake of brevity, we focus on LPLs in the paper. The CRPS are available from the first author upon request.}
The existing literature highlights the necessity to exploit large information sets \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see, for instance,][]{banbura2010large, carriero2015bayesian, giannone2015prior, koop2013forecasting}. Building on this evidence, we apply our techniques to a VAR model that features $m=165$ macroeconomic and financial variables. Out of these, we select three traditional target variables, namely output (GDPC1), consumer price inflation (CPIAUCSL) and the Federal Funds Rate (FEDFUNDS). Since the computational burden of using non-conjugate priors increases dramatically with model size, we do not use the SSVS prior for the large dataset.
Apart from this large-scale VAR (L-VAR), we investigate how our techniques perform across different model sizes and dimension-reduction techniques. These competing models range from small and medium-scale VARs to dynamic factor models in the spirit of bernanke2005measuring. These competing approaches are:\footnote{In Appendix (ref) we show a detailed list of the variables included along the transformation codes.}
All models feature $p=5$ lags of the endogenous variables.
Since we use VAR models that do not only differ in the size of their information sets, but also how these information is used during estimation, careful choice of the prior hyperparameters is necessary. Using the prior outlined above, we need to set $\theta_1$ as well as the sparsification parameters $\lambda$ and $\varpi$.
To set the hyperparameters of the Minnesota prior, we follow two different routes. The first (and simplest) way, is to set $\theta_1$ such that shrinkage increases with the size of the information set \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see, for instance,][]{banbura2010large, koop2013forecasting} and select $\theta_1$ on a grid of potential values.
The second approach is based on optimizing the marginal likelihood (which is available in closed form) a priori \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see][]{carriero2015bayesian}. One problem with this strategy, however, is that the marginal likelihood might be ill-behaved, which renders optimization difficult.\footnote{To circumvent this issue banbura2010large and koop2013forecasting, for example, define a training sample which serves the purpose to calibrate $\theta_1$ by minimizing the distance of the mean square error (MSE) of a large-scale model and a three-variable VAR estimated with OLS. Intuitively, this strategy implies that large dimensional models are shrunk to a larger degree than smaller-scale models \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{de2008forecasting}.}
Optimizing the marginal likelihood in huge models is often unfeasible due to numerical reasons. As a solution, we assess the sensitivity of the forecasts with respect to three choices of the hyperparameters $\theta_1 = \{0.025, 0.05, 0.075\}$. These three values all lead to an informative prior but the weight placed on prior information ranges from being large ($\theta_1 = 0.025$) to being more moderate ($\theta_1 = 0.075$). Using these values enables us to assess how shrinkage and sparsification interact. For example, if $\theta_1$is set to $0.025$ it is very likely that the SAVS estimator will lead to a sparse model while a shrinkage parameter $\theta_1=0.075$ allows for larger elements in $\bm a$ and thus a more dense model under the SAVS estimator.
For the small- and medium-scale models and the FA-VAR specification, we define a large grid of values for $\theta_1 \in \{0.01, 0.025, 0.050, 0.075, 0.10, 0.125, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50, 0.75, 1, 2, 5\}$. Using this grid, we seek the value of $\theta_1$ that maximizes the marginal likelihood. Over the hold-out sample, this procedure yields an average of $\theta_1 = 0.457$ for the small-scale model, $\theta_1 = 0.392$ for the FA-VAR specification and $\theta_1 = 0.149$ for the medium-scale model. These values indicate that the larger the model becomes, the more weight needs to be placed on the prior. Similar to carriero2015bayesian, we find that the hyperparameters tend to display little variation over the hold-out period. For example, in the case of the medium-scale model, we find that $\theta_1$ ranges from $0.125$ to $0.15$.
Finally, we investigate how forecasting performance changes for different values of $\lambda$ and $\varpi$, again, using a grid of candidate values. More precisely, we set $\lambda \in \{0.01, 0.1, 0.5, 1\}$ and $\varpi = \lambda/10$.
In this sub-section, we start by considering point forecasting accuracy of the different models and by comparing sparse with non-sparse models. Table (ref) depicts the relative RMSEs to a small-scale VAR with a Minnesota prior (henceforth called the benchmark model) and without sparsification for the three target variables. The asterisks indicate statistical significance for each model relative to the benchmark as measured by the Diebold-Mariano (DM) test \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{diebold1995comparing}. The numbers in parentheses refer to the ranking of the three best specifications using the procedure outlined in hansen2011model.
We start by considering the average model performance (in terms of computing average RMSEs across the three focus variables) for the one-step-ahead forecasts. In general, there exists no single superior modeling approach that outperforms its competitors in a statistically significant manner. However, when it comes to lowest RMSEs, our results suggest that the most accurate one-step-ahead forecasts can be found for smaller-sized models with sparsification (with the S-VAR with $\lambda = 0.5$ yielding the lowest average RMSEs). Irrespective of $\lambda$, small-scale models seem to perform well and are the only ones that outperform the benchmark Minnesota BVAR.
To analyze point forecasts in more detail, we now discuss the one-step-ahead forecasting performance across the three focus variables. For output and one-step-ahead forecasts, we observe that most large models display relative RMSEs below or close to one, suggesting that their predictions tend to improve forecast accuracy compared to the ones obtained from using the simple small-scale Bayesian VAR. When we compare sparsified large-scale BVARs with the non-sparsified counterparts we find limited evidence that sparsification improves point forecasts (but it also never substantially hurts predictive accuracy). Nevertheless, we find that sparsified large-scale BVARs with $\theta_1 = \{0.05, 0.075\}$ and $\lambda = 0.01$ improve upon the benchmark by around eleven to twelve percent in RMSE terms.
Turning to the medium-scale models, we find more pronounced improvements relative to the benchmark model and that using SAVS improves predictive accuracy. These performance gains depend on the specific value of $\lambda$. For example, setting $\lambda = 0.01$ yields a model that improves upon all competitors and produces statistically significant better forecasts than the non-sparsified benchmark model (at the five percent significance level). The second best performing model is the medium-scale VAR that sets $\lambda = 0.1$ (with differences between $\lambda = 0.01$ and $\lambda = 0.1$ being small). Finally, if $\lambda$ is set too large (i.e. equal to $0.5$ or $1$) we see that forecast quality tends to converge towards the benchmark model.
Although the M-VAR with a Minnesota prior tends to outperform the FA-VAR (irrespective of the choice of $\lambda$), a reverse picture emerges when considering the SSVS prior. With the SSVS prior a FA-VAR outperforms its medium-scale counterpart. For the SSVS model we would expect such an outcome, since the forecast performance of the spike-and-slab prior commonly deteriorates with model size \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see, for instance,][]{koop2013forecasting}.
Considering the one-step-ahead forecasts of inflation, it appears that none of the larger-scale competitors is capable of improving upon the benchmark small-scale VAR. Notice, however, that for the medium-sized model and the FA-VAR specification, introducing sparsity through our SAVS estimator sometimes yields more accurate short-run predictions and tends to outperform not only the non-sparsified model, but also the SSVS prior. In fact, a FA-VAR and M-VAR equipped with an SSVS prior produces quite poor one-quarter-ahead inflation forecasts. For the S-VAR, some sparsification slightly improves forecast accuracy. Setting $\lambda = 0.01$ and thus introducing a small penalty on non-zero elements in $\hat{\bm a}^*$ yields the most precise inflation forecasts. In general, and similarly to output, we find relative RMSEs close to unity (and in fact often slightly exceeding unity). This is consistent with the literature on inflation and mainly driven by the high persistence of inflation.
The one-quarter-ahead interest rate forecasts reveals two insights. First, we find that medium- to large-sized models yield forecasts that tend to be more accurate than the ones obtained from the benchmark VAR. This result is particularly pronounced for the moderately-sized VAR and $\lambda = 0.5$. Second, comparing sparse with non-sparse models shows that for interest rates, introducing sparsity pays off markedly. This is even visible for the small-scale model and $\lambda=1$ (ranked second), with RMSEs being over $26$ percent lower as compared to the shrinkage-only case. The largest gains, as expected, can be obtained in situations with increased model size. For the medium and large models, the increases in accuracy after sparsification are substantial and statistically significant. And these display little differences across the values of $\lambda$. We conjecture that some of the increases in predictive performance are driven by the zero lower bound, which is prominently featured in our hold-out sample. Methods that shrink but do not sparsify yield forecasts that are non-zero. This implies larger forecasts errors as compared to models that set coefficients in the interest rate equation to zero and thus predict that short-term interest rates to be zero in the next period.\footnote{A substantially lower RMSE error of the non-sparsified large-scale model with a relatively tight Minnesota prior ($\theta_1 = 0.025$) compared to the non-sparsified L-VAR with $\theta =\{0.05, 0.075\}$ underpins this observation.} This is the main source of the strong accuracy gains obtained by applying the SAVS step.
One important limitation of our forecasting exercise is that we rely on using final vintage data. This essentially implies that, for each point in the training sample, we use information that is not available at the time the forecast is produced. This effect is most pronounced in the case of the one-step-ahead forecast. Once we increase the forecast horizon, the corresponding information advantage becomes smaller (since we use information up to time $T$ to predict $\bm y_{T+h}$ for $h > 1$). This is because data revisions after a few quarters are typically very close to final vintage data.
The discussion above suggests that studying higher order forecasts in more detail is worthwhile. Next, we consider four and eight-step-ahead predictions. In summary, the results differ from the ones discussed above. While we find again mixed evidence that sparsification improves point predictions, one striking difference between one- and multi-step-ahead forecasts is that more information seems to be beneficial for accurately forecasting multiple periods ahead. Specifically, when interest centers on predicting one-year or two-year-ahead, using more information often improves predictive performance. For the one-year-ahead horizon, our findings suggest that a non-sparsified L-VAR with a tight Minnesota prior ($\theta_1 = 0.025$) yields the most accurate point forecasts. When two-year-ahead forecasts are taken under consideration, the sparse medium-scale model (with $\lambda = 0.5)$ dominates all competitors.
For output we observe small gains for the larger models. For one-year-ahead predictions, the single best performing model is the L-VAR with $\theta_1 = 0.075$ and $\lambda = 0.01$. But with close to two percent in RMSE terms, these gains are muted and not statistically significant. In general we find that the benchmark is difficult to beat and no model is capable of significantly outperforming it. This does not carry over to two-year-ahead forecasts. Here we find at least some models significantly improving (at the ten percent level) upon the benchmark by around $3.4$ to $3.2$ percent. The best performing model is the L-VAR with $\theta_1 = 0.025$ and $\lambda = 1$, closely followed by two other sparsified L-VARs (with $\theta_1 = 0.025$ and $\lambda = 0.5$ and with $\theta_1 = 0.075$ and $\lambda = 0.01$). For output, it is worth noting that we find larger gains from sparsification when two-year-ahead forecasts are being considered. In most cases, the sparsified models improve upon their dense counterparts.
Analyzing one-year-ahead inflation forecasts, again, provides limited evidence that SAVS improves forecast accuracy. While there are many cases where a given model can be slightly improved by adding the sparsification step, these differences appear to be insignificant and small. When we consider two-year-ahead inflation predictions, we find slightly more predictive evidence for sparsity (with the medium-scale VAR with $\lambda = 0.5$ performing best). For this forecast horizon, we generally observe that sparsification often yields accuracy gains which become more pronounced in larger models.
For the interest rate, the forecast gains of sparsification slightly diminish for higher-order forecasts. Generally speaking, there are little differences between all specifications, simply due to the fact that the forecasts of every stationary autoregressive model, irrespective of its size and degree of shrinkage, converge to the unconditional mean of the endogenous variables.
Overall, there is only mixed evidence that point forecasts can be improved by adding the SAVS step. There are cases where gains from SAVS become more pronounced (such as longer-run forecasts for output and inflation or short-term interest rate forecasts) but there also exist several instances where using a sparse model slightly hurts forecast accuracy. One important thing to note is that sparsification only very rarely harms predictive accuracy in a statistically significant manner.
To further substantiate this observation we consider the model confidence set (MCS) procedure of hansen2011model and implemented by bernardi2018model to obtain a measure for model uncertainty. (ref) shows the cardinality of the superior model set across variables and forecast horizons. Moreover, the table shows the number of sparsified models included in the MCS. Considering a total number of $33$ models (including the benchmark) and using a mean squared error loss function, the MCS procedure suggests a quite large superior model set, including around $26$ to $33$ models. These numbers strongly depend on the variable and forecast horizon considered. In general, only few specifications are eliminated by the MCS procedure with no single class of models significantly outperforming its alternatives. It is worth emphasizing that in most instances, the share of models using SAVS is high.
Notice, however, that focusing exclusively on point forecasts ignores what can be considered the main advantage of SAVS: the corresponding reduction in estimation uncertainty and the potentially positive effect on the full predictive distribution. This theme will be the subject of the next section.
When focusing on point forecasting accuracy we necessarily disregard any information on how well the different models capture higher order moments of the predictive distribution. However, the discussion in (ref) suggests that applying SAVS might improve density forecasts by zeroing out irrelevant predictors throughout MCMC sampling, eventually exerting a positive effect on the full predictive distribution by reducing the predictive variance.
Table (ref) shows differences in LPLs relative to the small-scale Minnesota VAR. Apart from focusing on variable-specific relative LPLs, the table also shows joint LPLs over the three target variables. These serve as a general measure on how well some approach performs in forecasting output, inflation and interest rates jointly. Numbers greater than zero imply that a given model improves upon the benchmark while negative values suggest that the benchmark yields more precise density predictions. Moreover, we compute a two-sided amisano2007comparing test for each specification with a null hypothesis of equal average LPLs relative to the benchmark and report the statistical significance based on the obtained $p$-values.\footnote{Note that the test statistics must be interpreted carefully since we use a recursive forecast design.}
We start with the joint predictive LPLs and the one-step-ahead horizon shown in Table (ref). It is worth noting that these numbers can be interpreted as a training sample marginal likelihood \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{geweke2010comparing}. For this measure, we find that the best performing model is the moderately-sized sparse VAR with $\lambda=0.1$, suggesting that using SAVS appreciably improves density predictions.
The large VARs are generally outperformed by smaller-sized ones according to one-step-ahead joint LPLs. For these models, we also find that sparsification hurts density performance, a finding that is absent when we consider the other models. For instance, in the case of the small and medium-scale VARs, we find statistical significant accuracy gains from using the additional SAVS step. This pattern can also be found in the case of the FA-VAR.
At the one-step-ahead horizon, varying $\lambda$ between $0.01$ and $0.1$ yields a similar forecasting performance for most models considered. If $\lambda$ is set too large, forecasting accuracy drops markedly. This can be traced back to the fact that a larger penalty leads to an overly sparse model and this, in turn, decreases the predictive variance too much. The resulting predictive distribution is too narrow and this makes capturing outliers increasingly difficult.
To dig into the sources on why some models work well while others perform poorly when all three variables are jointly considered, we now consider marginal LPLs. Considering the density forecasting performance for output, we observe that the model which does well for the joint predictive BF also excels (i.e. the medium-size VAR with $\lambda =0.01$). However, considering the large VARs for output alone reveals that they are also highly competitive and close to the single best performing specification (as long as $\lambda$ is set not too large) and that sparse models often outperform the non-sparse competitor.
Turning to inflation forecasts highlights that they are the main driver of the bad performance of most large-scale models. Irrespective of the choice of $\lambda$ and $\theta_1$, the small-scale VAR outperforms each specification and it seems that using SAVS only reduces predictive performance for inflation in large datasets. This is result is well-known in central bank practice, where these specifications are commonly outperformed by small and medium-sized models that include only a selected set of endogenous variables \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[e.g.][]{giannone2015prior}.
For interest rates, the story found for inflation is reversed. We find that the best models by large margins feature large information sets and are sparse. The main reason behind this strong performance is similar to the case raised in the case of point forecasts. A setup that sets most coefficients to zero yields a predictive distribution for the interest rate which is strongly centered on zero. During the period of the zero lower bound, the corresponding predictive distribution will feature a mean/median close to zero with a rather small variance and this yields large predictive gains in terms of LPLs.
After discussing all three marginal LPLs we find that large models yield competitive output and very strong interest rate forecasts while they poorly predict inflation. Hence, if little weight is placed on inflation predictions the large VARs appear to be competitive.
Inspecting multi-steps-ahead density forecasting performance yields similar insights to the one-step-ahead case. Considering relative joint LPLs shows that sparse medium-scale models perform best. The accuracy gains from using SAVS are especially pronounced for this model size. Interestingly, and in contrast to the one-step-ahead case, we find that larger values of $\lambda$ (i.e. $\lambda = 1$ for one-year-ahead and $\lambda = 0.5$ for two-year-ahead forecasts) translate into the largest gains. Again, we find that large VARs are generally beaten by medium-sized models but, as opposed to one-quarter-ahead forecasts, we find that small values of $\lambda$ improve upon the corresponding non-sparse counterpart.
Zooming into variable-specific performance highlights that large and sparse models yield precise density forecasts of GDP that are always better than the non-sparse variant of the model under scrutiny. These gains are often substantial and sometimes significant at the ten percent level. For multi-step-ahead GDP forecasts we even find that predictive performance increases in lockstep with model size. This pattern is quite consistent with one exception. The FA-VAR shows the weakest performance among all models used. This indicates that when the researcher wishes to produce multi-step-ahead forecasts, additional information that might be ignored by introducing a factor structure seems to be important.
In terms of inflation we again see that large models do not perform well, yielding the most inaccurate forecasts across models. For multi-step-ahead inflation predictions we moreover find little evidence that applying SAVS helps to improve forecasts since in most cases, the non-sparse model performs better than the sparse variant.
Short-term interest rates are again most precisely predicted using large models with sparsification. In both cases (one and two-year-ahead) we find the same ranking according to the MCS procedure, namely that the large VAR with $\theta_1 = 0.075$ and $\lambda = 1$ performs best. Lower values of $\lambda$ yield highly competitive interest rate forecasts with little differences across the different values of $\lambda$.
Similarly to the discussion of the point forecasts we investigate these statements using MCS. (ref) shows the number of models included in the MCS when we use the negative LPL as a loss function. Starting with an initial set of $33$ models we find that, depending on the variable and forecast horizon, a large number of competitors is eliminated. The majority of surviving models is comprised of sparse specifications. While the reduction of models in the case of joint LPLs is moderate to large, we find much smaller MCSs when we focus on variable-specific forecasting performance. Especially for output and interest rates, we find that the MCS consists of a handful models. In fact, for two-year-ahead predictions, all except for a single model are removed to form the MCS. In both cases, this specification is the large VAR with $\lambda = 1$ and $\theta_1 = 0.075$. For the remaining cases, the top three of these models can be read of from (ref) and always includes several sparsified variants. For inflation, the MCS corroborate the findings discussed above. The corresponding set of superior models remains elevated and only a relatively small number of models is consequently eliminated.
This discussion shows that using SAVS typically improves GDP and interest rate forecasts in large models while it is only of limited help when predicting inflation. One caveat of considering LPLs is that they are computed over the full hold-out period and thus provide a measure of average forecast quality. But it could be the case that some models perform well during selected parts of the hold-out whereas other models excel in other periods. To investigate whether this is the case we now consider density forecasting performance over time and across two distinct periods: the period before the global financial crisis (up to 2008:Q1) and the period afterwards.
Before discussing forecasting performance over time we provide guidance on how to select the penalty parameter $\lambda$ (and thus $\varpi$). As opposed to the simulation exercise, we find that forecast accuracy as measured by LPLs differs across specific values of $\lambda$ and forecast horizons. In general, we can recommend setting $\lambda$ to a rather small value (i.e. to $0.01$ or $0.1$) if interest centers on one-step-ahead predictions. This holds for most model sizes. By contrast, if higher order forecasts are of interest, we can recommend setting $\lambda$ to a larger value (i.e. $0.5$ or $1$). With very few exceptions, these are the values that yield the highest LPLs across models.
(ref) shows the evolution of cumulative joint LPLs relative to the benchmark model for one-quarter-, one-year-, and two-year-ahead predictions over time. The figure includes standard BVARs without SAVS (dashed lines) and BVARs post-processed with an additional SAVS step (solid lines) for all considered information sets. We focus on sparse models with the $\lambda$ that maximizes the LPLs at the end of the hold-out within each model class (see (ref)). Moreover, for the large-scale BVAR we depict the evolution of BFs for the different values of $\theta_1 \in \{0.025, 0.05, 0.075\}$.
Considering the period before the financial crisis clearly shows that using SAVS often yields more precise forecast distributions. This finding is especially pronounced for higher order forecasts, where we see strong and sustained gains over the period up to the financial crisis. The performance of the sparsified large models with $\theta_1 = \{0.075, 0.05\}$ and the medium-scale BVAR with SAVS appear to be the best specifications during that time span.
When we focus attention on the global financial crisis we observe a pronounced decline in model evidence for the large model with SAVS (right panels in (ref)). During the crisis, evidence in favour of sparsification also slightly decreases for the medium-scale model. This can be seen by comparing the solid and dashed red lines in (ref). While the non-sparse, medium-sized BVAR is outperformed in the run-up to the crisis, model evidence marginally supports the non-sparse variant during the recession. The main reason why predictive accuracy of sparsified models deteriorates in turbulent periods is that the forecast error variance becomes too small and large shocks become increasingly unlikely under the predictive distribution. In such a situation, dense models often feature a larger variance because of many small but non-zero $v_{ij}'$s in (ref). This makes capturing outliers (or rapid shifts) easier and thus improves LPLs.
After the financial crisis we see that applying SAVS helps for some models (especially the medium-scale VAR). However, predictive evidence in favor of SAVS declines with the forecast horizon. This pattern is the opposite of the one observed before the financial crisis. But note that this is mostly driven by the dismal performance during the recession. For one- and two-year-ahead forecasts we observe that the slopes of most LPL curves relative to sparse models are steep, indicating a period-by-period outperformance vis-\'{a}-vis to the benchmark model. Nevertheless, while we find that the medium-scale sparse VAR improves upon all benchmarks by the end of the hold-out, this does not carry over to two-year-ahead forecasts.
This discussion highlights that during tranquil periods, which we define to be characterized by small but frequent shocks, applying the SAVS step substantially improves predictive performance. This is even more pronounced if we consider higher order forecasts. In recessions, by contrast, predictions from sparse models tend to be too conservative. And this is deleterious for density forecast performance.
In this section we investigate how sparsification impacts the calibration of predictive densities. Therefore we follow giordani2010forecasting and clark2011real in analyzing normalized forecast errors that are obtained by applying the probability integral transform (PIT).\footnote{For more details on how the normalized forecast errors are computed, see clark2011real.} If a given model is well calibrated, these normalized forecast errors should be standard normally distributed. Deviations from the standard normal distribution provide information along which dimension the model might be miss-specified.
In this discussion we focus on one-year-ahead predictive densities for the reasons outlined above. The normalized forecast errors for one-quarter- and two-year-ahead forecasts look similar, displaying no discernible differences in qualitative terms.\footnote{The normalized forecast errors for the other forecast horizons are available from the first author upon request.}
In (ref) and (ref) we compare the one-year-ahead normalized forecast errors for the setup proposed in bhattacharya2018signal (with $\lambda = 1$) to the non-sparsified BVAR model across information sets and target variables. Since this purely visual analysis might be misleading, we add a legend to each panel that provides information on whether departures from standard normality are statistically significant. The null hypothesis is that the normalized forecast errors are zero mean, feature a variance of one and no serial correlation \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{berkowitz2001testing}.\footnote{Following clark2011real, the $p$-value of a zero mean is computed with a Newey-West variance with five lags, the $p$-value of a unit variance is obtained by regressing the squared normalized forecast errors on an intercept and using a Newey-West variance with three lags, and the $p$-value of no autocorrelation is computed with an AR($1$) model, featuring an unconditional mean and heteroskedasticity-robust standard errors.}
Considering each of the three focus variables separately provides some interesting insights. Before focusing on sparsified models, we analyze the non-sparsified counterparts because these are the models that are then sparsified.
For the standard BVAR forecasts, the mean is very often not significantly different from zero. With regards to the variance we find that normalized forecast errors display variances well below one for both GDP and the interest rate, while displaying a particularly high variance for inflation forecasts. A variance of forecast errors below one indicates that for many periods in the hold-out, the predictive distribution is too spread out. As discussed above, this may be at least in part attributed to the imprecisely estimated small parameters. By contrast, a variance above one would indicate that the predictive density tends to be too tight. Moreover, we observe a high autocorrelation of forecast errors for the interest rate. This feature is also commonly found in the literature \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see, for example][]{clark2011real}.
If we apply the additional SAVS step the properties of the normalized forecast errors improve. In the case of GDP, sparsifying the BVAR typically increases the variances at a small cost of a slightly more negative mean (which is not significantly different from zero in almost all instances). In particular, sparsified large-scale models attain a variance close to one (see, e.g., (ref) L-VAR with $\theta_1 = 0.05$). The variances of the normalized forecast errors are much closer to unity (and in fact do not differ from one in a statistical significant manner).
In the case of inflation, and especially during the financial crisis, sparsification does not seem to help. The non-sparsified BVAR estimates tend to produce already too tight predictive distributions. These become even narrower after applying sparsification. It is precisely this feature that hurts forecasting performance after the large drop in inflation in $2008$:Q$4$. In $2008$:Q$4$, models are generally overly confident and, in the following period ($2009$:Q$1$), somewhat too pessimistic due to the autoregressive nature of the models. Additionally, note that forecast errors are particularly high for large-scale models.
For the interest rate, BVAR models produce normalized forecast errors with variances well below unity. Using large VARs coupled with a relative loose prior and adding SAVS improves this. More precisely, for the sparsified large VAR with $\theta_1=0.075$ we find that the variance increases from $0.16$ to approximately $0.5$. While this is still below the unit variance we would expect under a well calibrated model we consider this a substantial improvement at the costs of slightly higher autocorrelation in the errors.
In general, the PITs (for the one-year-ahead forecasts) suggest that using SAVS has the potential to improve model calibration for some models and most variables. In the case of output and interest rates, we can improve model calibration by adding SAVS. A recommendation for practitioners might be that if the variance of normalized forecast errors of a BVAR is well below one, SAVS can substantially improve forecast density calibration (see, e.g., output and interest rate forecasts of large VARs). Conversely, if the predictive variance of a BVAR is already too tight, characterized by normalized forecast errors that feature a variance above one, SAVS tends to hurt predictive accuracy (see the case of inflation).
This paper proposes methods to shrink-and-sparsify VAR models with conjugate priors. The main feature of our SAVS approach is that we post-process each draw from the joint posterior by solving an optimization problem to search for a sparse coefficient vector. Without destroying the conjugacy of the model, this approach allows for different predictors across the equations in the VAR. And, instead of pushing coefficients close to zero, our approach introduces exact zeros, removing the lower bound on accuracy one can achieve under a popular shrinkage prior in the Minnesota tradition. Since the error covariance matrix in large VARs also features a large number of coefficients, we adapt techniques from the literature on graphical models to obtain a sparse estimate of the variance-covariance matrix of the system. Using a synthetic and real-data application, we illustrate the merits of combining shrinkage and sparsification in large multivariate models.
{\fontfamily{ptm}\selectfont\setstretch{0.85} \addcontentsline{toc}{section}{References} }