# Comparing Stochastic Volatility Specifications for Large Bayesian VARs

Joshua C.C. Chan
Purdue University
I would like to thank Roberto Casarin, Todd Clark and William McCausland for
their insightful comments and constructive suggestions. This paper has also
benefited from many helpful discussions with conference and seminar
participants at the 11th European Central Bank Conference on Forecasting
Techniques, Ca’ Foscari University of Venice, the University of Montreal and
the Federal Reserve Bank of Kansas City. All remaining errors are, of course,
my own.

(This version: August 2022
First version: March 2021 )

###### Abstract

Large Bayesian vector autoregressions with various forms of stochastic
volatility have become increasingly popular in empirical macroeconomics. One
main difficulty for practitioners is to choose the most suitable stochastic
volatility specification for their particular application. We develop Bayesian
model comparison methods—based on marginal likelihood estimators that combine
conditional Monte Carlo and adaptive importance sampling—to choose among a
variety of stochastic volatility specifications. The proposed methods can also
be used to select an appropriate shrinkage prior on the VAR coefficients,
which is a critical component for avoiding over-fitting in high-dimensional
settings. Using US quarterly data of different dimensions, we find that both
the Cholesky stochastic volatility and factor stochastic volatility outperform
the common stochastic volatility specification. Their superior performance,
however, can mostly be attributed to the more flexible priors that accommodate
cross-variable shrinkage.

Keywords: large vector autoregression, marginal likelihood, Bayesian model
comparison, stochastic volatility, shrinkage prior

##  1 Introduction

Large Bayesian vector autoregressions (VARs) are now widely used for empirical
macroeconomic analysis and forecasting thanks to the seminal work of Banbura,
Giannone, and Reichlin (2010).111The literature on using large Bayesian VARs
for structural analysis and forecasting is rapidly expanding. Early
applications include Carriero, Kapetanios, and Marcellino (2009), Koop (2013),
Koop and Korobilis (2013), Banbura, Giannone, Modugno, and Reichlin (2013) and
Carriero, Clark, and Marcellino (2015). Since it is well established that
time-varying volatility is vitally important for small VARs, it is expected to
be even more so for large systems.222For small VARs with only a few variables,
papers such as Clark (2011), D’Agostino, Gambetti, and Giannone (2013), Clark
and Ravazzolo (2014), Cross and Poon (2016) and Chan and Eisenstat (2018)
highlight the importance of time-varying volatility for model-fit and
forecasting. For large VARs with more than a dozen variables, Carriero, Clark,
and Marcellino (2016) present evidence that the data overwhelmingly favors a
VAR with stochastic volatility; Koop and Korobilis (2013), Carriero, Clark,
and Marcellino (2019) and Chan, Eisenstat, Hou, and Koop (2020) show that VARs
with stochastic volatility forecast better than their homoskedastic
counterparts, especially for density forecasts. Consequently, there has been a
lot of recent research devoted to designing stochastic volatility
specifications suitable for large systems. Prominent examples include the
common stochastic volatility models (Carriero, Clark, and Marcellino, 2016;
Chan, 2020), the Cholesky stochastic volatility models (Cogley and Sargent,
2005; Carriero, Clark, and Marcellino, 2019) and the factor stochastic
volatility models (Pitt and Shephard, 1999b; Chib, Nardari, and Shephard,
2006; Kastner, 2019). Since these stochastic volatility models are widely
different and the choice among these alternatives involves important trade-
offs—e.g., flexibility versus speed of estimation—one main issue facing
practitioners is the lack of tools to compare these high-dimensional, non-
linear and non-nested models.

Of course, the natural Bayesian model comparison criterion is the marginal
likelihood (or the marginal data density), and in principle it can be used to
select among these stochastic volatility models. In practice, however,
computing the marginal likelihood for high-dimensional VARs with stochastic
volatility is hardly trivial due to the large number of VAR coefficients and
the latent state variables (e.g., stochastic volatility and latent factors).
We tackle this obstacle by developing new methods to estimate the marginal
likelihood of large Bayesian VARs with a variety of stochastic volatility
specifications.

More specifically, we combine two popular variance reduction techniques,
namely, conditional Monte Carlo and adaptive importance sampling, to construct
our marginal likelihood estimators. We first analytically integrate out the
large number of VAR coefficients—i.e., we derive an analytical expression of
the likelihood unconditional on the VAR coefficients that can be evaluated
quickly. In the second step, we construct an adaptive importance sampling
estimator—obtained by minimizing the Kullback-Leibler divergence to the ideal
zero-variance importance sampling density—to integrate out the log-volatility
via Monte Carlo. By carefully combining these two ways of integration
(analytical and Monte Carlo integration), we are able to efficiently evaluate
the marginal likelihood of a variety of popular stochastic volatility models
for large Bayesian VARs.

Compared to earlier marginal likelihood estimators for Bayesian VARs with
stochastic volatility, such as Chan and Eisenstat (2018) and Chan (2020), the
new method offers two main advantages. First, it analytically integrates out
the large number of VAR coefficients. As such, it reduces the variance of the
estimator by eliminating the portion contributed by the VAR coefficients. This
reduction is expected to be substantial in large VARs. Second, earlier
marginal likelihood estimators are based on _local_ approximations of the
joint distribution of the log-volatility, such as a second-order Taylor
expansion of the log target density around the mode. Although these
approximations are guaranteed to approximate the target density well around
the neighborhood of the point of expansion, their accuracy typically
deteriorates rapidly away from the approximation point. In contrast, the new
method is based on a _global_ approximation that incorporates information from
the entire support of the target distribution. This is done by solving an
optimization problem to locate the closest density to the target posterior
distribution—measured by the Kullback-Leibler divergence—within a class of
multivariate Gaussian distributions.

In addition to comparing different stochastic volatility specifications, the
proposed method can also be used to select an appropriate shrinkage prior on
the VAR coefficients. Since even small VARs have a large number of parameters,
shrinkage priors are essential to avoid over-fitting in high-dimensional
settings. The most prominent example of these shrinkage priors is the
Minnesota prior first introduced by Doan, Litterman, and Sims (1984) and
Litterman (1986), not long after the seminal work on VARs by Sims (1980).
There are now a wide range of more flexible variants (see, e.g., Kadiyala and
Karlsson, 1993, 1997; Giannone, Lenza, and Primiceri, 2015; Chan, 2021), and
choosing among them for a particular application has become a practical issue
for empirical economists. In particular, we focus on Minnesota priors with two
potentially useful features: 1) allowing the overall shrinkage hyperparameters
to be estimated from the data rather than fixing them at some subjectively
elicited values; and 2) cross-variable shrinkage, i.e., the idea that
coefficients on ‘other’ lags should be shrunk to zero more aggressively than
those on ‘own’ lags. The proposed marginal likelihood estimators can provide a
way to compare shrinkage priors with and without these features.

Through a series of Monte Carlo experiments, we demonstrate that the proposed
estimators work well in practice. In particular, we show that one can
correctly distinguish the three different stochastic volatility
specifications: the common stochastic volatility (VAR-CSV), the Cholesky
stochastic volatility (VAR-SV) and the factor stochastic volatility (VAR-FSV).
In addition, the proposed marginal likelihood estimators can also be used to
identify the correct number of factors in factor stochastic volatility models.

In an empirical application using US quarterly data, we compare the three
stochastic volatility specifications in fitting datasets of different
dimensions (7, 15 and 30 variables). The model comparison results show that
the data overwhelmingly prefers VAR-SV and VAR-FSV over the more restrictive
VAR-CSV for all model dimensions. We also find strong evidence in favor of the
two aforementioned features of some Minnesota priors: both cross-variable
shrinkage and a data-based approach to determine the overall shrinkage
strength are empirically important. In fact, when we turn off the cross-
variable shrinkage in the priors of VAR-SV and VAR-FSV, they perform similarly
as VAR-CSV, suggesting that the superior performance of the former two models
can mostly be attributed to the more flexible priors that accommodate cross-
variable shrinkage. These results thus illustrate that in high-dimensional
settings, choosing a flexible shrinkage prior is as important as selecting a
flexible stochastic volatility specification. Lastly, we also demonstrate how
the proposed methodology can be used to improve other marginal likelihood
estimators and model selection criteria such as the deviance information
criterion.

The rest of the paper is organized as follows. We first outline in Section 2
the three stochastic volatility specifications designed for large BayBunesian
VARs, followed by an overview of various data-driven Minnesota priors. Then,
Section 3 describes the two components of the proposed marginal likelihood
estimators: adaptive importance sampling and conditional Monte Carlo. The
methodology is illustrated via a concrete example of estimating the marginal
likelihood of the common stochastic volatility model. We then conduct a series
of Monte Carlo experiments in Section 4 to assess how the proposed marginal
likelihood estimators perform in selecting the correct data generating
process. It is followed by an empirical application in Section 5, where we
compare the three stochastic volatility specifications in the context of
Bayesian VARs of different dimensions. Section 6 discusses various extensions
of using the proposed methods to improve other marginal likelihood estimators
and alternative model comparison methods. There we also consider a stochastic
volatility model with an outlier component. Lastly, Section 7 concludes and
briefly discusses some future research directions.

##  2 Stochastic Volatility Models for Large VARs

In this section we first describe a few recently proposed stochastic
volatility specifications designed for large Bayesian VARs. We then outline a
few data-driven Minnesota priors that are particularly useful for these high-
dimensional models.

###  2.1 Common Stochastic Volatility

Let $\mathbf{y}_{t}=(y_{1,t},\ldots,y_{n,t})^{\prime}$ be an $n\times 1$
vector of variables that is observed over the periods $t=1,\ldots,T.$ The
first specification is the common stochastic volatility model introduced in
Carriero, Clark, and Marcellino (2016). The conditional mean equation is the
standard reduced-form VAR with $p$ lags:

$\mathbf{y}_{t}=\mathbf{a}_{0}+\mathbf{A}_{1}\mathbf{y}_{t-1}+\cdots+\mathbf{A}_{p}\mathbf{y}_{t-p}+\boldsymbol{\varepsilon}_{t},$
(1)

where $\mathbf{a}_{0}$ is an $n\times 1$ vector of intercepts and
$\mathbf{A}_{1},\ldots,\mathbf{A}_{p}$ are all $n\times n$ coefficient
matrices. To allow for heteroskedastic errors, the covariance matrix of the
innovation $\boldsymbol{\varepsilon}_{t}$ is scaled by a common, time-varying
factor that can be interpreted as the overall macroeconomic volatility:

$\boldsymbol{\varepsilon}_{t}\sim\mathcal{N}(\mathbf{0},\text{e}^{h_{t}}\boldsymbol{\Sigma}).$
(2)

The log-volatility $h_{t}$ in turn follows a stationary AR(1) process:

$h_{t}=\phi h_{t-1}+u_{t}^{h},\quad u_{t}^{h}\sim\mathcal{N}(0,\sigma^{2}),$
(3)

for $t=2,\ldots,T$, where $|\phi|<1$ and the initial condition is specified as
$h_{1}\sim\mathcal{N}(0,\sigma^{2}/(1-\phi^{2}))$. Note that the unconditional
mean of the AR(1) process is assumed to be zero for identification. We refer
to this common stochastic volatility model as VAR-CSV.

One drawback of this volatility specification is that it appears to be
restrictive. For example, all variances are scaled by a single factor and,
consequently, they are always proportional to each other. Nevertheless, there
is empirical evidence that the error variances of macroeconomic variables tend
to move closely together (see, e.g., Carriero, Clark, and Marcellino, 2016;
Chan, Eisenstat, and Strachan, 2020), and a common stochastic volatility is a
parsimonious way to model this empirical feature.

Estimating large VARs is in general computationally intensive because of the
large number of VAR coefficients
$\mathbf{A}=(\mathbf{a}_{0},\mathbf{A}_{1},\ldots,\mathbf{A}_{p})^{\prime}$.
One main advantage of the common stochastic volatility specification is
that—if a natural conjugate prior on $(\mathbf{A},\boldsymbol{\Sigma})$ is
used—it leads to many useful analytical results that make estimation fast. In
particular, there are efficient algorithms to generate the large number of VAR
coefficients. Moreover, as demonstrated in Chan (2020), similar computational
gains can be achieved for a much wider class of VARs with non-Gaussian,
heteroskedastic and serially dependent innovations. Recent empirical
applications using this common stochastic volatility and its variants include
Mumtaz (2016), Mumtaz and Theodoridis (2017), Götz and Hauzenberger (2018),
Poon (2018), Louzis (2019), LeSage and Hendrikz (2019), Zhang and Nguyen
(2020) and Hartwig (2021).

###  2.2 Cholesky Stochastic Volatility

A more flexible way to model multivariate heteroscedasticity and time-varying
covariances is to incorporate multiple stochastic volatility processes, as
first proposed in Cogley and Sargent (2005). Specifically, consider the same
reduced-form VAR in (1), but the innovation $\boldsymbol{\varepsilon}_{t}$ has
a more flexible time-varying covariance matrix modeled as follows:

$\boldsymbol{\varepsilon}_{t}\sim\mathcal{N}(\mathbf{0},\boldsymbol{\Sigma}_{t}),\quad\boldsymbol{\Sigma}_{t}^{-1}=\mathbf{B}_{0}^{\prime}\mathbf{D}_{t}^{-1}\mathbf{B}_{0},$
(4)

where $\mathbf{B}_{0}$ is an $n\times n$ lower triangular matrix with ones on
the diagonal and
$\mathbf{D}_{t}=\operatorname{diag}(\text{e}^{h_{1,t}},\ldots,\text{e}^{h_{n,t}})$.
The law of motion for each element of
$\mathbf{h}_{t}=(h_{1,t},\ldots,h_{n,t})^{\prime}$ is specified as an
independent autoregressive process:

$h_{i,t}=\mu_{i}+\phi_{i}(h_{i,t-1}-\mu_{i})+u_{i,t}^{h},\quad
u_{i,t}^{h}\sim\mathcal{N}(0,\sigma_{i}^{2})$ (5)

for $t=2,\ldots,T$, where the initial condition is specified as
$h_{i,1}\sim\mathcal{N}(\mu_{i},\sigma_{i}^{2}/(1-\phi_{i}^{2}))$. This
stochastic volatility specification is sometimes called the Cholesky
stochastic volatility. We refer to this stochastic volatility model as VAR-SV.

This model can be estimated using the algorithm described in Del Negro and
Primiceri (2015). However, when the number of variables is large, the
conventional way of drawing all VAR coefficients jointly becomes excessively
computationally intensive. To tackle this computational problem, Carriero,
Clark, and Marcellino (2019) introduce a blocking scheme that makes it
possible to estimate the reduced-form VAR equation by equation that
drastically reduces the computational time.333As pointed out in Bognanni
(2022), the algorithm in Carriero, Clark, and Marcellino (2019) to sample the
VAR coefficients equation by equation can only be viewed as an approximation,
as it ignores certain integrating constants of the distributions of VAR
coefficients in other equations. Carriero, Chan, Clark, and Marcellino (2022)
provide an alternative algorithm to draw the VAR coefficients equation by
equation that has the same computational complexity of $\mathcal{O}(n^{4})$.

In contrast to the common stochastic volatility model, VAR-SV is more flexible
in that it contains $n$ stochastic volatility processes, which can accommodate
more complex co-volatility patterns. But this comes at a cost of more
intensive posterior computations: the complexity of estimating VAR-SV is
$\mathcal{O}(n^{4})$ compared to $\mathcal{O}(n^{3})$ for VAR-CSV (when a
natural conjugate prior is used). Recent empirical applications using this
Cholesky stochastic volatility in the context of large Bayesian VARs include
Banbura and van Vlodrop (2018), Bianchi, Guidolin, and Ravazzolo (2018), Huber
and Feldkircher (2019), Cross, Hou, and Poon (2019), Baumeister, Korobilis,
and Lee (2020), Koop, McIntyre, Mitchell, and Poon (2020), Tallman and Zaman
(2020), Zens, Böck, and Zörner (2020) and Chan (2021).

###  2.3 Factor Stochastic Volatility

The third stochastic volatility specification that is suitable for large
systems belongs to the class of factor stochastic volatility models (Pitt and
Shephard, 1999a; Chib, Nardari, and Shephard, 2006; Kastner, 2019). More
specifically, consider the same reduced-form VAR in (1), but the innovation is
instead decomposed as:

$\boldsymbol{\varepsilon}_{t}=\mathbf{L}\mathbf{f}_{t}+\mathbf{u}_{t},$

where $\mathbf{f}_{t}=(f_{1,t},\ldots,f_{r,t})^{\prime}$ is a $r\times 1$
vector of latent factors and $\mathbf{L}$ is the associated $n\times r$ factor
loading matrix. For identification purposes we assume that $\mathbf{L}$ is a
lower triangular matrix with ones on the main diagonal.444One main reason to
impose these identification restrictions on $\mathbf{L}$ is to avoid the so-
called label switching problem, which makes it a bit more difficult to
construct appropriate importance sampling densities for the log-volatility.
One way to handle this issue is to first postprocess the posterior draws to
sort them into the correct categories using, e.g., the approach in Kaufmann
and Schumacher (2019). Then, one can use these sorted MCMC draws to construct
the importance sampling densities discussed below. Furthermore, to ensure one
can separately identify the common and the idiosyncratic components, we adopt
a sufficient condition in Anderson and Rubin (1956) that requires
$r\leqslant(n-1)/2$.555In this paper we follow the common approach of fixing
the number of factors $r$ in the estimation and selecting $r$ via the marginal
likelihood. An alternative approach is to allow the inclusion of potentially
infinitely many factors, but to use a shrinkage prior on the factor loadings
so that they are increasingly shrunk to 0 as the column index increases
(Bhattacharya and Dunson, 2011; Kastner, 2019). This latter approach
circumvents the need to select the number of factors, but it makes the
interpretation of the factors more difficult.

The disturbances $\mathbf{u}_{t}$ and the latent factors $\mathbf{f}_{t}$ are
assumed to be independent at all leads and lags. Furthermore, they are jointly
Gaussian:

$\begin{pmatrix}\mathbf{u}_{t}\\\
\mathbf{f}_{t}\end{pmatrix}\sim\mathcal{N}\left(\begin{pmatrix}\mathbf{0}\\\
\mathbf{0}\end{pmatrix},\begin{pmatrix}\mathbf{D}_{t}&\mathbf{0}\\\
\mathbf{0}&\mathbf{G}_{t}\end{pmatrix}\right),$ (6)

where
$\mathbf{D}_{t}=\text{diag}(\text{e}^{h_{1,t}},\ldots,\text{e}^{h_{n,t}})$ and
$\mathbf{G}_{t}=\text{diag}(\text{e}^{h_{n+1,t}},\ldots,\text{e}^{h_{n+r,t}})$
are diagonal matrices. Here the correlations among the elements of the
innovation $\boldsymbol{\varepsilon}_{t}$ are induced by the latent factors.
In typical applications, a small number of factors would be sufficient to
capture the time-varying covariance structure even when $n$ is large.

Next, for each $i=1,\ldots,n+r$, the evolution of the log-volatility is
modeled as:

$h_{i,t}=\mu_{i}+\phi_{i}(h_{i,t}-\mu_{i})+u_{i,t}^{h},\quad
u_{i,t}^{h}\sim\mathcal{N}(0,\sigma_{i}^{2})$ (7)

for $t=2,\ldots,T$. The initial state is assumed to follow the stationary
distribution:
$h_{i,1}\sim\mathcal{N}(\mu_{i},\sigma_{i}^{2}/(1-\phi_{i}^{2}))$. We refer to
this factor stochastic volatility model as VAR-FSV.

In a sense VAR-FSV is even more flexible than VAR-SV as the former contains
$n+r$ stochastic volatility processes compared to $n$ in the latter. In terms
of estimation, both are more computationally intensive to estimate relative to
VAR-CSV, but they can still be fitted reasonably quickly even when $n$ is
large. In particular, given the latent factors VAR-FSV becomes $n$ unrelated
regressions. Consequently, the $n$-equation system can be estimated equation
by equation. While various factor stochastic volatility specifications are
widely applied in financial applications (see, e.g., Aguilar and West, 2000;
Jin, Maheu, and Yang, 2019; Hautsch and Voigt, 2019; Li and Scharth, 2020;
McCausland, Miller, and Pelletier, 2020), they are not yet commonly used in
conjunction with a large VAR (with the notable exception of Kastner and Huber,
2020).

###  2.4 Hierarchical Minnesota Priors

Next, we briefly describe the prior specifications for the above three
stochastic volatility models; the details of the priors are given in Appendix
A. In general, we assume exactly the same priors on the common parameters
across models. When this is not possible, we use similar priors so that they
are comparable across models. Below we focus on the shrinkage priors on the
VAR coefficients as they are critical in the current context.

More specifically, we consider Minnesota-type priors on the VAR coefficients.
This family of priors was first developed by Doan, Litterman, and Sims (1984)
and Litterman (1986), and a number of more flexible variants have been
developed since. We adopt a version in which the overall shrinkage
hyperparameters are treated as parameters to be estimated from the data, in
the spirit of Giannone, Lenza, and Primiceri (2015). Compared to the
conventional version in which the hyperparameters are fixed at some subjective
values, this data-based variant is often found to fit the data substantially
better and have better out-of-sample forecast performance (see, e.g.,
Carriero, Clark, and Marcellino, 2015; Chan, Jacobi, and Zhu, 2020).

First, for the VAR-CSV we consider the natural conjugate prior on
$(\mathbf{A},\boldsymbol{\Sigma})$ that depends on a shrinkage hyperparameter
$\kappa$:

$\boldsymbol{\Sigma}\sim\mathcal{IW}(\nu_{0},\mathbf{S}_{0}),\quad(\text{vec}(\mathbf{A})\,|\,\boldsymbol{\Sigma},\kappa)\sim\mathcal{N}(\text{vec}(\mathbf{A}_{0}),\boldsymbol{\Sigma}\otimes\mathbf{V}_{\mathbf{A}}),$

where the prior hyperparameters $\text{vec}(\mathbf{A}_{0})$ and
$\mathbf{V}_{\mathbf{A}}$ are chosen in the spirit of the Minnesota prior.
More specifically, since in the empirical applications we will be using data
in growth rates, we set $\text{vec}(\mathbf{A}_{0})=\mathbf{0}$ to shrink all
the coefficients to 0. For the covariance matrix $\mathbf{V}_{\mathbf{A}}$, it
is assumed to be diagonal and it depends on a single hyperparameter $\kappa$
that controls the overall shrinkage strength (see Appendix A for details). As
mentioned above, we treat $\kappa$ as an unknown parameter to be estimated.
This natural conjugate prior substantially speeds up posterior simulations due
to the availability of various analytical results. One major drawback of this
natural conjugate prior, however, is that it cannot accommodate cross-variable
shrinkage—i.e., the prior belief that coefficients on ‘own’ lags are on
average larger than ‘other’ lags.

Next, for the VAR-SV, we assume that the intercept and the VAR coefficients
are a priori independent across equations, and each vector
$\boldsymbol{\alpha}_{i}=(a_{i,0},\mathbf{A}_{i,1},\ldots,\mathbf{A}_{i,p})^{\prime}$,
for $i=1,\ldots,n$, where $\mathbf{A}_{i,j}$ denotes the $i$-th row of
$\mathbf{A}_{j}$, has a hierarchical normal prior:
$\boldsymbol{\alpha}_{i}\sim\mathcal{N}(\boldsymbol{\alpha}_{0,i},\mathbf{V}_{\boldsymbol{\alpha}_{i}})$.
Here the prior mean $\boldsymbol{\alpha}_{0,i}$ is set to be zero to shrink
the VAR coefficients to zero. The prior covariance matrix
$\mathbf{V}_{\boldsymbol{\alpha}_{i}}$ is assumed to be diagonal and it
depends on two hyperparameters: $\kappa_{1}$ and $\kappa_{2}$. The
hyperparameter $\kappa_{1}$ controls the overall shrinkage strength for VAR
coefficients on own lags, whereas $\kappa_{2}$ controls those on lags of other
variables. Allowing the two shrinkage parameters to be different across the
two groups of VAR coefficients is empirically important. For example,
Carriero, Clark, and Marcellino (2015) and Chan (2021) find empirical evidence
in support of cross-variable shrinkage. In addition, following recent works
such as Giannone, Lenza, and Primiceri (2015) and Amir-Ahmadi, Matthes, and
Wang (2020) that demonstrate the benefits of tuning these hyperparameters in a
data-based manner relative to fixing them at specific values, we treat both
$\kappa_{1}$ and $\kappa_{2}$ to be parameters to be estimated.

Let $\boldsymbol{\beta}_{i}$ denote the free elements in the $i$-th row of the
impact matrix $\mathbf{B}_{0}$ for $i=2,\ldots,n$. That is,
$\boldsymbol{\beta}_{i}=(B_{0,i,1},\ldots,B_{0,i,i-1})^{\prime}$, where
$B_{0,i,j}$ is the $(i,j)$ element of $\mathbf{B}_{0}$. We assume that
$\boldsymbol{\beta}_{i}$ has a hierarchical normal prior:
$\boldsymbol{\beta}_{i}\sim\mathcal{N}(\boldsymbol{\beta}_{0,i},\mathbf{V}_{\boldsymbol{\beta}_{i}})$.
We set the the prior mean $\boldsymbol{\beta}_{0,i}$ to be zero so that the
free elements in the impact matrix are shrunk to zero. The prior covariance
matrix $\mathbf{V}_{\boldsymbol{\beta}_{i}}$ is set to be diagonal, and it
depends on a hyperparameter $\kappa_{3}$, which controls the overall shrinkage
strength. We again treat $\kappa_{3}$ as an unknown parameter to be estimated.
Allowing this shrinkage parameter to be determined by the data is found to
substantially improve forecasts (see, e.g., Chan, 2021).

Finally, for the VAR-FSV we assume a similar normal prior on
$\boldsymbol{\alpha}_{i}=(a_{i,0},\mathbf{A}_{i,1},\ldots,\mathbf{A}_{i,p})^{\prime}$,
the intercept and the VAR coefficients in the $i$-th equation:
$\boldsymbol{\alpha}_{i}\sim\mathcal{N}(\boldsymbol{\alpha}_{0,i},\mathbf{V}_{\boldsymbol{\alpha}_{i}})$.
Again we set $\boldsymbol{\alpha}_{0,i}=\mathbf{0}$ and assume
$\mathbf{V}_{\boldsymbol{\alpha}_{i}}$ to be diagonal and depends on two
hyperparameters: one controls the overall shrinkage strength on coefficients
on own lags and the other on coefficients on other lags. Details of these
Minnesota priors and priors on other parameters are given in Appendix A.

##  3 Marginal Likelihood Estimation: Adaptive Importance Sampling and
Conditional Monte Carlo

In this section we first give an overview of the marginal likelihood and its
role in Bayesian model comparison. We then describe the two components of the
proposed marginal likelihood estimators: adaptive importance sampling and
conditional Monte Carlo. Specifically, conditional Monte Carlo is used to
integrate out the large number of VAR coefficients analytically, and adaptive
importance sampling is used to bias the distributions of the remaining
parameters to oversample the ‘important’ regions of the parameter space that
contribute most to the integral (and this bias is later corrected in the
estimation). Finally, to illustrate the methodology, we provide the details of
estimating the marginal likelihood of the VAR-CSV by combining adaptive
importance sampling and conditional Monte Carlo.

###  3.1 Overview of the Marginal Likelihood

Suppose we want to compare a collection of models $\\{M_{1},\ldots,M_{K}\\}$,
where each model $M_{k}$ is formally defined by a likelihood function
$p(\mathbf{y}\,|\,\boldsymbol{\theta}_{k},M_{k})$ and a prior on the model-
specific parameter vector $\boldsymbol{\theta}_{k}$ denoted by
$p(\boldsymbol{\theta}_{k}\,|\,M_{k})$. A natural Bayesian model comparison
criterion is the Bayes factor in favor of $M_{i}$ against $M_{j}$, defined as

$\text{BF}_{ij}=\frac{p(\mathbf{y}\,|\,M_{i})}{p(\mathbf{y}\,|\,M_{j})},$

where

$p(\mathbf{y}\,|\,M_{k})=\int
p(\mathbf{y}\,|\,\boldsymbol{\theta}_{k},M_{k})p(\boldsymbol{\theta}_{k}\,|\,M_{k})\text{d}\boldsymbol{\theta}_{k}$

is the _marginal likelihood_ under model $M_{k}$, $k=i,j.$ One can show that
if the data is generated from model $M_{i}$, the Bayes factor would on average
pick the correct model over some distinct model, say, $M_{j}$. To see that, we
compute the expected log Bayes factor in favor of model $M_{i}$ with respect
to the distribution $p(\mathbf{y}\,|\,M_{i})$:

$\mathbb{E}_{\mathbf{y}}\left[\log\frac{p(\mathbf{y}\,|\,M_{i})}{p(\mathbf{y}\,|\,M_{j})}\right]=\int
p(\mathbf{y}\,|\,M_{i})\,\log\frac{p(\mathbf{y}\,|\,M_{i})}{p(\mathbf{y}\,|\,M_{j})}\text{d}\mathbf{y}.$

This expression is the Kullback-Leibler divergence from
$p(\mathbf{y}\,|\,M_{i})$ to $p(\mathbf{y}\,|\,M_{j})$, and it is strictly
positive unless $p(\mathbf{y}\,|\,M_{i})=p(\mathbf{y}\,|\,M_{j})$, in which
case it is zero. In other words, the marginal likelihood of $M_{i}$ is larger
than that of $M_{j}$ on average.

Furthermore, the Bayes factor is related to the posterior odds ratio between
the two models as follows:

$\frac{\mathbb{P}(M_{i}\,|\,\mathbf{y})}{\mathbb{P}(M_{j}\,|\,\mathbf{y})}=\frac{\mathbb{P}(M_{i})}{\mathbb{P}(M_{j})}\times\text{BF}_{ij},$

where $\mathbb{P}(M_{i})/\mathbb{P}(M_{j})$ is the prior odds ratio. It
follows that if both models are equally probable a priori, i.e.,
$p(M_{i})=p(M_{j})$, the posterior odds ratio between the two models is then
equal to the Bayes factor or ratio of marginal likelihoods. In that case, if,
say, $\text{BF}_{ij}=100$, then model $M_{i}$ is 100 times more likely than
model $M_{j}$ given the data. For a more detailed discussion of the Bayes
factor and its role in Bayesian model comparison, see Koop (2003) or Chan,
Koop, Poirier, and Tobias (2019). From here onwards we suppress the model
indicator; for example we denote the marginal likelihood and the likelihood by
$p(\mathbf{y})$ and $p(\mathbf{y}\,|\,\boldsymbol{\theta})$, respectively.

###  3.2 Adaptive Importance Sampling

Next, we give an overview of an adaptive importance sampling approach called
the improved cross-entropy method for estimating the marginal likelihood. It
is based on the idea of biasing the sampling distribution in such a way that
more ‘important values’ are generated in the simulation. The sample is then
weighted to correct for the use of a different distribution to give an
unbiased estimator. In our context of estimating the marginal likelihood
$p(\mathbf{y})$, consider the following importance sampling estimator:

$\widehat{p(\mathbf{y})}_{\rm
IS}=\frac{1}{R}\sum_{r=1}^{R}\frac{p(\mathbf{y}\,|\,\boldsymbol{\theta}^{(r)})p(\boldsymbol{\theta}^{(r)})}{g(\boldsymbol{\theta}^{(r)})},$
(8)

where $\boldsymbol{\theta}^{(1)},\ldots,\boldsymbol{\theta}^{(R)}$ are
independent draws obtained from the importance sampling density $g(\cdot)$
that dominates $p(\mathbf{y}\,|\,\cdot)p(\cdot)$—i.e.,
$g(\boldsymbol{\theta})=0\Rightarrow
p(\mathbf{y}\,|\,\boldsymbol{\theta})p(\boldsymbol{\theta})=0$.

The importance sampling estimator in (8) is unbiased and simulation consistent
for any $g$ that dominates $p(\mathbf{y}\,|\,\cdot)p(\cdot)$. But its
performance depends critically on the choice of $g$. Below we follow Chan and
Eisenstat (2015) to use an improved version of the classic cross-entropy
method to construct $g$ optimally. The original cross-entropy method was
developed for rare-event simulation by Rubinstein (1997, 1999) using a multi-
level procedure to obtain the optimal importance sampling density (see also
Rubinstein and Kroese, 2004, for a book-length treatment). Chan and Eisenstat
(2015) show that this optimal importance sampling density can be obtained more
accurately in one step using MCMC methods. This adaptive importance sampling
estimator is then used to estimate the marginal likelihood.

To motivate the cross-entropy method, first note that there exists an
importance sampling density that gives a zero-variance estimator of the
marginal likelihood. In particular, it is easy to verify that if we use the
posterior distribution as the importance sampling density, i.e.,
$g^{*}(\boldsymbol{\theta})=p(\boldsymbol{\theta}\,|\,\mathbf{y})=p(\mathbf{y}\,|\,\boldsymbol{\theta})p(\boldsymbol{\theta})/p(\mathbf{y})$,
then the associated importance sampling estimator (8) has zero variance. Of
course in practice $g^{*}$ cannot be used as its normalizing constant is
exactly the marginal likelihood, the unknown quantity we wish to estimate.
Nevertheless, this provides a signpost to obtain an optimal importance
sampling density. Intuitively, if we choose an importance sampling density $g$
that is ‘close enough’ to $g^{*}$ so that both behave similarly, the resulting
importance sampling estimator should have reasonable accuracy. Hence, our goal
is to locate a convenient density that is in a well-defined sense ‘close’ to
$g^{*}$.

To that end, consider a parametric family
$\mathcal{F}=\\{f(\boldsymbol{\theta};\mathbf{v})\\}$ indexed by the parameter
vector $\mathbf{v}$ within which we locate the optimal importance sampling
density. We find the density
$f(\boldsymbol{\theta};\mathbf{v}^{*})\in\mathcal{F}$ such that it is the
‘closest’ to $g^{*}$. One convenient measure of closeness between densities is
the _Kullback-Leibler divergence_ or the _cross-entropy distance_.
Specifically, let $f_{1}$ and $f_{2}$ be two probability density functions.
Then, the cross-entropy distance from $f_{1}$ to $f_{2}$ is defined as:

$\mathcal{D}(f_{1},f_{2})=\int
f_{1}(\mathbf{x})\log\frac{f_{1}(\mathbf{x})}{f_{2}(\mathbf{x})}\text{d}\mathbf{x}.$

Given this measure, we can then locate the density
$f(\cdot;\mathbf{v})\in\mathcal{F}$ such that
$\mathcal{D}(g^{*},f(\cdot;\mathbf{v}))$ is minimized:
$\mathbf{v}_{\text{ce}}^{*}=\mathop{\rm
argmin}_{\mathbf{v}}\mathcal{D}(g^{*},f(\cdot;\mathbf{v}))$. It can be shown
that solving the CE minimization problem is equivalent to finding

$\mathbf{v}^{*}_{\text{ce}}=\mathop{\rm argmax}_{\mathbf{v}}\int
p(\mathbf{y}\,|\,\boldsymbol{\theta})p(\boldsymbol{\theta})\log
f(\boldsymbol{\theta};\mathbf{v})\text{d}\boldsymbol{\theta}.$

This optimization problem is often difficult to solve analytically. Instead,
we consider its stochastic counterpart:

$\widehat{\mathbf{v}}^{*}_{\text{ce}}=\mathop{\rm
argmax}_{\mathbf{v}}\frac{1}{M}\sum_{m=1}^{M}\log
f(\boldsymbol{\theta}_{m};\mathbf{v}),$ (9)

where $\boldsymbol{\theta}^{(1)},\ldots,\boldsymbol{\theta}^{(M)}$ are draws
from the posterior distribution $p(\boldsymbol{\theta}\,|\,\mathbf{y})\propto
p(\mathbf{y}\,|\,\boldsymbol{\theta})p(\boldsymbol{\theta})$.

In other words, $\widehat{\mathbf{v}}^{*}_{\text{ce}}$ is exactly the maximum
likelihood estimate for $\mathbf{v}$ if we treat
$f(\boldsymbol{\theta};\mathbf{v})$ as the likelihood function with parameter
vector $\mathbf{v}$ and
$\boldsymbol{\theta}^{(1)},\ldots,\boldsymbol{\theta}^{(M)}$ an observed
sample. Since finding the maximum likelihood estimator is a standard problem,
solving (9) is typically easy. In particular, analytical solutions to (9) can
be found explicitly for the exponential family (e.g., Rubinstein and Kroese,
2004, p. 70).

The parametric family $\mathcal{F}$ is often chosen so that each member
$f(\boldsymbol{\theta};\mathbf{v})$ is a product of densities, e.g.,
$f(\boldsymbol{\theta};\mathbf{v})=f(\boldsymbol{\theta}_{1};\mathbf{v}_{1})\times\cdots\times
f(\boldsymbol{\theta}_{B};\mathbf{v}_{B})$, where
$\boldsymbol{\theta}=(\boldsymbol{\theta}_{1},\ldots,\boldsymbol{\theta}_{B})$
and $\mathbf{v}=(\mathbf{v}_{1},\ldots,\mathbf{v}_{B})$. In that case, one can
reduce the possibly high-dimensional maximization problem (9) into $B$ low-
dimensional problems, which can then be readily solved. Once the optimal
density is obtained, we then set
$g(\cdot)=f(\cdot;\widehat{\mathbf{v}}^{*}_{\text{ce}})$ and use it to
construct the importance sampling estimator in (8).

###  3.3 Conditional Monte Carlo

Recall that the marginal likelihood can be written as
$p(\mathbf{y})=\mathbb{E}p(\mathbf{y}\,|\,\boldsymbol{\theta})$, where
$\mathbf{y}$ is fixed and the expectation is taken with respect to the prior
distribution $p(\boldsymbol{\theta})$. Suppose that there is a random vector
$\boldsymbol{\psi}$ with density function $p(\boldsymbol{\psi})$ such that the
conditional expectation
$\mathbb{E}[p(\mathbf{y}\,|\,\boldsymbol{\theta})\,|\,\boldsymbol{\psi}]$ can
be computed analytically. By the law of iterated expectation, we have

$p(\mathbf{y})=\mathbb{E}[\mathbb{E}[p(\mathbf{y}\,|\,\boldsymbol{\theta})\,|\,\boldsymbol{\psi}]],$

where the outer expectation is taken with respect to $p(\boldsymbol{\psi})$
and the inner expectation is taken with respect to the conditional
distribution $p(\boldsymbol{\theta}\,|\,\boldsymbol{\psi})$. Hence,
$\mathbb{E}[p(\mathbf{y}\,|\,\boldsymbol{\theta})\,|\,\boldsymbol{\psi}]$ is
an unbiased estimator of $p(\mathbf{y})$. In addition, by the variance
decomposition formula, we obtain

$\text{Var}(p(\mathbf{y}\,|\,\boldsymbol{\theta}))=\mathbb{E}[\text{Var}(p(\mathbf{y}\,|\,\boldsymbol{\theta})\,|\,\boldsymbol{\psi})]+\text{Var}(\mathbb{E}[p(\mathbf{y}\,|\,\boldsymbol{\theta})\,|\,\boldsymbol{\psi}]).$

Since
$\mathbb{E}[\text{Var}(p(\mathbf{y}\,|\,\boldsymbol{\theta})\,|\,\boldsymbol{\psi})]>0$
if the density $p(\boldsymbol{\theta}\,|\,\boldsymbol{\psi})$ is non-
degenerate, we have

$\text{Var}(\mathbb{E}[p(\mathbf{y}\,|\,\boldsymbol{\theta})\,|\,\boldsymbol{\psi}])<\text{Var}(p(\mathbf{y}\,|\,\boldsymbol{\theta})).$

Therefore, the conditional Monte Carlo estimator
$\mathbb{E}[p(\mathbf{y}\,|\,\boldsymbol{\theta})\,|\,\boldsymbol{\psi}]$
_always_ provides variance reduction. The degree of variance reduction, of
course, depends on the dimension and variability of $\boldsymbol{\psi}$.
Loosely speaking, to maximize variance reduction, one feasible approach is to
find the ‘smallest’ subset $\boldsymbol{\psi}$ of $\boldsymbol{\theta}$ such
that $\mathbb{E}[p(\mathbf{y}\,|\,\boldsymbol{\theta})\,|\,\boldsymbol{\psi}]$
is computable.

In our VARs the parameter vector $\boldsymbol{\theta}$ contains the VAR
coefficients and other latent variables such as the stochastic volatility and
dynamic factors. To have the ‘smallest’ $\boldsymbol{\psi}$ is equivalent to
integrating out as many elements in $\boldsymbol{\theta}$ as possible. In
particular, since there are a large number of VAR coefficients in our setting,
ideally we would like to remove their contribution to the variability in the
marginal likelihood estimation. Fortunately, for all the VARs we consider, we
are able to integrate them out analytically and obtain an analytical
expression of
$\mathbb{E}[p(\mathbf{y}\,|\,\boldsymbol{\theta})\,|\,\boldsymbol{\psi}]$. We
will give an explicit example in the following section.

###  3.4 An Example: Marginal Likelihood of VAR-CSV

As an illustration, we now provide the details of estimating the marginal
likelihood of the VAR-CSV using the proposed method by combining adaptive
importance sampling and conditional Monte Carlo. Other VARs with stochastic
volatility can be handled similarly and the technical details are given in
Appendix A.

Let
$\mathbf{x}_{t}^{\prime}=(1,\mathbf{y}_{t-1}^{\prime},\ldots,\mathbf{y}_{t-p}^{\prime})$
be a $1\times k$ vector of an intercept and lags with $k=1+np$. Then, stacking
the observations over $t=1,\ldots,T$, we can rewrite (1)–(2) as:

$\mathbf{Y}=\mathbf{X}\mathbf{A}+\boldsymbol{\varepsilon},\quad\text{vec}(\boldsymbol{\varepsilon})\sim\mathcal{N}(\mathbf{0},\boldsymbol{\Sigma}\otimes\mathbf{D}_{\mathbf{h}}),$

where
$\mathbf{D}_{\mathbf{h}}=\text{diag}(\text{e}^{h_{1}},\ldots,\text{e}^{h_{T}})$,
$\mathbf{A}=(\mathbf{a}_{0},\mathbf{A}_{1},\ldots,\mathbf{A}_{p})^{\prime}$ is
$k\times n$, and the matrices $\mathbf{Y}$, $\mathbf{X}$ and
$\boldsymbol{\varepsilon}$ are respectively of dimensions $T\times n$,
$T\times k$ and $T\times n$. In addition, the natural conjugate prior on
$(\mathbf{A},\boldsymbol{\Sigma})$ depends on a single hyperparameter $\kappa$
that controls the overall shrinkage strength, which is treated as an unknown
parameter. Finally, the log-volatility $h_{t}$ follows the AR(1) process in
(3) with AR coefficient $\phi$ and error variance $\sigma^{2}$.

Next, we define an appropriate conditional Monte Carlo estimator, which will
then be combined with adaptive importance sampling. Using the notations in
Section 3.3, the set of model parameters and latent variables is
$\boldsymbol{\theta}=\\{\mathbf{A},\boldsymbol{\Sigma},\phi,\sigma^{2},\kappa,\mathbf{h}\\}$,
and the conditional likelihood (i.e., conditional on the latent variables
$\mathbf{h}$) is given by
$p(\mathbf{y}\,|\,\boldsymbol{\theta})=p(\mathbf{Y}\,|\,\mathbf{h},\mathbf{A},\boldsymbol{\Sigma},\phi,\sigma^{2},\kappa)=p(\mathbf{Y}\,|\,\mathbf{h},\mathbf{A},\boldsymbol{\Sigma})$.
When $n$ is large, the dimensions of both $\mathbf{A}$ and
$\boldsymbol{\Sigma}$ are large, and we expect them to contribute the most to
the variability of
$p(\mathbf{Y}\,|\,\mathbf{h},\mathbf{A},\boldsymbol{\Sigma})$. Hence, we aim
to integrate them out. To that end, we let
$\boldsymbol{\psi}=\\{\mathbf{h},\phi,\kappa\\}$ and consider the conditional
expectation
$\mathbb{E}[p(\mathbf{y}\,|\,\boldsymbol{\theta})\,|\,\boldsymbol{\psi}]=p(\mathbf{Y}\,|\,\mathbf{h},\kappa)$,
which can be computed analytically:

$\displaystyle p(\mathbf{Y}\,|\,\mathbf{h},\kappa)$ $\displaystyle=\int
p(\mathbf{Y}\,|\,\mathbf{h},\mathbf{A},\boldsymbol{\Sigma})p(\mathbf{A},\boldsymbol{\Sigma}\,|\,\kappa)\text{d}(\mathbf{A},\boldsymbol{\Sigma})$
$\displaystyle=\pi^{-\frac{Tn}{2}}\text{e}^{-\frac{n}{2}\mathbf{1}_{T}^{\prime}\mathbf{h}}|\mathbf{V}_{\mathbf{A}}|^{-\frac{n}{2}}|\mathbf{K}_{\mathbf{A}}|^{-\frac{n}{2}}\frac{\Gamma_{n}\left(\frac{\nu_{0}+T}{2}\right)|\mathbf{S}_{0}|^{\frac{\nu_{0}}{2}}}{\Gamma_{n}\left(\frac{\nu_{0}}{2}\right)|\widehat{\mathbf{S}}|^{\frac{\nu_{0}+T}{2}}},$
(10)

where $\Gamma_{n}(\cdot)$ is the multivariate gamma function and

$\displaystyle\mathbf{K}_{\mathbf{A}}$
$\displaystyle=\mathbf{V}_{\mathbf{A}}^{-1}+\mathbf{X}^{\prime}\mathbf{D}^{-1}_{\mathbf{h}}\mathbf{X},\quad\widehat{\mathbf{A}}=\mathbf{K}_{\mathbf{A}}^{-1}(\mathbf{V}_{\mathbf{A}}^{-1}\mathbf{A}_{0}+\mathbf{X}^{\prime}\mathbf{D}^{-1}_{\mathbf{h}}\mathbf{Y}),$
$\displaystyle\widehat{\mathbf{S}}$
$\displaystyle=\mathbf{S}_{0}+\mathbf{A}_{0}^{\prime}\mathbf{V}_{\mathbf{A}}^{-1}\mathbf{A}_{0}+\mathbf{Y}^{\prime}\mathbf{D}^{-1}_{\mathbf{h}}\mathbf{Y}-\widehat{\mathbf{A}}^{\prime}\mathbf{K}_{\mathbf{A}}\widehat{\mathbf{A}}.$

Note that the overall shrinkage parameter $\kappa$ appears in the prior
covariance matrix $\mathbf{V}_{\mathbf{A}}$. The details of the derivations
are given in Appendix A. We note that (10) can be calculated quickly without
inverting any large matrices; see, e.g., Chan (2020) for various computational
shortcuts.

Now, the marginal likelihood can be written in terms of
$p(\mathbf{Y}\,|\,\mathbf{h},\kappa)$ as follows:

$p(\mathbf{y})=\int
p(\mathbf{Y}\,|\,\mathbf{h},\kappa)p(\mathbf{h}\,|\,\phi,\sigma^{2})p(\phi,\sigma^{2},\kappa)\text{d}(\mathbf{h},\phi,\sigma^{2},\kappa),$

where $p(\mathbf{h}\,|\,\phi,\sigma^{2})$ is a Gaussian density implied by the
state equation (3) and $p(\phi,\sigma^{2},\kappa)$ is the prior for
$(\phi,\sigma^{2},\kappa)$. Hence, the corresponding conditional Monte Carlo
estimator is given by

$\widehat{p(\mathbf{y})}_{\rm
CMC}=\frac{1}{M}\sum_{m=1}^{M}p(\mathbf{Y}\,|\,\mathbf{h}^{(m)},\kappa^{(m)}),$

where $\mathbf{h}^{(1)},\ldots,\mathbf{h}^{(M)}$ are drawn from the prior
$p(\mathbf{h})=\int
p(\mathbf{h}\,|\,\phi,\sigma^{2})p(\phi,\sigma^{2})\text{d}(\phi,\sigma^{2})$
and $\kappa^{(1)},\ldots,\kappa^{(M)}$ are from the prior $p(\kappa)$.

The conditional Monte Carlo estimator $\widehat{p(\mathbf{y})}_{\rm CMC}$
alone does not work well in practice, because most draws from the prior
distributions $p(\mathbf{h})$ and $p(\kappa)$ do not coincide with the high-
density region of $p(\mathbf{Y}\,|\,\mathbf{h},\kappa)$ (as a function of
$\mathbf{h}$ and $\kappa$). To tackle this problem, next we combine the
conditional Monte Carlo approach with importance sampling, where the
importance sampling density is obtained using the improved cross-entropy
method. Chan and Eisenstat (2015) first propose using the improved cross-
entropy method to estimate the marginal likelihood, but the models considered
there are low dimensional. In contrast, for our settings with a large number
of latent variables, additional consideration on parsimonious parameterization
is required for obtaining appropriate importance sampling densities.

More specifically, we consider the parametric family

$\mathcal{F}=\\{f_{\mathcal{N}}(\mathbf{h};\widehat{\mathbf{h}},\widehat{\mathbf{K}}_{\mathbf{h}}^{-1})f_{\mathcal{N}}(\phi;\widehat{\phi},\widehat{K}_{\phi}^{-1})f_{\mathcal{G}}(\kappa;\widehat{\nu}_{\kappa},\widehat{S}_{\kappa})\\},$

where $f_{\mathcal{N}}$ and $f_{\mathcal{G}}$ denote Gaussian and gamma
densities, respectively. Moreover, $\widehat{\mathbf{h}}$ and
$\widehat{\mathbf{K}}_{\mathbf{h}}$ are respectively the mean vector and
precision matrix—i.e., inverse covariance matrix—of the $T$-variate Gaussian
density of $\mathbf{h}$; $\widehat{\phi}$ and $\widehat{K}_{\phi}$ are the
mean and precision of the univariate Gaussian density of $\phi$; and
$\widehat{\nu}_{\kappa}$ and $\widehat{S}_{\kappa}$ are respectively the shape
and rate of the gamma importance density for $\kappa$.666Note that
$\sigma^{2}$ only appears in the prior of $\mathbf{h}$, and one can integrate
it out analytically: $p(\mathbf{h})=\int
p(\mathbf{h}\,|\,\phi,\sigma^{2})p(\phi,\sigma^{2})\text{d}(\phi,\sigma^{2})=\int
p(\mathbf{h}\,|\,\phi)p(\phi)\text{d}\phi$, where $p(\mathbf{h}\,|\,\phi)$ has
an analytical expression, which is given in Appendix A. Hence, there is no
need to simulate $\sigma^{2}$. Now, we aim to choose these parameters so that
the associated member in $\mathcal{F}$ is the closest in cross-entropy
distance to the theoretical zero-variance importance sampling density
$p(\mathbf{h},\phi,\kappa\,|\,\mathbf{y})\propto
p(\mathbf{Y}\,|\,\mathbf{h},\kappa)p(\mathbf{h}\,|\,\phi)p(\phi,\kappa)$,
which is simply the marginal posterior density. Draws from this marginal
distribution can be obtained using the posterior sampler described in Appendix
A—i.e., we obtain posterior draws from the full posterior
$p(\mathbf{A},\boldsymbol{\Sigma},\mathbf{h},\phi,\sigma^{2},\kappa\,|\,\mathbf{y})$
and keep only the draws of $\mathbf{h},\phi,$ and $\kappa$.

Now, given the posterior draws
$(\mathbf{h}^{(1)},\phi^{(1)},\kappa^{(1)}),\ldots,(\mathbf{h}^{(M)},\phi^{(M)},\kappa^{(M)})$
and the parametric family $\mathcal{F}$, the optimization problem in (9) can
be divided into 3 lower-dimensional problems: 1) obtain $\widehat{\mathbf{h}}$
and $\widehat{\mathbf{K}}_{\mathbf{h}}$ given the posterior draws of
$\mathbf{h}$; 2) obtain $\widehat{\phi}$ and $\widehat{K}_{\phi}$ given the
posterior draws of $\phi$; and 3) obtain $\widehat{\nu}_{\kappa}$ and
$\widehat{S}_{\kappa}$ given the posterior draws of $\kappa$. The latter two
steps are easy as both densities are univariate, and they can be solved using
similar methods as described in Chan and Eisenstat (2015). Below we focus on
the first step.

If we do not impose any restrictions on the Gaussian density
$f_{\mathcal{N}}(\mathbf{h};\widehat{\mathbf{h}},\widehat{\mathbf{K}}_{\mathbf{h}}^{-1})$,
in principle we can obtain $\widehat{\mathbf{h}}$ and
$\widehat{\mathbf{K}}_{\mathbf{h}}$ analytically given the posterior draws
$\mathbf{h}^{(1)},\ldots,\mathbf{h}^{(M)}$. However, there are two related
issues with this approach, First, if unrestricted, the precision matrix
$\widehat{\mathbf{K}}_{\mathbf{h}}$ is a full, $T\times T$ matrix. Evaluating
and sampling from this Gaussian density would be time-consuming. Second, since
$\widehat{\mathbf{K}}_{\mathbf{h}}$ is a symmetric but otherwise unrestricted
matrix, there are $T(T+1)/2$ parameters to be estimated. Consequently, one
would require a large number of posterior draws to ensure that
$\widehat{\mathbf{h}}$ and $\widehat{\mathbf{K}}_{\mathbf{h}}$ are accurately
estimated.

In view of these potential difficulties, we consider instead a restricted
family of Gaussian densities that exploits the time-series structure.
Specifically, this family is parameterized by the parameters $\rho$,
$\mathbf{a}=(a_{1},\ldots,a_{T})^{\prime}$ and
$\mathbf{b}=(b_{1},\ldots,b_{T})^{\prime}$ as follows. First, $h_{1}$ has the
marginal Gaussian distribution $h_{1}\sim\mathcal{N}(a_{1},b_{1})$. For
$t=2,\ldots,T,$

$h_{t}=a_{t}+\rho h_{t-1}+\eta_{t},\quad\eta_{t}\sim\mathcal{N}(0,b_{t}).$
(11)

In other words, the joint distribution of $\mathbf{h}$ is implied by an AR(1)
process with time-varying intercepts and variances. It is easy to see that
this parametric family includes the prior density of $\mathbf{h}$ implied by
the state equation (3) as a member.777We also investigate other Gaussian
families such as the joint distributions implied by the AR(2), MA(1) and
ARMA(1,1,) processes. None of them perform noticeably better than the simple
AR(1) process in (11).

To facilitate computations, we vectorize the AR(1) process and write

$\mathbf{H}_{\rho}\mathbf{h}=\mathbf{a}+\boldsymbol{\eta},\quad\boldsymbol{\eta}\sim\mathcal{N}(\mathbf{0},\mathbf{B}),$

where $\mathbf{B}=\text{diag}(b_{1},\ldots,b_{T})$ and

$\mathbf{H}_{\rho}=\begin{pmatrix}1&0&\cdots&0\\\ -\rho&1&\cdots&0\\\
\vdots&\ddots&\ddots&\vdots\\\ 0&\cdots&-\rho&1\end{pmatrix}.$

Since $|\mathbf{H}_{\rho}|=1$ for any $\rho$, $\mathbf{H}_{\rho}$ is
invertible. Then, the Gaussian distribution implied by the AR(1) process in
(11) has the form
$\mathbf{h}\sim\mathcal{N}(\mathbf{H}_{\rho}^{-1}\mathbf{a},(\mathbf{H}_{\rho}^{\prime}\mathbf{B}^{-1}\mathbf{H}_{\rho})^{-1})$.
Note that the number of parameters here is only 2$T$ \+ 1 instead of
unrestricted case of $T+T(T+1)/2$.

Next, given this parametric family and posterior draws
$\mathbf{h}^{(1)},\ldots,\mathbf{h}^{(M)}$, we can solve the maximization
problem (9) in two steps. First, note that given $\rho$, we can obtain the
maximizer
$\widehat{\mathbf{a}}=(\widehat{a}_{1},\ldots,\widehat{a}_{T})^{\prime}$ and
$\widehat{\mathbf{b}}=(\widehat{b}_{1},\ldots,\widehat{b}_{T})^{\prime}$
analytically. More specifically, by maximizing the log-likelihood

$\ell(\rho,\mathbf{a},\mathbf{b})=-\frac{TM}{2}\log(2\pi)-\frac{M}{2}\log|\mathbf{B}|-\frac{1}{2}\sum_{m=1}^{M}(\mathbf{H}_{\rho}\mathbf{h}^{(m)}-\mathbf{a})^{\prime}\mathbf{B}^{-1}(\mathbf{H}_{\rho}\mathbf{h}^{(m)}-\mathbf{a})$

with respective to $\mathbf{a}$ and $\mathbf{b}$, we obtain the maximizer
$\widehat{\mathbf{a}}=\frac{1}{M}\sum_{m=1}^{M}\mathbf{H}_{\rho}\mathbf{h}^{(m)}$
and
$\widehat{b}_{t}=\frac{1}{M}\sum_{m=1}^{M}(h^{(m)}_{t}-\widehat{a}_{t})^{2},t=1,\ldots,T$.

Second, given the analytical solution $\widehat{\mathbf{a}}$ and
$\widehat{\mathbf{b}}$—which are functions of $\rho$—we can readily evaluate
the one-dimensional concentrated log-likelihood
$\ell(\rho,\widehat{\mathbf{a}},\widehat{\mathbf{b}})$. Then, $\widehat{\rho}$
can be obtained numerically by maximizing
$\ell(\rho,\widehat{\mathbf{a}},\widehat{\mathbf{b}})$ with respect to $\rho$.
Finally, we use
$f_{\mathcal{N}}(\mathbf{h};\mathbf{H}_{\widehat{\rho}}^{-1}\widehat{\mathbf{a}},(\mathbf{H}_{\widehat{\rho}}^{\prime}\widehat{\mathbf{B}}^{-1}\mathbf{H}_{\widehat{\rho}})^{-1})$,
where
$\widehat{\mathbf{B}}=\text{diag}(\widehat{b}_{1},\ldots,\widehat{b}_{T})$, as
the importance sampling density for $\mathbf{h}$.

Chan and Eisenstat (2018) also consider a Gaussian importance sampling density
for the log-volatility $\mathbf{h}$. However, there the approximation is based
on a second-order Taylor expansion at the mode. Hence, it is a local
approximation that might not be close to the target density at points away
from the mode. In contrast, our proposed Gaussian importance sampling density
is a global approximation that takes into account of the whole support of the
distribution, and is therefore expected to behave more similarly to the target
distribution $p(\mathbf{h}\,|\,\mathbf{y})$ over its entire support.

##  4 Monte Carlo Experiments

In this section we conduct a series of Monte Carlo experiments to assess how
the proposed marginal likelihood estimators perform in selecting the correct
data generating process. More specifically, we focus on the following three
questions. First, can we distinguish the three stochastic volatility models
described in Section 2? Second, can we discriminate between models with time-
varying volatility against homoskedastic errors? Finally, for factor
stochastic volatility models, can we identify the correct number of factors?

###  4.1 Can We Distinguish the Stochastic Volatility Specifications?

In the first part of this Monte Carlo exercise, we address the question of
whether one can distinguish the three stochastic volatility models—namely, the
common stochastic volatility model (VAR-CSV), the Cholesky stochastic
volatility model (VAR-SV) and the factor stochastic volatility model (VAR-
FSV)—in the context of large Bayesian VARs. To that end, we generate 100
datasets from each of the three models. Each dataset consists of $n=10$
variables, $T=400$ observations and $p=2$ lags. Given each dataset, we compute
the log marginal likelihoods of the same three models using the proposed
method described in Section 3.

For VAR-CSV, we generate the intercepts from $\mathcal{U}(-10,10)$. The
diagonal elements of the first VAR coefficient matrix are iid
$\mathcal{U}(-0.2,0.4)$ and the off-diagonal elements are
$\mathcal{U}(-0.2,0.2)$; all elements of the second VAR coefficient matrix are
iid $\mathcal{N}(0,0.05^{2})$. The error covariance matrix
$\boldsymbol{\Sigma}$ is generated from the inverse-Wistart distribution
$\mathcal{IW}(n+5,0.7\times\mathbf{I}_{n}+0.3\times\mathbf{1}_{n}\mathbf{1}_{n}^{\prime})$,
where $\mathbf{1}_{n}$ is an $n\times 1$ column of ones. For the log-
volatility process, we set $\phi=0.98$ and $\sigma^{2}=0.1$. For VAR-SV, the
VAR coefficients are generated as before, and the free elements of the impact
matrix are generated independently from the $\mathcal{N}(0,0.5^{2})$
distribution. For the log-volatility process, we set $\mu_{i}=-1$,
$\phi_{i}=0.98$ and $\sigma^{2}_{i}=0.1$ for $i=1,\ldots,n$. Finally, for VAR-
FSV, we set the number of factor $r$ to be 3. The VAR coefficients are
generated the same way as other models. For the log-volatility process, we set
$\mu_{i}=-1$, $\phi_{i}=0.98$ and $\sigma^{2}_{i}=0.1$ for $i=1,\ldots,n$, and
$\mu_{n+j}=0$, $\phi_{n+j}=0.98$ and $\sigma^{2}_{n+j}=0.1$ for
$j=1,\ldots,r$. That is, the log stochastic volatility processes associated
with the factors have a larger mean, but are otherwise the same as the
idiosyncratic stochastic volatility processes.888Since the main goal of the
Monte Carlo experiments is to assess if one can distinguish the stochastic
volatility specifications, we aim to use similar priors across the models so
as to minimize their impacts. In particular, for the hierarchical Minnesota
priors we fix all shrinkage hyperparameters to be
$\kappa=\kappa_{1}=\kappa_{2}=0.2^{2}$. That is, we turn off the cross-
variable shrinkage feature of the Minnesota priors under VAR-SV and VAR-FSV,
so that they are comparable to the symmetric Minnesota prior under VAR-CSV. In
the first experiment, we generate 100 datasets from VAR-CSV as described
above. For each dataset, we then compute the log marginal likelihoods of VAR-
SV and VAR-FSV relative to that of the true model VAR-CSV. Specifically, we
subtract the latter log marginal likelihood from the log marginal likelihoods
of VAR-SV and VAR-FSV. The results are reported in Figure 1.

Since a model is preferred by the data if it has a larger log marginal
likelihood value, a difference that is negative indicates that the correct
model is favored. It is clear from the histograms that for all the datasets
the correct model VAR-CSV compares favorably to the other two stochastic
volatility specifications, often by a large margin.

Figure 1: Histograms of log marginal likelihoods under VAR-SV (left panel) and
VAR-FSV (right panel) relative to the true model (VAR-CSV). A negative value
indicates that the correct model is favored.

Next, we generate 100 datasets from VAR-SV. For each dataset, we then compute
the log marginal likelihoods of VAR-CSV and VAR-FSV relative to that of the
true model. The results are reported in Figure 2. Again, the model comparison
result shows that, for all datasets, the correct model VAR-SV is
overwhelmingly favored compared to the other two stochastic volatility
specifications.

Figure 2: Histograms of log marginal likelihoods under VAR-CSV (left panel)
and VAR-FSV (right panel) relative to the true model (VAR-SV). A negative
value indicates that the correct model is favored. Figure 3: Histograms of log
marginal likelihoods under VAR-CSV (left panel) and VAR-SV (right panel)
relative to the true model (VAR-FSV). A negative value indicates that the
correct model is favored.

Lastly, we generate 100 datasets from VAR-FSV, and for each dataset we compute
the log marginal likelihoods of VAR-CSV and VAR-SV relative to that of VAR-
FSV. The results are reported in Figure 3. Once again, these results show that
the proposed marginal likelihood estimators can select the correct model for
all the simulated datasets.

Appendix C provides additional Monte Carlo results to further assess the
performance of the proposed marginal likelihood estimators in higher
dimensions ($n=20$) and noisier VAR equations (magnitudes of the error
covariance matrices increase 10-fold). For these additional settings, the
correct model is again overwhelmingly favored compared to the other two
stochastic volatility models. All in all, this series of Monte Carlo
experiments show that using the proposed marginal likelihood estimators one
can clearly distinguish the three stochastic volatility specifications, even
for a moderate number of variables and sample size.

###  4.2 Can We Discriminate between Models with Time-Varying Volatility
against Homoskedastic Errors?

Since VARs with stochastic volatility are by design more flexible than
conventional VARs with homoskedastic errors, one concern of using these more
flexible models is that they might overfit the data. In this Monte Carlo
experiment we investigate if the proposed marginal likelihood estimators can
distinguish models with and without stochastic volatility. More specifically,
we generate 100 datasets from a standard VAR with homoskedasic errors. For
each dataset, we then compute the log marginal likelihoods of the three
stochastic volatility models: VAR-CSV, VAR-SV and VAR-FSV.999We also generate
data from the stochastic volatility models and compare their marginal
likelihoods with a homoskedastic VAR. The results overwhelmingly favor the
stochastic volatility models as the homoskedastic VAR does not fit the data
well at all. For space constraint, however, we do not report these model
comparison results. The differences in log marginal likelihoods relative to
the homoskedastic VAR are reported in Figure 4.

Recall that a model is preferred by the data if it has a larger log marginal
likelihood value. Hence, a negative value indicates that the correct
homoskedastic VAR is selected. Overall, the histograms show that for almost
all datasets the correct homoskedastic model compares favorably to the three
stochastic volatility models. This Monte Carlo experiment highlights that the
marginal likelihood does not always favor the model with the best model-fit;
it also has a built-in penalty for model complexity. And it is only when the
gain in model-fit outweighs the additional model complexity does the more
complex model have a larger marginal likelihood value.

Figure 4: Histograms of log marginal likelihoods under VAR-CSV (left panel),
VAR-SV (middle panel) and VAR-FSV (right panel) relative to the true model
(homoskedastic VAR). A negative value indicates that the correct model is
favored.

It is also interesting to note that VAR-CSV generally performs better than the
other two stochastic volatility models. This is perhaps not surprising, as the
VAR-CSV is the most parsimonious extension of the homoskedastic VAR—it has
only one stochastic volatility process $h_{t}$, and when $h_{t}$ is
identically zero it replicates the homoskedastic VAR. Nevertheless, even for
this closely related extension, the marginal likelihood criterion clearly
indicates that the stochastic volatility process is spurious and it tends to
favor the homoskedastic model.

###  4.3 Can We Identify the Correct Number of Factors?

One important specification choice for factor stochastic volatility models is
to select the number of factors. This represents the trade-off between a more
parsimonious model with fewer stochastic volatility factors versus a more
complex model with more factors but better model-fit. In this Monte Carlo
experiment, we investigate if the proposed marginal likelihood estimator for
the VAR-FSV can pick the correct number of factors.

To that end, we generate 100 datasets from VAR-FSV with $r=3$ factors. Then,
for each dataset, we compute the log marginal likelihood of VAR-FSV models
with $r=2,3$ and 4 factors. This set of number of factors is chosen to shed
light on the effects of under-fitting versus over-fitting. We report the log
marginal likelihoods of the 2- and 4-factor models relative to the true
3-factor model in Figure 5. The results clearly show that the proposed method
is able to identify the correct number of factors. In particular, for all
datasets the 3-factor model outperforms the more parsimonious 2-factor model
and the more flexible 4-factor model.

Figure 5: Histograms of log marginal likelihoods under the 2-factor model
(left panel) and the 4-factor model (right panel) relative to the true
3-factor model. A negative value indicates that the correct model is favored.

The 2-factor model under-fits the data and performs much worse than the
correct 3-factor model. On the other hand, while the 4-factor model is able to
replicate features of the 3-factor model, it also includes spurious features
that tend to over-fit the data. Consequently, the 3-factor model is strongly
preferred by the data over the two alternatives. However, the impacts of
under- and over-fitting are not symmetric. In particular, under-fitting
receives a much heavier penalty than over-fitting, as illustrated by the much
lager differences (in magnitude) between the 2- and 3-factor models than those
between the 3- and 4-factor models. Regardless of these differences, we
conclude that the proposed method is able to select the correct number of
factors.

##  5 Empirical Application

In this empirical application we compare the three stochastic volatility
specifications—i.e., the common stochastic volatility model (VAR-CSV), the
Cholesky stochastic volatility model (VAR-SV) and the factor stochastic
volatility model (VAR-FSV)—in the context of Bayesian VARs of different
dimensions. We use a dataset that consists of 30 US quarterly variables with a
sample period from 1959Q1 to 2019Q4. It is constructed from the FRED-QD
database at the Federal Reserve Bank of St. Louis as described in McCracken
and Ng (2020). The dataset contains a range of widely used macroeconomic and
financial variables, including Real GDP and its components, various measures
of inflation, labor market variables and interest rates. They are transformed
to stationarity, typically to annualized growth rates. We consider VARs that
are small ($n=7$), medium ($n=15$) and large ($n=30$). The complete list of
variables for each dimension and how they are transformed is given in Appendix
B.

In contrast to the Monte Carlo study where the shrinkage hyperparameters in
the hierarchical Minnesota priors are fixed at certain subjective values, here
we treat them as unknown parameters to be estimated. This is motivated by a
few recent papers that find this data-based approach outperforms subjective
prior elicitation in both in-sample model fit and out-of-sample forecast
performance (see, e.g., Carriero, Clark, and Marcellino, 2015; Giannone,
Lenza, and Primiceri, 2015; Amir-Ahmadi, Matthes, and Wang, 2020). For easy
comparison, we set the lag length to be $p=4$ for all VARs.101010The lag
length can of course be chosen by comparing the marginal likelihood. In a
preliminary analysis, we find that a lag length of $p=4$ is generally
sufficient. In addition, with our flexible hierarchical Minnesota priors on
the VAR coefficients, adding more lags than necessary does not substantially
affect model performance. For instance, for all of the three stochastic
volatility specifications, the best models for $n=7$ have $p=3$ lags. Adding
one more lag reduces the log marginal likelihood by only about 1-2 for all
models. We compute the log marginal likelihoods of the VARs using the proposed
hybrid algorithm described in Section 3. The importance sampling density for
each model is constructed using the cross-entropy method with 20,000 posterior
draws after a burn-in period of 1,000. Then, the log marginal likelihood
estimate is computed using an importance sample of size 10,000.

###  5.1 Comparing Stochastic Volatility Specifications

We first report the main model comparison results on comparing the three
stochastic volatility specifications: VAR-CSV, VAR-SV and VAR-FSV. As a
benchmark, we also include the standard homoskedastic VAR with the natural
conjugate prior (its marginal likelihood is available analytically). In
addition to the widely different likelihoods implied by these stochastic
volatility specifications, they also differ in the flexibility of the
shrinkage priors employed. More specifically, both the standard VAR and VAR-
CSV use the natural conjugate prior, which has only one hyperparameter that
controls the shrinkage strength of all VAR coefficients. In contrast, both the
priors under VAR-SV and VAR-FSV can accommodate cross-variable shrinkage. That
is, there are two separate shrinkage hyperparameters, one controls the
shrinkage strength on own lags, whereas the other controls that of lags of
other variables. In what follows we focus on the overall ranking of the
models; we will investigate the role of the priors in the next section.

Table 1 reports the log marginal likelihood estimates of the four VARs across
the three model dimensions ($n=7,15,30$). First, it is immediately apparent
that all three stochastic volatility models perform substantially better than
the standard homoskedastic VAR. For example, the log marginal likelihood
difference between VAR-CSV and VAR is about 760 for $n=30$, highlighting
overwhelming evidence in favor of the heteroskedastic model. This finding is
in line with the growing body of evidence that shows the importance of time-
varying volatility in modeling both small and large macroeconomic datasets.

Table 1: Log marginal likelihood estimates (numerical standard errors) of a standard homoskedastic VAR, VAR-CSV, VAR-SV and VAR-FSV. | VAR | VAR-CSV | VAR-SV | VAR-FSV
---|---|---|---|---
|  |  |  | $k=1$ | $k=2$ | $k=3$
$n=7$ |  $-$2,583.4 |  $-$2,410 |  $-$2,315 |  $-$2,318 |  $-$2,320 |  $-$2,331
| - | (0.1) | (0.3) | (0.3) | (0.7) | (0.7)
|  |  |  | $k=3$ | $k=4$ | $k=5$
$n=15$ |  $-$6,918.8 |  $-$6,618 |  $-$6,443 |  $-$6,468 |  $-$6,454 |  $-$6,485
| - | (0.1) | (0.4) | (0.7) | (0.8) | (1.3)
|  |  |  | $k=9$ | $k=10$ | $k=11$
$n=30$ |  $-$12,783.4 |  $-$12,024 |  $-$11,382 |  $-$11,571 |  $-$11,567 |  $-$11,608
| - | (0.1) | (0.6) | (1.2) | (1.8) | (1.8)

Second, among the three stochastic volatility specifications, the data
overwhelmingly prefers VAR-SV and VAR-FSV over the more restrictive VAR-CSV
for all model dimensions, due to a combination of the more flexible
likelihoods and priors. For example, the log marginal likelihood differences
between VAR-SV and VAR-CSV are 95 for $n=7$, 175 for $n=15$ and 642 for
$n=30$. Third, for all three model dimensions, VAR-SV is the most favored
stochastic volatility specification, though VAR-FSV comes in close second.
Finally, we note that the optimal number of factors changes across the model
dimension. It is perhaps not surprising that more factors are needed to model
the more complex error covariance structure as the model dimension increases.
For instance, for $n=7$ the 1-factor model performs the best, whereas one
needs 10 factors when $n=30$.

To corroborate these model comparison results, we next perform a recursive
out-of-sample forecasting exercise using an evaluation period from 1975Q1 to
the end of the sample. More specifically, Table 2 reports the joint density
forecasting results (as measured by the sum of log predictive likelihoods) for
the four VARs across the three model dimensions (for VAR-FSV, we fix $k=4$).
It is clear from the table that one can draw mostly similar conclusions from
these forecasting results. In particular, all three stochastic volatility
models perform substantially better than the homoskedastic VAR for both 1- and
4-step-ahead forecast horizons and for all model dimensions. This again
highlights the importance of allowing for some form of time-varying
volatility. In addition, for 1-step-ahead, both VAR-SV and VAR-FSV forecast
better than the more restrictive VAR-CSV for all model dimensions. However,
when the forecast horizon increases, VAR-CSV tends to perform better for
higher model dimensions.

Table 2: Sum of log predictive likelihoods of a standard homoskedastic VAR,
VAR-CSV, VAR-SV and VAR-FSV. A larger value indicates better forecasting
performance.

| 1-step-ahead | 4-step-ahead
---|---|---
| VAR | VAR-CSV | VAR-SV | VAR-FSV | VAR | VAR-CSV | VAR-SV | VAR-FSV
$n=7$ |  $-$1,858 |  $-$1,719 |  $-$1,625 |  $-$1,613 |  $-$2,573 |  $-$2,321 |  $-$2,311 |  $-$2,312
$n=15$ |  $-$5,070 |  $-$4,764 |  $-$4,612 |  $-$4,604 |  $-$6,243 |  $-$5,631 |  $-$5,821 |  $-$5,772
$n=30$ |  $-$11,271 |  $-$8,579 |  $-$8,035 |  $-$8,424 |  $-$15,715 |  $-$10,030 |  $-$11,739 |  $-$11,348

###  5.2 Comparing Shrinkage Priors

In this section we compare different types of Minnesota priors for each of the
three stochastic volatility specifications. In particular, we investigate the
potential benefits of allowing for cross-variable shrinkage and selecting the
overall shrinkage hyperparameters in a data-driven manner. To that end, we
consider two useful benchmarks. First, for VAR-SV and VAR-FSV we consider the
special case where $\kappa_{1}=\kappa_{2}$—i.e., we turn off the cross-
variable shrinkage and require the shrinkage hyperparameters on own and other
lags to be the same. We refer to this version as the symmetric prior. The
second benchmark is a set of subjectively chosen hyperparameter values that
apply cross-variable shrinkage. In particular, we follow Carriero, Clark, and
Marcellino (2015) and use the values $\kappa_{1}=0.04$ and
$\kappa_{2}=0.0016$. This second benchmark is referred to as the subjective
prior. Finally, our baseline prior, where $\kappa_{1}$ and $\kappa_{2}$ are
estimated from the data and could potentially be different, is referred to as
the asymmetric prior.

To fix ideas, we focus on VARs with $n=15$ variables. Table 3 reports the log
marginal likelihood estimates of the three stochastic volatility
specifications with different shrinkage priors. First, for both VAR-SV and
VAR-FSV, the asymmetric prior significantly outperforms the symmetric version
that requires $\kappa_{1}=\kappa_{2}$. This result suggests that it is
beneficial to shrink the coefficients on own lags differently than those on
lags of other variables. This makes intuitive sense as one would expect that,
on average, a variable’s own lags would be more important for its future
evolution than lags of other variables. By relaxing the restriction that
$\kappa_{1}=\kappa_{2}$, the log marginal likelihood values of VAR-SV and VAR-
FSV increase by 177 and 204, respectively.

Table 3: Log marginal likelihood estimates (numerical standard errors) of the stochastic volatility specifications with different shrinkage priors for $n=15$. | VAR-CSV | VAR-SV | VAR-FSV ($k=4$)
---|---|---|---
Subjective prior |  $-$6,702 |  $-$6,474 |  $-$6,491
| (0.1) | (0.6) | (0.9)
Symmetric prior |  $-$6,618 |  $-$6,620 |  $-$6,658
| (0.1) | (0.4) | (1.1)
Asymmetric prior | - |  $-$6,443 |  $-$6,454
|  | (0.7) | (0.8)

In addition, there are also substantial benefits of allowing the shrinkage
hyperparameters to be estimated instead of fixing them subjectively. For
example, the log marginal likelihood value of VAR-CSV with the symmetric prior
is 84 larger than that of the subjective prior; the log marginal likelihood
value of VAR-SV with the asymmetric prior is 31 larger than that of the
subjective prior. These results suggest that while those widely-used
subjective hyperparameter values seem to work well for some datasets, they
might not be suitable for others that contain different variables and span
different sample periods.

We reported in last section that the data overwhelmingly preferred the more
flexible VAR-SV and VAR-FSV over VAR-CSV. But it was unclear whether it was
due to the more flexible likelihoods or priors (recall that VAR-CSV can only
accommodate the symmetric prior). The result here suggests that the superior
performance of VAR-SV and VAR-FSV can be attributed to the more flexible
priors rather than the more flexible stochastic volatility specifications. In
particular, both VAR-SV and VAR-FSV with the symmetric prior actually perform
slightly worse than VAR-CSV (respectively, $-$6,620, $-$6,658 and $-$6,618 for
the three models). Only when the asymmetric prior is used do VAR-SV and VAR-
FSV outperform VAR-CSV. This stark result shows that choosing a suitable
shrinkage prior in high-dimensional settings is as important as, if not more
important than, selecting a suitable stochastic volatility specification.

Next, Table 4 reports the posterior estimates of the shrinkage hyperparameters
under the three stochastic volatility specifications. First, under VAR-CSV,
the posterior mean of $\kappa$ varies across the model dimensions, from 0.32
for $n=7$ to 0.22 for $n=15$ and 0.1 for $n=30$. This finding highlights the
fact that any particular subjectively elicited value—e.g., the widely-used
value of 0.04 for the natural conjugate prior; see Carriero, Clark, and
Marcellino (2016) and Chan (2020)—is unlikely to be suitable for all datasets
of different dimensions. In addition, the estimates also suggest that when the
model dimension increases, more aggressive shrinkage is needed, confirming the
findings in, e.g., Banbura, Giannone, and Reichlin (2010).

Second, it is clear that the estimates of $\kappa_{1}$ are orders of magnitude
larger than those of $\kappa_{2}$ under VAR-SV and VAR-FSV, the two models
that accommodate cross-variable shrinkage. For example, for VAR-FSV with
$n=15$, the estimates of $\kappa_{1}$ and $\kappa_{2}$ are, respectively, 0.17
and 0.0018, a difference of 2 orders of magnitude. These estimates confirm the
model comparison results above and suggest that the data strongly prefers
shrinking the coefficients on lags of other variables much more aggressively
to zero than those on own lags. Again, this is intuitive as one would expect
that own lags are more informative, on average, than lags of other variables
in forecasting a variable’s future evolution. In addition, these estimates are
also rather different from those widely-used subjectively chosen values of
$\kappa_{1}=0.04$ and $\kappa_{2}=0.0016$, reinforcing the conclusion that any
particular set of fixed hyperparameter values is unlikely to be suitable for
all datasets of different dimensions.

Table 4: Posterior means (standard deviations) of the shrinkage hyperparameters in VAR-CSV, VAR-SV and VAR-FSV. | VAR-CSV | VAR-SV | VAR-FSV
---|---|---|---
| $\kappa$ | $\kappa_{1}$ | $\kappa_{2}$ | $\kappa_{1}$ | $\kappa_{2}$
$n=7$ | 0.32 | 0.23 | 0.0032 | 0.24 | 0.0037
| (0.065) | (0.051) | (0.0014) | (0.054) | (0.0015)
$n=15$ | 0.22 | 0.16 | 0.0022 | 0.17 | 0.0018
| (0.027) | (0.030) | (0.0006) | (0.032) | (0.0004)
$n=30$ | 0.10 | 0.15 | 0.00002 | 0.15 | 0.00003
| (0.008) | (0.020) | (0.00001) | (0.020) | (0.00001)

All in all, our results confirm the substantial benefits of allowing for
cross-variable shrinkage and selecting the shrinkage hyperparameters using a a
data-based approach. They also highlight that in high-dimensional settings,
choosing a flexible shrinkage prior is as important as selecting a flexible
stochastic volatility specification.

##  6 Extensions

The proposed marginal likelihood estimators are based on two key ingredients:
a conditional Monte Carlo component to integrate out the large number of VAR
coefficients analytically and an adaptive importance sampling method to
integrate out the log-volatility via Monte Carlo simulation. These two ideas
can be more widely applied to other marginal likelihood estimators or
alternative model comparison methods to make them viable in high-dimensional
settings, as we demonstrate below. In addition, we further illustrate how the
proposed methods can be used to compute the marginal likelihood of a
stochastic volatility model with an outlier component.

###  6.1 Improving the Modified Harmonic Mean Estimator

We first illustrate how one can use the conditional Monte Carlo and adaptive
importance sampling methods to improve alternative estimators of the marginal
likelihood. To fix ideas, we focus on the modified harmonic mean estimator
proposed by Gelfand and Dey (1994); similar improvements can be made for other
estimators such as Chib’s method (Chib, 1995; Chib and Jeliazkov, 2001) and
sequential Monte Carlo (e.g., Bognanni and Zito, 2020).

The theoretical justification of the modified harmonic mean estimator is based
on the following observation: for any density function $f$ with support
contained in the support of the posterior density, we have

$\mathbb{E}\left(\frac{f(\boldsymbol{\theta})}{p(\boldsymbol{\theta})p(\mathbf{y}\,|\,\boldsymbol{\theta})}\;\bigg{|}\;\mathbf{y}\right)=\int\frac{f(\boldsymbol{\theta})}{p(\boldsymbol{\theta})p(\mathbf{y}\,|\,\boldsymbol{\theta})}\frac{p(\boldsymbol{\theta})p(\mathbf{y}\,|\,\boldsymbol{\theta})}{p(\mathbf{y})}\,\text{d}\boldsymbol{\theta}=p(\mathbf{y})^{-1},$
(12)

where the expectation is taken with respect to the posterior distribution
$p(\boldsymbol{\theta}\,|\,\mathbf{y})=p(\boldsymbol{\theta})p(\mathbf{y}\,|\,\boldsymbol{\theta})/p(\mathbf{y})$.
It follows that one can estimate the marginal likelihood $p(\mathbf{y})$ using
the following estimator:

$\text{GD}=\left\\{\frac{1}{R}\sum_{r=1}^{R}\frac{f(\boldsymbol{\theta}^{(r)})}{p(\boldsymbol{\theta}^{(r)})p(\mathbf{y}\,|\,\boldsymbol{\theta}^{(r)})}\right\\}^{-1},$
(13)

where $\boldsymbol{\theta}^{(1)},\ldots,\boldsymbol{\theta}^{(R)}$ are
posterior draws. This modified harmonic mean is simulation-consistent—it
converges to $p(\mathbf{y})$ in probability as $R$ tends to infinity—but it is
not unbiased in general.

Even though the estimator in (13) is simulation-consistent for any density $f$
with support contained in the support of
$p(\boldsymbol{\theta}\,|\,\mathbf{y})$, its finite-sample properties depend
critically on the choice of $f$. Geweke (1999) shows that if $f$ has tails
lighter than those of $p(\boldsymbol{\theta}\,|\,\mathbf{y})$, the estimator
has a finite variance. He further recommends a Gaussian approximation of
$p(\boldsymbol{\theta}\,|\,\mathbf{y})$ with tail truncations determined by
asymptotic arguments. More specifically, let $\widehat{\boldsymbol{\theta}}$
and $\mathbf{Q}_{\boldsymbol{\theta}}$ denote the posterior mean and
covariance matrix respectively. Then, $f$ is set to be the
$\mathcal{N}(\widehat{\boldsymbol{\theta}},\mathbf{Q}_{\boldsymbol{\theta}})$
density truncated within the region

$\\{\boldsymbol{\theta}\in\mathbb{R}^{m}:(\boldsymbol{\theta}-\widehat{\boldsymbol{\theta}})^{\prime}\mathbf{Q}_{\boldsymbol{\theta}}^{-1}(\boldsymbol{\theta}-\widehat{\boldsymbol{\theta}})<\chi^{2}_{\alpha,m}\\},$

where $\chi^{2}_{\alpha,m}$ is the $(1-\alpha)$ quantile of the $\chi^{2}_{m}$
distribution and $m$ is the dimension of $\boldsymbol{\theta}$. The GD
estimator in (13) with the truncated Gaussian tuning function is commonly-used
in economics; examples include Sims, Waggoner, and Zha (2008) and Justiniano
and Primiceri (2008).

One potential difficulty of implementing this approach in high-dimensional
setting is that the prescribed Gaussian density has a large number of
parameters to be estimated. As an example, for the VAR-SV model with
$T=300,n=20,$ and $p=4$, the Gaussian density is 2,133-dimensional and has
over 2.2 million parameters. As such, it is impractical to estimate them with
any reasonable accuracy using MCMC output.

Below we outline a few ways to make the GD estimator viable in high-
dimensional settings using the proposed adaptive importance sampling and
conditional Monte Carlo methods. First, we consider a tuning function of the
form
$f(\boldsymbol{\theta})=f_{\mathbf{h}}(\mathbf{h})f_{\boldsymbol{\theta}_{-\mathbf{h}}}(\boldsymbol{\theta}_{-\mathbf{h}})$,
where $\boldsymbol{\theta}_{-\mathbf{h}}$ denotes the set of model parameters
and latent variables excluding the log-volatility $\mathbf{h}$. Then, we use
the adaptive importance sampling method described in Section 3.2 to obtain
$f_{\mathbf{h}}$, whereas $f_{\boldsymbol{\theta}_{-\mathbf{h}}}$ is a
truncated Gaussian distribution constructed as described above. This approach
reduces the number of parameters of $f$ in two ways: 1) by factoring $f$ into
two lower-dimensional Gaussian densities, it reduces the number of covariance
parameters; 2) as described in Section 3.2, the way $f_{\mathbf{h}}$ is
parameterized reduces the number of parameters from $(T^{2}+3T)/2$ to
$2T+1$.111111As discussed in Section 3.4, the choice of the restricted family
of Gaussian densities exploits the time-series structure of the problem. In
particular, the log-volatility is assumed to follow an AR(1) process, which
implies the precision or inverse covariance matrix of the joint distribution
is tridiagonal. This strikes the balance between flexibility and parsimony.
However, this approach alone is unlikely to work well as there remains a large
number of parameters to be estimated. Using the same example of the VAR-SV
with $T=300,n=20,$ and $p=4$, the tuning function has still over 1.6 million
parameters.

Next, we incorporate the proposed conditional Monte Carlo estimator to improve
the GD estimator. Using the notations in Section 3.3, partition
$\boldsymbol{\theta}$ into two blocks, $\boldsymbol{\gamma}$ and
$\boldsymbol{\psi}$, and suppose we can compute
$g(\mathbf{y}\,|\,\boldsymbol{\psi})\equiv\mathbb{E}[p(\mathbf{y}\,|\,\boldsymbol{\theta})\,|\,\boldsymbol{\psi}]=\int
p(\mathbf{y}\,|\,\boldsymbol{\gamma},\boldsymbol{\psi})p(\boldsymbol{\gamma}\,|\,\boldsymbol{\psi})\text{d}\boldsymbol{\gamma}$
analytically, where $p(\boldsymbol{\gamma}\,|\,\boldsymbol{\psi})$ is the
prior of $\boldsymbol{\gamma}$ that can potentially depend on
$\boldsymbol{\psi}$. Since the marginal likelihood can be expressed as
$p(\mathbf{y})=\int
g(\mathbf{y}\,|\,\boldsymbol{\psi})p(\boldsymbol{\psi})\text{d}\boldsymbol{\psi}$,
a similar derivation as (12) suggests that we can estimate $p(\mathbf{y})$
using the following estimator:

$\text{GD}_{\boldsymbol{\psi}}=\left\\{\frac{1}{R}\sum_{r=1}^{R}\frac{f_{\boldsymbol{\psi}}(\boldsymbol{\psi}^{(r)})}{p(\boldsymbol{\psi}^{(r)})g(\mathbf{y}\,|\,\boldsymbol{\psi}^{(r)})}\right\\}^{-1},$
(14)

where $\boldsymbol{\psi}^{(1)},\ldots,\boldsymbol{\psi}^{(R)}$ are posterior
draws from
$p(\boldsymbol{\psi}\,|\,\mathbf{y})=g(\mathbf{y}\,|\,\boldsymbol{\psi})p(\boldsymbol{\psi})/p(\mathbf{y})$.
These draws can be obtained using a standard MCMC algorithm: first obtain
posterior draws of the pairs
$(\boldsymbol{\gamma}^{(1)},\boldsymbol{\psi}^{(1)}),\ldots,(\boldsymbol{\gamma}^{(R)},\boldsymbol{\psi}^{(R)})$
and then discard the first components. One can use the truncated Gaussian
density described above as the tuning function $f_{\boldsymbol{\psi}}$. Using
the same VAR-SV example, the number of parameters of the tuning function here
is drastically reduced to about 46,000. Even so, estimating these parameters
with sufficient accuracy using MCMC output is still computationally intensive.

Finally, one can combine both approaches to further improve the performance of
the estimator. Specifically, consider the estimator in (14) and a tuning
function of the form
$f_{\boldsymbol{\psi}}(\boldsymbol{\psi})=f_{\mathbf{h}}(\mathbf{h})f_{\boldsymbol{\psi}_{-\mathbf{h}}}(\boldsymbol{\psi}_{-\mathbf{h}})$,
where $f_{\mathbf{h}}$ is obtained using the adaptive importance sampling
method and $f_{\boldsymbol{\psi}_{-\mathbf{h}}}$ is a truncated Gaussian
density. For the same VAR-SV example, the number of parameters for this
combined approach is only 610, and they can be estimated with adequate
accuracy using MCMC output.

To illustrate how these improved estimators perform empirically, we revisit
the model comparison exercise in Section 5 and compute the marginal likelihood
of VAR-CSV with $n=15$ using two versions of the GD estimator. Recall that the
set of model parameters and latent variables for VAR-CSV is
$\boldsymbol{\theta}=\\{\mathbf{A},\boldsymbol{\Sigma},\phi,\sigma^{2},\kappa,\mathbf{h}\\}$.
We first consider the estimator in (14), where we integrate out
$\boldsymbol{\gamma}=\\{\mathbf{A},\boldsymbol{\Sigma},\sigma^{2}\\}$
analytically and set $f_{\boldsymbol{\psi}}$ to be a truncated Gaussian
density with $\boldsymbol{\psi}=\\{\phi,\kappa,\mathbf{h}\\}$. Second, we use
the same estimator in (14), but with a different tuning function: we set
$f_{\boldsymbol{\psi}}(\boldsymbol{\psi})=f_{\mathbf{h}}(\mathbf{h})f_{\boldsymbol{\psi}_{-\mathbf{h}}}(\boldsymbol{\psi}_{-\mathbf{h}})$,
where $f_{\mathbf{h}}$ is obtained using the adaptive importance sampling
method and $f_{\boldsymbol{\psi}_{-\mathbf{h}}}$ is a truncated bivariate
Gaussian density with $\boldsymbol{\psi}_{-\mathbf{h}}=\\{\phi,\kappa\\}$. For
each estimator, we obtain 100 marginal likelihood estimates, where each
estimate is based on 50,000 posterior draws. To summarize the results we
report the boxplots of the log marginal likelihood estimates in Figure 6. The
middle line of each box denotes the median, while the lower and upper lines
represent, respectively, the 25- and the 75-percentiles. The whiskers extend
to the maximum and minimum. For comparison we also include the estimates from
the proposed method.121212The original GD estimator in (13) and the version
with only the adaptive importance sampling method have too many parameters,
and they are unstable in this high-dimensional setting.

Figure 6: Boxplots of the log marginal likelihoods for VAR-CSV with $n=15$
estimated using three estimators: the GD estimator in which
$\boldsymbol{\gamma}=\\{\mathbf{A},\boldsymbol{\Sigma},\sigma^{2}\\}$ is
integrated out analytically (GD-1), the GD estimator in which
$\boldsymbol{\gamma}$ is integrated out analytically and the adaptive
importance sampling method is used to obtain $f_{\mathbf{h}}$ (GD-2) and the
proposed method.

While the original GD estimator in (13) has too many parameters to provide
stable estimates, the version in which
$\boldsymbol{\gamma}=\\{\mathbf{A},\boldsymbol{\Sigma},\sigma^{2}\\}$ is
integrated out analytically (GD-1) performs reasonably well. In particular,
the variance of the estimator is relatively small and most of the marginal
likelihood estimates are between $-6,630$ and $-6,620$. Nevertheless, the
distribution of the estimator has a heavier left tail and it seems to
underestimate the value of the marginal likelihood. In contrast, when both the
conditional Monte Carlo and adaptive importance sampling methods are used to
improve the GD estimator (GD-2), its variance is substantially reduced and its
performance is virtually the same as the proposed estimator that is based on
the cross-entropy method. In particular, all estimates of both estimators are
around $-6,618$.

###  6.2 The Observe-Data Deviance Information Criterion

Another commonly-used model comparison criterion is the deviance information
criterion (DIC) introduced in the seminal paper by Spiegelhalter, Best,
Carlin, and van der Linde (2002). Below we discuss how one can use the
proposed conditional Monte Carlo and adaptive importance sampling methods to
efficiently compute the DIC for high-dimensional stochastic volatility models.

The DIC is based on the concept of deviance, which is defined as

$D(\boldsymbol{\theta})=-2\log p(\mathbf{y}\,|\,\boldsymbol{\theta}),$

where $p(\mathbf{y}\,|\,\boldsymbol{\theta})$ is the (observed-data)
likelihood.131313Here we distinguish between the conditional likelihood, i.e.,
the data density conditional on the latent variables, and the observed-data
likelihood (data density marginal of the latent variables). While the DIC can
be defined in terms of the conditional likelihood (Celeux, Forbes, Robert, and
Titterington, 2006), a few recent papers have pointed out that this version is
problematic for theoretical and practical reasons (see, e.g., Chan and Grant,
2016; Li, Yu, and Zeng, 2020). We therefore focus on the original DIC based on
the observed-data likelihood. Next, the effective number of parameters,
$p_{D}$, of the model is defined to be

$p_{D}=\overline{D(\boldsymbol{\theta})}-D(\widetilde{\boldsymbol{\theta}}),$

where
$\overline{D(\boldsymbol{\theta})}=-2\mathbb{E}_{\boldsymbol{\theta}}[\log
p(\mathbf{y}\,|\,\boldsymbol{\theta})\,|\,\mathbf{y}]$ is the posterior mean
deviance and $\widetilde{\boldsymbol{\theta}}$ is an estimate of
$\boldsymbol{\theta}$, which is typically taken as the posterior mean or mode.
Finally, the deviance information criterion is the sum of the posterior mean
deviance, which can be viewed as a Bayesian measure of model fit, and the
effective number of parameters that measures model complexity:

$\text{DIC}=\overline{D(\boldsymbol{\theta})}+p_{D}.$

The DIC may therefore be interpreted as a trade-off between model fit and
model complexity. Given a set of competing models, the preferred model is the
one with the minimum DIC value.

The main difficulty in computing the DIC for high-dimensional stochastic
volatility models is that the likelihood is not available in closed-form and
it is typically computationally intensive to evaluate. We show how one can use
the conditional Monte Carlo and adaptive importance sampling methods developed
in Section 3 to overcome this difficulty. To fix ideas, we focus on the
likelihood evaluation of the VAR-FSV model; likelihood evaluation for other
stochastic volatility models can be done similarly.

The latent variables in the VAR-FSV model are the factors and the log-
volatility. We aim to integrate out the factors analytically and the log-
volatility using importance sampling. To that end, we stack the observations
$\mathbf{y}=(\mathbf{y}_{1}^{\prime},\ldots,\mathbf{y}_{T}^{\prime})^{\prime}$
over $t=1,\ldots,T$ and write the VAR-FSV model as

$\mathbf{y}=\mathbf{X}\boldsymbol{\alpha}+(\mathbf{I}_{T}\otimes\mathbf{L})\mathbf{f},\quad\mathbf{u}^{y}\sim\mathcal{N}(\mathbf{0},\mathbf{D}),\quad\mathbf{f}\sim\mathcal{N}(\mathbf{0},\mathbf{G}),$

where $\mathbf{X}$ is an appropriately defined design matrix consisting of
intercepts and lagged values,
$\mathbf{D}=\text{diag}(\mathbf{D}_{1},\ldots,\mathbf{D}_{T})$ and
$\mathbf{G}=\text{diag}(\mathbf{G}_{1},\ldots,\mathbf{G}_{T})$ with
$\mathbf{D}_{t}=\text{diag}(\text{e}^{h_{1,t}},\ldots,\text{e}^{h_{n,t}})$ and
$\mathbf{G}_{t}=\text{diag}(\text{e}^{h_{n+1,t}},\ldots,\text{e}^{h_{n+r,t}})$.
Hence, the distribution of $\mathbf{y}$ marginal of the factors $\mathbf{f}$
is
$(\mathbf{y}\,|\,\boldsymbol{\alpha},\mathbf{l},\mathbf{h})\sim\mathcal{N}(\mathbf{X}\boldsymbol{\alpha},\mathbf{S}_{y})$
with density function

$p(\mathbf{y}\,|\,\boldsymbol{\alpha},\mathbf{l},\mathbf{h})=(2\pi)^{-\frac{Tn}{2}}|\mathbf{S}_{y}|^{-\frac{1}{2}}\text{e}^{-\frac{1}{2}(\mathbf{y}-\mathbf{X}\boldsymbol{\alpha})^{\prime}\mathbf{S}_{y}^{-1}(\mathbf{y}-\mathbf{X}\boldsymbol{\alpha})},$

where
$\mathbf{S}_{y}=\mathbf{D}+(\mathbf{I}_{T}\otimes\mathbf{L})\mathbf{G}(\mathbf{I}_{T}\otimes\mathbf{L}^{\prime})$.
Finally, the likelihood function (marginal of both $\mathbf{f}$ and
$\mathbf{h}$ but conditional on the time-invariant model parameters) can be
evaluated via:

$\displaystyle
p(\mathbf{y}\,|\,\boldsymbol{\alpha},\mathbf{l},\boldsymbol{\mu},\boldsymbol{\phi},\boldsymbol{\sigma}^{2})=\int
p(\mathbf{y}\,|\,\boldsymbol{\alpha},\mathbf{l},\mathbf{h})\prod_{j=1}^{n+r}p(\mathbf{h}_{j,\boldsymbol{\cdot}}\,|\,\mu_{j},\phi_{j},\sigma_{j}^{2})\text{d}\mathbf{h},$

where $\mathbf{h}_{j,\boldsymbol{\cdot}}=(h_{j,1},\ldots,h_{j,T})^{\prime}$
and $p(\mathbf{h}_{j,\boldsymbol{\cdot}}\,|\,\mu_{j},\phi_{j},\sigma_{j}^{2})$
is a Gaussian density that is available analytically. To implement the
importance sampling approach described in Section 3.2, we consider the
parametric family

$\mathcal{F}=\left\\{\prod_{j=1}^{n+r}f_{\mathcal{N}}(\mathbf{h}_{j,\boldsymbol{\cdot}};\widehat{\mathbf{h}}_{j,\boldsymbol{\cdot}},\widehat{\mathbf{K}}_{\mathbf{h}_{j,\boldsymbol{\cdot}}}^{-1})\right\\}.$

The parameters of the importance sampling densities are chosen by solving the
maximization problem in (9). In particular, all the $T$-variate Gaussian
densities are obtained using the procedure described in Section 3.4.

To illustrate the methods, we compute the observed-data DIC of the 3
stochastic volatility models with $n=7$, and the results are reported in Table
5. In contrast to the marginal likelihood results, VAR-CSV is the most favored
model according to the observed-data DIC. This is due to a combination of
better model-fit (as shown by a lower posterior mean deviance value) and lower
model complexity (in terms of a lower effective number of parameters). Despite
being commonly used in the literature, rigorous justification of the DIC as a
model selection criterion has only been recently investigated (Li, Yu, and
Zeng, 2021). More research to better understand the theoretical and empirical
differences between the DIC and the marginal likelihood is therefore needed.

Table 5: The observed-data DIC, the posterior mean deviance $\overline{D(\boldsymbol{\theta})}$ and the effective number of parameters $p_{D}$ of VAR-CSV, VAR-SV and VAR-FSV for $n=7$. | VAR-CSV | VAR-SV | VAR-FSV
---|---|---|---
DIC | 3,367 | 3,440 | 3,517
$\overline{D(\boldsymbol{\theta})}$ | 3,181 | 3,195 | 3,239
$p_{D}$ | 186 | 245 | 278

###  6.3 A Stochastic Volatility Model with an Outlier Component

There has been a lot of recent work to develop VARs that can handle COVID-19
outliers; examples include Hartwig (2021), Ng (2021), Lenza and Primiceri
(2022) and Carriero, Clark, Marcellino, and Mertens (2022). In particular,
Carriero, Clark, Marcellino, and Mertens (2022) propose various stochastic
volatility specifications that have an explicit component to model infrequent
volatility outliers. Below we consider one of these stochastic volatility
models. It extends the Cholesky stochastic volatility by incorporating a
discrete mixture representation proposed in Stock and Watson (2016) for
handling outliers in unobserved component models. More specifically, with the
same reduced-form VAR in (1), the innovation $\boldsymbol{\varepsilon}_{t}$ is
now distributed as:

$\boldsymbol{\varepsilon}_{t}\sim\mathcal{N}(\mathbf{0},\boldsymbol{\Sigma}_{t}),\quad\boldsymbol{\Sigma}_{t}^{-1}=\mathbf{B}_{0}^{\prime}\mathbf{D}_{t}^{-1}\mathbf{B}_{0}o_{t}^{-2}.$
(15)

The new addition is the outlier parameter $o_{t}$, which is modeled to have a
mixture distribution that distinguishes between regular observations with
$o_{t}=1$ and outliers for which $o_{t}>2$. More precisely, with probability
$1-p_{o}$, $o_{t}=1$; and with probability $p_{o}$, $o_{t}$ follows a uniform
distribution on $(2,20)$: $o_{t}\sim\mathcal{U}(2,20)$. The outlier
probability $p_{o}$ is assumed to have a beta prior, where its hyperparameters
are calibrated so that the mean outlier frequency is once every 4 years in
quarterly data. We refer to this model as VAR-SVO.

Following Stock and Watson (2016) and Carriero, Clark, Marcellino, and Mertens
(2022), the distribution of $o_{t}$ is discretized using a grid. Consequently,
$o_{t}$ follows a discrete distribution that can be easily handled. We refer
the readers to Carriero, Clark, Marcellino, and Mertens (2022) for estimation
details. We fit VAR-SVO using the datasets ($n=7,15,30$) in Section 5 and
compute the corresponding marginal likelihood values. The results are reported
in Table 6. For convenience, we also reproduce the marginal likelihood values
of other models (the value of the best VAR-FSV for each dimension is
reported).

Table 6: Log marginal likelihood estimates (numerical standard errors) of a standard homoskedastic VAR, VAR-CSV, VAR-SV, VAR-FSV and VAR-SVO. | VAR | VAR-CSV | VAR-SV | VAR-FSV | VAR-SVO
---|---|---|---|---|---
$n=7$ |  $-$2,583.4 |  $-$2,410 |  $-$2,315 |  $-$2,318 |  $-$2,306
| – | (0.1) | (0.3) | (0.4) | (0.6)
$n=15$ |  $-$6,918.8 |  $-$6,618 |  $-$6,443 |  $-$6,454 |  $-$6,386
| – | (0.1) | (0.7) | (0.8) | (0.8)
$n=30$ |  $-$12,783.4 |  $-$12,024 |  $-$11,382 |  $-$11,567 |  $-$11,255
| – | (0.1) | (1.2) | (1.8) | (1.9)

For all model dimensions, VAR-SVO is the most favored specification by the
data. These results suggest that there are infrequent, extreme observations
even before the COVID-19 pandemic. Moreover, they also highlight that an
explicit outlier component is useful for handling these apparent outliers,
even in the presence of time-varying volatility.

##  7 Concluding Remarks and Future Research

As large Bayesian VARs are now widely used in empirical applications, choosing
the most suitable stochastic volatility specification and shrinkage priors for
a particular dataset has become an increasingly pertinent problem. We took a
first step to address this issue by developing Bayesian model comparison
methods to select among a variety of stochastic volatility specifications and
shrinkage priors in the context of large VARs. We demonstrated via a series of
Monte Carlo experiments that the proposed method worked well—particularly it
could be used to discriminate VARs with different stochastic volatility
specifications.

Using US datasets of different dimensions, we showed that the data strongly
preferred the Cholesky stochastic volatility, whereas the factor stochastic
volatility was also competitive. This finding thus suggests that more future
research on factor stochastic volatility models would be fruitful, given that
they are not as commonly used in empirical macroeconomics. Our results also
confirmed the vital role of flexible shrinkage priors: both cross-variable
shrinkage and a data-based approach to determine the overall shrinkage
strength were empirically important.

In future work, it would be useful to extend the proposed estimators to
compare large time-varying parameter VARs. Existing evidence seems to suggest
that in a large VAR, only a few of the coefficients are time-varying. Such a
model comparison method would be helpful for comparing different dynamic
shrinkage priors, and thus provide useful guidelines for empirical
researchers.

## Appendix A: Estimation Details

In this appendix we provide the details of the priors and the estimation of
the Bayesian VARs with three different stochastic volatility specifications.
We also discuss the details of the adaptive importance sampling approach for
estimating the marginal likelihoods for these models.

### Common Stochastic Volatility

We first outline the estimation of the common stochastic volatility model
given in (1)–(3). To that end, let
$\mathbf{x}_{t}^{\prime}=(1,\mathbf{y}_{t-1}^{\prime},\ldots,\mathbf{y}_{t-p}^{\prime})$
be a $1\times k$ vector of an intercept and lags with $k=1+np$. Then, stacking
the observations over $t=1,\ldots,T$, we have

$\mathbf{Y}=\mathbf{X}\mathbf{A}+\boldsymbol{\varepsilon},$

where
$\mathbf{A}=(\mathbf{a}_{0},\mathbf{A}_{1},\ldots,\mathbf{A}_{p})^{\prime}$ is
$k\times n$, and the matrices $\mathbf{Y}$, $\mathbf{X}$ and
$\boldsymbol{\varepsilon}$ are respectively of dimensions $T\times n$,
$T\times k$ and $T\times n$. Under the common stochastic volatility model, the
innovations are distributed as
$\text{vec}(\boldsymbol{\varepsilon})\sim\mathcal{N}(\mathbf{0},\boldsymbol{\Sigma}\otimes\mathbf{D}_{\mathbf{h}})$,
where
$\mathbf{D}_{\mathbf{h}}=\text{diag}(\text{e}^{h_{1}},\ldots,\text{e}^{h_{T}})$.
The log-volatility $h_{t}$ evolves as a stationary AR(1) process given in (3),
which for convenience we reproduce below:

$h_{t}=\phi h_{t-1}+u_{t}^{h},\quad u_{t}^{h}\sim\mathcal{N}(0,\sigma^{2}),$

for $t=2,\ldots,T$, where the process is initialized as
$h_{1}\sim\mathcal{N}(0,\sigma^{2}/(1-\phi^{2}))$.

Next, we describe the priors on the model parameters
$\mathbf{A},\boldsymbol{\Sigma},\kappa,\phi$ and $\sigma^{2}$. First, we
consider a natural conjugate prior on
$(\mathbf{A},\boldsymbol{\Sigma}\,|\,\kappa)$:

$\boldsymbol{\Sigma}\sim\mathcal{IW}(\nu_{0},\mathbf{S}_{0}),\quad(\text{vec}(\mathbf{A})\,|\,\boldsymbol{\Sigma},\kappa)\sim\mathcal{N}(\text{vec}(\mathbf{A}_{0}),\boldsymbol{\Sigma}\otimes\mathbf{V}_{\mathbf{A}}),$

where the prior hyperparameters $\text{vec}(\mathbf{A}_{0})$ and
$\mathbf{V}_{\mathbf{A}}$ are chosen in the spirit of the Minnesota prior.
More specifically, we set $\mathbf{A}_{0}=\mathbf{0}$ and the covariance
matrix $\mathbf{V}_{\mathbf{A}}$ is assumed to be diagonal with diagonal
elements $v_{\mathbf{A},ii}=\kappa/(l^{2}\widehat{s}_{r})$ for a coefficient
associated to the $l$-th lag of variable $r$ and $v_{\mathbf{A},ii}=100$ for
an intercept, where $\widehat{s}_{r}$ is the sample variance of an AR(4) model
for the variable $r$. Note that here a single hyperparameter $\kappa$ controls
the overall shrinkage strength and this prior does not distinguish ‘own’
versus ‘other’ lags. Here we treat $\kappa$ to be an unknown parameter with a
hierarchical gamma prior: $\kappa\sim\mathcal{G}(c_{1},c_{2})$. Finally, for
the prior distributions of $\phi$ and $\sigma^{2}$, they are respectively
truncated normal and inverse-gamma:
$\phi\sim\mathcal{N}(\phi_{0},V_{\phi})1(|\phi|<1)$ and
$\sigma^{2}\sim\mathcal{IG}(\nu_{\sigma^{2}},S_{\sigma^{2}})$, where
$1(\cdot)$ denotes the indicator function.

With the priors on the model parameters specified above, one can obtain
posterior draws by sequentially sampling from:

  1. 1.

$p(\mathbf{A},\boldsymbol{\Sigma}\,|\,\mathbf{Y},\mathbf{h},\phi,\sigma^{2},\kappa)=p(\mathbf{A},\boldsymbol{\Sigma}\,|\,\mathbf{Y},\mathbf{h},\kappa)$;

  2. 2.

$p(\mathbf{h}\,|\,\mathbf{Y},\mathbf{A},\boldsymbol{\Sigma},\phi,\sigma^{2},\kappa)=p(\mathbf{h}\,|\,\mathbf{Y},\mathbf{A},\boldsymbol{\Sigma},\phi,\sigma^{2})$;

  3. 3.

$p(\phi\,|\,\mathbf{Y},\mathbf{A},\boldsymbol{\Sigma},\mathbf{h},\sigma^{2},\kappa)=p(\phi\,|\,\mathbf{h},\sigma^{2})$;

  4. 4.

$p(\sigma^{2}\,|\,\mathbf{Y},\mathbf{A},\boldsymbol{\Sigma},\mathbf{h},\phi,\kappa)=p(\sigma^{2}\,|\,\mathbf{h},\phi)$;

  5. 5.

$p(\kappa\,|\,\mathbf{Y},\mathbf{A},\boldsymbol{\Sigma},\mathbf{h},\phi,\sigma^{2})=p(\kappa\,|\,\mathbf{A},\boldsymbol{\Sigma})$.

Step 1: we use the results proved in Chan (2020) that
$(\mathbf{A},\boldsymbol{\Sigma}\,|\,\mathbf{Y},\mathbf{h},\kappa)$ has a
normal-inverse-Wishart distribution. More specifically, let

$\displaystyle\mathbf{K}_{\mathbf{A}}$
$\displaystyle=\mathbf{V}_{\mathbf{A}}^{-1}+\mathbf{X}^{\prime}\mathbf{D}^{-1}_{\mathbf{h}}\mathbf{X},$
$\displaystyle\widehat{\mathbf{A}}$
$\displaystyle=\mathbf{K}_{\mathbf{A}}^{-1}(\mathbf{V}_{\mathbf{A}}^{-1}\mathbf{A}_{0}+\mathbf{X}^{\prime}\mathbf{D}^{-1}_{\mathbf{h}}\mathbf{Y}),$
$\displaystyle\widehat{\mathbf{S}}$
$\displaystyle=\mathbf{S}_{0}+\mathbf{A}_{0}^{\prime}\mathbf{V}_{\mathbf{A}}^{-1}\mathbf{A}_{0}+\mathbf{Y}^{\prime}\mathbf{D}^{-1}_{\mathbf{h}}\mathbf{Y}-\widehat{\mathbf{A}}^{\prime}\mathbf{K}_{\mathbf{A}}\widehat{\mathbf{A}}.$

Then $(\mathbf{A},\boldsymbol{\Sigma}\,|\,\mathbf{Y},\mathbf{h},\kappa)$ has a
normal-inverse-Wishart distribution with parameters $\nu_{0}+T$,
$\widehat{\mathbf{S}}$, $\widehat{\mathbf{A}}$ and
$\mathbf{K}_{\mathbf{A}}^{-1}$. We can sample
$(\mathbf{A},\boldsymbol{\Sigma}\,|\,\mathbf{Y},\mathbf{h},\kappa)$ in two
steps. First, sample $\boldsymbol{\Sigma}$ marginally from the inverse-Wishart
distribution
$(\boldsymbol{\Sigma}\,|\,\mathbf{Y},\mathbf{h})\sim\mathcal{IW}(\nu_{0}+T,\widehat{\mathbf{S}})$.
Then, given the sampled $\boldsymbol{\Sigma}$, obtain $\mathbf{A}$ from the
normal distribution:

$(\text{vec}(\mathbf{A})\,|\,\mathbf{Y},\boldsymbol{\Sigma},\mathbf{h})\sim\mathcal{N}(\text{vec}(\widehat{\mathbf{A}}),\boldsymbol{\Sigma}\otimes\mathbf{K}_{\mathbf{A}}^{-1}).$

We note that one can sample from this normal distribution efficiently without
explicitly computing the inverse $\mathbf{K}_{\mathbf{A}}^{-1}$; we refer the
readers to Chan (2020) for computational details.

Step 2: note that

$p(\mathbf{h}\,|\,\mathbf{Y},\mathbf{A},\boldsymbol{\Sigma},\phi,\sigma^{2})\propto
p(\mathbf{h}\,|\,\phi,\sigma^{2})\prod_{t=1}^{T}p(\mathbf{y}_{t}\,|\,\mathbf{A},\boldsymbol{\Sigma},h_{t}),$

where $p(\mathbf{h}\,|\,\phi,\sigma^{2})$ is a Gaussian density implied by the
state equation (3). The log conditional likelihood
$p(\mathbf{y}_{t}\,|\,\mathbf{A},\boldsymbol{\Sigma},h_{t})$ has the following
explicit expression:

$\log
p(\mathbf{y}_{t}\,|\,\mathbf{A},\boldsymbol{\Sigma},h_{t})=c_{t}-\frac{n}{2}h_{t}-\frac{1}{2}\text{e}^{-h_{t}}\boldsymbol{\varepsilon}_{t}^{\prime}\boldsymbol{\Sigma}^{-1}\boldsymbol{\varepsilon}_{t},$

where $c_{t}$ is a constant independent of $h_{t}$. It is easy to check that

$\displaystyle\frac{\partial}{\partial h_{t}}\log
p(\mathbf{y}_{t}\,|\,\mathbf{A},\boldsymbol{\Sigma},h_{t})$
$\displaystyle=-\frac{n}{2}+\frac{1}{2}\text{e}^{-h_{t}}\boldsymbol{\varepsilon}_{t}^{\prime}\boldsymbol{\Sigma}^{-1}\boldsymbol{\varepsilon}_{t},$
$\displaystyle\frac{\partial^{2}}{\partial h_{t}^{2}}\log
p(\mathbf{y}_{t}\,|\,\mathbf{A},\boldsymbol{\Sigma},h_{t})$
$\displaystyle=-\frac{1}{2}\text{e}^{-h_{t}}\boldsymbol{\varepsilon}_{t}^{\prime}\boldsymbol{\Sigma}^{-1}\boldsymbol{\varepsilon}_{t}<0.$

Given the above first and second derivatives of the log conditional
likelihood, one can use the Newton-Raphson algorithm to obtain the mode of
$\log
p(\mathbf{h}\,|\,\mathbf{Y},\mathbf{A},\boldsymbol{\Sigma},\phi,\sigma^{2})$
and compute the negative Hessian evaluated at the mode, which are denoted as
$\widehat{\mathbf{h}}$ and $\mathbf{K}_{\mathbf{h}}$, respectively. Since the
Hessian is negative definite everywhere, $\mathbf{K}_{\mathbf{h}}$ is positive
definite. Next, using
$\mathcal{N}(\widehat{\mathbf{h}},\mathbf{K}_{\mathbf{h}}^{-1})$ as a proposal
distribution, one can sample $\mathbf{h}$ directly using an acceptance-
rejection Metropolis-Hastings step. We refer the readers to Chan (2020) for
details.

Step 3: the conditional density of $\phi$ is of the form

$p(\phi\,|\,\mathbf{h},\sigma^{2})\propto
p(\phi)g(\phi)\text{e}^{-\frac{1}{2\sigma^{2}}\sum_{t=2}^{T}(h_{t}-\phi
h_{t-1})^{2}},$

where
$g(\phi)=(1-\phi^{2})^{\frac{1}{2}}\text{e}^{-\frac{1}{2\sigma^{2}}(1-\phi^{2})h_{1}^{2}}$
and $p(\phi)$ is the truncated normal prior. This density is nonstandard, but
we can sample $\phi$ via an independence-chain Metropolis-Hastings step using
the proposal distribution
$\mathcal{N}(\widehat{\phi},K_{\phi}^{-1})1(|\phi|<1)$, where

$K_{\phi}=V_{\phi}^{-1}+\frac{1}{\sigma^{2}}\sum_{t=2}^{T}h_{t-1}^{2}\quad\widehat{\phi}=K_{\phi}^{-1}\left(V_{\phi}^{-1}\phi_{0}+\frac{1}{\sigma^{2}}\sum_{t=2}^{T}h_{t}h_{t-1}\right).$

Then, given the current draw $\phi$, a proposal $\phi^{*}$ is accepted with
probability $\min(1,g(\phi^{*})/g(\phi))$; otherwise the Markov chain stays at
the current state $\phi$.

Step 4: $\sigma^{2}$ can be sampled easily as it has the following inverse-
gamma distribution:

$(\sigma^{2}\,|\,\mathbf{h},\phi)\sim\mathcal{IG}(\nu_{\sigma^{2}}+T/2,\widetilde{S}_{\sigma^{2}}),$

where
$\widetilde{S}_{\sigma^{2}}=S_{\sigma^{2}}+\left[(1-\phi^{2})h_{1}^{2}+\sum_{t=2}^{T}(h_{t}-\phi
h_{t-1})^{2}\right]/2$.

Step 5: first note that $\kappa$ only appears in its gamma prior
$\kappa\sim\mathcal{G}(c_{1},c_{2})$ and $\mathbf{V}_{\mathbf{A}}$, the prior
covariance matrix of $\mathbf{A}$. Recall that $\mathbf{V}_{\mathbf{A}}$ is a
$k\times k$ diagonal matrix in which the first element—corresponding to the
prior variance of the intercept—does not involve $\kappa$. More explicitly,
for $i=2,\ldots,k$, write the $i$-th diagonal element of
$\mathbf{V}_{\mathbf{A}}$ as $v_{\mathbf{A},ii}=\kappa C_{i}$ for some
constant $C_{i}$. Then, we have

$\displaystyle p(\kappa\,|\,\mathbf{A},\boldsymbol{\Sigma})$
$\displaystyle\propto\kappa^{c_{1}-1}\text{e}^{-c_{2}\kappa}\times|\mathbf{V}_{\mathbf{A}}|^{-\frac{n}{2}}\text{e}^{-\frac{1}{2}\text{tr}\left(\boldsymbol{\Sigma}^{-1}(\mathbf{A}-\mathbf{A}_{0})^{\prime}\mathbf{V}_{\mathbf{A}}^{-1}(\mathbf{A}-\mathbf{A}_{0})\right)}$
$\displaystyle\propto\kappa^{c_{1}-\frac{(k-1)n}{2}-1}\text{e}^{-c_{2}\kappa}\text{e}^{-\frac{1}{2}\text{tr}\left(\mathbf{V}_{\mathbf{A}}^{-1}(\mathbf{A}-\mathbf{A}_{0})\boldsymbol{\Sigma}^{-1}(\mathbf{A}-\mathbf{A}_{0})^{\prime}\right)}$
$\displaystyle\propto\kappa^{c_{1}-\frac{(k-1)n}{2}-1}\text{e}^{-\frac{1}{2}\left(2c_{2}\kappa+\kappa^{-1}\sum_{i=2}^{k}Q_{i}/C_{i}\right)},$

where $Q_{i}$ is the $i$-th diagonal element of
$\mathbf{Q}=(\mathbf{A}-\mathbf{A}_{0})\boldsymbol{\Sigma}^{-1}(\mathbf{A}-\mathbf{A}_{0})^{\prime}$.
Note that this is the kernel of the generalized inverse Gaussian distribution
$\mathcal{GIG}\left(c_{1}-\frac{(k-1)n}{2},2c_{2},\sum_{i=2}^{k}Q_{i}/C_{i}\right)$.
Draws from the generalized inverse Gaussian distribution can be obtained using
the algorithm in Devroye (2014).

Next, we derive the expression of $p(\mathbf{Y}\,|\,\mathbf{h},\kappa)$ given
in (10). First let $k_{1}$ denote the normalizing constant of the normal-
inverse-Wishart prior:
$k_{1}=(2\pi)^{-\frac{nk}{2}}2^{-\frac{n\nu_{0}}{2}}|\mathbf{V}_{\mathbf{A}}|^{-\frac{n}{2}}\Gamma_{n}(\nu_{0}/2)^{-1}|\mathbf{S}_{0}|^{\frac{\nu_{0}}{2}}$.
Then, by direct computation:

$\displaystyle p(\mathbf{Y}\,|\,\mathbf{h},\kappa)$ $\displaystyle=\int
p(\mathbf{Y}\,|\,\mathbf{A},\boldsymbol{\Sigma},\mathbf{h})p(\mathbf{A},\boldsymbol{\Sigma}\,|\,\kappa)\text{d}(\mathbf{A},\boldsymbol{\Sigma})$
$\displaystyle=\int(2\pi)^{-\frac{Tn}{2}}|\boldsymbol{\Sigma}|^{-\frac{T}{2}}\text{e}^{-\frac{n}{2}\mathbf{1}_{T}^{\prime}\mathbf{h}}\text{e}^{-\frac{1}{2}\text{tr}\left(\boldsymbol{\Sigma}^{-1}(\mathbf{Y}-\mathbf{X}\mathbf{A})^{\prime}\mathbf{D}^{-1}_{\mathbf{h}}(\mathbf{Y}-\mathbf{X}\mathbf{A})\right)}$
$\displaystyle\quad\times
k_{1}|\boldsymbol{\Sigma}|^{-\frac{\nu_{0}+n+k+1}{2}}\text{e}^{-\frac{1}{2}\text{tr}\left(\boldsymbol{\Sigma}^{-1}\mathbf{S}_{0}\right)}\text{e}^{-\frac{1}{2}\text{tr}\left(\boldsymbol{\Sigma}^{-1}(\mathbf{A}-\mathbf{A}_{0})^{\prime}\mathbf{V}_{\mathbf{A}}^{-1}(\mathbf{A}-\mathbf{A}_{0})\right)}\text{d}(\mathbf{A},\boldsymbol{\Sigma})$
$\displaystyle=k_{1}(2\pi)^{-\frac{Tn}{2}}\text{e}^{-\frac{n}{2}\mathbf{1}_{T}^{\prime}\mathbf{h}}\int|\boldsymbol{\Sigma}|^{-\frac{\nu_{0}+T+n+k+1}{2}}\text{e}^{-\frac{1}{2}\text{tr}\left(\boldsymbol{\Sigma}^{-1}\widehat{\mathbf{S}}\right)}\text{e}^{-\frac{1}{2}\text{tr}\left(\boldsymbol{\Sigma}^{-1}(\mathbf{A}-\widehat{\mathbf{A}})^{\prime}\mathbf{K}_{\mathbf{A}}^{-1}(\mathbf{A}-\widehat{\mathbf{A}})\right)}\text{d}(\mathbf{A},\boldsymbol{\Sigma})$
$\displaystyle=\pi^{-\frac{Tn}{2}}\text{e}^{-\frac{n}{2}\mathbf{1}_{T}^{\prime}\mathbf{h}}|\mathbf{V}_{\mathbf{A}}|^{-\frac{n}{2}}|\mathbf{K}_{\mathbf{A}}|^{-\frac{n}{2}}\frac{\Gamma_{n}\left(\frac{\nu_{0}+T}{2}\right)|\mathbf{S}_{0}|^{\frac{\nu_{0}}{2}}}{\Gamma_{n}\left(\frac{\nu_{0}}{2}\right)|\widehat{\mathbf{S}}|^{\frac{\nu_{0}+T}{2}}},$

where the shrinkage hyperparameter $\kappa$ appears in
$\mathbf{V}_{\mathbf{A}}$, and the last equality holds because

$\displaystyle\int|\boldsymbol{\Sigma}|^{-\frac{\nu_{0}+T+n+k+1}{2}}$
$\displaystyle\text{e}^{-\frac{1}{2}\text{tr}\left(\boldsymbol{\Sigma}^{-1}\widehat{\mathbf{S}}\right)}\text{e}^{-\frac{1}{2}\text{tr}\left(\boldsymbol{\Sigma}^{-1}(\mathbf{A}-\widehat{\mathbf{A}})^{\prime}\mathbf{K}_{\mathbf{A}}^{-1}(\mathbf{A}-\widehat{\mathbf{A}})\right)}\text{d}(\mathbf{A},\boldsymbol{\Sigma})$
$\displaystyle=(2\pi)^{\frac{nk}{2}}2^{\frac{n(\nu_{0}+T)}{2}}|\mathbf{K}_{\mathbf{A}}^{-1}|^{\frac{n}{2}}\Gamma_{n}\left(\frac{\nu_{0}+T}{2}\right)|\widehat{\mathbf{S}}|^{-\frac{\nu_{0}+T}{2}}.$

Finally, the marginal density $p(\mathbf{h}\,|\,\phi)$ has the following
analytical expression:

$\displaystyle p(\mathbf{h}\,|\,\phi)$ $\displaystyle=\int
p(\mathbf{h}\,|\,\phi,\sigma^{2})p(\sigma^{2})\text{d}\sigma^{2}$
$\displaystyle=(2\pi)^{-\frac{T}{2}}|1-\phi^{2}|^{\frac{1}{2}}S_{\sigma^{2}}^{\nu_{\sigma^{2}}}\Gamma(\nu_{\sigma^{2}})^{-1}\int(\sigma^{2})^{-(\nu_{\sigma^{2}}+\frac{T}{2}+1)}\text{e}^{-\frac{1}{\sigma^{2}}\widetilde{S}_{\sigma^{2}}}\text{d}\sigma^{2}$
$\displaystyle=(2\pi)^{-\frac{T}{2}}|1-\phi^{2}|^{\frac{1}{2}}\frac{\Gamma\left(\nu_{\sigma^{2}}+\frac{T}{2}\right)S_{\sigma^{2}}^{\nu_{\sigma^{2}}}}{\Gamma(\nu_{\sigma^{2}})\widetilde{S}_{\sigma^{2}}^{\nu_{\sigma^{2}}+\frac{T}{2}}},$

where $\Gamma(\cdot)$ is the gamma function and $\widetilde{S}_{\sigma^{2}}$
is defined in Step 4 of the MCMC algorithm.

### Cholesky Stochastic Volatility

Next, we outline the estimation of the VAR-SV model:

$\mathbf{y}_{t}=\mathbf{a}_{0}+\mathbf{A}_{1}\mathbf{y}_{t-1}+\cdots+\mathbf{A}_{p}\mathbf{y}_{t-p}+\boldsymbol{\varepsilon}_{t},\quad\boldsymbol{\varepsilon}_{t}\sim\mathcal{N}(\mathbf{0},\boldsymbol{\Sigma}_{t}),$
(16)

where
$\boldsymbol{\Sigma}_{t}^{-1}=\mathbf{B}_{0}^{\prime}\mathbf{D}_{t}^{-1}\mathbf{B}_{0}$,
$\mathbf{D}_{t}=\operatorname{diag}(\text{e}^{h_{1,t}},\ldots,\text{e}^{h_{n,t}})$
and $\mathbf{B}_{0}$ is an $n\times n$ lower triangular matrix with ones on
the diagonal. Each element of
$\mathbf{h}_{t}=(h_{1,t},\ldots,h_{n,t})^{\prime}$ follows an autoregressive
process:

$h_{i,t}=\mu_{i}+\phi_{i}(h_{i,t-1}-\mu_{i})+u_{i,t}^{h},\quad
u_{i,t}^{h}\sim\mathcal{N}(0,\sigma_{i}^{2})$

for $t=2,\ldots,T$, and the initial condition is specified as
$h_{i,1}\sim\mathcal{N}(\mu_{i},\sigma^{2}/(1-\phi_{i}^{2}))$.

We aim to estimate the system equation by equation, which substantially speeds
up the computation time. Recall that $\boldsymbol{\alpha}_{i}$ denotes the
$k\times 1$ vector that consists of the intercept and VAR coefficients in the
$i$-th equation, and $\boldsymbol{\beta}_{i}$ represents the $(i-1)\times 1$
vector of free elements in the $i$-th row of the impact matrix
$\mathbf{B}_{0}$. Then, the parameters for the $i$-th equation are
$\boldsymbol{\alpha}_{i},\boldsymbol{\beta}_{i},\mu_{i},\phi_{i},$ and
$\sigma_{i}^{2}$. We assume the following independent prior distributions on
the parameters:

$\begin{split}(\boldsymbol{\alpha}_{i}\,|\,\boldsymbol{\kappa})&\sim\mathcal{N}(\boldsymbol{\alpha}_{0,i},\mathbf{V}_{\boldsymbol{\alpha}_{i}}),\;(\boldsymbol{\beta}_{i}\,|\,\boldsymbol{\kappa})\sim\mathcal{N}(\boldsymbol{\beta}_{0,i},\mathbf{V}_{\boldsymbol{\beta}_{i}}),\\\
\mu_{i}&\sim\mathcal{N}(\mu_{0,i},V_{\mu_{i}}),\;\phi_{i}\sim\mathcal{N}(\phi_{0,i},V_{\phi_{i}})1(|\phi_{i}|<1),\;\sigma_{i}^{2}\sim\mathcal{IG}(\nu_{i},S_{i}),\end{split}$
(17)

where $\boldsymbol{\kappa}$ is a vector of hyperparameters that is described
in more detail below. The prior mean vector $\boldsymbol{\alpha}_{0,i}$ and
the prior covariance matrix $\mathbf{V}_{\boldsymbol{\alpha}_{i}}$ are
selected to mimic the Minnesota prior. More specifically, we set
$\boldsymbol{\alpha}_{0,i}=\mathbf{0}$ to shrink the VAR coefficients to zero.
For $\mathbf{V}_{\boldsymbol{\alpha}_{i}}$, we assume it to be diagonal with
the $k$-th diagonal element $V_{\boldsymbol{\alpha}_{i},kk}$ set to be:

$V_{\boldsymbol{\alpha}_{i},kk}=\left\\{\begin{array}[]{ll}\frac{\kappa_{1}}{l^{2}},&\text{for
the coefficient on the $l$-th lag of variable }i,\\\
\frac{\kappa_{2}s_{i}^{2}}{l^{2}s_{j}^{2}},&\text{for the coefficient on the
$l$-th lag of variable }j,j\neq i,\\\ 100s_{i}^{2},&\text{for the
intercept},\\\ \end{array}\right.$

where $s_{r}^{2}$ denotes the sample variance of the residuals from an AR(4)
model for the variable $r$ for $r=1,\ldots,n$. In addition, we set the prior
mean vector $\boldsymbol{\beta}_{0,i}$ to be zero to shrink the impact matrix
to the identity matrix. In addition, the prior covariance matrix
$\mathbf{V}_{\boldsymbol{\beta}_{i}}$ is assumed to be diagonal, where the
$j$-th diagonal element is set to be $\kappa_{3}s_{i}^{2}/s_{j}^{2}$. Finally,
we treat the shrinkage hyperparameters
$\boldsymbol{\kappa}=(\kappa_{1},\kappa_{2},\kappa_{3})^{\prime}$ as unknown
parameters to be estimated with hierarchical gamma priors
$\kappa_{i}\sim\mathcal{G}(c_{j,1},c_{j,2}),j=1,2,3.$

Let $\mathbf{y}_{i,\boldsymbol{\cdot}}=(y_{i,1},\ldots,y_{i,T})^{\prime}$
denote the vector of observed values for the $i$-th variable for
$i=1,\ldots,n$. We similarly define
$\mathbf{h}_{i,\boldsymbol{\cdot}}=(h_{i,1},\ldots,h_{i,T})^{\prime}$ . Next,
stack
$\mathbf{y}=(\mathbf{y}_{1,\boldsymbol{\cdot}}^{\prime},\ldots,\mathbf{y}_{n,\boldsymbol{\cdot}}^{\prime})^{\prime},\mathbf{h}=(\mathbf{h}_{1,\boldsymbol{\cdot}}^{\prime},\ldots,\mathbf{h}_{n,\boldsymbol{\cdot}}^{\prime})^{\prime}$
and
$\boldsymbol{\alpha}=(\boldsymbol{\alpha}_{1}^{\prime},\ldots,\boldsymbol{\alpha}_{n}^{\prime})^{\prime}$;
similarly define $\boldsymbol{\beta},\boldsymbol{\mu},\boldsymbol{\phi}$ and
$\boldsymbol{\sigma}^{2}$. Then, posterior draws can be obtained by sampling
sequentially from:

  1. 1.

$p(\boldsymbol{\alpha}\,|\,\mathbf{y},\boldsymbol{\beta},\mathbf{h},\boldsymbol{\mu},\boldsymbol{\phi},\boldsymbol{\sigma}^{2},\boldsymbol{\kappa})$;

  2. 2.

$p(\boldsymbol{\beta}\,|\,\mathbf{y},\boldsymbol{\alpha},\mathbf{h},\boldsymbol{\mu},\boldsymbol{\phi},\boldsymbol{\sigma}^{2},\boldsymbol{\kappa})=\prod_{i=2}^{n}p(\boldsymbol{\beta}_{i}\,|\,\mathbf{y},\boldsymbol{\alpha},\mathbf{h}_{i,\boldsymbol{\cdot}},\boldsymbol{\kappa})$;

  3. 3.

$p(\mathbf{h}\,|\,\mathbf{y},\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\mu},\boldsymbol{\phi},\boldsymbol{\sigma}^{2},\boldsymbol{\kappa})=\prod_{i=1}^{n}p(\mathbf{h}_{i,\boldsymbol{\cdot}}\,|\,\mathbf{y},\boldsymbol{\alpha},\boldsymbol{\beta},\mu_{i},\phi_{i},\sigma^{2}_{i})$;

  4. 4.

$p(\boldsymbol{\mu}\,|\,\mathbf{y},\boldsymbol{\alpha},\boldsymbol{\beta},\mathbf{h},\boldsymbol{\phi},\boldsymbol{\sigma}^{2},\boldsymbol{\kappa})=\prod_{i=1}^{n}p(\mu_{i}\,|\,\mathbf{h}_{i,\boldsymbol{\cdot}},\phi_{i},\sigma^{2}_{i})$;

  5. 5.

$p(\boldsymbol{\phi}\,|\,\mathbf{y},\boldsymbol{\alpha},\boldsymbol{\beta},\mathbf{h},\boldsymbol{\mu},\boldsymbol{\sigma}^{2},\boldsymbol{\kappa})=\prod_{i=1}^{n}p(\phi_{i}\,|\,\mathbf{h}_{i,\boldsymbol{\cdot}},\mu_{i},\sigma^{2}_{i})$;

  6. 6.

$p(\boldsymbol{\sigma}^{2}\,|\,\mathbf{y},\boldsymbol{\alpha},\boldsymbol{\beta},\mathbf{h},\boldsymbol{\mu},\boldsymbol{\phi},\boldsymbol{\kappa})=\prod_{i=1}^{n}p(\sigma_{i}^{2}\,|\,\mathbf{h}_{i,\boldsymbol{\cdot}},\mu_{i},\phi_{i})$;

  7. 7.

$p(\boldsymbol{\kappa}\,|\,\mathbf{y},\boldsymbol{\alpha},\boldsymbol{\beta},\mathbf{h},\boldsymbol{\mu},\boldsymbol{\phi},\boldsymbol{\sigma}^{2})=p(\boldsymbol{\kappa}\,|\,\boldsymbol{\alpha},\boldsymbol{\beta}).$

Step 1: we sample the reduced-form VAR coefficients equation by equation using
the triangular algorithm proposed in Carriero, Chan, Clark, and Marcellino
(2022). To that end, define $\mathbf{A}_{i=0}$ to be a $k\times n$ matrix that
has exactly the same elements as
$\mathbf{A}=(\mathbf{a}_{0},\mathbf{A}_{1},\ldots,\mathbf{A}_{p})^{\prime}$
except for the $i$-th column, which is set to be zero, i.e.,
$\mathbf{A}_{i=0}=(\boldsymbol{\alpha}_{1},\ldots,\boldsymbol{\alpha}_{i-1},\mathbf{0},\boldsymbol{\alpha}_{i+1},\ldots,\boldsymbol{\alpha}_{n}).$
Then, we can rewrite (16) as

$\mathbf{B}_{0}(\mathbf{y}_{t}-\mathbf{A}_{i=0}^{\prime}\mathbf{x}_{t})=\mathbf{B}_{0,1:n,i}\boldsymbol{\alpha}_{i}^{\prime}\mathbf{x}_{t}+\widetilde{\boldsymbol{\varepsilon}}_{t},\quad\widetilde{\boldsymbol{\varepsilon}}_{t}\equiv\mathbf{B}_{0}\boldsymbol{\varepsilon}_{t}\sim\mathcal{N}(\mathbf{0},\mathbf{D}_{t}),$

where
$\mathbf{x}_{t}=(1,\mathbf{y}_{t-1}^{\prime},\ldots,\mathbf{y}_{t-p}^{\prime})^{\prime}$
and $\mathbf{B}_{0,1:n,i}$ is the $i$-th column of $\mathbf{B}_{0}$.
Transposing and stacking the above equation over $t=1,\ldots,T,$, we can write
the system more compactly as

$(\mathbf{Y}-\mathbf{X}\mathbf{A}_{i=0})\mathbf{B}_{0}^{\prime}=\mathbf{X}\boldsymbol{\alpha}_{i}\mathbf{B}_{0,1:n,i}^{\prime}+\widetilde{\boldsymbol{\varepsilon}},$

where the matrices $\mathbf{Y}$, $\mathbf{X}$ and
$\widetilde{\boldsymbol{\varepsilon}}$ are respectively of dimensions $T\times
n$, $T\times k$ and $T\times n$. Vectoring the above system, we obtain

$\mathbf{z}^{i}=\mathbf{W}^{i}\boldsymbol{\alpha}_{i}+\text{vec}(\widetilde{\boldsymbol{\varepsilon}}),\quad\text{vec}(\widetilde{\boldsymbol{\varepsilon}})\sim\mathcal{N}(\mathbf{0},\mathbf{D}),$

where
$\mathbf{z}^{i}_{t}=\text{vec}((\mathbf{Y}-\mathbf{X}\mathbf{A}_{i=0})\mathbf{B}_{0}^{\prime})$,
$\mathbf{W}^{i}=\mathbf{B}_{0,1:n,i}\otimes\mathbf{X}$ and
$\mathbf{D}=\text{diag}(\mathbf{D}_{1,\boldsymbol{\cdot}},\ldots,\mathbf{D}_{n,\boldsymbol{\cdot}})$
with
$\mathbf{D}_{i,\boldsymbol{\cdot}}=\text{diag}(\text{e}^{h_{i,1}},\ldots,\text{e}^{h_{i,T}}).$
Given the conditionally Gaussian prior on $\boldsymbol{\alpha}_{i}$ in (17)
(conditional on $\boldsymbol{\kappa}$), it follows from standard linear
regression results that

$(\boldsymbol{\alpha}_{i}\,|\,\mathbf{y},\boldsymbol{\beta},\boldsymbol{\alpha}_{-i},\mathbf{h},\boldsymbol{\mu},\boldsymbol{\phi},\boldsymbol{\sigma}^{2},\boldsymbol{\kappa})\sim\mathcal{N}(\widehat{\boldsymbol{\alpha}}_{i},\mathbf{K}_{\boldsymbol{\alpha}_{i}}^{-1}),$

where
$\boldsymbol{\alpha}_{-i}=(\boldsymbol{\alpha}_{1}^{\prime},\ldots,\boldsymbol{\alpha}_{i-1}^{\prime},\boldsymbol{\alpha}_{i+1}^{\prime},\ldots,\boldsymbol{\alpha}_{n}^{\prime})^{\prime}$,

$\mathbf{K}_{\boldsymbol{\alpha}_{i}}=\mathbf{V}_{\boldsymbol{\alpha}_{i}}^{-1}+\mathbf{W}^{i^{\prime}}\mathbf{D}^{-1}\mathbf{W}^{i},\quad\widehat{\boldsymbol{\alpha}}_{i}=\mathbf{K}_{\boldsymbol{\alpha}_{i}}^{-1}\left(\mathbf{V}_{\boldsymbol{\alpha}_{i}}^{-1}\boldsymbol{\alpha}_{0,i}+\mathbf{W}^{i^{\prime}}\mathbf{D}^{-1}\mathbf{z}^{i}\right).$

Some of the above computations can be simplified by noting that
$\mathbf{B}_{0}$ is a lower triangular matrix. In addition, random draws from
the high-dimensional
$\mathcal{N}(\widehat{\boldsymbol{\alpha}}_{i},\mathbf{K}_{\boldsymbol{\alpha}_{i}}^{-1})$
distribution can be obtained efficiently without inverting any large matrices;
see, e.g., Chan (2021) for computational details.

Step 2: to sample $\boldsymbol{\beta}$, the lower triangular elements of
$\mathbf{B}_{0}$, first note that given $\mathbf{y}$ and
$\boldsymbol{\alpha}$, the innovations
$\boldsymbol{\varepsilon}_{1},\ldots,\boldsymbol{\varepsilon}_{T}$ can be
computed using (16). Next, we rewrite the model as a system of $(n-1)$
regressions in which $\varepsilon_{i,t}$ is regressed on the negative values
of $\varepsilon_{1,t},\ldots,\varepsilon_{i-1,t}$ for $i=2,\ldots,n$, and
$B_{0,i,1},\ldots,B_{0,i,i-1}$ are the corresponding regression coefficients.
More specifically, note that

$\displaystyle\mathbf{B}_{0}\boldsymbol{\varepsilon}_{t}$
$\displaystyle=\begin{pmatrix}1&0&0&\cdots&0\\\ B_{0,2,1}&1&0&\cdots&0\\\
B_{0,3,1}&B_{0,3,2}&1&\cdots&\vdots\\\ \vdots&\vdots&&\ddots&\vdots\\\
B_{0,n,1}&B_{0,n,2}&B_{0,n,3}&\cdots&1\end{pmatrix}\begin{pmatrix}\varepsilon_{1,t}\\\
\varepsilon_{2,t}\\\ \varepsilon_{3,t}\\\ \vdots\\\
\varepsilon_{n,t}\end{pmatrix}=\begin{pmatrix}\varepsilon_{1,t}\\\
\varepsilon_{2,t}+B_{0,2,1}\varepsilon_{1,t}\\\
\varepsilon_{3,t}+B_{0,3,1}\varepsilon_{1,t}+B_{0,3,2}\varepsilon_{2,t}\\\
\vdots\\\
\varepsilon_{n,t}+\sum_{j=1}^{n-1}B_{0,n,j}\varepsilon_{j,t}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\varepsilon_{1,t}\\\ \varepsilon_{2,t}\\\
\varepsilon_{3,t}\\\ \vdots\\\
\varepsilon_{n,t}\end{pmatrix}-\begin{pmatrix}0&0&0&0&0&\cdots&\cdots&0\\\
-\varepsilon_{1,t}&0&0&0&0&\cdots&&\vdots\\\
0&-\varepsilon_{1,t}&-\varepsilon_{2,t}&0&0&\cdots&&0\\\
\vdots&&&\ddots&\ddots&&\cdots&0\\\
0&\cdots&0&\cdots&0&-\varepsilon_{1,t}&\cdots&-\varepsilon_{t,n-1}\end{pmatrix}\begin{pmatrix}B_{0,2,1}\\\
B_{0,3,1}\\\ B_{0,3,2}\\\ \vdots\\\ B_{0,n,n-1}\end{pmatrix}.$

Or more succinctly,

$\mathbf{B}_{0}\boldsymbol{\varepsilon}_{t}=\boldsymbol{\varepsilon}_{t}-\mathbf{E}_{t}\boldsymbol{\beta}.$

Since $|\boldsymbol{\Sigma}_{t}|=|\mathbf{D}_{t}|$, we can rewrite the
likelihood implied by (16) as

$\displaystyle
p(\mathbf{y}\,|\,\boldsymbol{\alpha},\boldsymbol{\beta},\mathbf{h})$
$\displaystyle\propto\left(\prod_{t=1}^{T}|\mathbf{D}_{t}|^{-\frac{1}{2}}\right)\exp\left(-\frac{1}{2}\sum_{t=1}^{T}\boldsymbol{\varepsilon}_{t}^{\prime}(\mathbf{B}_{0}^{\prime}\mathbf{D}_{t}^{-1}\mathbf{B}_{0})\boldsymbol{\varepsilon}_{t}\right)$
$\displaystyle=\left(\prod_{t=1}^{T}|\mathbf{D}_{t}|^{-\frac{1}{2}}\right)\exp\left(-\frac{1}{2}\sum_{t=1}^{T}(\boldsymbol{\varepsilon}_{t}-\mathbf{E}_{t}\boldsymbol{\beta})^{\prime}\mathbf{D}_{t}^{-1}(\boldsymbol{\varepsilon}_{t}-\mathbf{E}_{t}\boldsymbol{\beta})\right).$

That is, the likelihood is the same as that implied by the regression

$\boldsymbol{\varepsilon}_{t}=\mathbf{E}_{t}\boldsymbol{\beta}+\boldsymbol{\eta}_{t},$
(18)

where $\boldsymbol{\eta}_{t}\sim\mathcal{N}(\mathbf{0},\mathbf{D}_{t}).$ In
addition, since the covariance matrix $\mathbf{D}_{t}$ is diagonal and
$\mathbf{E}_{t}$ is block-diagonal, we can sample
$\boldsymbol{\beta}_{2},\ldots,\boldsymbol{\beta}_{n}$ equation by equation to
improve computational efficiency.

To that end, we stack the $i$-th equation, $i=2,\ldots,n$, over $t=1,\ldots,T$
and obtain

$\boldsymbol{\varepsilon}_{i,\boldsymbol{\cdot}}=\mathbf{E}_{i,\boldsymbol{\cdot}}\boldsymbol{\beta}_{i}+\boldsymbol{\eta}_{i,\boldsymbol{\cdot}},$

where
$\boldsymbol{\varepsilon}_{i,\boldsymbol{\cdot}}=(\varepsilon_{i,1},\ldots,\varepsilon_{i,T})^{\prime}$,
$\mathbf{E}_{i,\boldsymbol{\cdot}}=(\boldsymbol{\varepsilon}_{1,\boldsymbol{\cdot}},\ldots,\boldsymbol{\varepsilon}_{i-1,\boldsymbol{\cdot}})$
and
$\boldsymbol{\eta}_{i,\boldsymbol{\cdot}}\sim\mathcal{N}(\mathbf{0},\mathbf{D}_{i,\boldsymbol{\cdot}})$
with
$\mathbf{D}_{i,\boldsymbol{\cdot}}=\operatorname{diag}(\text{e}^{h_{i,1}},\ldots,\text{e}^{h_{i,T}})$.
Finally, given the hierarchical Gaussian prior
$(\boldsymbol{\beta}_{i}\,|\,\boldsymbol{\kappa})\sim\mathcal{N}(\boldsymbol{\beta}_{0,i},\mathbf{V}_{\boldsymbol{\beta}_{i}})$,
it follows that

$(\boldsymbol{\beta}_{i}\,|\,\mathbf{y},\boldsymbol{\alpha},\mathbf{h}_{i,\boldsymbol{\cdot}},\boldsymbol{\kappa})\sim\mathcal{N}(\widehat{\boldsymbol{\beta}}_{i},\mathbf{K}_{\boldsymbol{\beta}_{i}}^{-1}),$

where