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.
106,261 characters · 20 sections · 113 citation commands
Forecasting macroeconomic data with Bayesian VARs: Sparse or dense? It depends!
\def\spacingset#1{ {#1}} \spacingset{1}
\if00 {
Variants of this paper were presented at the 12th European Seminar on Bayesian Econometrics (ESOBE 2022), the 5th Vienna Workshop on High-Dimensional Time Series in Macroeconomics and Finance 2022, the 2022 ISBA World Meeting, the 16th International Conference on Computational and Financial Econometrics 2022, the Bayesian Young Statisticians Meeting (BAYSM 2023), and the 44th International Symposium on Forecasting (ISF 2024). The authors thank all participants, in particular Sylvia Fr\"uhwirth-Schnatter, Neil Shephard, Massimiliano Marcellino, Florian Huber, Sylvia Kaufmann, Hedibert Lopes, Mike West, as well as anonymous referees and the associate editor for crucially valuable comments and suggestions.
The authors acknowledge funding from the Austrian Science Fund (FWF) for the project “High-dimensional statistical learning: New methods to advance economic and sustainability policy” (ZK 35), jointly carried out by the University of Klagenfurt, Paris Lodron University Salzburg, TU Wien, and the Austrian Institute of Economic Research (WIFO).}} \affil[\empty]{Department of Statistics, University of Klagenfurt, Austria}
} \fi
\if10 {
} \fi
{\it Keywords:} density forecasting, hierarchical priors, illusion of sparsity, (semi-)global-local shrinkage, stochastic volatility
\onehalfspacing
The recent literature suggests that predictions of macroeconomic variables benefit from exploiting large data sets. Especially vector autoregressions (VARs), where the number of free parameters is relatively large compared to the limited length of macroeconomic time series, are prone to overfitting. Bayesian methods have shown to be an effective regularization technique in order to reduce estimation uncertainty by imposing additional structure on the model litterman_forecasting_1986, sims_bayesian_1998, koop_forecasting_2013, giannone_prior_2015, chan_minnesota-type_2021, griffin_luo_2022.
Prior elicitation for VARs is a long-standing issue which goes back at least to the Minnesota prior in litterman_forecasting_1986. It stipulates that lagged coefficients of the respective dependent variables are shrunken less than lagged coefficients of other variables, and that coefficients are generally shrunken stronger with an increasing number of lags. The Minnesota prior, therefore, is a classic example of a prior based on domain knowledge. On the other side of the spectrum, hierarchical priors such as stochastic search variable selection (SSVS) priors or global-local (GL) priors -- where the latter most recently attract more and more attention in combination with large VARs follett_achieving_2019, huber_adaptive_2019, kastner_sparse_2020 -- offer off-the-shelf shrinkage, since they are neither problem-specific nor designed for VARs per se.
The first contribution of this paper is the introduction of the class of semi-global priors, which aims at combining the merits of domain knowledge with off-the-shelf shrinkage. In short, this framework is characterized by shrinkage parameters acting on pre-specified subgroups of the parameter space, i.e., on the semi-global level. It can easily be applied to already existing priors, e.g., it can be applied to any GL prior or to all kinds of SSVS priors. The framework is very flexible and nests other extension to GL priors, like the extensions to the normal-gamma prior in huber_adaptive_2019 and chan_minnesota-type_2021, as special cases. We propose a grouping of the coefficients which mimics some features of the Minnesota prior, namely the discrimination between lags in general, and within lags the discrimination between own-lag and cross-lag coefficients.
Second, we provide a concise comparison of shrinkage priors that are commonly used in the VAR literature, but are scattered across different works and consequently lack a systematic comparison. Here, our attempt is to shed more light on the “prior zoo” by an in-depth comparison focusing on the two most important features of shrinkage priors:
This in turn gives guidance for practitioners in meaningfully selecting hyperparameters. The priors under scrutiny are: The Minnesota prior, a semi-hierarchical version of the Minnesota prior (SHM), the stochastic search variable selection (SSVS) prior as in george_bayesian_2008, the Horseshoe (HS) prior carvalhi2010, the normal-gamma (NG) prior griffin2010, the Dirichlet-Laplace (DL) prior bhattacharya_dirichletlaplace_2015 and the $R^2$-induced Dirichlet decomposition (R2D2) prior zhang_bayesian_2020. Concerning the latter and to the best of our knowledge, we are the first who investigate the R2D2 prior in the context of VARs.
Third, we contribute to the ongoing “Illusion of Sparsity” (IoS) debate cross_macroeconomic_2020, fava_illusion_2020, giannone_economic_2021. Loosely speaking, within the Bayesian framework, there are two contrasting approaches in order to reduce parameter uncertainty in high dimensions. On the one hand, there are sparse modeling techniques that aim at selecting a small set of important predictors. On the other hand, dense modeling techniques support the assumption that all possible explanatory variables might be important; their individual impact for prediction, however, is expected to be small. giannone_economic_2021 design one specific prior -- which is basically a discrete mixture prior with point mass at zero -- in order to detect whether specific datasets are best summarized by means of many equally important covariates or rather by a small subset of covariates. Analyzing the posterior results, they conclude that macroeconomic data is rather dense. fava_illusion_2020 show that with little changes to the prior used in giannone_economic_2021, i.e., using fatter tailed distributions, the posterior results appear to be sparser. \footnote{{
We note at this point that sparsity in the VAR coefficients might imply unrealistic constraints from a purely economic point of view. However, as detailed in simsMacroeconomicsReality1980, these can nevertheless be useful for forecasting purposes, as restrictions can reduce forecast errors even if they are false.}}
In light of these results, it becomes clear that the level of posterior sparsity can be very sensitive to prior assumptions. Thus, we augment the narrative by quantifying the sparseness of our considered priors before looking at the data. By applying the Hoyer sparseness measure hoyer_non-negative_2004 to simulated data from the prior distributions, we demonstrate that they can be ordered from sparse to dense. Sparsity is best expressed by GL priors; denseness is best expressed by Minnesota(-type) priors. In addition, we find that prior sparseness translates to posterior sparseness. In other words, the prior and the posterior sparsity ranking is similar.
Other than fava_illusion_2020 and giannone_economic_2021, who investigate a linear model with a single response variable, the flexible framework of VARs allows us to analyze results for different response variables with only one data set. In an extensive simulation study based on various data generating processes (DGPs), we test the validity of our claims. In sparse DGPs, the variants of GL priors have the highest concentration around the true model parameters. In purely dense DGPs, the Minnesota priors seem to be superior to all competing models.
The gold standard of model evaluation in economic and financial applications is comparing forecasting performance: A model is considered good if it performs well in predicting the future. In our empirical application, we adopt a variant of the quarterly data set of the US economy proposed by stock_watson_2012 and provided by mccracken_fred-qd_2020. Overall, the semi-global versions of GL priors, with specific shrinkage for own-lag and cross-lag coefficients, achieve the highest support from the data in terms of one-step-ahead predictive likelihoods. Nevertheless, we find heterogeneity in model performance over time. The heterogeneous results demonstrate that there is no prior, and hence no sparse or dense modeling approach, that performs best during the whole evaluation period and for all variables. To combine the merits of different priors, we also discuss the possibility of dynamic model averaging as proposed in raftery_online_2010 and the loss discounting framework proposed in bernaciakLossDiscountingFramework2024.
Overall, this paper provides a broad overview of specification choices with respect to Bayesian VARs in their reduced forms featuring stochastic volatility (SV). Although all of our considered priors have originally been proposed or put forward to the VAR framework using the reduced form, lately, many authors opt for using the structural form cross_macroeconomic_2020, chan_minnesota-type_2021. Placing priors on the structural form coefficients is alluring because then one can employ a faster Markov chain Monte Carlo (MCMC) algorithm. Nevertheless, placing the prior on the structural coefficients constitutes a different model where posteriors are (potentially substantially) different as well. Beyond that, there is early empirical evidence that modeling reduced-form coefficients leads to better out-of-sample forecasting performance bernardi2022. Modeling reduced-form coefficients, however, comes at the cost of a higher computational burden. To render MCMC computation feasible, we apply the corrected triangular algorithm as in carriero_corrigendum_2021.\footnote{carriero_corrigendum_2021 present a correction to the broadly applied algorithm based on equation per equation estimation put forward in carriero_large_2019.}
{ cross_macroeconomic_2020 investigate and compare a similar set of prior distributions for VARs. We differ in three important aspects. First, we demonstrate how forecasting performance of GL priors can be improved with the semi-global shrinkage approach. Second, the aforementioned paper places the priors on the structural form VAR coefficients. Third, cross_macroeconomic_2020 evaluate forecasting performance only on a three-variate subset of considered variables, whereas we consider also the full-system predictive performance in addition to further robustness checks.}
The remainder of the paper is organized as follows: Section (ref) describes the econometric framework and lays out the various prior specifications, while Section (ref) compares the shrinkage behavior and sparsity of the priors. Section (ref) presents the results of an extensive simulation study considering different time series length within sparse and dense data-generating scenarios. In Section (ref), we apply the different model setups to data of the US economy. After inspecting the posterior distributions, we perform a forecasting exercise to assess the predictive performance of the considered models. Finally, Section (ref) concludes the article.
For $t=1,\dots,T$, let $\bm{y}_t$ denote an $M$-dimensional column vector containing observations on $M$ time series variables. In a $\mathrm{VAR}$ model of order $p$, $\mathrm{VAR}(p)$, $\bm{y}_t$ is determined by\footnote{For simplicity of exposition we omit the intercept in the following (which we nonetheless include in the empirical application).}
where $\bm{A}_j$ is an $M \times M$ matrix of coefficients, the lag $p$ a positive integer, and $\bm{\varepsilon}_t$ an $M$-dimensional vector of errors.
We assume the distribution of the errors to be multivariate Gaussian with time varying variance-covariance matrix, i.e., $\bm{\varepsilon}_t \sim \bm{N(0,\Sigma_t)}$, and follow cogley_drifts_2005 in applying the decomposition in the form of
where $\bm{U}^{-1}$ is an upper triangular matrix with ones on the diagonal, and $\bm{D}_t$ is a diagonal matrix. We denote the free off-diagonal elements in $\bm{U}$ as $\bm{u} = (u_{12}, u_{13}, \dots, u_{(M-1)M})^\prime$. The transformed errors $\bm{\xi}_t = \bm{U}^\prime\bm{\varepsilon}_t$ then have a diagonal variance-covariance matrix $\bm{D}_t$. Borrowing from jacquier_bayesian_1994 and kim_stochastic_1998, the orthogonalized errors are assumed to follow univariate SV models: For $i=1,\dots,M$,
where $\bm{\epsilon}_t$ and $\bm{\eta}_t$ are assumed to be i.i.d.\ $N(\bm{0, I}_M)$. Hence, the $ii$th element of $\bm{D}_t$ is $d_{ii,t}=\exp(h_{it})$. The log-variance process $\bm{h_{i}}=(h_{i1},\dots,h_{iT})^\prime$ is initialized by $h_{i0} \sim N\left(\mu_i, \frac{\sigma_i^2}{1-\rho_i^2}\right)$. Here, $\mu_i$ is the level of log-variance, $\rho_i \in (-1,1)$ the persistence of log-variance, and $\sigma_i$ the volatility of log-variance.
To facilitate prior implementation, it proves to be convenient to rewrite the model in matrix form. Define a $K \times 1$ vector of predictors $\bm{x}_t = (\bm{y}^\prime_{t-1}, \dots, \bm{y}^\prime_{t-p})^\prime$ and a $K \times M$ matrix of coefficients $\bm{\Phi}=(\bm{A}_1, \dots, \bm{A}_p)^\prime $, where $K = mp$ is the number of vectorautoregressive coefficients per equation. Then the VAR can be written as
where $\bm{Y} = (\bm{y}, \dots, \bm{y}_T)^\prime$, $\bm{X}=(\bm{x}_1,\dots,\bm{x}_T)^\prime$, and $\bm{E}=(\bm{\varepsilon}_1,\dots,\bm{\varepsilon}_T)^\prime$. $\bm{Y}$ and $\bm{E}$ are $T \times M$ matrices and $\bm{X}$ is a $T \times K$ matrix. Further, let $\bm{\phi} = \mathrm{vec}(\bm{\Phi})=(\phi_1,\dots, \phi_n)^\prime$, where $n=KM$ denotes the number of autoregressive coefficients.
As our approach to inference is Bayesian, we have to specify prior distributions. The generic prior for the coefficient vector is multivariate normal: $\bm{\phi} \sim \bm{N(0, \underline{V})}$. For growth rates and/or approximately stationary transformed data it is common to center the prior at zero george_bayesian_2008,koop_bayesian_2010,cross_macroeconomic_2020, kastner_sparse_2020, whereas for data in levels often the prior mean of own-lag coefficients in the first lag is set to one litterman_forecasting_1986, sims_bayesian_1998. Moreover, $\underline{\bm{V}}$ is a diagonal $n \times n$ matrix with diagonal elements $v_1, \dots, v_n$. The shrinkage priors under scrutiny which are to be detailed in Sections (ref) to (ref) distinguish themselves in their treatment of $\bm{\underline{V}}$.
In view of the next sections, we want to highlight that the prior for the VAR coefficients does not depend on $\bm{\Sigma}_t$ and hence neither on $\bm{U}$. This implies, that all considerations in Sections (ref) to (ref) and Sections (ref) to (ref) also apply to other specifications of the error covariance matrix, e.g.\ the order-invariant factor stochastic volatility specification without restrictions on the factor loadings proposed in kastner_sparse_2020, the order invariant stochastic volatility specification in chan_large_2021 or the order invariant specification based on the eigendecomposition of the error covariance matrix proposed in WuKoop2022.
{ Instead of modelling the reduced-form coefficients, one could also model the coefficients $\bm B=\bm{\Phi U}$ of the VAR in structural form: $\bm{y}_t^\prime \bm{U} = \bm{x}_t^\prime \bm{B} + \bm{\xi}_t^\prime$, where $\bm B = \bm \Phi U$ and $\bm \xi_t = \bm \epsilon_t \bm U \sim N(\bm 0, \bm D_t)$. Often, this is done to speed up computations, as it allows for embarrassingly parallel “equation per equation” estimation in the spirit of carriero_large_2019. Discussing whether it is better to shrink the reduced-form or the structural-form coefficients is out of scope of this paper. However, note that prior distributions are not translation invariant per se, and any comparison amongst the different settings must be treated with great care (cf.\ Appendix (ref) for a demonstration). Performance of state-of-the art shrinkage priors in the structural setting are discussed in cross_macroeconomic_2020 and chan_minnesota-type_2021.}
\paragraph{Traditional Minnesota prior} The prior proposed in litterman_forecasting_1986 -- we refer to as traditional Minnesota prior (MP_LIT) in the following -- is mainly characterized by two assumptions based on domain knowledge: First, own-lags are assumed to account for most of the variation of a given variable. Second, recent lags are assumed to be more important in predicting current values than distant lags. Hence, $\underline{\bm{V}}$ is structured in a way, such that the { diagonal elements of the blocks in $\bm{\Phi}$ corresponding to $\bm A_j$, $j=1,\dots,p$ (the entries of the own-lag coefficients), are shrunken less than the off-diagonal elements.} And, coefficients associated with more recent lags are shrunken less than the ones associated with more distant lags. We follow koop_bayesian_2010 in setting the notation: Denote $\mathbf{\underline{V}}_i$ the block of $\mathbf{\underline{V}}$ that corresponds to the $K$ coefficients in the $i$th equation, and let $\mathbf{\underline{V}}_{i,jj}$ be its diagonal elements. The diagonal elements are set to
where $\hat{\sigma}^2_{i}$ is the OLS variance of a univariate AR(6) model of the $i$th variable. We set $\lambda_1 > \lambda_2$ to regularize cross-lag coefficients more heavily. The term $r^2$ in the denominator automatically imposes more shrinkage on the coefficients towards their prior mean as lag length increases. The term $\frac{\hat{\sigma}^2_{i}}{\hat{\sigma}^2_{j}}$ adjusts not only for different scales in the data, it is also intended to account for different scales of the responses of one economic variable to another litterman_forecasting_1986.
\paragraph{Semi-hierarchical Minnesota prior} The semi-hierarchical Minnesota (SHM) prior is a generalization of the traditional Minnesota prior. Instead of fixing the shrinkage parameters $\lambda_1$ and $\lambda_2$, SHM treats them as unknown quantities and learns them in a data-based fashion. That is, the strong assumption that $\lambda_1 > \lambda_2$ is now replaced with the weaker assumption that own-lags and cross-lags might account for different amounts in the variation of a given variable. We follow huber_adaptive_2019 in imposing gamma priors on $\lambda_1$ and $\lambda_2$, i.e., $ \lambda_i \sim G(c_i,d_i) \text{ for } i=1,2, $ where $G(\alpha,\beta)$ denotes the probability density function of the gamma distribution with shape $\alpha$ and rate $\beta$.
As default choice, we set $c_i=d_i=0.01$ for $i=1,2$. The resulting gamma distribution has expectation 1 and variance 100. Moreover, the density diverges to infinity as $\lambda_j \rightarrow 0$. As most mass of this distribution is at zero, the prior can heavily shrink the parameters if necessary, though as soon as there are some moderate to large nonzero coefficients (in absolute values), the shrinkage will be only moderate to weak. Since the prior variances of cross-lag coefficients still have to be scaled with $\frac{\hat{\sigma}^2_i}{\hat{\sigma}^2_j}$, it is not a fully hierarchical prior, and we refer to it as the semi-hierarchical Minnesota prior.
Other than the Minnesota priors, global-local (GL) shrinkage priors are not built upon domain knowledge and are neither specifically designed for VARs nor for forecasting macroeconomic data. They are very popular in sparse and high-dimensional applications, and therefore also well suited for large VARs follett_achieving_2019, huber_adaptive_2019, kastner_sparse_2020. Following polson_shrink_2010, GL priors can be written as continuous scale mixture distributions. Let $K(\delta)$ denote a symmetric unimodal density with variance $\delta$, then a typical GL prior is of the following form:
where $\zeta$ represents the global shrinkage and $\vartheta_{i}$ the local shrinkage. While the global parameter determines the overall shrinkage, the local parameters act to detect the relevant signals. Hence, $g$ is usually a density with substantial mass at zero. In case of $\zeta$ being very small, $f$ must be a heavy-tailed density such that $\vartheta_i$ can override the effect of $\zeta$ for a nonzero coefficient. GL priors distinguish themselves in their choices for $f$ and $g$, which will be discussed in the next section.
In this section, we aim at building a framework that allows to put more structure on off-the-shelf shrinkage priors, but with as little restrictions as possible.
Even though many applications have shown that a simple prior like the traditional Minnesota prior -- where all hyperparameters are treated as known and therefore are not estimated from the data in a hierarchical fashion -- can be very competitive, it lacks flexibility. First, due to the fact that own-lag shrinkage and cross-lag shrinkage are fixed a priori, the traditional Minnesota prior is not (automatically) adaptable. Especially in higher dimensions, well-established default values are difficult to determine, which in turn makes selecting them hard to justify. This issue is partly dealt with by considering semi-hierarchical versions of the prior (cf.\ Section (ref)). Second, it relies on tuning the prior variances with OLS estimates in order to adapt for the different scales of the endogenous variables. Third, shrinkage of higher order lags is deterministic, i.e., prior variances are additionally scaled by $r^x$ where $r$ denotes the lag and usually $x=2$. Theoretically, GL priors could remedy those shortcomings with their local scales, which individually adjust the regularization for each coefficient. A problem arises, however, when the noise level and the absolute level of the coefficients are different in different parts of the parameter space. Then, the adjustment is probably not strong enough. In what follows, we introduce a simple framework that tackles this problem. It only needs little input from the researcher and does not require any pre-tuning of hyperparameters.
Let $\mathcal{A}_j$ denote the generic index set that labels the coefficients of the $j$th group in $\bm{\phi}$ (e.g., the first group could be the own-lag coefficients associated with the first lag, the second group could be the cross-lag coefficients associated with the first lag, etc.), and let $n_j$ denote the number of elements in $\mathcal{A}_j$. Then, as a semi-global local prior with $k$ groups, we define the following hierarchical representation:
The only additional input required is the partitioning of $\bm{\Phi}$ into $k$ subgroups. We propose a partitioning that mimics two features of the Minnesota prior, namely the distinction between own-lag and cross-lag coefficients and the distinction between lags in general: In each lag, the diagonal elements (the own-lags) and the off-diagonal elements (the cross-lags) constitute separate groups, which makes $2p$ groups in total. { The vanilla global-local prior, i.e.\ no partitioning of $\bm{\Phi}$, is nested as the special case where $k=1$.}
The semi-global framework nests other modifications to global-local priors within the VAR framework as special cases. E.g., huber_adaptive_2019 consider both, equation-specific and covariate-specific shrinkage. The first means that the covariates of each equation form $M$ separate groups (individual degree of shrinkage for each equation) and the latter that the specific covariates across all equations form $K$ separate groups (individual degree of shrinkage for each covariate across all equations). The multiplicative lag-wise specification considered in the same paper, however, is not nested in the semi-global framework, since we assume independence across the $k$ groups. The adaptive Minne\-sota prior in chan_minnesota-type_2021 could also be cast in the semi-global framework. It is the unification of the Minnesota prior with a global-local prior, namely the normal-gamma prior griffin2010. That is, there are two semi-global shrinkage parameters: One for own-lags and one for cross-lags. Additionally, each coefficient is assigned a local scale. The prior variances are rounded off with the scaling term $\frac{\hat{\sigma}^2_{i}}{r^2\hat{\sigma}^2_j}$ of the Minnesota prior in analogy to Eq. (ref). Besides using quantities estimated from the data to tune the prior variances, the adaptive Minnesota prior differs from our approach as it is a prior for the structural coefficients. { Another class of priors, which could be extended to semi-global priors are the asymmetric conjugate priors for homoskedastic VARs proposed in chan2022. Similar as with the adaptive Minnesota prior, the asymmetric conjugate priors exploit the structural form of the VAR.}
In the remainder, we briefly describe several state-of-the art off-the-shelf shrinkage priors and how they can (in one case even have to) be modified, such that the pooling strategy introduced by the semi-global framework materializes. \paragraph{DL prior} The Dirichlet-Laplace (DL) prior, introduced in bhattacharya_dirichletlaplace_2015 and put forward to the VAR framework in kastner_sparse_2020, takes the following hierarchical form:
where $DE(\delta)$ denotes the double exponential (sometimes called Laplace) distribution with variance $2\delta^2$, and $Dir(\alpha_1,\dots,\alpha_K)$ denotes the Dirichlet distribution on the open $K-1$ simplex with concentration parameters $\bm{\alpha}=(\alpha_1,\dots,\alpha_K)$. Letting $\sqrt{\vartheta_i/2} = \varrho_i \omega_j$ for $i \in \mathcal{A}_j$ and $j=1,\dots,k$, it can be shown that $\sqrt{\vartheta_i/2} \sim G(a^\gamma,1/2)$ independently for $i \in \mathcal{A}_j$, $j=1,\dots,k$ zhou_carin_2012. Using this result, it becomes clear that the proposed semi-global pooling strategy is not compatible with the vanilla DL prior, since there is no information sharing on the (semi-)global level. Therefore, to make things work, we replace the global shape parameter $a^\gamma$ with the semi-global shape parameter $a_j^\gamma$, $j=1,\dots,k$. Applying the Lemma, the semi-global-local DL prior has the following form:
where $Var(\phi_i|\vartheta_i,\zeta_j)=\vartheta_i \zeta_j$. Posterior results can be very sensitive to the choice of $a^\gamma_j$, as it is the parameter that controls the strength of regularization (cf.\ Section (ref)). Hence, we place i.i.d.\ discrete hyperpriors $Pr(a_j^\gamma = \tilde{a}^\gamma_r) = p^\gamma_r$ on $a^\gamma_j$ for $j=1,\dots,k$, where $\bm{\tilde{a}}^\gamma=(\tilde{a}^\gamma_1,\dots,\tilde{a}^\gamma_R)^\prime$ is the vector of strictly positive support points. Updates of $a_j^\gamma$ are straightforward, as the full conditional posterior is again discrete on the same support points.
Last but not least, the DL prior can be expressed as a Gaussian scale mixture by introducing the auxiliary scaling parameter $\psi_i \sim Exp(1/2)$, s.t.\ $\phi_i|\vartheta_i,\zeta_j \sim DE(\sqrt{\vartheta_i\zeta_j/2}) \Longleftrightarrow \phi_i|\psi_i,\vartheta_i,\zeta_j \sim N(0,\psi_i\vartheta_i\zeta_j/2)$. Consequently, for DL, $v_i = \psi_i\vartheta_i\zeta_j/2$. \paragraph{NG prior} The normal-gamma prior proposed by griffin2010 is nowadays well-established in different kinds of multivariate time-series modelling setups huber_adaptive_2019, kastner_sparse_2019, bitto2019achieving. It takes the following form:
i.e., for NG, $v_i=\vartheta_i \zeta_j$. As with DL, we place on $a^\delta_j$, $j=1,\dots,k$, i.i.d.\ discrete hyperpriors $Pr(a_j^\delta = \tilde{a}^\delta_r) = p^\delta_r$, where $\bm{\tilde{a}}^\delta=(\tilde{a}^\delta_1,\dots,\tilde{a}^\delta_R)^\prime$ is the vector of strictly positive support points. \paragraph{R2D2 prior} Only recently, zhang_bayesian_2020 proposed the $R^2$-induced Dirichlet decomposition (R2D2) prior. It is based on a beta prior on $R^2 \coloneqq \frac{W}{W + 1}$, where $W$ denotes the sum of all prior variances. The induced prior then takes the following form:
Apart from a higher level of hierarchy, it has a similar structure to the DL prior. A major difference is that, for R2D2, the prior variance follows a gamma distribution, whereas for DL, the prior scale follows a gamma distribution. Using the same result as before, conditional on $\xi_j$, $\varpi_i=\varrho_i\omega_j \sim G(a_\pi, \xi_j)$ independently. Exploiting the scaling property of the gamma distribution, the R2D2 prior can be cast in the form of Eq.\ (ref):
where $\zeta_j^{-1}=\frac{2\xi_j}{a_\pi}$ and $\vartheta_i = \frac{\varpi_i}{\zeta_j}$. By introducing auxiliary scaling parameters $\psi_i \sim Exp(1/2)$, the R2D2 prior can be written as a Gaussian scale mixture, s.t.\ $\phi_i|\vartheta_i,\zeta_j \sim DE\left(\sqrt{\vartheta_i\zeta_j/2}\right) \Longleftrightarrow \phi_i|\psi_i,\vartheta_i,\zeta_j \sim N\left(0,\psi_i\vartheta_i\zeta_j/2\right)$, i.e.\ for R2D2, $v_i=\psi_i\vartheta_i\zeta_j/2$.
The R2D2 prior and the NG prior are closely related. For R2D2, the variance of a double exponential distribution follows a gamma distribution, whereas for NG, the variance of a normal distribution follows a gamma distribution. The special case of NG with $c=\frac{a^\delta_j}{2}$, hence, could be seen as an $R^2$-induced Dirichlet decomposition prior with a normal kernel.
As with DL and NG, we place on $a^{\pi}_j$, $j=1,\dots,k$, i.i.d.\ discrete hyperpriors $Pr(a_j^\pi = \tilde{a}^\pi_r) = p^\pi_r$, where $\bm{\tilde{a}}^\pi=(\tilde{a}^\pi_1,\dots,\tilde{a}^\pi_R)^\prime$ is the vector of strictly positive support points. \paragraph{Horseshoe} The horseshoe (HS) prior proposed in carvalhi2010 is given by
where $\mathcal{C}^+(\cdot,\cdot)$ denotes the half-Cauchy distribution. Other than NG and R2D2, HS places the hyperprior distributions on the standard deviation, rather than on the variance. Compared to all other priors discussed so far, it has the advantage that no tuning parameter has to be specified by the user. \paragraph{Stochastic search variable selection} Another popular prior for BVARs is the stochastic search variable selection prior george_bayesian_2008, which is a discrete mixture of two normal distributions. Although this prior is neither a global-local shrinkage prior nor a continuous scale mixture prior, it can be modified in order to comply with the semi-global-local framework. In hierarchical form and adapted to the semi-global-local framework, the prior is given through
where $\gamma_i$ is an auxiliary dummy variable, $\tau_{0i}$ should be small (close to zero) and it must hold that $\tau_{0i} \ll \tau_{1i}$. Intuitively speaking, if $\phi_i$ is assigned to the spike variance $\tau^2_{0i}$, it is virtually constrained to zero. If $\phi_i$ is assigned to the slab variance $\tau^2_{1i}$, it is shrunken depending on how $\tau^2_{1i}$ is chosen. Further, we assume that the prior inclusion probability $\underline{p}_j$ might differ between the $k$ groups. In order to learn the parameter in a data-based fashion, we impose i.i.d\ beta priors on each $\underline{p}_j$, $j=1,\dots,k$: $\underline{p}_j \sim Beta(s_1,s_2)$. Note, that for prediction purposes and for inference on the joint posterior $p(\phi_i,\gamma_i|\bullet)$, learning $\underline{p}_j$ in a data-based fashion only has an effect if there is information sharing between coefficients. In the case where each $\phi_i$ constitutes a separate group, i.e., each $\phi_i$ is assigned an independent $\underline{p}_i$ cross_macroeconomic_2020, it is straightforward to show that the full conditional posterior $p(\gamma_i|\bullet)$ only depends on the prior expectation of $\underline{p}_i$.
To complete the models, we have to specify the priors on the decomposed variance-covariance matrix $\bm{\Sigma}_t = \bm{U}^{\prime -1}\bm{D}_t\bm{U}^{-1}$. Note that the precision matrix -- the inverse of the variance-covariance matrix -- can be written as $\bm{\Sigma}_t^{-1} = \bm{U D}_t^{-1} \bm{U}^\prime$. It is well-known that the precision matrix should be a sparse matrix, as zero entries on the off-diagonals imply conditional independence among the respective equations west_bayesian_2020. Hence, for $\bm{u}$, we choose the HS prior: $u_{ij} \sim N(0,\vartheta_{ij}^u\zeta^u)$, $\sqrt{\vartheta_{ij}^u} \sim \mathcal{C}^+(0,1)$, $\sqrt{\zeta^u} \sim \mathcal{C}^+(0,1)$, $i=1,\dots,M-1$ and $j=2,\dots M$.\footnote{In Section (ref), we provide some robustness checks regarding different priors for $\bm{u}$.} Following kim_stochastic_1998, for $i=1, \dots, M$, we choose a normal prior for the level of the log variance, $\mu_i \sim N(0, 100^2)$, and for the persistence parameter, $\rho_i$, a beta distribution on $\frac{\rho_i + 1}{2} \sim Beta(20, 1.5)$. Finally, we follow kastner_ancillarity-sufficiency_2014 by imposing a gamma prior on the variance of the log variance, $\sigma_i^2 \sim G(1/2, 1/2)$.
{ Regarding our prior choices for the SV parameters and related to the finding in rossiniLossbasedPriorDegrees2024, that forecasting performance of homoskedastic VARs in conjunction with a Wishart prior on $\bm \Sigma^{-1}$ can be negatively influenced by choosing an unfavorable default value for the degrees of freedom hyperparameter, we acknowledge that choosing different priors for the vector of SV parameters $(\mu_i, \rho_i, \sigma_i^2)$ can have a non-negligible effect on the posterior distribution of those parameters. However, the impact on the estimated log-variances $h_{it}$, which are most relevant for forecasting, is typically small.}
Efficient computer implementations of all posterior samplers discussed in this paper are conveniently bundled into the R R_core package bayesianVARs, interfacing to C/C++ via Rcpp rcpp and RcppArmadillo rcpparmadillo for increased computational efficiency. The package is available under the General Public License (GPL $\ge$ 3) from the Comprehensive R Archive Network (CRAN) at \url{https://cran.r-project.org/package=bayesianVARs}. We thus refrain from describing the various posterior samplers in detail. Instead, we provide their main building blocks, the necessary conditionals used for Gibbs sampling, in Appendix (ref).
{ Concerning the scalability of our approach, we present some rough estimates of computation time. The estimates are based on an Intel\textsuperscript{\textregistered} Core\textsuperscript{\texttrademark} i7-10610U processor using one core. For the dimensions of our empirical application, i.e., $M=21$, $T=198$ and $p = 1,\dots,5$, generating 100 draws from the posterior distribution takes approximately 1.9, 5.7, 11.4, 19.1, and 28.6 seconds, respectively. Since the complexity of the underlying correct triangular algorithm carriero_corrigendum_2021 is $O(M^4)$, the estimation of larger datasets up to $M\approx50$ is surely possible. For even higher dimensions, approximate estimation techniques or a different structure of the error variance-covariance matrix could be considered. Instead of using the correct triangular algorithm, one could use the equation-per-equation algorithm described in carriero_large_2019 as an approximation. In that paper, a VAR with $p=13$ lags is successfully estimated using a time series with $M=113$ and $T=648$, which is similar in dimensionality to the application in banbura_large_2010. The approach of the latter, however, builds on a substantially different model, namely a homoskedastic and conjugate VAR. Variational inference as proposed in bernardi2022, another approximate estimation technique, can also reduce estimation time of VARs which are similar in spirit to those presented in this paper. Concerning the error structure, instead of using the Cholesky stochastic volatility specification in (ref) and (ref), one could use the the factor stochastic volatility (FSV) specification proposed in kastner_sparse_2020. Their variant is also compatible with the priors for the VAR coefficients described in Sections (ref) to (ref) and Sections (ref) to (ref). By conditioning on the factors and their loadings, this approach makes conditional equation-per-equation estimation straightforward. In situations where $K \gg T$, equation-per-equation estimation can be further combined with the algorithm for high-dimensional regressions proposed in bhattacharyaFastSamplingGaussian2016. This then has complexity $O(pM^2T^2)$ and kastner_sparse_2020, e.g., demonstrate that MCMC inference for a VAR with $p=5$ lags using a time series with $M=215$ and $T=228$ becomes feasible.}
{ In Section (ref), we analyse both the tail behavior and the concentration at zero of the various shrinkage priors. In Section (ref), we characterize the sparsity of the priors.}
The limiting behavior of univariate global-local prior densities is typically analyzed conditionally on the global scale carvalhi2010, griffin2010, i.e.\ marginalized over the local scales only. Since the global parameters often introduce dependence among all $\phi_i$, $i \in \mathcal{A}_j$, $j=1,\dots,k$, we follow this line of thinking and conduct the analysis in Section (ref) in terms of univariate conditionally independent prior densities $p(\phi_i|\cdot)$, where $\cdot$ stands for all global hyperparameters.
For the proofs of Theorem (ref)(ref)((ref)) and the first part of (ref)(ref)((ref)), we refer to zhang_bayesian_2020, and for (ref)(ref)((ref)) and (ref)(ref)((ref)), we refer to carvalhi2010. The remaining proofs can be found in Appendix (ref).
Some noteworthy remarks follow. All global-local priors share a pole at $\phi_i=0$, at least for certain hyperparameter settings, which is important for handling vectors with many zero elements. The Minnesota priors and SSVS do not diverge to infinity at zero, which could be seen as a disadvantage in sparse settings. For $a_j^\gamma<1$, $a^\delta_j<\frac{1}{2}$, and $a^{\pi}_j<\frac{1}{2}$, DL, NG and R2D2, respectively, diverge to infinity with polynomial speed; this is faster than HS. Concerning the tails: HS has the heaviest (Cauchy-like) tails. It is followed by DL and R2D2, respectively, which are still heavy-tailed in the sense that they have heavier tails than the exponential distribution. NG is lighter-tailed than DL and R2D2, but still heavier-tailed than the Minnesota priors and SSVS, which have Gaussian-like tails.\footnote{A marginal note regarding R2D2: zhang_bayesian_2020 claim that R2D2 is polynomial in both regions at zero and in the tails, which seems contradictory to Theorem (ref)(ref)((ref)). For the avoidance of doubt, we want to clarify that the univariate density of R2D2 is polynomial in the tails only if the global scale $\zeta_j$ is integrated out. This is somewhat misleading, since for all other priors under scrutiny in both papers, i.e. here and in zhang_bayesian_2020, the global scale is fixed.}
Figure (ref) depicts the univariate densities conditional on the global scales and the prior densities for the global scales separately. For reasons of comparability, we show the densities for $\sqrt{\zeta}_j$ for all priors. It can easily be verified that for NG and R2D2, the gamma hyperpriors on $\zeta_j$ imply generalized gamma hyperpriors on $\sqrt{\zeta}_j$. Clearly, the division of tasks among global and local scales varies among different GL priors. Whereas the global scale of HS has considerable mass at zero, the global scale of NG and R2D2 is bounded away from zero. The global scale of DL is a Dirac-mass at one. That is, for DL, NG, and R2D2, the task of the local scales is twofold: In addition to detecting signals, they must be able to overwhelm the global scale in shrinking the zero coefficients to zero.
In order to describe the implied level of sparsity of the different priors, we now analyze their multivariate distributions. We begin by characterizing the sparsity of the priors qualitatively and continue by quantifying it via a sparseness measure.
Figure (ref) depicts three-dimensional scatter plots of prior simulations. The star-like shapes of GL priors (here the HS prior) have most of its mass where all coefficients are almost zero or where at most one is substantial different from zero, indicating sparseness. By contrast, the ball-like shape of SHM, which does not allow for outliers, indicates denseness. SSVS is somewhere in the middle of the contrasting approaches. The figure also demonstrates the effect of structured shrinkage. Assume $\phi_1, \phi_2 \in \mathcal{A}_1$ and $\phi_3 \in \mathcal{A}_2$, then for all priors there is relatively more mass where $\phi_1=\phi_2=0 \And \phi_3 \neq 0$ or else $\phi_3=0 \And \phi_1\neq 0 \And \phi_2\neq 0$, i.e., the structuring introduces discrimination between the groups. For GL priors, this implies discrimination between groups, while being sparse within groups. On the other hand, SHM only discriminates between own-lag and cross-lag coefficients. Within own-lag and within cross-lag coefficients it remains dense, which explains the spinning-top shape.
We quantify the sparseness of each prior by means of the Hoyer sparseness measure hoyer_non-negative_2004. The Hoyer measure for the vector $\bm{\phi}_j=(\phi_i)_{i \in \mathcal{A}_j}$,
is a normalized measure where $H=1$ indicates maximum sparseness, i.e., only one single nonzero value, and $H=0$ indicates maximum denseness, i.e., the absolute values of all components are equal. More importantly, it can express sparseness when all elements of $\bm{\phi}_j$ are nonzero, but when almost all are very close to zero and only very few are large (in absolute values). This is important, as we analyze continuous prior distributions which never strictly exclude single coefficients.
We estimate the Hoyer measure through Monte Carlo simulation. Table (ref) reports the arithmetic means stemming from a Monte Carlo simulation with $10,000$ iterations and simulated vectors of length $n=1000$. Interestingly, the simulation suggests that for MP_LIT/SHM, i.e., all priors of the form $\bm{\phi} \sim N(\bm{0},\bm{I} c)$, $c \sim f$, regardless the specific choice of $f$, the expectation of the Hoyer measure converges to $0.21$ for $n=1000$. This behavior of the measure is convenient, since sparsity should not depend on the overall noise level. For the remaining priors, results can be sensitive to the choice of hyperparameters; thus, for each prior, we consider two different scenarios. In scenario A, the hyperparameters for DL, R2D2, and NG are selected to ensure that the concentration at zero is comparable to the one of HS. In scenario B, we select the hyperparameters such that the concentration at zero is higher compared to scenario A: For DL we set $a^\gamma=0.001$, for R2D2 we set $a^{\pi}=0.0005$ and for NG we set $a^{\delta}=0.0005$, which results in the same behavior around the origin for those three priors. Our numerical experiments suggest, that the sparseness of DL, R2D2, and NG does not depend upon choices for the global hyperparameters. For SSVS, we set $\tau_0=0.01$ and $\tau_1=100$ in both scenarios. We set $\underline{p}=0.5$ in scenario A and $\underline{p}=0.01$ in scenario B.
The exercise corroborates our preceding considerations. From sparse to dense, the priors can be ranked in the following order: DL, R2D2, NG, SSVS, HS, and at some distance SHM/MP_LIT. Moreover, the sparseness of DL, R2D2, NG, and SSVS is very sensitive to the choices of the parameters that control the concentration at zero, namely $a^\gamma, a^\pi$, $a^\delta$, and $\underline{p}$, respectively.
This section aims at comparing the performances of the various priors on simulated data. Before that, the merits of the semi-global framework are illustrated on synthetic data.
Consider the following example of a VAR(1): Assume that all diagonal elements of $\bm{\Phi}$ are nonzero, whereas there is only one nonzero off-diagonal element. That is, taken as a whole, $\bm{\Phi}$ is not sparse. Figure (ref) demonstrates that a vanilla global-local prior -- here the HS prior -- overshrinks the diagonal elements and does not shrink the off-diagonal elements strong enough. The structured semi-global local version, however, adapts the degree of shrinkage in the two different subspaces of $\bm{\Phi}$.
In this section, we compare the performances of the different priors on simulated data. The structure of the simulation study is borrowed from kastner_sparse_2020. We consider dense and sparse data generating processes (DGPs), as well as different dimensions in $T \in \lbrace 50, 100, 200 \rbrace$. Dimensionality is $M=20$ throughout. In each setting, we sample 20 replicates of the DGP. In all scenarios, we distinguish between own- and cross-lags. The nonzero coefficients of the DGPs are drawn from Gaussians with mean $\mu_i$ and standard deviation $\sigma_i$ for $i \in \{ol,cl\}$, where $ol$ and $cl$ denote own-lag and cross-lag, respectively. In both scenarios, the probability of own-lag coefficients to be nonzero is 0.8, i.e., own-lags are always considered to be dense. The probability of cross-lag coefficients to be nonzero is 0.01 in the sparse, 0.1 in the medium, and 0.8 in the dense scenario. Further, the scenarios differ regarding the signal-to-noise ratio. In the sparse scenario, we set $\mu_{ol} = \sigma_{ol}=\mu_{cl} = \sigma_{cl}=0.3$. In the medium scenario, we set $\mu_{ol} = \sigma_{ol}=0.15$ and $\mu_{cl} = \sigma_{cl}=0.1$. In the dense scenario, we set $\mu_{ol} = \sigma_{ol}=0.15$ and $\mu_{cl} =\sigma_{cl}=0.01$. To simulate stable VAR processes, we control that none of the eigenvalues of the companion form have modulus greater than one lutkepohl_new_2005.
Elements in $\bm{u}$ are nonzero with probability 0.1. Furthermore, we set $\mu_{u} = \sigma_{u} = 0.001$. The AR(1) processes driving the log-variances of the orthogonalized errors have mean $\mu_{\sigma i} = -10$, persistencies $\rho_{\sigma i}$ in the range $[0.85,0.98]$ and standard deviations $\sigma_{\sigma i}$ in the range $[0.1,0.3]$.
The root mean squared error (RMSE) between the posterior mean of the coefficients and the true parameter values is a common measure to evaluate the performances:
where $\phi_i^{(r)}$ denotes the $r$-th posterior sample and $\hat{\phi}_i = \frac{1}{R} \sum_{r=1}^R \phi_i^{(r)}$ is the posterior mean of the $i$-th VAR coefficient.
The following prior specifications are considered: DL, NG, and R2D2 denote the vanilla global-local implementations with the concentration parameters being fixed. In DL, $a^\gamma=\frac{1}{K}$, the value proposed in kastner_sparse_2020. In NG and R2D2, $a^\delta=a^\pi=\frac{1}{2K}$, i.e.\ DL, NG, and R2D2 are specified to have the same concentration at zero. Further, in R2D2, $b^\pi=0.5$, the value proposed in zhang_bayesian_2020. In NG, $b^\delta=0.5$ and $c=\frac{a^\delta}{2}$, such that NG and R2D2 differ only in their respective kernels: Whereas NG defines hyperpriors for the variance of a normal distribution, R2D2 defines hyperpriors for the variance of a double-exponential distribution. { In DL$_a$, NG$_a$, and R2D2$_a$, the subscript “a” indicates that the following discrete prior is placed on the corresponding concentration parameter:} The support points of the discrete priors $\tilde{\bm{a}}^\gamma=\tilde{\bm{a}}^\delta=\tilde{\bm{a}}^\pi=\tilde{\bm{a}}=(\tilde{a}_1=\frac{1}{1000},\dots,\tilde{a}_{1000}=1)^\prime$ are equally spaced, and the probabilities are proportional to the density of an exponential distribution with rate $1$. In MP_LIT, $\lambda_1=0.16$ and $\lambda_2=0.004$. In SHM, $c_i=d_i=0.01$ for $i=1,2$. For all SSVS priors, we closely follow the semi-automatic approach in george_bayesian_2008 by setting $\tau_{0i} = \frac{1}{100} \sqrt{\widehat{var}(\phi_i)}$ and $\tau_{1i}=100 \sqrt{\widehat{var}(\phi_i)}$. Here, $\widehat{var}(\phi_i)$ denotes the variance of the posterior distribution resulting from a conjugate flat normal-Wishart prior, which is available in closed form. george_bayesian_2008 use the variance of the ordinary least squares (OLS) estimator. However, in situations where the number of covariates per equation exceeds the number of observations, the OLS estimator might not exist. { For SSVS, the prior inclusion probability is fixed at $\underline{p}=0.5$, whereas for SSVS$_p$, the subscript “$p$” indicates that a $\text{Beta}(1,1)$ is placed on $\underline{p}$}. A superscript “$^*$” following a prior indicates a semi-global modification, e.g., the semi-global-local HS prior is denoted by HS$^*$. {
Sub- and superscripts also can be combined: For instance, NG$_a^*$ stands for the semi-global-local NG prior where the concentration parameters $a^\delta_j$, $j=1,\dots,2p$, are treated hierarchically, whereas NG$^*$ stands for the semi-global-local NG prior where the concentration parameters $a^\delta_j=\frac{1}{2K}$, $j=1,\dots,2p$, are fixed}.
Table (ref) reveals the results of the simulation study. HS$^*$ has the lowest RMSE in six out of nine scenarios: It is best in all sparse scenarios, in the medium scenarios with at least 100 observations and in the dense scenarios with only 50 observations. In the dense scenarios with at least 100 observations, SHM has the lowest scores. Concerning GL priors, the semi-global versions perform better than the vanilla versions in all scenarios. In the sparse and the medium scenarios, GL priors in general, i.e., the vanilla and the semi-global versions, tend to perform better than the Minnesota priors (with one notable exception). In the dense scenarios, SHM is overall best, closely followed by HS$^*$. SHM is considerably better than MP_LIT in the dense scenarios with at least 100 observations, whereas in the remaining scenarios their performances in terms of RMSEs is comparable. Turning to the various specifications of SSVS, we tend to observe competitive performance only in the sparse scenarios.
In Section (ref), we shortly summarize the data set, present the model specifications and the forecasting design. In Section (ref), we inspect the posterior distributions of VAR coefficients arising from different prior distributions. In Section (ref), we assess out-of-sample forecasting accuracy of the different models. In Section (ref), we demonstrate the merits of combining forecasts of different models in a dynamic fashion via dynamic model averaging (DMA).
The aim of the empirical application is to forecast economically relevant US time series. We use the quarterly data set provided by mccracken_fred-qd_2020, which is based on the well-known stock_watson_2012 data set. The sample period ranges from 1959:Q4 to 2020:Q1. All in all, we include $M=21$ quarterly time series with the intention to cover the most important segments of the US economy. The complete list of all variables used in the following illustration is provided in Table (ref) (Appendix (ref)).
{ Similar to huber_adaptive_2019, cross_macroeconomic_2020, and chan_minnesota-type_2021, most variables enter the models in log differences to be interpreted as growth rates, except interest rates, which are already defined in rates and hence taken in levels. In our experience, this choice results in more stable predictive densities compared to data taken in (log-)levels. Another stream of literature deals explicitly with the problem of potentially non-stationary data by eliciting priors accordingly sims_bayesian_1998,villaniSteadystatePriorsVector2009, giannonePriorsLongRun2019. It has to be acknowledged that results about the performances of the priors might be different when considering different transformations. A preliminary analysis available upon request, however, suggests that many of the results are qualitatively similar when the data is taken in log-levels.}
A very important and often underemphasized choice when modeling macroeconomic time series with VARs is the number of lags of endogenous variables included as predictors. In the Bayesian framework, fixing $p$ at a certain number implies prior information that higher lags do not carry any important information. Hence, litterman_forecasting_1986 suggests estimating as many lags as is computationally feasible, in which the prior has to be designed such that irrelevant coefficients get shrunken to zero. However, the more lags we include, the severer the problem of overfitting becomes, and the harder it will be to regularize the parameter space. Not to mention the increase of computational burden. Typical choices for quarterly data in similar dimensions are either four or five lags chan_minnesota-type_2021,cross_macroeconomic_2020, huber_adaptive_2019, giannone_prior_2015. We estimate VAR($p$)s for $p\in\lbrace1,2,3,4,5\rbrace$ to investigate to what degree longer lags improve/worsen forecasting performance.
The forecasting performance is evaluated by means of a recursive pseudo out-of-sample forecasting exercise. Based on the initial estimation window ranging from 1959:Q4 to 1980:Q1, one-step-ahead and { four-steps-ahead} predictive densities are evaluated. After that, the estimation window is expanded by one quarter and the model is re-estimated. For each estimation step, we run ten independent MCMC chains. Per chain, we keep 15,000 draws after a burn-in of 5,000 iterations. This procedure is consequently repeated until the end of the sample 2020:Q1 is reached.
The prior specifications we consider in the empirical application are the same as in Section (ref). { In what is to follow, we also include an intercept term. The prior for the intercept is independent normal centered at zero with variance 1000 for all models.}
Before presenting the results of the forecasting exercise, it is worth analyzing the posterior distributions arising from the different priors. Figure (ref) depicts posterior medians and posterior interquartile ranges of different VAR(2) models.\footnote{Although we find signals with higher lag orders, the out-of-sample results indicate that two lags usually are enough (cf.\ Section (ref)).} All models, even the vanilla SSVS and GL priors, detect many signals associated with the first own-lag and little signals in the reaming part of the coefficient space. Concerning the cross-lag coefficients, the posteriors arising from the Minnesota priors show relatively high uncertainty, whereas the remaining priors only find sparse signals. The differences between the vanilla and the semi-global versions of SSVS and GL priors are hardly detectable by visual exploration only.
Table (ref) displays the posterior mean of the Hoyer measure.\footnote{We compute the Hoyer measure for every single coefficient posterior draw, which then provides us with the posterior distribution of the measure.} In the first lag, the posteriors under all priors are dense, whereas for the Minnesota priors own-lags in general are dense. Moreover, the Minnesota priors are always denser than the remaining priors for all considered lag-lengths with respect to both own-lag and cross-lag coefficients. Last but not least, especially on the off-diagonals, the posteriors under all priors except SHM and MP_LIT can get extremely sparse. The heterogeneous results regarding sparsity show that posterior inference is strongly affected by the prior assumptions, which stresses the importance of careful prior elicitation. { This corroborates findings in jarocinski2017: In order to answer the question whether one variable is relevant for another variable, the authors derive a closed-form expression for the posterior probability of Granger noncausality and a generalization thereof in a homoskedastic VAR with a conjugate prior. They find that the posterior probability for zero restrictions increases as the prior gets looser. Since prior choices matter,} we think it can be delusive to draw general conclusions about sparseness mainly by analyzing posterior quantities. fava_illusion_2020 elaborate on how such posterior results often just spread illusions.
Comparing the sparseness between GL priors, we observe that the own-lag coefficients of the semi-global modifications are denser compared to their vanilla GL counterparts. HS$^*$ discriminates the strongest between own-lags and cross-lags in terms of sparsity, especially in more distant lags. In this context, it is worth analyzing the posterior distributions of the (semi-)global hyperparameters. Figure (ref) depicts, exemplary for a VAR(4), box plots of the posterior distribution of $\sqrt{\zeta}_j$ for HS$^*$, NG$^*_a$, and R2D2$^*_a$, and the posterior distribution of $a_j^\gamma$ for DL$^*_a$, $a_j^\delta$ for NG$^*_a$, $a_j^\pi$ for R2D2$^*_a$, $j \in \{\text{lag}: 1,\dots,4\} \times \{\text{own-lag, cross-lag}\}$. For all models there is a clear discrimination in all lags between own-lag and cross-lag, i.e., shrinkage is always stronger for cross-lag coefficients than for own-lag coefficients in all lags. The strength of discrimination between own-lag and cross-lag decreases with higher lags, and is the strongest for HS. Last but not least, the posterior distributions of the semi-global hyperparameters concentrate around smaller values for more distant lags than for closer lags, i.e., shrinkage increases with increasing lags.
{ Following geweke_comparing_2010 we evaluate forecasting performance in terms of log predictive likelihoods, treating the period ranging from 1959:Q4 to 1980:Q1 as training sample (initial estimation window). Evaluation of the predictive densities begins in 1980:Q2 or 1981:Q1 for one-step-ahead forecasts and four-steps-ahead forecasts, respectively. The estimation window is recursively expanded by one observation until the end of the evaluation period 2020:Q1 is reached.}
Table (ref) depicts sums of one-step- { and four-steps}-ahead log predictive likelihoods relative to the traditional Minnesota prior with four lags (MP_LIT(4)). { We start by inspecting the one-step-ahead scores.} The most important findings are the following: The overall best performance is achieved by the semi-global modification of HS with 4 lags (HS$^*$(4)). In almost all cases, the semi-global modifications clearly improve forecasting performance compared to their vanilla GL counterparts. Moreover, the semi-global modifications usually perform better than the Minnesota priors, whereas the vanilla GL priors do not. { SHM priors perform well in general.} Comparing the scores of SHM and MP_LIT reveals that learning the shrinkage parameters of the Minnesota prior from the data always improves forecasting performance, which corroborates findings in giannone_prior_2015 and cross_macroeconomic_2020. Among all considered priors, SSVS priors achieve the lowest scores.
Concerning the performance of GL priors compared to Minnesota priors: The vanilla GL priors do not produce more accurate forecasts than SHM. This contradicts results in huber_adaptive_2019 at first sight. The weaker performance of SHM in huber_adaptive_2019 can be explained (at least to some extent) by the vague prior imposed on $\bm{u}$, a normal prior with zero mean and variance 10. The fact that the prior choice for $\bm{u}$ indeed plays a crucial role in terms of forecasting performance is demonstrated in Section (ref). Interestingly, among the vanilla GL priors, HS has the lowest scores: For all orders, it has even lower scores than MP_LIT(4) { concering one-step-ahead forecasts}. The story changes radically, though, when looking at the semi-global modifications. For $p \in \{3,4,5\}$, HS$^*$ performs best. DL$_a^*$, R2D2$_a^*$, and NG$_a^*$ with $p \in \{2,3,4\}$ also perform considerably better than the best Minnesota prior, SHM(3).
{ Regarding four-steps-ahead predictions, we again observe that the semi-global modifications usually perform better than the vanilla implementations. Here, the overall best performance is achieved by NG$^*_a$(3). What stands out compared to one-step-ahead forecasts is that even the vanilla GL priors are doing better than the Minnesota priors. In other words, at longer forecasting horizons there is stronger evidence in favor of sparsity inducing priors.}
Turning the focus to lag-length, we find that for the data set at hand, competitive forecasts require at least two lags of the endogenous variables as predictors. { Most priors perform best with $p=2$ or $p=3$ lags. HS$^*$ and MP_LIT perform best with $p=4$ lags.} Interestingly, { concerning one-step-ahead forecasts, the performances of HS, DL, DL$_a$, DL$_a^*$, NG, NG$_a$, NG$_a^*$, R2D2, R2D2$_a$, and R2D2$_a^*$ deteriorate with $p>3$. In this respect, HS$^*$ and SHM -- the only priors without a clear drop in scores for $p=5$ -- are the most robust priors}.
The forecasting exercise also reveals both potentials and limitations of the hierarchical treatment of hyperparameters. The Minnesota prior clearly benefits from treating both shrinkage parameters as random variables. For GL priors, the story is more complex. On the one hand, DL, NG, and R2D2 in their vanilla implementations do not improve by placing priors on the concentration parameters $a^\gamma$, $a^\delta$, and $a^\pi$, respectively. On the other hand, the semi-global modifications of both NG and R2D2 unfold their full potential only if both $a^\delta_j$ and $a^\pi_j$, $j=1,\dots,2p$, are equipped with hyperpriors. This finding confirms our conjecture that GL priors have limitations when the noise level and the sparsity differs in different regions of the parameter space.
Figure (ref) depicts how the forecasting scores evolve over time for an illustrative selection of different priors. It reveals that the { time between the early 2000s recession and the financial crisis 2007/2008 marks a turning point} in the evolution of the scores. { Especially during the time of the zero lower bound (ZLB), which we define as the period after the financial crisis when the FFR is below 20 basis points mavroeidisIdentificationZeroLower2021, } the sparse models -- the semi-global modifications of the GL priors -- accumulate much higher scores than the denser SHM. { The sparsity inducing priors seem to capture the ZLB considerably well, by virtually zeroing almost all cross-lag coefficients in the equations of the interest rates.} Before the crisis, however, SHM is very competitive. The only sparse prior en par with SHM before the crisis is HS*. Overall, there is no single model that performs best over the whole evaluation period. { The findings suggest that there may be no general answer to the question whether sparse or dense modeling techniques perform better in forecasting economic data.}
{ In real applications, there might be interest in forecasts of only a small set of variables banbura_large_2010, feldkircherSophisticatedSmallSimple2024. Figure (ref) depicts cumulative log predictive likelihoods where instead of the full 21-dimensional predictive density the joint predictive density of three focal variables, namely GDP, CPI and FFR, is evaluated. Similar to before, we observe that SHM dominates until the mid 2000s, and the semi-global-local shrinkage priors since the mid 2000s.}
The previous sections highlighted the heterogeneity of model performance over time. The question arises whether combining weighted forecasts can improve forecasting performance.
Dynamic model averaging raftery_online_2010,koop_forecasting_2012,onorante_dynamic_2016 is a straightforwardly implementable averaging scheme, where the weights are sequentially updated as new observations become available, depending on the forecasting performance up to this point in time. A model with strong support from the data for a given period will receive a higher weight for the following period. By contrast, a model performing relatively worse will receive a lower weight.
In a nutshell, DMA works as follows. Denote $PL_{t|t-1,i}$ the one-step-ahead PL in $t$ for model $i$ within model space $M$. The predicted weight $\omega_{t|t-1,i}$ associated with model $i$ is computed as follows,
where $0 \leq \alpha \leq 1$ denotes a discount factor. The model updating equation then is
To shed more light on the discount factor: Standard (static) Bayesian model averaging is achieved by setting $\alpha=1$, as $\omega_{t|t-1,i}$ would be proportional to the (training-sample) marginal likelihood of model $i$ using data through the time $t-1$. We follow raftery_online_2010 in specifying $\alpha = 0.99$, indicating that $PL_{t-1|t-2}$ receives $99\%$ as much weight as $PL_{t|t-1}$.
The sparse models that we consider for DMA are DL$_a^*$, HS$^*$, NG$_a^*$, and R2D2$_a^*$. MP_LIT and SHM, by contrast, are representing the dense models. Concerning lag-length, we average over VARs of order $p \in \{1,2,3,4,5\}$. Last but not least, we initialize the weights to be equally distributed over all considered models for the first prediction.
Turning towards the results, we first consider the one-step-ahead predictive density of all 21 variables. Figure (ref) shows the evolution of model weights and cumulative log predictive likelihoods of the considered models relative to DMA over the whole evaluation period. Until the end of the 2000s, SHM with $p \in \{2,3,4,5\}$ and HS$^*$ with $p=4$ are dominating alternately. From { 2009 until the end of the evaluation period, HS$^*$ with $p=4$ receives the majority of the weights, the remaining semi-global-local priors have some weights, and SHM has quasi zero weights.} Analyzing the cumulative log predictive likelihoods, it turns out that DMA is practically never worse than any single model specification.
DMA can also improve forecasts if only a subset of variables is considered. Figure (ref) depicts DMA model weights and cumulative log predictive likelihoods relative to DMA when the predictive density of GDP, CPI and FFR is taken into account, only. By and large, the evolution of the predictive scores is similar to when all variables are considered. However, the analysis of model weights reveals that the distribution of weights over the models is now more dispersed. DL$_a^*$(2), which is practically not apparent when evaluating the full 21-dimensional predictive density, now receives the largest weights after 2010.
{ The result of DMA, of course, depends on the choice of $\alpha$. To address this issue, bernaciakLossDiscountingFramework2024 introduce the loss discounting framework (LDF) which extends DMA to general loss functions and more general discounting dynamics. We consider the following variant of LDF: Let $S_\alpha=\{\alpha_1,\dots,\alpha_s\}$ denote a set of different discount factors. In a first layer, we apply DMA with discount factor $\alpha_i$ for $i=1,\dots,s$ to our model pool. This yields $s$ dynamically averaged meta-models. In the second layer, we apply DMA with discount factor $\alpha_i$ for $i=1,\dots,s$ to the $s$ meta-models. This process can be repeated arbitrarily. In the final meta-model layer a single discount factor needs to be chosen. Adding more layers has a diminishing effect on the sequence of predictive distributions and therefore on the final result. It can be shown that the weights converge as the number of layers approaches infinity. In this respect, the final result is least affected by the choice of the discount factor in the last layer when a large number of layers $L$ is chosen. Concretely, we set $S_\alpha=\{0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,0.95,0.99,1\}$ and $l=4$, which appears to be sufficient in our case in the sense that the exact choice of the discount factor in the final layer shows negligible impact. According to Figure (ref), LDF outperforms all single model specifications and DMA throughout the whole evaluation period with respect to the full 21-dimensional predictive density. Focusing on the joint predictive distribution of the three focal variables, LDF is en par with DMA and not worse than any single model in the long run. }
To sum up, model performance is not only heterogeneous over time, but also across evaluated variables. Model averaging is a very effective tool in combining the merits of various modeling techniques, in contrast to the apparently unrewarding search for the single model that can do them all.
{ In this section, we investigate the sensitivity of results with respect to different prior choices for $\bm \Sigma_t$. The main focus is on different prior choices for $\bm u$. Additionally, we briefly discuss a homoskedastic specification of the variance-covariance matrix.}
The elements in $\bm u$ are often interpreted as contemporaneous coefficients. Hence, an obvious choice for $\bm u$ would be to use a prior of the same family as the prior on $\bm \phi$, e.g., DL/DL$_a^*$ on $\bm \phi$ is combined with DL on $\bm u$. For DL on $\bm u$, { for simplicity we specify $a^\gamma=\frac{1}{n_u}$, although this parameter could be learned from the data as well. Here, $n_u$ is the number of elements in $\bm u$.} For NG and R2D2 on $\bm u$ we set $a^\delta = a^\pi = \frac{1}{2 n_u}$, such that for DL, NG, and R2D2, the concentration at zero is matched. { For SSVS on $\bm u$ we specify $\tau_{0i} = \frac{1}{1000}$ and $\tau_{1i} = 1.$} As a second alternative, we investigate a relatively uninformative normal prior with zero mean and variance 10 for each element in $\bm{u}$, independently, which we refer to as FLAT prior. As a third option, we investigate a global(-but-not-local) shrinkage (GS) prior. More specifically, the prior is conditionally normal with zero mean and a global gamma hyperprior on the variance parameter: $u_{ij} | \lambda_g \sim N(0,\lambda_g)$, $\lambda_g \sim G(0.01,0.01)$. GS is similar in spirit to SHM. In combination with SHM on $\bm \phi$, then the contemporaneous coefficients $\bm u$ represents the third distinct group besides own-lag and cross-lag coefficients.
Table (ref) shows sums of one-step-ahead log predictive likelihoods relative to a VAR with $p=4$ lags with MP_LIT on $\bm \phi$ and HS on $\bm u$ for an illustrative selection of VARs. GS or FLAT on $\bm{u}$ clearly deteriorates forecasting performance of all models compared to HS used in the main analysis. Interestingly, GS is even worse than FLAT, which indicates that there are few strong signals that GS simply overshrinks. GL priors and SSVS on $\bm u$ are serious contenders to HS. Overall, the results highlight the importance of a sparsity-inducing prior for $\bm u$. The sensitivity of forecasting performance with respect to the prior choice for $\bm u$ strengthen our view to use the same prior on $\bm u$ for all models in the main analysis and also could explain to some extent the relatively weak performance of SHM in huber_adaptive_2019, where it is combined with the FLAT prior on $\bm{u}$.
{ Concerning the role of SV in forecasting macroeconomic data, we also consider homoskedastic VARs in a preliminary exercise. We assume that the diagonal matrix $\bm D_t=\bm D$ is constant for $t=1,\dots,T$. In that case, we specify independent inverse gamma priors with both shape and scale equalling 0.01 for the diagonal elements in $\bm D$. The prior on $\bm u$ is HS as described in Section (ref). In terms of forecasting performance, the results (available upon request) are overwhelmingly clear in favor of stochastic volatility. Therefore, we do not continue this route of research any further. More comprehensive treatments of this matter with similar findings can be found in koopLargeTimevaryingParameter2013,clark2015, carriero_large_2019, among others.}
Careful prior elicitation in Bayesian vector autoregressions goes back at least to litterman_forecasting_1986 and has been refined ever since. Our main contribution is a discussion of modern variants of prior choices, both established and novel, and a side-by-side comparison in terms of their structural (sparsity) behavior. In doing so, we borrow both from well-established (Minnesota-type) structural knowledge as well as from more automatic (shrinkage-type) approaches, thereby aiming at combining the best of both worlds. As a result of this aim, we propose what we term semi-global-local priors. We compare the various specifications in terms of their in- and out-of-sample performance on simulated data, as well as on a practically interesting (and notoriously hard) macroeconomic prediction problem. { Overall, we find that the semi-global-local priors outperform the vanilla global-local priors and are strong competitors to the semi-hierarchical Minnesota prior, which itself already constitutes an improvement to the traditional Minnesota prior.} Finally, we investigate how dynamic model averaging can help to shed light on the temporal patterns representing historical regime shifts. It becomes clear that no single model dominates the others; much rather, we find that some models are sometimes better than others, and dynamically combining them is fruitful.
While trying to comprehensively and integratively shed light on the Illusion of Sparsity debate, we also want to point out that many questions remain unanswered. Issues left for future research include an in-depth analysis where forecasting performance between modeling reduced-form and modeling structural-form coefficients is compared side by side. Moreover, a systematic analysis of the various different approaches to modeling reduced-form VARs with time-varying variance-covariance matrices (e.g.\ Cholesky-SV vs.\ Factor-SV) is still missing. { A first notable attempt in this direction is undertaken by chanComparingStochasticVolatility2023. As mentioned by an anonymous referee, the semi-global approach could be useful for putting structure on the variance-covariance matrix as well, e.g., by imposing group-specific shrinkage according to industry and economic sectors kastner_sparse_2019.} {
Furthermore, while our forecasting exercise takes into account the entire predictive distribution, we do not specifically focus on the assessment of tail risk. Threshold- and quantile-weighted scoring rules as discussed in gneiting2011 or the asymmetric continuous probability score (ACPS) proposed in iacopiniProperScoringRules2023 could serve as starting points in this regard. } Lastly, we believe that this work could be expanded to cater for VARs with time-varying regression parameters bitto2019achieving, huber2019should, knaussfs, feldkircherSophisticatedSmallSimple2024 where the issue of choosing good shrinkage priors is exacerbated in the sense that in addition to shrinkage towards zero, shrinkage towards constancy needs to be dealt with.