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.
80,594 characters · 13 sections · 67 citation commands
Variational inference for large Bayesian vector autoregressions
\def\spacingset#1{ {#1}} \spacingset{1}
\thispagestyle{empty}
\centerline{\bf Abstract} We propose a novel variational Bayes approach to estimate high-dimensional vector autoregression (VAR) models with hierarchical shrinkage priors. Our approach does not rely on a conventional structural VAR representation of the parameter space for posterior inference. Instead, we elicit hierarchical shrinkage priors directly on the matrix of regression coefficients so that (1) the prior structure directly maps into posterior inference on the reduced-form transition matrix, and (2) posterior estimates are more robust to variables permutation. An extensive simulation study provides evidence that our approach compares favourably against existing linear and non-linear Markov Chain Monte Carlo and variational Bayes methods. We investigate both the statistical and economic value of the forecasts from our variational inference approach within the context of a mean-variance investor allocating her wealth in a large set of different industry portfolios. The results show that more accurate estimates translate into substantial statistical and economic out-of-sample gains. The results hold across different hierarchical shrinkage priors and model dimensions.
Keywords: Bayesian methods, variational inference, hierarchical shrinkage prior, high-dimensional models, vector autoregressions, industry returns predictability.
JEL codes: C11, C32, C55, C53, G11
\doublespacing
\pagenumbering{arabic}
\textcolor{black}{Hierarchical shrinkage priors have been shown to represent an effective regularization technique when estimating large vector autoregression (VAR) models.} The use of these priors often relies on a Cholesky decomposition of the residuals covariance matrix so that a large system of equations is reduced to a sequence of univariate regressions. \textcolor{black}{This allows for more efficient computations as priors can be elicited on the structural VAR representation implied by the Cholesky factorization and posterior inference is carried out equation-by-equation.}
Such a conventional approach has two important implications for posterior inference: first, priors are not order-invariant, meaning that posterior inference is sensitive to permutations of the endogenous variables for a given prior specification. This is particularly relevant in high dimensions whereby logical orders of the endogenous variables might be unclear or a full search among all possible ordering combinations might be unfeasible (see, e.g., chan2021large). \textcolor{black}{Second, imposing a shrinkage prior on the structural VAR formulation does not necessarily help to pin down the significance of cross-correlations in the reduced-form VAR formulation. This is especially relevant in forecasting applications whereby the main objective is to accurately identify predictive relationships across variables, rather than to identify structural shocks.}
{\color{black}In this paper, we take a different approach towards posterior inference with hierarchical shrinkage priors in large VAR models. Specifically, we propose a novel variational Bayes estimation approach which allows for fast and accurate estimates of the reduced-form regression coefficients without leveraging on a structural VAR representation.} This allows us to elicit hierarchical shrinkage priors directly on the matrix of regression coefficients so that (1) the prior structure directly maps into the posterior inference of the reduced-form transition matrix, and (2) posterior estimates are more robust to variables permutation. \textcolor{black}{We also account for the effect of “exogenous” covariates and stochastic volatility in the residuals.}
\textcolor{black}{The key feature of our approach is that by abstracting from the linearity constraints implied by a structural VAR formulation, one can provide a more direct identification of the reduced-form regression parameters.} This could have important implications for forecasting within the context of weak predictability whereby the transition matrix and/or the coefficients on exogenous predictors are potentially sparse in nature (see, e.g., bianchi2023dynamic). The main advantage of our variational inference approach is that an accurate identification of the regression parameters does not translate into a higher computational cost compared to existing Bayesian estimation methods. This is particularly relevant in practice for recursive forecasting implementations with higher frequency data, such as portfolio returns.
We investigate the accuracy of the posterior estimates based on an extensive simulation study for different model dimensions and variables permutation. As benchmarks, we consider a variety of established estimation approaches developed for large Bayesian VAR models, such as the linearized MCMC proposed by chan2018bayesian,cross2020macroeconomic and its variational Bayes counterpart proposed by chan_yu2020,gefang2023forecasting. Both approaches are built upon a structural VAR formulation. \textcolor{black}{In addition, we compare our variational Bayes method against the MCMC approach developed by gruber2022forecasting, which is not constrained by a Cholesky factorization for parameters identification, similar to our approach.} We test each estimation method for different hierarchical priors, such as the adaptive-Lasso of Leng.2014, an adaptive version of the Normal-Gamma of griffin_brown.2010, and the Horseshoe of carvalho_etal.2010.
Overall, the simulation results show that our variational inference approach represents the best trade-off between estimation accuracy and computational efficiency. Specifically, posterior inference from our variational Bayes method is as accurate as non-linear MCMC methods (see, e.g., gruber2022forecasting) but is considerably more efficient. At the same time, our approach is as efficient as conventional MCMC and variational Bayes methods based on a structural VAR formulation, but is considerably more accurate and less sensitive to variables permutation.
\textcolor{black}{Our approach towards posterior inference in large VARs is guided by the principle that a more accurate identification of the reduced-form transition matrix should ultimately lead to better out-of-sample forecasts and financial decision making.} To test this assumption, we investigate both the statistical and economic value of the forecasts from our variational Bayes approach within the context of a mean-variance investor who allocates her wealth between an industry portfolio and a risk-free asset based on lagged cross-industry returns and a series of macroeconomic predictors.
Although the model is general and can be applied to any type of financial returns, as far as data are stationary, our focus on different industry portfolios is motivated by a keen interest from researchers (see, e.g., fama1997industry,hou2006industry) and practitioners alike. Indeed, the implications of industry returns predictability are arguably far from trivial. If all industries are unpredictable, then the market return, which is a weighted average of the industry portfolios, should also be unpredictable. As a result, the abundant evidence of aggregate market return predictability (see, e.g., rapach2013forecasting), implies that at least some industry portfolio return is predictable.
The main results show that our variational inference approach fares better than competing methods in terms of out-of-sample point and density forecasts. We show that more accurate forecasts translate into larger economic gains as measured by certainty equivalent returns spreads vis-\'{a}-vis a naive investor which take investment decisions based on sample estimates of the conditional mean and variance of the returns. This holds across different hierarchical prior specifications. Overall, the empirical results support our view that by a more accurate identification of weak correlations between predictors and portfolio returns, one can significantly improve -- both statistically and economically -- the out-of-sample performance of large-scale multivariate time-series models.
Our paper connects to a growing literature exploring the use of Bayesian methods to estimate high-dimensional VAR models with shrinkage priors. A non-exhaustive list of works on the topic contains \citet*{chan2018bayesian,carriero2019large,huber2019adaptive,chan_yu2020,cross2020macroeconomic,kastner2020sparse,chan2021large,chan2021minnesota,carriero2022corrigendum,gruber2022forecasting,gefang2023forecasting}, among others. \textcolor{black}{We contribute to this literature by providing a fast and accurate variational Bayes method which generalize posterior inference of quantities of interest by abstracting from a conventional structural VAR representation.}
A second strand of literature we contribute to is related to the predictability of stock returns. More specifically, we contribute to the ongoing struggle to understand the dynamics of risk premiums by looking at industry-based portfolios. As highlighted by lewellen2010skeptical, the time series variation of industry portfolios is particularly problematic to measure, since conventional risk factors do not seem to capture significant comovements and cross-signals which might improve out-of-sample predictability. Early exceptions are ferson1991variation,ferson1995arbitrage,ferson1999conditioning and Avramov:2004. We extend this literature by investigating the out-of-sample predictability of industry portfolios through the lens of a novel estimation method for large Bayesian VAR models.
Let $\mathbf{y}_t=\left(y_{1,t},\dots,y_{d,t}\right)^\intercal\in\mathbb{R}^d$ be a multivariate normal random variable and denote by $\mathbf{x}_t=\left(1,x_{1,t},\dots,x_{p,t}\right)^\intercal\in\mathbb{R}^{(p+1)}$ a vector of covariates at time $t$. {\color{black}A vector autoregressive model with exogenous covariates and stochastic volatility is defined in compact form as:}
with $\mathbf{z}_{t-1} = (\mathbf{y}_{t-1}^\intercal,\mathbf{x}_{t-1}^\intercal)^\intercal$ and $\boldsymbol{\Theta} = (\boldsymbol{\Phi},\boldsymbol{\Gamma})$ consistently partitioned, where $\boldsymbol{\Phi}\in\mathbb{R}^{d\times d}$ is the transition matrix containing the autoregression coefficients and $\boldsymbol{\Gamma}\in\mathbb{R}^{d\times (p+1)}$ is the matrix of regression parameters for the exogenous predictors. Here, $\mathbf{u}_t\in\mathbb{R}^d$ is a sequence of uncorrelated innovation terms such that $\mathbf{u}_{t-k}\perp \mathbf{u}_{t-j}$ $\forall k,j$ with $k\neq j$ and {\color{black}$\boldsymbol{\Omega}_t\in\mathbb{S}^d_{++}$ being a symmetric and positive-definite time-varying precision matrix. A modified Cholesky factorization of $\boldsymbol{\Omega}_t$ can be conveniently exploited to re-write the model in Eq.(ref) with orthogonal innovations rothman_etal.2010.
Let $\boldsymbol{\Omega}_t = \mathbf{L}^\intercal\mathbf{V}_t\mathbf{L}$, where $\mathbf{L}\in\mathbb{R}^{d\times d}$ is unit-lower-triangular and $\mathbf{V}_t\in\mathbb{S}_{++}^d$ is diagonal with time-varying elements $\mathbf{V}_t=\text{Diag}(\nu_{1,t},\ldots,\nu_{d,t})$ huber2019adaptive,gefang2023forecasting.} By multiplying both sides of Eq.(ref) by $\mathbf{L}=\mathbf{I}_d - \mathbf{B}$ one can obtain two alternative re-parametrizations of the same model:
where $\mathbf{A}=\mathbf{L}\boldsymbol{\Theta}$ and $\mathbf{B}$ has a strict-lower-triangular structure with elements $\beta_{j,k}=-l_{j,k}$ for $j=2,\ldots,d$ and $k=1,\ldots,j-1$. The key difference is that Eq.(ref) is non-linear in the parameters, while Eq.(ref) is linear. More importantly, Eq.(ref) is known as structural VAR representation, widely used in existing MCMC and variational Bayes estimations methods for high-dimensional VAR models (see, e.g., chan2018bayesian,chan_yu2020,gefang2023forecasting). Instead, Eq.(ref) is the reduced-form parametrization at the core of our variational inference approach. This has also been used within the context of MCMC for smaller dimensions (see, e.g., huber2019adaptive,gruber2022forecasting).
From Eq.(ref) one can obtain an equation-by-equation representation in which the $j$-th component of $\mathbf{y}_t$ becomes:
for all $j=1,\ldots,d$ and $t=1,\ldots,T$, where $\boldsymbol{\beta}_j\in\mathbb{R}^{j-1}$ is a row vector containing the non-null elements in the $j$-th row of $\mathbf{B}$, $\boldsymbol{\vartheta}_j$ and $\mathbf{a}_j$ denote the $j$-th row of $\boldsymbol{\Theta}$ and $\mathbf{A}$, respectively. For any $j=1,\dots,d$, let $\mathbf{r}_{j,t}=\mathbf{y}_t^j - \boldsymbol{\Theta}^j \mathbf{z}_{t-1}$ denotes the the vector of residuals up to the $(j-1)$-th regression, with $\mathbf{y}_t^j = (y_{1,t},\ldots,y_{j-1,t})^\intercal\in\mathbb{R}^{j-1}$ being the sub-vector of $\mathbf{y}_t$ collecting the variables up to the $(j-1)$-th and $\boldsymbol{\Theta}^j\in\mathbb{R}^{(j-1)\times d}$ is the sub-matrix containing the first $j-1$ rows of $\boldsymbol{\Theta}$. {\color{black}We follow gefang2023forecasting,chan_yu2020 and model the time variation in $\nu_{j,t}^{-1} = \exp\left(h_{j,t}\right)$ assuming a log-volatility process $h_{j,t}=h_{j,t-1}+e_{j,t}$ with $e_{j,t}\sim\mathsf{N}(0,\psi_j)$, where the initial state $h_{0,j}\sim\mathsf{N}(0,k_0\,\psi_j)$, $k_0\gg 0$, is unknown.}
\paragraph{A discussion on variables permutation.} \textcolor{black}{Existing Bayesian approaches for large VAR models often rely on the structural representation in Eq.(ref), and therefore consider the elements in $\mathbf{A}$ as the parameters of interest.} This has the key merit of simplifying the implementation of MCMC (see, e.g., chan2018bayesian) and variational Bayes algorithms (see, e.g., gefang2023forecasting). Under the re-parametrization $\mathbf{A}=\mathbf{L}\boldsymbol{\Theta}$, each element $\vartheta_{i,j}$ -- which denotes the $(i,j)$-entry of $\mathbf{\Theta}$ -- is a linear combination $\vartheta_{i,j}=a_{i,j}+\sum_{k=1}^{i-1}c_{i,k}a_{k,j}$, where $a_{i,j}$ and $c_{i,j}$ are the $(i,j)$-entry of $\mathbf{A}$ and $\mathbf{L}^{-1}$, respectively.
This raises two main issues: first, $a_{i,j}=0$ does not imply $\vartheta_{i,j}=0$, that is a shrinkage prior on $\mathbf{A}$ does not preserve the structure of $\boldsymbol{\Theta}$. Second, the estimate $\widehat{\boldsymbol{\Theta}}=\widehat{\mathbf{L}}^{-1}\widehat{\mathbf{A}}$ for a given prior is potentially highly sensitive to variables permutation due to its dependence on the Cholesky factorization (see gruber2022forecasting for a related discussion). Figure (ref) provides a visual representation of this argument by comparing the estimates obtained based on Eq.(ref) vs Eq.(ref), for two different permutations of $\mathbf{y}_t$.
The evidence confirms that the estimates based on the transformation $\widehat{\boldsymbol{\Theta}}=\widehat{\mathbf{L}}^{-1}\widehat{\mathbf{A}}$ clearly diverge from the true $\mathbf{\Theta}$. In addition, the posterior estimates are influenced by the variables permutation. \textcolor{black}{Instead, inference based on the representation in Eq.(ref) provides a more accurate identification of $\mathbf{\Theta}$ which is also less sensitive to variables permutation.} Before taking this intuition to task both in simulation and on actual forecasting, in the next Section we provide details of our variational Bayes inference approach.
A variational approach to Bayesian inference requires to minimize the Kullback-Leibler ($\mathit{KL}$) divergence between an approximating density $q(\boldsymbol{\xi})$ and the true posterior density $p(\boldsymbol{\xi}|\mathbf{y})$, where $\boldsymbol{\xi}$ denotes the set of parameters of interest. ormerod_wand.2010 show that minimizing the $\mathit{KL}$ divergence can be equivalently stated as the maximization of the “effective lower bound” (ELBO) denoted by $\underline{p}\left(\mathbf{y};q\right)$:
where $q^*(\boldsymbol{\xi})\in\mathcal{Q}$ represents the optimal variational density and $\mathcal{Q}$ is a space of density functions. \textcolor{black}{Depending on the assumption on $\mathcal{Q}$, one falls into different variational paradigms. For instance, given a partition of the parameters vector $\boldsymbol{\xi}=\{ \boldsymbol{\xi}_1,\dots,\boldsymbol{\xi}_p\}$, a mean-field variational Bayes (MFVB) approach assumes a factorization of the form $q(\boldsymbol{\xi})=\prod_{j=1}^p q_i(\boldsymbol{\xi}_j)$.} A closed form expression for each optimal variational density $q^\ast(\boldsymbol{\xi}_j)$ can be defined as:
where the expectation is taken with respect to the joint approximating density with the $j$-th element of the partition removed $q^\star(\boldsymbol{\xi}\setminus\boldsymbol{\xi}_j)$. This allows to implement an efficient iterative algorithm to estimate the optimal density $q^*(\boldsymbol{\xi})$, although some components $q^*(\boldsymbol{\xi}_j)$ may remain too complex to handle and further restrictions are needed. If we assume that $q^*(\boldsymbol{\xi}_j)$ belongs to a pre-specified parametric family of distributions, the MFVB outlined above is sometimes labelled as {\it semi-parametric} rohde2016semiparametric.
We present a factorization of the variational density $q(\boldsymbol{\xi})$ for the model outlined in Eq.(ref). As a benchmark, we consider a non-informative Normal prior for the regression coefficients. For each entry of $\boldsymbol{\Theta}$, let $\vartheta_{j,k} \sim \mathsf{N}(0,\upsilon)$, for $j=1,\dots,d$ and $k=1,\dots,d+p+1$. In addition, let $\psi_j \sim \mathsf{InvGa}(a_{\psi},b_{\psi})$ for $j=1,\dots,d$, and $\beta_{j,k} \sim \mathsf{N}(0,\tau)$, for $j=2,\dots,d$ and $k=1,\dots,j-1$. Here, $\mathsf{InvGa}(\cdot,\cdot)$ denotes the Inverse-Gamma distribution, and $a_\psi>0$, $b_\psi>0$, $\tau \gg 0$ and $\upsilon \gg 0$ are the related hyper-parameters. Let $\boldsymbol{\xi}=(\boldsymbol{\vartheta}^\intercal,\mathbf{h}^\intercal,\mbox{\boldmath $\psi$}^\intercal,\mbox{\boldmath $\beta$}^\intercal)^\intercal$ be the set of parameters of interest, the corresponding variational density can be factorised as $q(\boldsymbol{\xi}) = q(\boldsymbol{\vartheta})q(\mathbf{h})q(\boldsymbol{\psi})q(\boldsymbol{\beta})$, where:
For the ease of exposition, in the main text of the paper we summarize the optimal variatonal density for the main parameters of interest $\boldsymbol{\Theta}$, with both a baseline non-informative prior and three alternative hierarchical shrinkage priors. \textcolor{black}{The parameters and the full derivations of the optimal variational densities $q^*(\mathbf{h}_j) \equiv \mathsf{N}_{T+1}(\boldsymbol{\mu}_{q(h_j)}, \mathbf{\Sigma}_{q(h_j)})$, $q^*(\psi_j) \equiv \mathsf{InvGa}(a_{q(\psi_j)}, b_{q(\psi_j)})$, and $q^*(\boldsymbol{\beta}_j) \equiv \mathsf{N}_{j-1}(\boldsymbol{\mu}_{q(\beta_j)}, \boldsymbol{\Sigma}_{q(\beta_j)})$ for $j=1,\ldots,d$, are reported in Proposition (ref), (ref) and (ref) of Appendix (ref), respectively. Notice these optimal variational densities are invariant across different shrinkage prior specifications for $\boldsymbol{\Theta}$. We leave to Proposition (ref) in Appendix (ref) also the derivations for the constant volatility case with $\nu_{j,t}=\nu_j$ and $\nu_j \sim \mathsf{Ga}(a_{\nu},b_{\nu})$ for $j=1,\dots,d$, where $\mathsf{Ga}(\cdot,\cdot)$ denotes the gamma distribution, and $a_\nu>0$, $b_\nu>0$. For the interested reader, Appendix (ref) also provides the analytical form of the lower bound for each set of parameters.}
Proposition (ref) provides the optimal variational density for the $j$-th row of $\boldsymbol{\Theta}$ under the baseline Normal prior specification $\vartheta_{j,k} \sim \mathsf{N}(0,\upsilon)$. The proof and analytical derivations are available in Appendix (ref).
Notice that despite the multivariate model is reduced to a sequence of univariate regressions, the analytical form of the variational mean $\boldsymbol{\mu}_{q(\vartheta_j)}$ in Proposition (ref) depends on all the other rows through $\mbox{\boldmath $\mu$}_{q(\mathbf{\vartheta}_{-j})}$. As a result, the variational estimates of $\boldsymbol\vartheta_j$ explicitly depend on all of the other $\boldsymbol\vartheta_{-j}$. This addresses the issue in the MCMC algorithm of carriero2019large, which has been highlighted by bognanni2022comment and corrected by carriero2022corrigendum.
\paragraph{Bayesian adaptive-Lasso.} The Bayesian adaptive-Lasso of Leng.2014 extends the original work of Park_Casella.2008 by assuming a different shrinkage for each regression parameter based on a laplace distribution with an individual scaling parameter $\vartheta_{j,k}\vert\lambda_{j,k} \sim\mathsf{Lap}(\lambda_{j,k})$, for $j=1,\dots,d$ and $k=1,\dots,d+p+1$. The latter can be represented as a scale mixture of normals with an exponential mixing density, $\vartheta_{j,k}|\upsilon_{j,k} \sim \mathsf{N}(0,\upsilon_{j,k})$, $\upsilon_{j,k}|\lambda^2_{j,k} \sim \mathsf{Exp}(\lambda^2_{j,k}/2)$. The scaling parameters $\lambda^2_{j,k}$ are not fixed but inferred from the data by assuming a common hyper-prior distribution $\lambda^2_{j,k} \sim \mathsf{Ga}(h_1,h_2)$, where $h_1, h_2>0$.
Let $\boldsymbol{\xi}_{\text{L}}=(\boldsymbol{\xi}^\intercal, \boldsymbol{\upsilon}^\intercal,(\boldsymbol{\lambda}^2)^\intercal))^\intercal$ be the vector $\boldsymbol{\xi}$ augmented with the adaptive-Lasso prior parameters. The distribution $q(\boldsymbol{\xi}_{\text{L}})$ can be factorised as,
Proposition (ref) provides the optimal variational density for the $j$-th row of $\boldsymbol{\Theta}$ under Bayesian adaptive-Lasso prior specification $\vartheta_{j,k}|\upsilon_{j,k} \sim \mathsf{N}(0,\upsilon_{j,k})$, $\upsilon_{j,k}|\lambda^2_{j,k} \sim \mathsf{Exp}(\lambda^2_{j,k}/2)$, and $\lambda^2_{j,k} \sim \mathsf{Ga}(h_1,h_2)$. The proof and analytical derivations are available in Appendix (ref).
\paragraph{Adaptive Normal-Gamma.} We expand the original Normal-Gamma prior of griffin_brown.2010 by assuming that each regression coefficient has a different shrinkage parameter, similar to the adaptive-Lasso. The hierarchical specification requires that $\vartheta_{j,k}|\upsilon_{j,k} \sim \mathsf{N}(0,\upsilon_{j,k})$, and $\upsilon_{j,k}|\eta_j,\lambda_{j,k} \sim \mathsf{Ga}\left(\eta_j,\eta_j\lambda_{j,k}/2\right)$ for $j=1,\dots,d$ and $k=1,\dots,d+p+1$. Notice that by restricting $\eta_j=1$ one could obtain the adaptive-Lasso prior. Marginalization over the variance $\upsilon_{j,k}$ leads to $p(\vartheta_{j,k}\vert\eta_j,\lambda_{j,k})$ which corresponds to a Variance-Gamma distribution. The hyper-parameters $\eta_j$ and $\lambda_{j,k}$ are not fixed but are inferred from the data by assuming two common hyper-priors $\lambda_{j,k} \sim \mathsf{Ga}(h_1,h_2)$ and $\eta_j \sim \mathsf{Exp}(h_3)$, where $h_l>0$ for $l=1,2,3$.
Let $\boldsymbol{\xi}_{\text{NG}}=(\boldsymbol{\xi}^\intercal, \boldsymbol{\upsilon}^\intercal,\boldsymbol{\lambda}^\intercal,\boldsymbol{\eta}^\intercal)^\intercal$ be the vector $\boldsymbol{\xi}$ augmented with the parameters of the adaptive Normal-Gamma prior. The joint distribution $q(\boldsymbol{\xi}_{\text{NG}})$ can be factorised as,
Proposition (ref) provides the optimal variational density for the $j$-th row of $\boldsymbol{\Theta}$ under an adaptive Normal-Gamma specification $\upsilon_{j,k}|\eta_j,\lambda_{j,k} \sim \mathsf{Ga}\left(\eta_j,\eta_j\lambda_{j,k}/2\right)$, $\lambda_{j,k} \sim \mathsf{Ga}(h_1,h_2)$ and $\eta_j \sim \mathsf{Exp}(h_3)$. The proof and analytical derivations are available in Appendix (ref).
Notice that the optimal density for the parameter $\eta_j$ is not a known distribution function. Proposition (ref) in Appendix (ref) provides an analytical approximation of its moments so that the optimal density can be calculated via numerical integration. \paragraph{Horseshoe prior.} As a third hierarchical shrinkage prior we consider the Horseshoe prior as proposed by Carvalho.2009,carvalho_etal.2010. This is based on the hierarchical specification $\vartheta_{j,k}|\upsilon^2_{j,k}$, $\gamma^2 \sim \mathsf{N}(0,\gamma^2\upsilon^2_{j,k})$, $\gamma \sim \mathsf{C}^+(0,1)$, $\upsilon_{j,k}\sim \mathsf{C}^+(0,1)$, where $\mathsf{C}^+(0,1)$ denotes the standard half-Cauchy distribution with probability density function equal to $f(x) = 2/\{\pi(1 + x^2)\}\mathbbm{1}_{(0,\infty)}(x)$. The Horseshoe is a global-local prior that implies an aggressive shrinkage of weak signals without affecting the strong ones polson_etal.2011. We follow Wand.2011 and leverage on a scale mixture representation of the half-Cauchy distribution as,
where the local and global shrinkage parameters are $\upsilon^2_{j,k}$ and $\gamma^2$ respectively.
Let $\boldsymbol{\xi}_{\text{HS}}=(\boldsymbol{\xi}^\intercal,(\boldsymbol{\upsilon}^2)^\intercal,\gamma^2,\boldsymbol{\lambda}^\intercal,\eta)^\intercal$ be the vector $\boldsymbol{\xi}$ augmented with the parameters of the Horseshoe prior. The joint distribution $\boldsymbol{\xi}_{\text{HS}}$ can be factorized as,
Proposition (ref) provides the optimal variational density for the $j$-th row of $\boldsymbol{\Theta}$ under the Horseshoe prior outlined in Eq.(ref). The proof and analytical derivations are available in Appendix (ref).
In addition to computational tractability, shrinking rather than selecting is a defining feature of the hierarchical priors outlined in Section (ref). That is, posterior estimates of $\mathbf{\Theta}$ are non-sparse, and thus can not provide exact differentiation between significant vs non-significant predictors. The latter is particularly relevant since we ultimately want to assess the accuracy of our variational inference approach -- versus existing MCMC and variational Bayes algorithms -- in identifying the exact structure of $\mathbf{\Theta}$.
To address this issue, we build upon pallavi_battacharya2019savs and implement a Signal Adaptive Variable Selector (SAVS) algorithm to induce sparsity in $\widehat{\mathbf{\Theta}}$, conditional on a given prior. The SAVS is a post-processing algorithm which divides signals and nulls on the basis of the point estimates of the regression coefficients \citep*[see, e.g.,][]{hauzenberger2021combining}. Specifically, let $\widehat{\vartheta}_j$ the posterior estimate of $\vartheta_j$ and $\mathbf{z}_j$ the associated vector of covariates. If $|\widehat{\vartheta}_j|\,||\mathbf{z}_j||^2\leq|\widehat{\vartheta}_j|^{-2}$ we set $\widehat{\vartheta}_j=0$, where $||\cdot||$ denotes the euclidean norm.
{\color{black} The reason why we rely on the SAVS post-processing to induce sparsity in the posterior estimates is threefold. First, as highlighted by pallavi_battacharya2019savs, the SAVS represents an automatic procedure in which the sparsity-inducing property directly depends on the effectiveness of the shrinkage performed on $\widehat{\vartheta}_j$. This refers to the precision of the posterior mean estimates; that is, the more accurate is $\widehat{\vartheta}_j$, the more precise is the identification of the non-zero elements in $\mathbf{\Theta}$. Second, the SAVS is “agnostic” with respect to the shrinkage prior or estimation approach adopted, so it represents a natural tool to compare different estimation methods. Third, it is decision theoretically motivated as it grounds on the idea of minimizing the posterior expected loss \citep*[see, e.g.,][]{huber2021inducing}.
In addition to SAVS, we also expand on hahn2015decoupling (HC henceforth) and provide a multivariate extension to their least-angle regression which has originally been built for univariate regressions. Appendix (ref) provides the full derivation of our extended HC approach as well as a complete discussion of the drawbacks compared to SAVS. In addition, for the interested reader, Appendix (ref) provides a direct comparison between the SAVS and our multivariate extension to hahn2015decoupling based on simulated data (see also the discussion in Section (ref)).}
Consider the posterior distribution $p(\boldsymbol{\xi}|\mathbf{z}_{1:t})$ given the information set $\mathbf{z}_{1:t}=\left\{\mathbf{y}_{1:t},\mathbf{x}_{1:t}\right\}$ and the conditional likelihood $p(\mathbf{y}_{t+1}|\mathbf{z}_{t},\boldsymbol{\xi})$. A standard predictive density takes the form,
Given an optimal variational density $q^*(\boldsymbol{\xi})$ that approximates $p(\boldsymbol{\xi}|\mathbf{z}_{1:t})$, we follow gunawan_etal.2020 and obtain the variational predictive distribution
Although an analytical expression for Eq.(ref) is not available, a simulation-based estimator for $q(\mathbf{y}_{t+1}|\mathbf{z}_{1:t})$ can be obtained through Monte Carlo integration by averaging $p(\mathbf{y}_{t+1}|\mathbf{z}_{t},\boldsymbol{\xi}^{(i)})$ over the draws $\boldsymbol{\xi}^{(i)}\sim q^\ast(\boldsymbol{\xi})$, such that $\widehat{q}(\mathbf{y}_{t+1}|\mathbf{z}_{1:t}) = N^{-1}\sum_{i=1}^N p(\mathbf{y}_{t+1}|\mathbf{z}_{t},\boldsymbol{\xi}^{(i)})$. {\color{black}Notice that a complete characterization of the optimal variational predictive density entails $q^*(\mathbf{\Omega}_t)$ with $\mathbf{\Omega}_t=\mathbf{L}^\intercal\mathbf{V}_t\mathbf{L}$. Proposition (ref) shows that, conditional on $\mathbf{L}$ and $\mathbf{V}_t$, the optimal distribution of $\mathbf{\Omega}_t$ can be approximated by a $d$-dimensional Wishart distribution $\mathsf{Wishart}_d(\delta_t,\mathbf{H}_t)$, where $\delta_t$ and $\mathbf{H}_t$ are the degrees of freedom parameter and the scaling matrix, respectively.}
The complete proof is available in Appendix (ref) and is based on the Expectation Propagation (EP) approach proposed by minka.2001. In order to implement this approach, there is no need to know $q^*(\mathbf{\Omega}_t)$, but it is sufficient to be able to compute $\mathbb{E}_q(\mathbf{\Omega}_t)$. The latter can be reconstructed based on the optimal variational densities of the Cholesky factor $q^*(\mbox{\boldmath $\beta$})$ -- and therefore for $\mathbf{L}$ -- and of $q^*(\mathbf{V}_t)$. The simulation results in Appendix (ref) show that the proposed Wishart distribution provides an accurate approximation of $q^*(\mathbf{\Omega}_t)$ for both small and large dimensional models.
{\color{black}Based on Proposition (ref), we can further simplify Eq.(ref) by integrating $\mathbf{\Omega}_t$ such that:
where $h(\mathbf{y}_{t+1}|\mathbf{z}_t,\boldsymbol{\vartheta})$ denotes the probability density function of a multivariate Student-$t$ distribution $\mathsf{t}_v(\mathbf{m},\mathbf{S})$ with mean $\mathbf{m}=\mathbf{\Theta}\mathbf{z}_t$, scaling matrix $\mathbf{S}=(v\widehat{\mathbf{H}})^{-1}$, and degrees of freedom parameter $v=\widehat{\delta}-d+1$. As a result, the predictive distribution can be approximated by averaging the density of the multivariate Student-$t$ $h(\mathbf{y}_{t+1}|\mathbf{z}_t,\boldsymbol{\vartheta}^{(i)})$ over the draws $\boldsymbol{\vartheta}^{(i)}\sim q^\ast(\boldsymbol{\vartheta})$, for $i=1,\ldots,N$, such that $\widehat{q}(\mathbf{y}_{t+1}|\mathbf{z}_{1:t}) = N^{-1}\sum_{i=1}^N h(\mathbf{y}_{t+1}|\mathbf{z}_{t},\boldsymbol{\vartheta}^{(i)})$. This allows for a more efficient sampling from the predictive density.}
\textcolor{black}{Notice that the main advantage of the approximation obtained from Proposition 3.5 is to allow for a considerably faster computation of the variational predictive density, compared to using $q^*(\mathbf{L})$ and $q^*(\mathbf{V}_t)$ as stationary distributions to sample $\mbox{\boldmath $\Omega$}_t$, similar to an MCMC. This is because the scaling matrix of the Wishart distribution is available in closed form and the computation of degrees of freedom requires only a one-dimensional optimization. In Appendix (ref) we discuss a further simplification that minimizes the KL divergence between the multivariate Student-$t$ and a multivariate Normal distribution.}
In this section, we report the results of an extensive simulation study designed to compare the properties of our estimation approach against both MCMC and variational Bayes methods for large VAR models. To begin, we compare our {\tt VB} algorithm against the MCMC approach of chan2018bayesian,cross2020macroeconomic and the variational inference framework proposed by chan_yu2020,gefang2023forecasting. Both these approaches are built upon the structural VAR representation in Eq.(ref). {\color{black} Then, we also compare our {\tt VB} method against the MCMC approach developed by huber2019adaptive,gruber2022forecasting which is based upon a non-linear parametrization as in Eq.(ref), similar to our approach.}
{\color{black}For the sake of comparability with gruber2022forecasting,gefang2023forecasting, which do not consider the presence of exogenous predictors, we consider a standard VAR(1) as data generating process. Consistent with the empirical implementations, we set $T=360$ and $d=30,49$. The choice of $d$ is due to the two alternative industry classifications which are explored in the main empirical analysis.} We assume either a moderate -- $50\%$ of zeros -- or a high -- $90\%$ of zeros -- level of sparsity in the true matrix $\boldsymbol{\Theta}$. {\color{black} The latter is generated as follows: we fix to zero $s\cdot d^2$ entries at random, with $s=0.5,0.9$ and $d=30,49$, while the remaining non-zero coefficients are sampled from a mixture of two normal distributions with means equal to $\pm 0.08$ and standard deviation $0.1$. Appendix (ref) provides additional details on the data generating process and additional simulation results for $d=15$.}
As a measure of point estimation accuracy, we first look at the Frobenius norm $\Vert\boldsymbol{\Theta}-\widehat{\boldsymbol{\Theta}}\Vert_F$, which measures the difference between the true $\boldsymbol{\Theta}$ observed at each simulation and its estimate $\widehat{\boldsymbol{\Theta}}$. {\color{black}In addition, we compare the ability of each estimation method to identify the non-zero elements in the true $\boldsymbol{\Theta}$ based on the F1 score. The latter can be expressed as a function of counts of true positives ($tp$), false positives ($fp$) and false negatives ($fn$),
The F1 score takes value one if identification is perfect, i.e., no false positives and no false negatives, and zero if there are no true positives. We compute both measures of estimation accuracy on $N=100$ replications to compare each estimation method and prior specification. The estimates from the MCMC specifications are based on 5,000 posterior simulations, after discarding the first 5,000 as a burn-in sample.}
\paragraph{Point estimates.} Figure (ref) shows the box charts summarizing the Frobenius norm $\Vert\boldsymbol{\Theta}-\widehat{\boldsymbol{\Theta}}\Vert_F$ across $N=100$ replications. We label the linearized MCMC and variational methods with {\tt LMCMC} and {\tt LVB}, respectively, with {\tt MCMC} the non-linear method of gruber2022forecasting and with {\tt VB} our variational inference method, respectively. To increase readability, we separate the results by prior and color-code the four different estimation methods. For instance, for a given sub-plot we report the results for the Normal, adaptive-Lasso, adaptive Normal-Gamma and Horseshoe priors from the left to the right panel. Within each panel, the simulation results for the {\tt LMCMC}, {\tt LVB}, {\tt MCMC} and {\tt VB} estimates are reported in red, yellow, light-blue and green, respectively.
Beginning with the moderate sparsity case (top panels), the simulation results show that {\tt LMCMC} and {\tt LVB} approaches tend to perform equally across different shrinkage priors, with the only exception of the Normal-Gamma prior, in which {\tt LMCMC} slightly outperforms {\tt LVB}. However, the discrepancy between the two structural VAR representation methods tend to increase when sparsity becomes more pervasive (see bottom panels).
Overall, the simulation results support our view that, by eliciting shrinkage priors directly on $\boldsymbol{\Theta}$ -- as per the parametrization in Eq.(ref) -- the accuracy of the posterior estimates improves. The mean squared errors obtained from {\tt MCMC} and {\tt VB} are lower compared to both {\tt LMCMC} and {\tt LVB}. This holds for all priors and the model dimension. The accuracy with $d=30$ of the {\tt MCMC} and {\tt VB} is virtually the same. Yet, with $d=49$ our {\tt VB} produces slightly more accurate estimates than {\tt MCMC} for both the adaptive-Lasso and the Horseshoe prior.
\paragraph{Sparsity identification.} Figure (ref) shows the box charts of F1 scores across $N=100$ simulations. The labeling is the same as in Figure (ref). Both {\tt LMCMC} and {\tt LVB} produce a rather dismal identification of the non-zero elements in $\mathbf{\Theta}$ across prios and model dimensions. This is due to the fact that $\widehat{\boldsymbol{\Theta}}=\widehat{\mathbf{L}}^{-1}\widehat{\mathbf{A}}$ in Eq.(ref), so that a sparse estimate of $\widehat{\mathbf{A}}$ does not map into a sparse estimate of $\widehat{\boldsymbol{\Theta}}$, and therefore produces a lower accuracy in identifying the non-zero coefficients in the true $\boldsymbol{\Theta}$. As the level of sparsity increases, the divergence between $\mathbf{A}$ and $\boldsymbol{\Theta}$ increases.
Consistent with our argument in favor of the parametrization in Eq.(ref), both the {\tt MCMC} and {\tt VB} approaches produce a more accurate identification of the non-zero coefficients in $\mathbf{\Theta}$, as shown by the F1 score. The gap between {\tt LMCMC}, {\tt LVB} versus {\tt MCMC} and {\tt VB} becomes larger for higher levels of sparsity. This result holds across different hierarchical shrinkage priors and for different VAR dimensions. Yet, our {\tt VB} approach turns out to be more accurate than {\tt MCMC} under the adaptive-Lasso and Horseshoe priors for higher levels of sparsity.
{\color{black}As outlined in Section (ref), sparsity in the posterior estimates for $\widehat{\boldsymbol{\Theta}}$ for different hierarchical shrinkage priors is induced in the simulation results by using the SAVS algorithm of pallavi_battacharya2019savs. Appendix (ref) provides additional simulation results obtained by implementing a multivariate version of the post-processing method proposed by hahn2015decoupling as an alternative to the SAVS. A full derivation is provided in Appendix (ref). The F1 scores are largely the same across methods; in fact, the evidence is even more in favour of our {\tt VB}, compared to its {\tt MCMC} counterpart when using the extended hahn2015decoupling approach: our {\tt VB} is more accurate than {\tt MCMC} with a Normal-Gamma prior.}
\paragraph{Computational efficiency.} chan_yu2020 and gefang2023forecasting highlight that one of the main advantages of variational Bayes methods is computational efficiency. {\color{black}Figure (ref) reports the computational time -- expressed in a log-minute scale -- required by each estimation approach under different shrinkage priors. To highlight the performance for a given prior, we separate the results by estimation methods and color-code the four different shrinkage priors. For instance, for a given sub-plot, we report the results for the {\tt LMCMC}, {\tt LVB}, {\tt MCMC} and {\tt VB} estimates from left to right panel. Within each panel, the Normal, adaptive-Lasso, adaptive Normal-Gamma, and Horseshoe priors are colored in shades of gray from light (left) to dark (right) grey, respectively. To guarantee a more accurate comparability, we re-coded all competing methods in {\bf Rcpp} and use the same 2.5 GHz Intel Xeon W-2175 with 32GB of RAM for all implementations.}
{\color{black}The results highlight that our {\tt VB} approach has a clear computational advantage compared to both linear and non-linear MCMC methods. For instance, for $d=30$ our {\tt VB} is more than 100 times faster than the {\tt MCMC} of gruber2022forecasting and more than 10 times faster than the {\tt LMCMC} of cross2020macroeconomic, respectively. The gap in favour of our {\tt VB} method compared to both {\tt LMCMC} and {\tt MCMC} increases in larger dimensions; for $d=49$ the {\tt MCMC} approach takes almost 60 minutes, on average, to generate comparably accurate posterior estimates to our {\tt VB}, which instead takes approximately between 30 to 40 seconds, on average. Such efficiency gap between {\tt VB} and {\tt MCMC} has profound implications for a practical forecasting implementation, especially within the context of recursive predictions with higher frequency data such as stock returns (see Section (ref)). Perhaps not surprisingly, the {\tt LVB} approach of chan_yu2020,gefang2023forecasting is highly competitive in terms of computational efficiency. However, being built on a structural VAR formulation, we showed in Figures (ref) and (ref) that such computational efficiency comes at the cost of a lower estimation accuracy.}
\textcolor{black}{Appendix (ref) also provides a broader qualitative discussion on the computational costs of some of the existing MCMC approaches. Specifically, we review some of the results reported in the original papers and show that these largely align with our own findings. In addition, we also discuss some of the limitations of the non-linear {\tt MCMC} for the recursive forecasting implementation (see Section (ref) for more details).}
\paragraph{Robustness to variables permutation.} {\color{black}At the outset of the paper, we argue that a conventional structural VAR formulation potentially generates posterior estimates which are not permutation-invariant. That is, posterior estimates of $\mathbf{\Theta}$ are sensitive to the ordering imposed on the target variables $\mathbf{y}_t$, conditional on a given prior. To highlight this issue, in Appendix (ref), we report a set of additional simulation results for all estimation methods and shrinkage priors under variables permutation.
The results show that the accuracy of the posterior estimates from both {\tt LMCMC} and {\tt LVB} changes once the variables ordering is reversed (see Figure (ref)). This is especially clear for the Normal-Gamma and Horseshoe priors, and when the amount of zero coefficients in $\boldsymbol{\Theta}$ is more pervasive. On the other hand, the estimation accuracy of both the {\tt MCMC} approach of gruber2022forecasting and our {\tt VB} method does not substantially deteriorates by arbitrarily changing ordering of the target variables. Overall a substantially higher computational efficiency coupled with a comparable accuracy with complex MCMC, makes our {\tt VB} extremely competitive within the context of recursive forecasts with higher frequency data.}
We investigate both the statistical and economic value of our variational Bayes approach within the context of US industry returns predictability. \textcolor{black}{To expand the scope of the testing framework, we consider two alternative industry aggregations: $d=30$ industry portfolios from July 1926 to May 2020, and a larger cross section of $d=49$ industry portfolios from July 1969 to May 2020. The size of the cross sections change due to a different industry classification.} At the end of June of year $t$ each NYSE, AMEX, and NASDAQ stock is assigned to an industry portfolio based on its four-digit SIC code at that time. Thus, the returns on a given value-weighted portfolio are computed from July of $t$ to June of $t+1$. The sample periods cover major events, from the great depression to the Covid-19 outbreak.
In addition to cross-industry portfolio returns, we consider a variety of predictors, such as the returns on the market portfolio ({\tt mkt}), and the returns on four alternative long-short investment strategies based on market capitalization ({\tt smb}), book-to-market ratios ({\tt hml}), operating profitability ({\tt rmw}) and firm investments ({\tt cma}) (see fama2015five). We also consider a set of additional macroeconomic predictors from Goyal2008, such as the log price-dividend ratio ({\tt pd}), the difference between the long term yield on government bonds and the T-bill ({\tt term}), the BAA-AAA bond yields difference ({\tt credit}), the monthly log change in the CPI ({\tt infl}), the aggregate market book-to-market ratio ({\tt bm}), the net-equity issuing activity ({\tt ntis}) and the corporate bond returns ({\tt corpr}).
In order to highlight some of the main properties of different estimation methods, we first report the in-sample estimates of $\mathbf{\Theta}$ for the $d=49$ industry case across all priors. Figure (ref) compares $\widehat{\mathbf{\Theta}}$ based on the full sample obtained from the {\tt LMCMC} and the {\tt LVB} with constant volatility, and our {\tt VB} with and without stochastic volatility. Appendix (ref) reports the additional in-sample estimates for $d=30$ industry portfolios.
The in-sample estimates highlight three key results. First, there are visible differences across shrinkage priors. For instance, the Horseshoe tend to shrink parameters more aggressively towards zero so that $\widehat{\mathbf{\Theta}}$ is more sparse compared to, for e.g., the adaptive Normal-Gamma. Second, consistent with gefang2023forecasting, the estimates of the {\tt LMCMC} and {\tt LVB} tend to be closely related. Yet, these in-sample estimates are substantially different compared to our {\tt VB} approach. This is due to the re-parametrization $\widehat{\boldsymbol{\Theta}}=\widehat{\mathbf{L}}^{-1}\widehat{\mathbf{A}}$ in Eq.(ref); that is, the estimated $\widehat{\mathbf{A}}$ is not translation-invariant, unlike in our approach. Third, with the exception of the adaptive-Lasso prior, the estimates $\widehat{\boldsymbol{\Theta}}$ from {\tt VB} are remarkably stable between constant vs stochastic volatility specifications.
Intuitively, different estimates of $\mathbf{\Theta}$ should reflect in different conditional forecasts. To test this intuition we now compare the {\tt LMCMC}, {\tt LVB} and the {\tt VB} estimation approaches with and without stochastic volatility. For the sake of completeness, we also consider a series of univariate model specifications ({\tt U} henceforth), which corresponds to assuming conditional independence across industry portfolios. We consider a 360 months rolling window period for each model estimation; for instance for the 30-industry classification the out-of-sample period is from July 1957 to May 2020.
Notice that given the recursive nature of the empirical implementation we do not consider the {\tt MCMC} approach of gruber2022forecasting. This is because the computational cost would make such implementation prohibitive in practice, as discussed in the simulation study based on Figure (ref). For instance, on a 2.5 GHz Intel Xeon W-2175 with 32GB of RAM and 14 cores it would take $20\ \text{min}\times 767\ \text{forecasts} \times 4\ \text{priors}= 61,360$ minutes, or 42 days, to implement the {\tt MCMC} approach for recursive forecasting for the 30 industry portfolios with constant volatility. The computational cost would be even more prohibitive when adding stochastic volatility and/or for the 49 industry portfolios. Appendix (ref) provides an additional discussion on the computational costs of some of the existing MCMC approaches and the key relevance for a higher-frequency forecasting implementation such as ours.
\paragraph{Point forecasts.} We begin by inspecting the accuracy of point forecasts for each industry based on the out-of-sample predictive R squared (see, e.g., Goyal2008),
where $t_0$ is the date of the first prediction, $\overline{y}_{jt}$ is the naive forecast from the recursive mean -- using the same rolling window of observations -- and $\widehat{y}_{jt}\left(\mathcal{M}_s\right)$ is the conditional mean returns for industry $j=1,\ldots,d$ for a given model $\mathcal{M}_s$.
{\color{black}The left panels of Figure (ref) show the box charts with the distribution of the $R_{j,oos}^2$ across $j=1,\ldots,d$ industries. For a given sub-plot the results for the Normal, Bayesian Lasso, Normal-Gamma and Horseshoe priors are reported from the left to the right. Within each panel of a sub-plot, the forecasting results for the {\tt U}, {\tt LMCMC}, {\tt LVB}, and {\tt VB} estimates are color coded in orange, red, yellow, and green (from left to right), respectively. The vertical dashed line within each panel separates between constant and stochastic volatility specifications. Based on the same separation across methods and priors, the right panels of Figure (ref) report a breakdown of the industries for which the corresponding $R_{j,oos}^2\left(\mathcal{M}_s\right)>0$.}
The out-of-sample $R_{j,oos}^2\left(\mathcal{M}_s\right)$ tend to be mostly negative across estimation methods and shrinkage priors. This is consistent with the existing evidence on stock returns predictability: a simple naive forecast based on a rolling sample mean represents a challenging benchmark to beat (see, e.g., campbell2007). However, our variational inference approach substantially improves upon univariate regressions, as well as upon the {\tt LMCMC} and {\tt LVB} methods, which are both based on a structural VAR representation.
For instance, our {\tt VB} with stochastic volatility generates a positive $R_{j,oos}^2\left(\mathcal{M}_s\right)$ for more than half of the 30 industry portfolios based on the adaptive Normal-Gamma and the Horseshoe. This compares to 4 (adaptive Normal-Gamma) and 3 (Horseshoe) positive $R_{j,oos}^2\left(\mathcal{M}_s\right)$ obtained from {\tt LMCMC} with stochastic volatility. The gap further increases within the 49-industry classification; our {\tt VB} method is virtually the only approach that can systematically generate positive $R_{j,oos}^2\left(\mathcal{M}_s\right)$ across industries. Although concentrated on the Horseshoe prior, the out-performance of our method relative to both {\tt LMCMC} and {\tt VB} holds across different priors.
\paragraph{Density forecasts.} {\color{black}We follow fisher2020optimal and assess the accuracy of the density forecasts across priors and estimation methods based on the average log-score (ALS) differential with respect to a “no-predictability” benchmark,
where $\ln{S_{jt}}\left(\mathcal{M}_s\right)$ denotes the log-score at time $t$ for industry $j$ obtained by evaluating a Normal density with the conditional mean and variance forecast from the model $\mathcal{M}_s$. Consistent with the rationale of $R_{j,oos}^2\left(\mathcal{M}_s\right)$, the log-score for the no-predictability benchmark $\ln{\overline{S}_{j,t}}$ is constructed by evaluating a Normal density based on recursive mean and variance.}
Figure (ref) reports the results. The labeling is the same as in Figure (ref). \textcolor{black}{Not surprisingly, we find that by adding stochastic volatility the accuracy of density forecasts substantially improves across priors and estimation methods.} For instance, our {\tt VB} method with stochastic volatility generate positive log-score differentials for almost all of the portfolios for the 30 industry classification and for more than half of the 49 industry portfolios. Interestingly, when it comes to density forecasts rather than modeling expected returns, the gefang2023forecasting variational method built on a structural VAR representation performs on par with our {\tt VB} method. This is likely due to stochastic volatility alone, since our {\tt VB} still stands out within the constant volatility specifications. More generally, our {\tt VB} approach outperforms the competing estimation methods under all prior specifications.
\paragraph{Returns predictability over the business cycle.}
Existing literature suggests that expected returns are counter-cyclical and that returns predictability is more concentrated during period of economic contractions vs expansions (see, e.g., Goyal2010). Thus, we investigate if the forecasting performance of our modeling framework changes over the business cycle. More precisely, we split the data into recession and expansionary periods using the NBER dates of peaks and troughs. This information is considered {\it ex-post} and is not used at any time in the estimation and/or forecasting process. We compute the corresponding $R_{j,oos}^2\left(\mathcal{M}_s\right)$ for the recession periods only.
Figure (ref) reports the industries for which $R_{j,oos}^2\left(\mathcal{M}_s\right)>0$ for both the 30 (left panel) and the 49 (right panel) industry classification. The corresponding cross-sectional distribution of the $R_{j,oos}^2\left(\mathcal{M}_s\right)$ and the relative log-scores are reported in Appendix (ref). The labeling of Figure (ref) is the same as in Figure (ref). By comparing Figure (ref) with the results for the full sample, it suggests that the accuracy of the predictions substantially improves across methods and priors. Nevertheless, our {\tt VB} method outperforms the naive forecast from the rolling mean for a larger fraction of industry portfolios compared to other methods, in particular when stochastic volatility is considered. The difference between the recession and the full-sample performance persists when considering the 49 industry classification, especially for the adaptive Normal-Gamma and the Horseshoe prior.
A positive predictive performance does not necessarily translate into economic value. However, in practice an investor is obviously keenly interested in the economic value of returns predictability, perhaps even more than the statistical performance. Hence, it is of paramount importance to evaluate the extent to which apparent gains in predictive accuracy translates into better investment performances.
Following existing literature (see, e.g., Goyal2008,Goyal2010), we consider a representative investor with a single-period horizon and mean-variance preferences who allocates her wealth between an industry portfolio and a risk-free asset. {\color{black}Thus, the investor optimal allocation to stocks for period $t+1$ based on information at time $t$ is given by $w_{jt} = \frac{1}{\gamma}\frac{\widehat{y}_{jt}}{\widehat{\nu}^{-1}_{jt}}$, where $\widehat{y}_{jt}$ represents the returns conditional mean forecast for industry $j=1,\ldots,d$ and $\widehat{\nu}^{-1}_{jt}$ the corresponding volatility forecast at time $t$. We also constraint the weights for each of the industry to $-0.5 \le w_{jt} \le 1.5$ to prevent extreme short-sales and leverage positions. We assume a risk aversion coefficient of $\gamma=5$ (see, e.g., Dangl:Halling:2012).}
Figure (ref) reports the average utility gain -- in monthly % -- obtained by using a given forecast $\widehat{y}_{jt}$ instead of the recursive sample mean $\overline{y}_{jt}$. {\color{black}The average utility for a given model is calculated as $\widehat{u}_{j} = \overline{r}_{j} - 0.5\gamma\overline{\sigma}_{j}^2$ where $\overline{r}_{j}$ and $\overline{\sigma}_j^2$ represent the sample mean and variance, respectively, of the portfolio return $r_{jt+1} = w_{jt}y_{jt+1}$ realized over the forecasting period for the industry $j=1,\ldots,d$ under a given prior specification and estimation method. The utility gain is calculated by subtracting the average utility of a given model $\widehat{u}_{j}$ to the average utility obtained by using the naive forecast from the recursive mean and variance to calculate $w_{jt}$. A positive value for the utility gain indicates the fee that a risk-averse investor is willing to pay to access the investment strategy implied by $\mathcal{M}_s$.}
The economic value of each forecast largely confirms the same evidence offered by the out-of-sample statistical performance. From a pure economic standpoint, the forecast from a recursive mean are quite challenging to beat: we observe that the average utility gain is mostly negative, with the only exception of those provided by {\tt VB} under an Horseshoe prior specification. Economically, the results show that a representative investor with mean-variance utility is willing to pay, on average, a monthly fee of almost 15 basis points monthly to access the strategy based on our variational inference with stochastic volatility. In addition, the right panels of Figure (ref) show that the positive economic value obtained from our {\tt VB} is more broadly spread across industries compared to alternative methods. This holds especially for the 30 industry classification, but also applies to the more granular 49 industry classification.
\textcolor{black}{We propose a novel variational inference method for large Bayesian vector autoregressions (VAR) with exogenous predictors and stochastic volatility. Differently from most existing estimation methods for high-dimensional VAR models, our approach does not rely on a structural form representation. This allows a fast and accurate identification of the regression coefficients without leveraging on a standard Cholesky-based transformation of the parameter space. We show both in simulation and empirically that our estimation approach outperforms across different prior specifications, both statistically and economically, forecasts from existing benchmark estimation strategies, such as equivalent, non-linear MCMC algorithms (see, e.g., gruber2022forecasting) linearized MCMC (see, e.g., cross2020macroeconomic) and linearized variational inference methods (see, e.g., gefang2023forecasting).}
\vskip50pt \singlespacing {.05cm }
\onehalfspacing