Large Bayesian Tensor VARs with Stochastic Volatility

Large Bayesian Tensor VARs with Stochastic Volatility

\onehalfspace

\thispagestyle{empty}

Introduction

Large Bayesian vector autoregressions (BVARs) are now commonly used in macroeconomic forecasting and structural analysis, following early influential papers by \citet*{BGR10}, CKM09 and koop13. The dominant approach to tackle the challenge of parameter proliferation in large systems is to use shrinkage priors to regularize the variations in VAR coefficients; there is now an extensive literature on various shrinkage priors designed for BVARs.\footnote{The most widely-used shrinkage priors for BVARs are the family of Minnesota priors developed in a series of papers by \citet*{DLS84}, litterman86 and KK93, KK97. Recent additions to this family include \citet*{GLP15} and chan22. Another popular family is the adaptive hierarchical shrinkage priors that can be represented as scale mixtures of normals. Examples include the Bayesian Lasso PC08,KP19, the normal-gamma prior GB10,HF19, the horseshoe prior CPS10horseshoe,FY19 and the Dirichlet-Laplace prior DLP15,KH20. Naturally, one can also combine these adaptive hierarchical priors with the Minnesota priors, as proposed in chan21. HHK24 provide an excellent recent review on the state-of-the-art shrinkage priors developed for macroeconomic forecasting using BVARs.} In contrast, other dimension reduction techniques are relatively unexplored.

We investigate the usefulness of specifying a low-rank structure on the VAR coefficients for forecasting. More specifically, we follow the approach proposed by WZLL22 to treat the VAR coefficients collectively as a three-dimensional array or tensor: for a BVAR with $n$ endogenous variables and $p$ lags, we arrange the $n\times n$ coefficient matrices $\mathbf{A}_1, \ldots, \mathbf{A}_p$ in the third dimension to construct the third-order tensor $\mathcal{A}\in\mathbb{R}^{n\times n\times p}$. We then model $\mathcal{A}$ using a rank-$R$ CP decomposition. We call these BVARs constructed via the CP decomposition tensor VARs or TVARs.

Using this tensor decomposition, the number of free parameters is reduced from $n^2p$ to $(2n + p)R$. Since the number of free parameters under this tensor decomposition grows linearly in $n$, it is especially suitable for applications with a large number of variables. This approach is related to the reduced-rank VAR CKM11, which may be viewed as a special case in which the rank of $\mathcal{A}$ is reduced along one of the three possible dimensions.

Departing from the homoskedastic framework in WZLL22, we formulate the TVARs using a general stochastic volatility specification, which can represent a wide variety of multivariate stochastic volatility and COVID-19 outlier-augmented models commonly used for BVARs. This more general setup is motivated by the increasing recognition of the importance of allowing time-varying volatility for forecasting macroeconomic and financial variables, especially after the extreme economic turbulence triggered by the COVID-19 pandemic LP22,CCMM22.

We develop efficient estimation procedures for these TVARs with stochastic volatility. In particular, we propose two types of algorithms to sample the components of the CP decomposition of the VAR coefficient tensor $\mathcal{A}$: we can either sample each block of the components jointly as a matrix or sample each column of the component matrix separately. The latter is inspired by the equation-by-equation estimation approach designed for BVARs proposed in CCM19 and CCCM22, which can drastically reduce the computational burden when the number of endogenous variables $n$ is very large.

We illustrate the methodology using a forecasting exercise that involves 40 US quarterly variables, such as GDP, industrial production, labor market variables and a variety of inflation and interest rates. We consider various TVARs with different stochastic volatility specifications, and compare them to a standard BVAR with the Minnesota prior in an out-of-sample forecasting exercise. The results show that TVARs clearly outperform the standard BVAR, highlighting the usefulness of the low-rank specification for the VAR coefficient tensor. We also find that models with some form of time-varying volatility substantially forecast better than their homoskedastic counterparts. Interesting, the parsimonious common stochastic volatility model of CCM16 tends to outperform the more flexible stochastic volatility model of CS05, suggesting that it is foremost important to capture the strong comovements in the macroeconomic volatilities.

Our paper is closely related to the recent work by LG23, who also consider a BVAR with a tensor decomposition. While they focus on the Cholesky stochastic volatility model of CS05, we consider a more general setup that can represent a wide range of multivariate stochastic volatility models. Our paper is also related to the emerging literature on modeling multidimensional tensors, where the coefficient tensors are typically parameterized using CP or Tucker decompositions; see, for example, BCIK23 and WZL24.

The rest of this paper is organized as follows. Section (ref) first introduces a general framework of BVARs with a generic time-varying error covariance matrix and discusses how it can be used to represent a variety of popular stochastic volatility and COVID-19 outlier-augmented models. We then outline the proposed approach of arranging the VAR coefficient matrices as a third-order tensor in Section (ref). Next, we introduce the efficient sampling algorithms in Section (ref). Section (ref) considers a recursive out-of-sample forecasting application that involves 40 US macroeconomic variables. Lastly, Section (ref) concludes and briefly discusses some future research directions.

Large Bayesian VARs with Stochastic Volatility

Let $\mathbf{y}_t$ be an $n\times 1$ vector of endogenous variables at time $t$ for $t=1,\ldots, T.$ Consider the following VAR($p$): \[ \mathbf{y}_t = \mathbf{A}_1 \mathbf{y}_{t-1} + \cdots + \mathbf{A}_p\mathbf{y}_{t-p} + \mathbf{u}_t, \] where $\mathbf{A}_1, \ldots, \mathbf{A}_p$ are $n\times n$ coefficient matrices. We omit an intercept term for ease of exposition; an intercept or any exogenous variables can be added to the model with minor modifications. Let $\mathbf{A} = (\mathbf{A}_1, \ldots, \mathbf{A}_p)'$ be the $np\times n$ matrix of VAR coefficients and let $\mathbf{x}_t = (\mathbf{y}_{t-1}',\ldots, \mathbf{y}_{t-p}')'$ denote a vector of lag variables of dimension $np\times 1$. Then, stacking the observations over $t=1,\ldots, T$, we can rewrite the VAR more succinctly as

where the matrices $\mathbf{Y}$, $\mathbf{X}$ and $\mathbf{U}$ are, respectively, $T\times n$, $T\times np$ and $T\times n$. In a standard homoskedastic VAR, the reduced-form errors $\mathbf{u}_1,\ldots, \mathbf{u}_T$ are assumed to be independent and identically distributed (iid) as $\mathcal{N}(\mathbf{0}_n,\boldsymbol \Sigma)$, where $\mathbf{0}_n$ is an $n\times 1$ vector of zeros and $\boldsymbol \Sigma$ is an $n\times n$ covariance matrix. However, it is increasingly recognized that some form of time-varying volatility is needed in modeling typical macroeconomic time-series. Early influential papers such as CS05, Primiceri05 and SZ06 have highlighted the secular variations in volatility. There is now a large empirical literature that demonstrates the importance of allowing time-varying volatility in improving model-fit and forecasting performance in the context of Bayesian VARs; examples include clark11, \citet*{KK13}, \citet*{DGG13}, CR15, CP16 and CE18.

Below we outline a few stochastic volatility models suitable for large BVARs. In particular, the innovations of the VAR in (ref) are now distributed as

where $\boldsymbol \Sigma_t$ is a generic time-varying covariance matrix.

One of the first stochastic volatility models designed for large BVARs is the common stochastic volatility proposed in \citet*{CCM16}. Their model is motivated by the empirical observation that the estimated time-varying error variances of many macroeconomic variables have broadly similar low-frequency movements. A parsimonious way to model these comovements is to introduce a time-varying latent factor to scale the error covariance matrix via

where $\boldsymbol \Omega$ is a time-invariant covariance matrix. The log-volatility $h_t$ is modeled using a zero-mean stationary AR(1) process:

for $t=2,\ldots, T$, where $|\phi|<1$ and the process is initialized as $h_{1}\sim\mathcal{N}(0,\sigma_h^2/(1-\phi^2))$. The common stochastic volatility model specified in (ref)--(ref) can be extended to incorporate other useful features. For example, chan20 introduces a general framework that can accommodate heavy-tailed, heteroskedastic and serially dependent innovations.

Another widely-used stochastic volatility specification for BVARs is the Cholesky stochastic volatility model---based on the modified Cholesky decomposition of the covariance matrix---developed in CS05. More specifically, consider the decomposition of $\boldsymbol \Sigma_t $ via

where $\mathbf{B}_{0}$ is an $n \times n$ lower triangular matrix with ones on the diagonal and $\mathbf{D}_t = \text{diag}(\text{e}^{h_{1,t}}, \ldots, \text{e}^{h_{n,t}})$. Each element of the vector $\mathbf{h}_t = (h_{1,t}, \ldots, h_{n,t})'$ is modeled via an independent autoregressive process:

for $t=2,\ldots, T$, where $h_{i,1}$ is initialized as $h_{i,1}\sim\mathcal{N}(\mu_i,\sigma_i^2/(1-\phi_i^2))$ for $i=1,\ldots, n.$

This Cholesky stochastic volatility specification is more flexible than the common stochastic volatility model in (ref), since the former contains $n$ stochastic volatility processes and can accommodate more complex covolatility patterns. This flexibility, however, comes at a cost of higher model complexity. Whether this stochastic volatility specification forecasts better than alternatives in large systems is an empirical question.

Another difference between the Cholesky stochastic volatility and the common stochastic volatility is that the latter is order-invariant---i.e., parameter estimates are invariant to reordering the endogenous variables in $\mathbf{y}_t$---whereas the former is not. One practical implication is that forecasts from the Cholesky stochastic volatility model could differ substantially across different variable orderings, as documented in ARRS23 using a similar model of Primiceri05.

The reason why the Cholesky stochastic volatility is not order-invariant is partly due to the use of the lower triangular parameterization of $\mathbf{B}_0$ in (ref). Motivated by this simple observation, CKY24 extend the model by relaxing this lower triangular assumption and instead specify $\mathbf{B}_0$ to be any non-degenerate square matrix. They prove that the model is order invariant. Moreover, based on the results in BB22, $\mathbf{B}_0$ is also identified up to permutations and sign switches.

The extreme movements in many macroeconomic variables at the onset of the COVID-19 pandemic have motivated much recent work on modeling outliers in macroeconomic time-series. An example is the outlier-augmented stochastic volatility model proposed by CCMM22, which builds on the Cholesky stochastic volatility model and the discrete mixture representation for the innovations introduced in SW16. This outlier-augmented model can also be represented using (ref).

There are many other multivariate stochastic volatility models for $\boldsymbol \Sigma_t$. CM23 and chan24 provide two recent reviews on BVARs with a wide range of stochastic volatility and outlier-augmented specifications.

Tensor Decomposition of VAR Coefficients

For high-dimensional settings, there are two key challenges related to the proliferation of the VAR coefficients. First, the number of VAR coefficients increases quadratically in $n$. In many large-scale applications, there are far more VAR coefficients than the number of observations, which makes it necessary to regularize these VAR coefficients. Second, due to the proliferation of VAR coefficients, sampling them tends to be very computational intensive.

These two related challenges are typically tackled separately in the literature. For instance, Bayesian shrinkage priors are widely used to regularize the VAR coefficients. These include the family of Minnesota priors \citep*{DLS84, litterman86, KK93, KK97, GLP15,chan22} and various adaptive hierarchical shrinkage priors HF19,KP19,KH20,chan21. The computational challenge is addressed by developing efficient MCMC or variational methods to sample the large number of VAR coefficients CCM19,CCCM22,GKP23,BBB24. We instead take an alternative approach that tackles these two challenges simultaneously by imposing a low-rank structure on the VAR coefficients.

One possibility is the reduced-rank VAR CKM11, in which the matrix $\mathbf{A}' = (\mathbf{A}_1, \ldots, \mathbf{A}_p)$ is assumed to have a reduced rank $R<n$. That is, the dimension of the column space of the VAR coefficient matrices $\mathbf{A}_1, \ldots, \mathbf{A}_p$ has rank $R$. Alternatively, one could impose a low-rank structure on $(\mathbf{A}_1', \ldots, \mathbf{A}_p')$ or $(\text{vec}(\mathbf{A}_1), \ldots, \text{vec}(\mathbf{A}_p))$, whose ranks correspond to the dimensions of the row space and vectorized matrix space of the VAR coefficient matrices, respectively. Each of these options essentially reduces the dimension along one of the three different directions.

This motivates us to follow the approach in WZLL22 to treat the VAR coefficients collectively as a three-dimensional array or tensor. That is, we arrange the $n\times n$ coefficient matrices $\mathbf{A}_1, \ldots, \mathbf{A}_p$ in the third dimension to form the third-order tensor $\mathcal{A}\in\mathbb{R}^{n\times n\times p}$. With this tensor representation, we then use a rank-$R$ CP decomposition to construct $\mathcal{A}$:

where $\boldsymbol \theta_1^{(r)}, \boldsymbol \theta_2^{(r)}\in\mathbb{R}^{n}$ and $\boldsymbol \theta_3^{(r)} \in\mathbb{R}^{p}$. Each component $\boldsymbol \theta_1^{(r)}\circ \boldsymbol \theta_2^{(r)} \circ \boldsymbol \theta_3^{(r)}$ is a rank-1 third-order tensor whose $(i,j,k)$ element is the product of the $i$-th, $j$-th and $k$-th elements of $\boldsymbol \theta_1^{(r)}, \boldsymbol \theta_2^{(r)}$ and $\boldsymbol \theta_3^{(r)}$, respectively. Using this tensor decomposition with a small $R$, the number of parameters is reduced from $n^2p$ to $(2n+p)R$. We refer the readers to TB09 for a general introduction to tensors and their operations.

Bayesian Estimation

In this section we describe the priors and outline the posterior simulator. In particular, we focus on the sampling of the components $\boldsymbol \theta_1^{(r)}, \boldsymbol \theta_2^{(r)}$ and $\boldsymbol \theta_3^{(r)}, r=1,\ldots, R$. To that end, let $\boldsymbol \Theta_1 = (\boldsymbol \theta_1^{(1)}, \ldots, \boldsymbol \theta_1^{(R)})$ and similarly define $\boldsymbol \Theta_2$ and $\boldsymbol \Theta_3$. The dimensions of $\boldsymbol \Theta_1, \boldsymbol \Theta_2$ and $\boldsymbol \Theta_3$ are, respectively, $n\times R$, $n\times R$ and $p\times R$. For later reference, stack $\boldsymbol \theta_1 = \text{vec}(\boldsymbol \Theta_1)$, $\boldsymbol \theta_2 = \text{vec}(\boldsymbol \Theta_2')$ and $\boldsymbol \theta_3 = \text{vec}(\boldsymbol \Theta_3')$. For reasons that will become transparent later, note that $\boldsymbol \theta_1$ is constructed by stacking the columns of $\boldsymbol \Theta_1$, whereas $\boldsymbol \theta_2$ and $\boldsymbol \theta_3$ are formed by stacking the rows of $\boldsymbol \Theta_2$ and $\boldsymbol \Theta_3$, respectively.

Next, consider the following independent Gaussian priors on $\boldsymbol \theta_1, \boldsymbol \theta_2$ and $\boldsymbol \theta_3$:

Here we consider simple Gaussian priors, but any adaptive hierarchical shrinkage priors that have a conditionally Gaussian representation---such as the normal-gamma prior, the horseshoe prior or the Dirichlet-Laplace prior---can be used.

In what follows, we derive the full conditional posterior distributions of $(\boldsymbol \theta_1 \,|\, \mathbf{Y}, \boldsymbol \theta_{2}, \boldsymbol \theta_{3}, \boldsymbol \Sigma)$, $(\boldsymbol \theta_2 \,|\, \mathbf{Y}, \boldsymbol \theta_{1}, \boldsymbol \theta_3, \boldsymbol \Sigma)$ and $(\boldsymbol \theta_3 \,|\, \mathbf{Y}, \boldsymbol \theta_{1}, \boldsymbol \theta_2, \boldsymbol \Sigma)$, where $\mathbf{\boldsymbol \Sigma}=\text{diag}\left(\boldsymbol \Sigma_1, \ldots, \boldsymbol \Sigma_T\right)$.

Sampling of $\boldsymbol \theta_1$

For sampling $\boldsymbol \theta_1$, first let

Then, the mode-1 matricization of $\mathcal{A}$ can be written as TB09:

Furthermore, let $\mathcal{X}\in \mathbb{R}^{T\times n \times p}$ denote the third-order tensor constructed by stacking the $n\times p$ matrices $(\mathbf{y}_{t-1},\ldots, \mathbf{y}_{t-p}), t=1,\ldots, T,$ along the first dimension so that $\text{vec}(\mathcal{X}_{t,:,:}) = \mathbf{x}_t$. It is easy to verify that its mode-1 matricization is the $T\times np$ matrix $\mathbf{X}$, i.e., $ \mathcal{X}_{(1)}=\mathbf{X}$. Hence, combining (ref) and (ref), we have

It follows that \[ \text{vec}(\mathbf{Y}')=(\mathcal{X}_{(1)}\boldsymbol \Theta_{-1} \otimes \mathbf{I}_{n})\text{vec}(\boldsymbol \Theta_1)+\text{vec}(\mathbf{U}'), \] where $\text{vec}(\mathbf{U}') \sim \mathcal{N}(\mathbf{0}_{Tn}, \boldsymbol \Sigma)$ with $\mathbf{\boldsymbol \Sigma}=\text{diag}\left(\boldsymbol \Sigma_1, \ldots, \boldsymbol \Sigma_T\right)$. Given the Gaussian prior on $\boldsymbol \theta_1 = \text{vec}(\boldsymbol \Theta_1)$ specified in (ref), the full conditional posterior of $\boldsymbol \theta_1$ is given by \[ (\boldsymbol \theta_1 \,|\, \mathbf{Y}, \boldsymbol \theta_{2}, \boldsymbol \theta_{3}, \boldsymbol \Sigma) \sim \mathcal{N}\left(\widehat{\boldsymbol \theta}_1, \mathbf{K}_{\boldsymbol \theta_1}^{-1}\right), \] where

Since $\boldsymbol \theta_1$ is of length $nR$, sampling $\boldsymbol \theta_1$ generally involves $\mathcal{O}(n^3R^3)$ elementary operations. When both $n$ and $R$ are large, this sampling step could be computationally intensive. An alternative is to sample each $\boldsymbol \theta_1^{(r)}$ separately, $r=1,\ldots, R$. For a fixed $r$, this can be done by substituting \[ \boldsymbol \Theta_1\boldsymbol \Theta_{-1}' = \sum_{s=1}^R \boldsymbol \theta_1^{(s)}\boldsymbol \theta_{-1}^{(s)\prime}, \] where $\boldsymbol \theta_{-1}^{(s)} = \boldsymbol \theta_3^{(s)}\otimes \boldsymbol \theta_2^{(s)}$, into the transpose of (ref) to obtain \[ \mathbf{Y}' = \sum_{s=1}^R\boldsymbol \theta^{(s)}_{1}\boldsymbol \theta_{-1}^{(s)\prime}\mathcal{X}_{(1)}'+ \mathbf{U}'. \] Next, let $\mathbf{Y}_{1r} = \mathbf{Y}'-\sum_{s \neq r}\boldsymbol \theta_{1}^{(s)}\boldsymbol \theta_{-1}^{(s)\prime}\mathcal{X}_{(1)}'$ and vectorize the above equation, we have \[ \text{vec}(\mathbf{Y}_{1r})= (\mathcal{X}_{(1)}\boldsymbol \theta^{(r)}_{-1} \otimes \mathbf{I}_{n})\boldsymbol \theta_{1}^{(r)}+\text{vec}(\mathbf{U}'). \] If the marginal prior for $\boldsymbol \theta_1^{(r)}$ is \[ \boldsymbol \theta_1^{(r)} \sim \mathcal{N}\left(\boldsymbol \theta_{1,0}^{(r)}, \mathbf{V}_{\boldsymbol \theta_1^{(r)}}\right), \] the posterior distribution $(\boldsymbol \theta_1^{(r)} \,|\, \mathbf{Y}, \{\boldsymbol \theta_{1}^{(s)}\}_{s\neq r}, \boldsymbol \theta_{2}, \boldsymbol \theta_{3}, \boldsymbol \Sigma)$ has the form \[ (\boldsymbol \theta_1^{(r)} \,|\, \mathbf{Y}, \{\boldsymbol \theta_{1}^{(s)}\}_{s\neq r}, \boldsymbol \theta_{2}, \boldsymbol \theta_{3}, \boldsymbol \Sigma ) \sim \mathcal{N}\left(\widehat{\boldsymbol \theta}_1^{(r)}, \mathbf{K}_{\boldsymbol \theta_1^{(r)}}^{-1}\right), \] where

The drawback of sampling each $\boldsymbol \theta_1^{(r)}$ separately is that this tends to increase the autocorrelation of the constructed Markov chain. But this approach is substantially faster and remains computationally feasible even when both $n$ and $R$ are large.

Sampling of $\boldsymbol \theta_2$ and $\boldsymbol \theta_3$

Next, we derive the conditional distribution of $(\boldsymbol \theta_2 \,|\, \mathbf{Y}, \boldsymbol \theta_{1}, \boldsymbol \theta_3, \boldsymbol \Sigma)$ and show that it is Gaussian. To start, we aim to write the VAR in (ref) as a linear regression in $\boldsymbol \theta_2=\text{vec}(\boldsymbol \Theta_2')$.

We first introduce some notations. Let $\mathbf{e}_i^p$ be the $i$-th column of $\mathbf{I}_p$. Let $\mathbf{P}$ denote the $np \times np$ commutation matrix so that $\mathbf{P}'\text{vec}(\mathbf{Z}) = \text{vec}(\mathbf{Z}')$ for any $n\times p$ matrix $\mathbf{Z}$. Explicitly, $\mathbf{P}$ can be constructed by setting the $(k,l)$ element to be 1, i.e., $\mathbf{P}_{k,l} = 1$, if there exist $i$ and $j$ such that $k = (i-1)n+j$ and $l=(j-1)p + i$; otherwise, set $\mathbf{P}_{k,l}=0$.

The proof of this proposition is given in Appendix A. Now, given the representation in (ref) and the Gaussian prior on $\boldsymbol \theta_2$ specified in (ref), by standard linear regression results, one can verify that the full conditional posterior of $\boldsymbol \theta_2$ is given by \[ (\boldsymbol \theta_2 \,|\, \mathbf{Y}, \boldsymbol \theta_{1}, \boldsymbol \theta_3, \boldsymbol \Sigma) \sim \mathcal{N}\left(\widehat{\boldsymbol \theta}_2, \mathbf{K}_{\boldsymbol \theta_2}^{-1}\right), \] where

Therefore, one can sample $\boldsymbol \theta_2$ in one block. When $n$ and $R$ are large, it might only be feasible to sample each $\boldsymbol \theta_2^{(r)}$ at a time for $r=1,\ldots, R$. We provide the details of this alternative approach in Appendix C.

Likewise, the conditional distribution of $(\boldsymbol \theta_3 \,|\, \mathbf{Y}, \boldsymbol \theta_{1}, \boldsymbol \theta_2, \boldsymbol \Sigma)$ can be shown to be Gaussian. More specifically, we first write the VAR in (ref) as a linear regression in $\boldsymbol \theta_3=\text{vec}(\boldsymbol \Theta_3')$. To that end, let $\mathbf{e}_i^n$ denote the $i$-th column of $\mathbf{I}_n$, let $\mathcal{X}_{(3)}$ represent the mode-3 matricization of $\mathcal{X}$, and define $\mathbf{Q}$ to be the $n^2 \times n^2$ commutation matrix so that $\mathbf{Q}'\text{vec}(\mathbf{Z}) = \text{vec}(\mathbf{Z}')$ for any $n\times n$ matrix $\mathbf{Z}$. Explicitly, $\mathbf{Q}$ can be constructed by setting the $\mathbf{Q}_{k,l} = 1$ if there exist $i$ and $j$ such that $k = (i-1)n+j$ and $l = (j-1)n + i$; otherwise, set $\mathbf{Q}_{k,l} = 0$.

Given the representation in (ref), it is easy to verify that \[ (\boldsymbol \theta_3 \,|\, \mathbf{Y}, \boldsymbol \theta_{1}, \boldsymbol \theta_2, \boldsymbol \Sigma) \sim \mathcal{N}\left(\widehat{\boldsymbol \theta}_3, \mathbf{K}_{\boldsymbol \theta_3}^{-1}\right), \] where

Alternatively, one can sample each $\boldsymbol \theta_3^{(r)}$ at a time for $r=1,\ldots, R$ when $n$ and $R$ are large as before.

Finally, sampling the time-varying error covariance matrices $\boldsymbol \Sigma_1,\ldots, \boldsymbol \Sigma_T$ naturally depends on the stochastic volatility specification used. For a wide variety of stochastic volatility specifications commonly employed in applied work, such as the common stochastic volatility and Cholesky stochastic volatility discussed earlier, efficient algorithms are available to sample the latent variables and model parameters. We refer the readers to chan23JE for more details.

An Empirical Application

We conduct an out-of-sample forecasting exercise to evaluate the performance of the proposed tensor VARs compared to a standard benchmark. More specifically, we construct a dataset of 40 quarterly macroeconomic and financial variables from the FRED-QD database MN20. The sample spans from 1969Q1 to 2024Q1. The data include key macroeconomic variables such as GDP, inflation rates, labor market variables and various interest rates. We refer the readers to Appendix B for the detailed description of the time series and their transformations.

We consider three tensor VARs: a homoskedastic tensor VAR (TVAR), TVARs with the common volatility (TVAR-CSV) and the Cholesky stochastic volatility (TVAR-SV). For now, we set the rank of the VAR coefficient tensor $\mathcal{A}$ to be $R=1$ for all TVARs. As a benchmark, we also include a standard Bayesian VAR (BVAR) with the Minnesota prior (implemented as the natural conjugate prior). The evaluation period of the forecasting exercise begins in 2010Q1 and ends in 2024Q1.

To assess the performance of jointly forecasting all $n=40$ time series, we calculate the average log predictive likelihoods for each model over one- and four-quarter-ahead forecast horizons. The results are presented in Table (ref). Higher values of log predictive likelihoods signify better forecast performance.

A few observations can be drawn from these forecasting results. Firstly, comparing the two homoskedastic models, TVAR and BVAR, it is clear that the former substantially outperforms the latter, suggesting that the specification of a low rank structure on the VAR coefficient tensor $\mathcal{A}$ is more appropriate. Secondly, allowing some form of time-varying volatility clearly improves forecast performance. For example, extending a homoskedastic TVAR to a version with the common stochastic volatility increases the average log predictive likelihood from $-8$ to about $-3$ for one-step-ahead density forecasts. This finding is line with the large body of empirical evidence that demonstrates the importance of allowing time-varying volatility in macroeconomic forecasting clark11, DGG13, chan23JE. Finally, among the two TVARs with stochastic volatility, the version with the parsimonious common stochastic volatility outperforms the one with the more flexible Cholesky stochastic volatility, indicating strong comovements in the macroeconomic volatilities.

Next, we look at the point forecast performance of the TVARs for individual time series. In particular, Table (ref) reports the root mean squared forecast errors (RMSFEs) of a few key macroeconomic variables over the evaluation period relative to the benchmark BVAR. Values less than one indicate better forecast performance than the benchmark.

The results show that the TVARs tend to outperform the benchmark BVAR for both one- and four-quarter-ahead horizons. For example, the RMSFE of the homoskedastic TVAR for forecasting GDP is only 65% of that of the benchmark. Interestingly, allowing time-varying volatility does not appear to improve point forecasts. In fact, the more flexible TVAR-SV often provides inferior point forecasts compared to the homoskedastic TVAR (though the performance of the more parsimonious TVAR-CSV is virtually identical to TVAR).

So far we have set $R$, the rank of the VAR coefficient tensor $\mathcal{A}$, to be 1 for all TVARs. Naturally, one might wonder how this choice affects the density and point forecast performance. To investigate the impact of the choice of $R$, we present in Table (ref) the one-quarter-ahead average log predictive likelihoods of the three TVARs for $R=1,3,5,10$.

It is interesting to note that the choice of $R$ does not seem to have a large impact on the joint density forecast performance across the three TVARs. In particular, it does not affect the relative ranking of the three models; the two TVARs with time-varying volatility outperform the homoskedastic version regardless of the rank $R$. Overall, low-rank TVARs tend to work as well as, if not better than, TVARs with $R=10$.

Table (ref) reports the one-quarter-ahead RMSFEs of the three TVARs for forecasting a few key macroeconomic variables. To facilitate comparison, the results for each TVAR are relative to the corresponding TVAR with $R=1$. With the exception of CPI inflation, increasing the rank $R$ does not appear to substantially improve point forecast performance.

Overall, these forecasting results demonstrate the benefits of the proposed approach of specifying a low rank structure on the VAR coefficient tensor. In addition, our results highlight the importance of accommodating time-varying volatility in forecasting macroeconomic variables. Notably, in our forecasting exercise, the more parsimonious TVAR-CSV tends to forecast better than the more flexible TVAR-SV.

Concluding Remarks and Future Research

We have developed Bayesian tensor VARs in which the VAR coefficients are arranged as a third-order tensor and parameterized using a rank-$R$ CP decomposition. We then introduced efficient sampling algorithms to simulate the components of the tensor decomposition. Through a forecasting exercise, we showed that these TVARs outperform the standard BVAR with the Minnesota prior.

For future research, it would be useful to extend these TVARs to allow for time-varying coefficients in the mean equations. This tensor framework is especially suitable for developing time variation in the VAR coefficients, as the number of free parameters grows only linear in $n$. Of course, for large systems one might need additional shrinkage. In those cases, one can consider either the static shrinkage approach in chan23 or the dynamic shrinkage approaches proposed in KK18 and HKO19.