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.
79,109 characters · 13 sections · 61 citation commands
Sparse Bayesian Vector Autoregressions in Huge Dimensions
\onehalfspacing
Keywords: Efficient MCMC; Shrinkage; Factor stochastic volatility; Dirichlet-Laplace prior; Normal-gamma prior; Minnesota prior.
Acknowledgments: The authors acknowledge funding from the Austrian Science Fund (FWF) for the project “High-dimensional statistical learning: New methods to advance economic and sustainability policies” (ZK 35), jointly carried out by WU Vienna University of Economics and Business, Paris Lodron University Salzburg, TU Wien, and the Austrian Institute of Economic Research (WIFO).
Previous research has identified two important features that macroeconometric models should possess: the ability to exploit high dimensional information sets banbura2010, stock2011dynamic, roc-mca:dyn, koop2019bayesian and the possibility to capture non-linear features of the underlying time series Cogley2002, primiceri2005time,clark2011real, Clark2015, bitto2015achieving, huber2019should. While the literature suggests several paths to estimate large models, the majority of such approaches implies that once non-linearities are taken into account, analytical solutions are no longer available and the computational burden becomes prohibitive.\footnote{Two recent exceptions are Koop2013185 and carriero2015common.} This implies that high dimensional non-linear models can practically be estimated only under strong (and often unrealistic) restrictions on the dynamics of the model. However, especially in forecasting applications or in structural analysis, successful models should generally be able to exploit lots of information and also control for breaks in the autoregressive parameters or, more importantly, changes in the volatility of economic shocks primiceri2005time, Sims2006a, koop2009evolution.
Two reasons limit the use of large (or even huge) non-linear models. The first reason is statistical. Since the number of parameters in a standard vector autoregression rises quadratically with the number of time series included and commonly used macroeconomic time series are rather short, in-sample overfitting turns out to be a serious issue. As a solution, the Bayesian literature on VAR modeling Doan1984,Litterman1986,Sims1998, george2008bayesian, banbura2010,clark2011real,koop2013forecasting,Clark2015, korobilis2016, fol-yu:ach, huber2016adaptive, ankargren2019flexible suggests shrinkage priors that push the parameter space towards some stylized prior model like a multivariate random walk. On the other hand, ahe-etal:spa suggest to view VARs as graphical models and perform model selection drawing from the literature on sparse directed acyclic graphs. This typically leads to much improved forecasting properties and more meaningful structural inference. Moreover, the majority of the literature on Bayesian VARs imposes conjugate priors on the autoregressive parameters, allowing for analytical posterior solutions and thus avoiding simulation based techniques like Markov chain Monte Carlo (MCMC). Frequentist approaches often consider multi-step approaches dav-etal:spa.
The second reason is computational. Non-linear Bayesian models typically have to be estimated by means of MCMC, and computational intensity increases vastly when the number of component series becomes large. This increase stems from the fact that standard algorithms for multivariate regression models call for the inversion of large covariance matrices. Especially for sizable systems, this can quickly turn prohibitive since the inverse of the posterior variance-covariance matrix on the coefficients has to be computed for each sweep of the MCMC algorithm. For natural conjugate models, this step can be vastly simplified because the likelihood possesses a convenient Kronecker structure, implying that all equations in the VAR feature the same set of explanatory variables. This speeds up computation by large margins but restricts the flexibility of the model. carriero2015common, for instance, exploit this fact and introduce a simplified stochastic volatility specification. Another strand of the literature augment each equation of the VAR by including the residuals of the preceding equations carriero2019 which also provides significant improvements in terms of computational speed. Finally, in a recent contribution, koop2019bayesian reduce the dimensionality of the problem at hand by randomly compressing the lagged endogenous variables in the VAR.
All papers mentioned hitherto focus on capturing cross-variable correlation in the conditional mean through the VAR part and the co-movement in volatilities is captured by a rich specification of the error variance primiceri2005time or by a single factor carriero2015common. Another strand of the literature, typically used in financial econometrics, utilizes factor models to provide a parsimonious representation of a covariance matrix, focusing exclusively on the second moment of the predictive density. For instance, pitt1999time and aguilar2000bayesian assume that the variance-covariance matrix of a broad panel of time series might be described by a lower dimensional matrix of latent factors featuring stochastic volatility and a variable-specific idiosyncratic stochastic volatility process.
The present paper combines the virtues of exploiting large information sets and allowing for movements in the error variance. The overfitting issue mentioned above is solved as follows. First, we use a Dirichlet-Laplace (DL) prior specification bhattacharya2015dirichlet on the VAR coefficients. This prior is a global-local shrinkage prior in the spirit of polson2010shrink that enables us to heavily shrink the parameter space but at the same time provides enough flexibility to allow for non-zero regression coefficients if necessary. Second, a factor stochastic volatility model on the VAR errors grants a parsimonious representation of the time-varying error variance-covariance matrix of the VAR. To deal with the computational complexity, we exploit the fact that, conditionally on the latent factors and their loadings, equation-by-equation estimation becomes possible within each MCMC iteration. Moreover, we apply recent advances for fast sampling from high dimensional multivariate Gaussian distributions bhattacharya2015fast that permit estimation of models with hundreds of thousands of autoregressive parameters and an error covariance matrix with tens of thousands of nontrivial time-varying elements on a quarterly US dataset in a reasonable amount of time. In a careful analysis, we show to what extent our proposed method improves upon a set of standard algorithms typically used to simulate from the joint posterior distribution of large dimensional Bayesian VARs.
We first assess the merits of our approach in an extensive simulation study based on a range of different data generating processes. Relative to a set of competing benchmark specifications we show that, in terms of point estimates, the proposed global-local shrinkage prior yields precise parameter estimates and successfully introduces shrinkage in the modeling framework, without overshrinking significant signals.
In an empirical application, we adopt a modified version of the quarterly dataset proposed by stock2011dynamic and mccracken2016fred. To illustrate the out-of-sample performance of our model, we forecast important economic indicators such as output, consumer price inflation and short-term interest rates, amongst others. The proposed model is benchmarked against several alternatives. Our findings suggests that it performs well in terms of one-step-ahead predictive likelihoods. In addition, investigating the time profile of the cumulative log predictive likelihood reveals that allowing for large information sets in combination with the factor structure especially pays off in times of economic stress.
The remainder of this paper is structured as follows. Section (ref) introduces the econometric framework. Section (ref) details the Bayesian estimation approach, including an elaborated account of the (shrinkage) prior setup adopted and the corresponding conditional posterior distributions. Section (ref) provides an analysis of the computational gains of our algorithm relative to a set of established algorithms. Section (ref) presents the results of an extensive simulation study comparing the performance of carefully selected shrinkage priors for different time series lengths and model dimensions within various (sparse and dense) data generating scenarios. Section (ref), after giving a brief overview of the dataset used along with the model specification, illustrates our modeling approach by fitting a single factor model to 215-dimensional quarterly US data. Moreover, we perform a forecasting exercise to assess the predictive performance of our approach and discuss the choice of the number of latent factors. Finally, Section (ref) concludes.
Suppose interest centers on modeling an $m \times 1$ vector of time series denoted by $\boldsymbol{y}_t$ with $t=1,\dots,T$. We assume that $\boldsymbol{y}_t$ follows a heteroscedastic VAR($p$) process,\footnote{For simplicity of exposition we omit the intercept term in the following discussion (which we nonetheless include in the empirical application).}
Each $\boldsymbol{A}_j~(j=1,\dots,p)$ is an $m \times m$ matrix of autoregressive coefficients. The error term is assumed to follow a multivariate Gaussian distribution with time-varying variance-covariance matrix $\boldsymbol{\Omega}_t$. To permit reliable and parsimonious estimation when $m$ is large, we decompose the residual covariance matrix into
where both $\boldsymbol{\Sigma}_t = \text{diag}(\sigma^2_{1t},\dots,\sigma^2_{mt})$ and $\boldsymbol{V}_t = \text{diag}(e^{h_{1t}},\dots,e^{h_{qt}})$ are diagonal matrices with dimension $m$ and $q$, respectively, and $\boldsymbol{\Lambda}$ denotes an $m \times q$ matrix of factor loadings with typical element $\lambda_{ij}$ ($i=1,\dots,m$; $j=1,\dots,q$). The logarithms of the diagonal elements of $\boldsymbol{\Sigma}_t$ and $\boldsymbol{V}_t$ follow AR($1$) processes,
To identify the scaling of the elements of $\boldsymbol{\Lambda}$, the process specified in ((ref)) is assumed to have mean zero while $\mu_{\sigma j}$ in ((ref)) is the unconditional mean of the log-elements of $\boldsymbol{\Sigma}_t$ to be estimated from the data kas-etal:eff. The parameters $\rho_{hj}$ and $\rho_{\sigma i}$ are a priori restricted to the interval $(-1,1)$ and denote the persistences of the latent log variances. The error terms $e_{hj, t}$ and $e_{\sigma i,t}$ constitute independent zero mean innovations with variances $\varsigma^2_{hj}$ and $\varsigma^2_{\sigma i}$, respectively. This specification implies that the volatilities are mean reverting and thus bounded in the limit.
This error structure is known as the factor stochastic volatility model pitt1999time, aguilar2000bayesian. It can be equivalently written by introducing $q$ conditionally independent latent factors $\boldsymbol{f}_t \sim \mathcal{N}_q(\boldsymbol{0}, \boldsymbol{V}_t)$ and rewriting the error term in ((ref)) as
Note that off-diagonal entries of $\boldsymbol{\Omega}_t$ exclusively stem from the volatilities of the $q$ factors while the diagonal entries of $\boldsymbol{\Omega}_t$ are allowed to feature idiosyncratic deviations driven by the elements of $\boldsymbol{\Sigma}_t$. This specification reduces the number of free elements in $\boldsymbol{\Omega}_t$ from $m(m+1)/2$ to $mq$, where the latter quantity is typically much smaller than the former. In addition, by conditioning on the latent factors, this representation enables us to derive an efficient Gibbs sampler that allows for conditional equation-by-equation estimation. As will be discussed in more detail in Section (ref), this constitutes a key feature for computationally feasible Bayesian inference when the dimensionality $m$ becomes large.
The model described by Eqs.\ ((ref)) to ((ref)) is related to several alternative specifications commonly used in the literature. For instance, assuming that $\boldsymbol{V}_t = \boldsymbol{I}$ and $\boldsymbol{\Sigma}_t \equiv \boldsymbol{\Sigma}$ for all $t$ leads to the specification adopted in stock2005understanding. Setting $q=1$ and $\boldsymbol{\Sigma}_t \equiv \boldsymbol{\Sigma}$ yields a specification that is similar to the one stipulated in carriero2015common, with the difference that our model imposes restrictions on the covariances whereas carriero2015common estimate a full (but constant) covariance matrix. In addition, our model implies that the stochastic volatility enters $\boldsymbol{\Omega}_t$ in an additive fashion.
Before proceeding to the next subsection it is worth summarizing the key features of the model given by Eqs.\ ((ref)) to ((ref)). First, we capture cross-variable movements in the conditional mean through the VAR block of the model and assume that co-movement in conditional variances is captured by a factor structure. Second, the model introduces stochastic volatility by assuming that a large panel of volatilities may be efficiently summarized through a set of latent heteroscedastic factors. This choice is more flexible than a single factor model for the volatility, effectively providing a parsimonious representation of $\boldsymbol{\Omega}_t$ that is flexible enough to replicate the dynamic behavior of the variances of a broad set of macroeconomic quantities.
Our approach to estimation and inference is Bayesian. This implies that after specifying a suitable prior distribution on the model parameters, we can combine this prior with the likelihood implied by the data and the model to obtain the corresponding posterior distribution.
For prior implementation, it proves to be convenient to define a $k \times 1$ vector of predictors $\boldsymbol{x}_t =( \boldsymbol{y}'_{t-1},\dots,\boldsymbol{y}_{t-p}')'$ and an $m \times k$ coefficient matrix $\boldsymbol{B}=(\boldsymbol{A}_1,\dots,\boldsymbol{A}_p)$ with $k=mp$ to rewrite the model in ((ref)) more compactly as $ \boldsymbol{y}_t = \boldsymbol{B} \boldsymbol{x}_t + \boldsymbol{\varepsilon}_t. $ Stacking the rows of $\boldsymbol{y}_t$, $\boldsymbol{x}_t$, and $\boldsymbol{\varepsilon_t}$ yields
where $\boldsymbol{Y} = (\boldsymbol{y}_1,\dots,\boldsymbol{y}_T)'$, $\boldsymbol{X} = (\boldsymbol{x}_1, \dots, \boldsymbol{x}_T)'$, and $\boldsymbol{E} = (\boldsymbol{\varepsilon}_1,\dots,\boldsymbol{\varepsilon}_T)'$ denote the corresponding full data matrices.
Typically, the matrix $\boldsymbol{B}$ is a sparse matrix with non-zero elements mainly located on the main diagonal of $\boldsymbol{A}_1$. In fact, existing priors in the Minnesota tradition tend to strongly push the system towards the prior model in high dimensions. However, especially in large models an extremely tight prior on $\boldsymbol{B}$ might lead to severe overshrinking, effectively zeroing out coefficients that might be important to explain $\boldsymbol{y}_t$. If the matrix $\boldsymbol{B}$ is characterized by a relatively low number of non-zero regression coefficients, a possible solution is a global-local shrinkage prior polson2010shrink.
A recent variant that falls within the class of global-local shrinkage priors is the Dirichlet-Laplace (DL) prior put forward in bhattacharya2015dirichlet. This prior possesses convenient shrinkage properties in the presence of a large degree of sparsity of the parameter vector $\boldsymbol{b}=\text{vec}(\boldsymbol{B})$. In what follows, we impose the DL prior on each of the $K=mk$ elements of $\boldsymbol{b}$, denoted as $b_j$ for $j=1,\dots,K$,
where $\mathcal{DE}$ denotes the double exponential (Laplace) and $\mathcal{E}$ the exponential distribution, $\psi_j$ is an auxiliary scaling parameters to achieve conditional normality, and the elements of $\boldsymbol{\vartheta}=(\vartheta_1,\dots,\vartheta_K)'$ are local auxiliary scaling parameters that are bounded to the $(K-1)$-dimensional simplex $\mathcal{S}^{K-1}=\{\boldsymbol{\vartheta}: \vartheta_j \ge 0, \sum_{j=1}^n \vartheta_j =1\}$. A natural prior choice for $\vartheta_j$ is the (symmetric) Dirichlet distribution with hyperparameter $a$, $ \vartheta_j\sim \mathcal{D}(a,\dots,a). $ In addition, $\zeta$ is a global shrinkage parameter that pushes all elements in $\boldsymbol{B}$ towards zero and exhibits an important role in determining the tail behavior of the marginal prior distribution on $b_j$, obtained after integrating out the $\vartheta_j$s. Thus, we follow bhattacharya2015dirichlet and adopt a fully Bayesian approach by specifying a gamma distributed prior on $ \zeta \sim \mathcal{G}(K a, 1/2). $ It is noteworthy that this prior setup has at least two convenient features that appear to be of prime importance for VAR modeling. First, it exerts a strong degree of shrinkage on all elements of $\boldsymbol{B}$ but still provides additional flexibility such that non-zero regression coefficients are permitted. This critical property is a feature which a large class of global-local shrinkage priors share griffin2010inference,carvalho2010horseshoe,polson2010shrink and has been recently adopted in a VAR framework by huber2016adaptive and within the general context of state space models by bitto2015achieving. Second, implementation is simple and requires relatively little additional input from the researcher. In fact, the prior heavily relies on a single structural hyperparameter that has to be specified with care, namely $a$.
The hyperparameter $a$ influences the empirical properties of the proposed shrinkage prior along several important dimensions. Smaller values of $a$ lead to heavy shrinkage on all elements of $\boldsymbol{B}$. To see this, note that lower values of $a$ imply that more prior mass is placed on small values of $\zeta$ a priori. Similarly, when $a$ is small, the Dirichlet prior places more mass on values of $\vartheta_j$ close to zero. Since lower values of $\zeta$ translate into thicker tails of the marginal prior on $b_j$, the specific choice of $a$ not only influences the overall degree of shrinkage but also the tail behavior of the prior. Letting $\tilde p$ denote the number of predictors, bhattacharya2015dirichlet show that if $a$ is specified as $\tilde p^{-(1+\Delta)}$ for any $\Delta>0$ to be small, the DL prior displays excellent posterior contraction rates, and pat-etal:pos discuss the shrinkage properties of the proposed prior within the context of factor models. In our application, $\tilde p=K$ (when considering the total number of predictors) or $\tilde p=k$ (when considering the number or predictors per equation).
For the factor loadings we independently use a standard normally distributed prior on each element $\lambda_{ij}\sim \mathcal{N}(0,1)$ for $i=1,\dots,m$ and $j=1,\dots,q$. In the empirical application (Section (ref)), we in addition consider the row-wise normal-gamma griffin2010inference shrinkage prior discussed in kas:spa, i.e.\ $\lambda_{ij}|\tau_{ij}^2 \sim \mathcal{N}(0, \tau_{ij}^2)$, $\tau_{ij}^2|\upsilon_{i}^2 \sim \mathcal{G}(a_\lambda, a_\lambda \upsilon_{ij}^2/2)$, $\lambda_{ij}^2 \sim \mathcal{G}(c_\lambda, d_\lambda).$ Furthermore, we impose a normally distributed prior on the mean of the log-volatility $\mu_{\sigma j} \sim \mathcal{N}(0,M_\mu)$ with $M_\mu$ denoting the prior variance, and the commonly employed Beta distributed prior on the transformed persistence parameter of the log-volatility $\frac{\rho_{sj}+1}{2} \sim \mathcal{B}(a_0, b_0)$ for $s \in \{h, \sigma\}$ and $a_0, b_0 \in \mathbb{R}^+$ to ensure stationarity. Finally, we use a restricted gamma prior on the innovation variances in Eqs.\ ((ref)) and ((ref)), $\varsigma_{sj}^2 \sim \mathcal{G}(\frac{1}{2}, \frac{1}{2 \xi})$. Here, $\xi$ is a hyperparameter used to control the tightness of the prior. This choice, motivated in fruhwirth2010stochastic implies that if the data is not informative on the degree of time variation of the log volatilities then we do not bound $\varsigma_{sj}^2$ artificially away from zero, effectively applying more shrinkage than the standard inverted gamma prior.
Conditional on the latent factors and the corresponding loadings, the model in ((ref)) can be cast as a system of $m$ unrelated regression models for the elements in $\boldsymbol{z}_t =\boldsymbol{y}_t-\boldsymbol{\Lambda} \boldsymbol{f}_t$, labeled $z_{it}$, with heteroscedastic errors,
Here we let $\boldsymbol{B}_{i \bullet}$ denote the $i$th row of $\boldsymbol{B}$ and $\eta_{it}$ is the $i$th element of $\boldsymbol{\eta}_t$. The corresponding posterior distribution of $\boldsymbol{B}'_{i\bullet}$ is $k$-variate Gaussian,
with $\bullet$ indicating that we condition on the remaining parameters and latent quantities of the model. The posterior variance and mean are given by
The diagonal prior covariance matrix of the coefficients related to the $i$th equation is given by $\boldsymbol{\Phi}_i$, the respective $k \times k$ diagonal submatrix of $\boldsymbol{\Phi}=\zeta \times \text{diag}(\psi_1 \vartheta_1^2 ,\dots,\psi_K \vartheta_K^2)$. Moreover, $\tilde{\boldsymbol{X}}_i$ is a $T\times k$ matrix with typical row $t$ given by $\boldsymbol{X}_t/\sigma_{it}$ and $\tilde{\boldsymbol{z}}_{i}$ is a $T$-dimensional vector with the $t$th element given by $z_{it}/\sigma_{it}$. This normalization renders ((ref)) conditionally homoscedastic with standard normally distributed white noise errors.
The full conditional posterior distribution of $\psi_j$ is inverse Gaussian,
The conditional posterior of the global shrinkage parameter $\zeta$ follows a generalized inverse Gaussian (GIG) distribution,
To draw from this distribution, we use the efficient algorithm of hoe-ley:gen. Moreover, we sample the scaling parameters $\vartheta_j$ by first sampling $L_j$ from $ L_j |\bullet \sim \mathcal{GIG}(a-1, 1, 2|b_j|), $ and then setting $ \vartheta_j = L_j/\sum_{i=1}^K L_i. $
The conditional posterior distributions of the factors are Gaussian and thus straightforward to draw from. The factor loadings are sampled using “deep interweaving” kas-etal:eff, and the parameters in ((ref)) and ((ref)) along the full histories of the latent log-volatilities are sampled as in kastner2014ancillarity using the R-packages factorstochvol r:fac and stochvol kas:dea.
Our MCMC algorithm iteratively draws from the conditional posterior distributions outlined above and discards the first $J$ draws as burn-in. In terms of computational requirements, the single most intensive step is the simulation from the joint posterior of the autoregressive coefficients in $\boldsymbol{B}$. Because this step is implemented on an equation-by-equation basis, speed improvements relative to the standard approach are already quite substantial. However, note that if $k$ is large (i.e.\ of the order of several thousands), even the commonly employed equation-by-equation sampling fails to deliver a sufficient amount of draws within a reasonable time window. Consequently, we outline an alternative algorithm to draw from a high-dimensional multivariate Gaussian distribution under a Bayesian prior that features a diagonal prior variance-covariance matrix in the upcoming section.
The typical approach to sampling from ((ref)) is based on the full system and simultaneously samples from the full conditional posterior of $\boldsymbol{B}$, implying that the corresponding posterior distribution is a $K$-dimensional Gaussian distribution with a $K \times K$ dimensional variance-covariance matrix. Under a non-conjugate prior, the computational difficulties arise from the need to invert the $K\times K$ variance-covariance matrix which requires operations of order $O(m^6 p^3)$ under Gaussian elimination.
If a conjugate prior in combination with a constant carriero2015common specification of $\boldsymbol{\Omega}_t$ is used, the corresponding variance-covariance features a Kronecker structure which is computationally cheaper to invert and scales better in large dimensions. Specifically, the manipulations of the corresponding covariance matrix are of order $O(m^3+k^3)$, a significant gain relatively to the standard approach. However, this comes at a cost since all equations have to feature the same set of variables, the prior on the VAR coefficients has to be symmetric and any stochastic volatility specification that preserves conjugacy is necessary overly simplistic.
By contrast, recent studies emphasize the computational gains that arise from utilizing a framework that is based on equation-by-equation estimation. carriero2019 and koop2019bayesian augment each equation of the system by either contemporaneous values of the endogenous variables of the preceding equations or the residuals from the previous equations. Here, our approach renders the equations of the system conditionally independent by conditioning on the factors. From a computational perspective, the differences between using a factor model to disentangle the equations and an approach based on augmenting specific equations by quantities that aim to approximate covariance parameters are negligible. If we sample from ((ref)) directly, the computations involved are of order $O(m k^3)=O(m^4 p^3)$. This already poses significant improvements relative to full system estimation.
One contribution of the present paper is the application of the algorithm proposed by bhattacharya2015fast and developed for univariate regression models under a global-local shrinkage prior. This algorithm is applied to each equation in the system and cycles through the following steps:
This algorithm outperforms all competing variants discussed previously in situations where $k \gg T$, a situation commonly encountered when dealing with large VAR models. In such cases, steps (1) to (4) can be carried out using $O(p m^2 T^2)$ floating point operations. In situations where $k \approx T$, the computational advantages relative to the standard equation-by-equation algorithm mentioned above are modest or even negative. However, note that the cost is quadratic in $m$ and linear in $p$ and thus scales much better when the number of endogenous variables and/or lags thereof is increased. More information on the empirical performance of our algorithm can be found in Section (ref).
This section aims at comparing the performance of the DL prior with a range of commonly used alternatives. We investigate sparse, intermediate, and dense data generating processes (DGPs) where $T \in \{50, 100, 150, 200, 250\}$ and $m \in \{10, 20, 50, 100\}$. The probability of an off-diagonal entry to be non-zero is $0.01$, $0.1$, and $0.8$ in each of the respective scenarios. In all scenarios, each intercept entry has a $0.1$ probability of being non-zero and all diagonal elements are non-zero with probability $0.8$. The non-zero elements are randomly generated from Gaussian distributions roughly tuned to yield stable VARs. More concretely, both the mean $\mu_I$ and the standard deviation $\sigma_I$ of the intercept are set to $0.01$, whereas mean and standard deviation of the diagonal ($D$) and the off-diagonal ($O$) elements are chosen as follows:
Concerning the errors, we use a single factor SV specification. The factor loadings are generated from $\mathcal{N}(0.001, 0.001^2)$ to roughly match the above scaling. The AR(1) processes driving the idiosyncratic log-variances are assumed to have mean $\mu_{\sigma i} = -12$ with persistences $\rho_{\sigma i}$ ranging from $0.85$ to $0.98$ and innovation standard deviations $\varsigma_{\sigma i}$ from $0.3$ to $0.1$. The process driving the factor log variance is assumed to be highly persistent with $\rho_{h1} = 0.99$ and $\varsigma_{h1} = 0.1$.
For each of the $60$ settings, we simulate $10$ data sets. For each of these, we run our MCMC algorithm to obtain $2000$ posterior draws after a burn-in of $1000$. Consequently, the posterior means are compared to the true values and root mean squared errors are computed. Finally, the median of each of these is reported in Table (ref). Alongside the DL prior with weak ($a_{DL} = a = 1/2$) and strong ($a_{DL} = a = 1/k$ and $a_{DL} = a = 1/K$) shrinkage, we also consider the NG prior with a single global shrinkage parameter huber2016adaptive and a standard conjugate Minnesota prior with a single shrinkage parameter $a_M$, implemented by using dummy observations. For the NG prior we specify the prior on the global shrinkage parameter to induce heavy shrinkage (by setting both hyperparameters of the gamma prior equal to $0.01$) and the prior controlling the excess kurtosis $a_{NG}$ is set equal to $1$, corresponding to the Bayesian Lasso par-cas:bay, and $a_{NG}=0.1$. The latter choice places significant prior mass around zero but at the same time leads to a heavy tailed marginal prior. Finally, we report RMSEs of the OLS estimator (if it exists).
As is to be expected, Table (ref) reveals strong to severe overfitting of OLS (corresponding to the posterior mode under a flat prior) which can be mitigated to a certain extent when the Minnesota prior with $a_M = 0.001$ is employed instead. Similarly, the DL prior with weak shrinkage ($a_{DL} = 1/2$) displays a tendency to overfit, in particular when $T$ is small. By contrast, the more aggressive DL and NG shrinkage priors show superior performance. Overall, DL($1/k$) and NG($0.1$) exhibit lowest RMSEs, where DL($1/k$) performs best in the sparse scenarios, NG($0.1$) performs best in the intermediate settings, and no clear winner is to be found in the dense context. Turning towards NG($1$) and DL($1/K$) we tend to observe acceptable but slightly inferior overall performance. The Minnesota prior with $a_M = 0.0001$ yields an extreme degree of shrinkage, translating into estimates of autoregressive coefficients that are very close to zero, irrespectively of the contribution from the likelihood. In that sense, it overshrinks most of the nonzero coefficients. Nevertheless, in scenarios with extremely low signal-to-noise ratios (such as the dense scenario with $T=50$ and $m=100$), this can be beneficial for the overall performance.
For further illustration, we showcase four exemplary scenarios in Figures (ref) to (ref) in the Appendix.
In Section (ref) we first summarize the data set adopted and present the model specification choices made. The section that follows (Section (ref)) estimates a simple one factor model to outline the virtues of our proposed framework. Section (ref) presents the main findings of our forecasting exercise and discusses the choice of the number of factors used for modeling the error covariance structure.
Aim of the empirical application is to forecast a set of key US macroeconomic quantities. To this end, we use the quarterly dataset provided by mccracken2016fred, a variant of the well-known stock2011dynamic dataset for the US.\footnote{In addition to quarterly observations, mccracken2016fred also provide a subset of the data which is observed monthly. Of course, our method is analogously applicable to higher frequency observations. However, given that the computational cost of the bhattacharya2015fast approach is quadratic in $T$, the run-time gains of their approach in comparison to equation-by-equation estimation is then smaller and can, depending on the number of lags, even become be negative.} The data spans the period ranging from 1959:Q1 to 2015:Q4. We include $m=215$ quarterly time series, capturing information on $14$ important segments of the economy and follow mccracken2016fred in transforming the data to be approximately stationary. Furthermore, we standardize each component series to have zero mean and variance one. In the empirical examples we include $p=1$ lags of the endogenous variables.\footnote{We have also experimented with higher lag orders and also found some evidence of signals at lag two for the dataset at hand; see Figures (ref) to (ref) in the Appendix for an illustration. However, out-of-sample predictive studies favored one lag only (cf.\ Section (ref)).} The hyperparameters are chosen as follows: $M_\mu = 10$, $a_0=20$, $b_0=1.5$, $\xi = 1$, $a_\lambda = 0.1$, $c_\lambda = d_\lambda = 1$.
To provide some intuition on how our modeling approach works in practice, we first estimate a simple one factor model (i.e.\ $q=1$) and investigate several features of our empirical model. In the next section we will perform an extensive forecasting exercise and discuss the optimal number of factors in terms of forecasting accuracy.
We start by inspecting the posterior distribution of $\boldsymbol{\Lambda}$ and assess what variables load heavily on the latent factor. It is worth emphasizing that most quantities\footnote{Hereby we refer to the one-step-ahead forecast error related to a given time series.} associated with real activity (i.e.\ industrial production and its components, GDP growth, employment measures) load heavily on the factor. Moreover, expectation measures, housing markets, equity prices and spreads also load heavily on the joint factor.
To assess whether spikes in the volatility associated with the factor coincide with major economic events, the bottom panel of Figure (ref) depicts the evolution of the posterior distribution of factor volatility over time. A few findings are worth mentioning. First, volatility spikes sharply during the midst of the 1970s, a period characterized by the first oil price shock and the bankruptcy of Franklin National Bank in 1974. After declining markedly during the second half of the 1970s, the shift in US monetary policy towards aggressively fighting inflation and the second oil price shock again translate into higher macroeconomic uncertainty. Note that from the mid 1980s onward, we observe a general decline in macroeconomic volatility that lasts until the beginning of the 1990s. There we observe a slight increase in volatility possibly caused by the events surrounding the first gulf war. The remaining years up to the beginning of the 2000s has been relatively unspectacular, with volatility levels being muted most of the time. In 2000/2001, volatility again increases due to the burst of the dot-com bubble and the 9/11 terrorist attacks. Finally, we observe marked spikes in volatility during recessionary episodes like the recent financial crisis in 2008.
Finally, we assess how well the DL prior with $a=1/k$ performs in shrinking the coefficients in $\boldsymbol{B}$ to zero. The top panel of Figure (ref) depicts a heatmap that gives a rough feeling on the size of each regression coefficient based on the posterior median of $\boldsymbol{B}$. The bottom panel of Figure (ref) depicts the posterior interquartile range, providing some evidence on posterior uncertainty.\footnote{Since the corresponding posterior distribution is quite heavy-tailed, using posterior standard deviations, while providing a qualitatively similar picture, tend to be slightly exaggerated.} The DL prior apparently succeeds in shrinking the vast majority of the approximately $50\,000$ coefficients towards zero. Even though not discussed in detail to conserve space, we note that at higher lag orders this very strong shrinkage effect is even more pronounced; see also Figures (ref) to (ref) in the Appendix.
The top panel of Figure (ref) displays the posterior median estimates when the shrinkage parameter $a$ is chosen to be $1/2$ bhattacharya2015dirichlet. While $a = 1/2$ appears to provide a fair amount of shrinkage in other applications, for our huge dimensional example this prior exerts only relatively little shrinkage and tends to lead to overfitting. The diagonal pattern in the first lag appears here as well, but there is a considerable amount of nonzero medians elsewhere. Correspondingly, the interquartile ranges visualized in the bottom panel of Figure (ref) are also very large compared to those obtained with $a = 1/k$.
Interestingly, for selected time series measuring inflation (both consumer and producer price inflation) we find that lags of monetary aggregates are allowed to load on the respective inflation series. This result points towards a big advantage of our proposed prior relative to standard VAR priors in the Minnesota tradition: while these priors have been shown to work relatively well in huge dimensions banbura2010, they also display a tendency to overshrink when the overall tightness of the prior is integrated out in a Bayesian framework, effectively pushing the posterior distribution of $\boldsymbol{B}$ towards the prior mean and thus ruling out patterns observed under the DL prior.
Inspection of the interquartile range also indicates that the proposed shrinkage prior succeeds in reducing posterior uncertainty markedly. Note that the pattern found for the posterior median of $\boldsymbol{B}$ can also be found in terms of the posterior dispersion. We again observe that the coefficients associated with the first, own lag of a given variable are allowed to be non-zero whereas in most other cases the associated posterior is strongly concentrated around zero.
We focus on forecasting gross domestic product (GDPC96), industrial production (INPRO), total nonfarm payroll (PAYEMS), civilian unemployment rate (UNRATE), new privately owned housing units started (HOUST), consumer price index inflation (CPIAUCSL), producer price index for finished goods inflation (PPIFGS), effective federal funds rate (FEDFUNDS), 10-year treasury constant maturity rate (GS10), U.S./U.K.\ exchange rate (EXUSUKx), and the S&P 500 (S.P.500). This choice includes the variables investigates by koop2019bayesian and some additional important macroeconomic indicators that are commonly monitored by practitioners, resulting in a total of eleven series.
To assess the forecasting performance of our model, we conduct a pseudo out-of-sample forecasting exercise with initial estimation sample ranging from 1959:Q3 to 1990:Q2. Based on this estimation period, we compute one-quarter-ahead predictive densities for the first period in the hold-out (i.e.\ 1990:Q3). After obtaining the corresponding predictive densities and evaluating the corresponding log predictive likelihoods, we expand the estimation period and re-estimate the model. This procedure is repeated $100$ times until the final point of the full sample is reached. The quarterly scores obtained this way are then accumulated.
Our model with $q \in \{0,1\dots,4\}$ factors is benchmarked against the prior model, a pure factor stochastic volatility (FSV) model with conditional mean equal to zero (i.e. $\boldsymbol{B}=\boldsymbol{0}_{m \times k})$. In what follows we label this specification FSV $0$. To assess the merits of the proposed shrinkage prior vis-\'{a}-vis a Minnesota prior and a NG shrinkage prior we also include the models described in Section (ref). Moreover, we include two models that impose the restriction that $\boldsymbol{A}_1=\boldsymbol{I}_m$ and $\boldsymbol{A}_1=0.8 \times \boldsymbol{I}_m$ while $\boldsymbol{A}_j$ for $j>1$ are set equal to zero matrices in both cases. The first model, labeled FSV $1$, assumes that the conditional mean of $\boldsymbol{y}_t$ follows a random walk process and the second specification, denoted as FSV $0.8$, imposes the restriction that the variables in $\boldsymbol{y}_t$ feature a rather strong degree of persistence but are stationary. The exercise serves to evaluate whether it pays off to impose a VAR structure on the first moment of the joint density of our data and to assess how many factors are needed to obtain precise multivariate density predictions for our eleven variables of interest.
Overall log predictive scores (LPSs) are summarized in Table (ref). An immediate finding is that ignoring the error covariance structure (using zero factors) produces rather inaccurate forecasts for all models considered. While a single factor model improves predictive accuracy by a large margin, allowing for more factors (i.e.\ even more flexible modeling of the covariance structure) further increases the forecasting performance. For this specific exercise, we identify two or three factors to be a reasonable choice for most models when the joint log predictive scores of the aforementioned variables are considered. We would like to stress that this choice critically depends on the number of variables we include in our prediction set. If we focus attention on the marginal predictive densities (i.e.\ the univariate predictive densities obtained after integrating out the remaining elements in $\boldsymbol{y}_t$) we find that fewer or even no factors receive more support (see Table (ref)), whereas in the case of higher dimensional prediction sets more than two factors lead to more accurate density predictions kas-etal:eff. As a general remark, we note that identifying the optimal number of factors in high-dimensional FSV models is a challenging problem in practice. Using the deviance information criterion chan2016fast may be an option but is likely to be unstable in very high dimensions. The approach adopted in the paper at hand, namely the decomposition of the marginal likelihood into predictive likelihoods geweke2010comparing tends to be more stable, in particular when interest ist placed on predicting subsets only. Moreover, it can be trivially parallelized, thus becoming computationally feasible on high performance computing infrastructures.
Considering forecasting accuracy across models reveals that our proposed VAR($1$)-FSV with a DL($1/k$) prior displays excellent forecasting capabilities, outperforming all competitors. Amongst the VAR($1$) models, DL($1/K$) and NG($0.1$) also do well, and the Bayesian Lasso (NG($1$)) as well as the Minnesota prior with medium shrinkage (Min($0.01$)) show decent performance. Clearly, DL($1/2$) overfits and Min($0.001$) overshrinks. Note that higher lag orders seem to rarely increase predictive accuracy. However, comparing the differences between the benchmark pure FSV models and the VAR-FSV models considered, we find that explicitly modeling the conditional mean improves the forecasting accuracy in practically all cases.
To investigate whether forecasting performance is homogeneous over time, Figure (ref) visualizes the cumulative LPSs relative to the zero-factor FSV model over time. The benefit of the flexible SV structure in the VAR residuals is particularly pronounced during the 2008 financial crisis which can be seen by comparing the solid lines to the broken lines. During this period, time-varying covariance modeling appears to be of great importance and the performance of models that ignore contemporaneous dependence deteriorates. This finding is in line with kas:spa who reports analogous results for US asset returns. The increase in predictive accuracy can be traced back to the fact that within an economic downturn, the correlation structure of our dataset changes markedly, with most indicators that measure real activity sharply declining in lockstep. A model that takes contemporaneous cross-variable linkages seriously is thus able to fully exploit such behavior which in turn improves predictions.
Up to this point, we focused exclusively on the joint performance of our model for the specific set of variables considered. To gain a deeper understanding on how our model performs for relevant selected quantities, Table (ref) displays marginal LPSs for the two most promising prior specifications with one, two, and five lags. The variables we consider are inflation (CPIAUCSL), short-term interest rates (FEDFUNDS), and output growth (GDPC96).
In contrast to the findings based on joint LPSs, we observe that models without a factor structure tend to perform better than models that set $q>0$, with the exception of interest rates where all models predict more or less equally badly. This finding corroborates our conjecture stated above, implying that if the set of focus variables is subsequently enlarged, more factors are necessary in order to obtain precise density predictions. Here, we only focus on marginal model performance, implying that for each variable, contemporaneous relations between the elements in $\bm{y}_t$ are integrated out. This, in turn, implies that the additional gain in model flexibility is offset by the comparatively larger number of parameters. Concerning the difference between VAR priors, it appears that NG slightly outperforms DL for inflation whereas DL is superior when it comes to predicting output growth.
Even though the efficient sampling schemes outlined in this paper help to overcome absolutely prohibitive computational burdens, the CPU time needed to perform fully Bayesian inference in a model of this size can still be considered substantial. In what follows we shed light on the estimation time required and how it is related to the length of the time series $T$, the lag length $p$ and to the number of latent factors $q \in \{0, 50\}$. Figure (ref) shows the time needed to perform a single draw from the joint posterior distribution of the $215 + 215^2p$ coefficients and their corresponding $2(215 + 215^2p) + 1$ auxiliary shrinkage quantities, the $qT$ factor realizations and the associated $215q$ loadings, alongside $(T+1)(215 + q)$ latent volatilities with their corresponding $645 + 2q$ parameters. This amounts to $166\,841$ random draws for the smallest model considered (one lag, no factors, $T=124$) and $776\,341$ random draws for the largest model ($5$ lags, $50$ factors, $T=224$) at each MCMC iteration.
As mentioned above, the computation time rises approximately linearly with the number of lags included. Dotted lines indicate the time in seconds needed to perform a single draw from a model with $50$ factors included while solid lines refer to the time needed to estimate a model without factors and a diagonal time-varying variance-covariance matrix $\boldsymbol{\Sigma}_t$. Interestingly, the additional complexity when moving from a model without factors to a highly parameterized model with $50$ factors appears to be negligible, increasing the time needed by a fraction of a second on average. The important role of the length of the sample can be seen by comparing the green, red and black lines. The time necessary to perform a simple MCMC draw quickly rises with the length of our sample, consistent with the statements made in Section (ref). This feature of our algorithm, however, is convenient especially when researchers are interested in combining many short time series or performing recursive forecasting based on a tiny initial estimation sample.
In this paper we propose an alternative route to estimate huge dimensional VAR models that allow for time-variation in the error variances. The Dirichlet-Laplace prior, a recent variant of a global-local shrinkage prior, enables us to heavily shrink the parameter space towards the prior model while providing enough flexibility that individual regression coefficients are allowed to be unrestricted. This prior setup alleviates overfitting issues generally associated with large VAR models. To cope with computational issues we assume that the one-step-ahead forecast errors of the VAR feature a factor stochastic volatility structure that enables us to perform equation-by-equation estimation, conditional on the loadings and the factors. Since posterior simulation of each equation's autoregressive parameters involves manipulating large matrices, we implement an alternative recent algorithm that improves upon existing methods by large margins, rendering a fully fledged Bayesian estimation of truly huge systems possible.
In an empirical application we first present various key features of our approach based on a single factor model. This single factor which summarizes the joint dynamics of the VAR errors can be interpreted as an uncertainty measure that closely tracks observed factors such as the volatility index. The question whether such a simplistic structure proves to be an adequate representation of the time-varying covariance matrix naturally arises and we thus provide a detailed forecasting exercise to evaluate the merits of our approach relative to the prior model and a set of competing models with a different number of latent factors in the errors.
Finally, three potential extensions are worth mentioning. First, given the fact that systematic and in-depth empirical comparisons of the various recently developed roads towards handling high-dimensional VARs with time-varying contemporaneous covariance in a Bayesian framework (VAR-FSV, VAR-Cholesky-SV, compressed VAR-SV, etc.) is still missing and it is not clear whether one of these models turns out to dominate the others for all points in time, one could consider to average/select dynamically. Second, note that it is trivial to relax the assumption of symmetry for the DL components. In the context of VARs, this might be of particular interest for distinguishing diagonal ($a_D$ large) from off-diagonal ($a_O$ small) elements in the spirit of the Minnesota prior or increasing the amount of shrinkage with increasing lag order huber2016adaptive. Third, we would like to stress that our approach could also be used to estimate huge dimensional time-varying parameter VAR models with stochastic volatility. To cope with the computational difficulties associated with the vast state space, a possible approach could be to rely on an additional layer of hierarchy that imposes a (dynamic) factor structure on the time-varying autoregressive coefficients in the spirit of eisenstat2018reducing and thus reduce the computational burden considerably.
\printbibliography