# Fast and Accurate Variational Inference for Large Bayesian VARs with
Stochastic Volatility

Joshua C.C. Chan
Purdue University Xuewen Yu
Purdue University

(This version: June 2022
First version: November 2020 )

###### Abstract

We propose a new variational approximation of the joint posterior distribution
of the log-volatility in the context of large Bayesian VARs. In contrast to
existing approaches that are based on _local_ approximations, the new proposal
provides a _global_ approximation that takes into account the entire support
of the joint distribution. In a Monte Carlo study we show that the new global
approximation is over an order of magnitude more accurate than existing
alternatives. We illustrate the proposed methodology with an application of a
96-variable VAR with stochastic volatility to measure global bank network
connectedness.

Keywords: large vector autoregression, stochastic volatility, Variational
Bayes, volatility network, connectedness

JEL classifications: C11, C32, C55, G21

##  1 Introduction

Since the influential work of Banbura, Giannone, and Reichlin (2010), large
Bayesian vector autoregressions (VARs) have been widely used to characterize
the comovements of a large number of macroeconomic and financial
variables.111Notable examples include Carriero, Kapetanios, and Marcellino
(2009), Koop (2013), Banbura, Giannone, Modugno, and Reichlin (2013),
Carriero, Clark, and Marcellino (2015), McCracken, Owyang, and Sekhposyan
(2015), Ellahie and Ricco (2017) and Morley and Wong (2020). Since a vast
empirical literature has demonstrated the importance of allowing for time-
varying volatility in small systems,222See, for example, Cogley and Sargent
(2005), Primiceri (2005), Clark (2011), D’Agostino, Gambetti, and Giannone
(2013), and Cross and Poon (2016). there is a lot of recent work that aims to
develop stochastic volatility specifications for large VARs, including Koop
and Korobilis (2013), Carriero, Clark, and Marcellino (2016, 2019), Chan
(2020a), Chan, Eisenstat, and Strachan (2020) and Kastner and Huber (2020).
Despite recent advances, estimating large VARs with flexible stochastic
volatility specifications using conventional Markov chain Monte Carlo (MCMC)
methods remains computationally intensive.

In view of the computational burden, some recent papers, such as Koop and
Korobilis (2018) and Gefang, Koop, and Poon (2019), have adopted an
alternative approach of using Variational Bayesian methods to approximate the
posterior distributions of large VARs with stochastic volatility. The main
advantage of these Variational Bayesian methods is that they are substantially
faster than MCMC, especially for high-dimensional models such as large VARs,
making estimation of very large systems possible. For example, fitting a
100-variable system takes only a few minutes compared to hours when MCMC is
used. However, they are approximate methods — as opposed to MCMC that can be
made arbitrarily accurate by increasing the simulation size — and the
approximation accuracy depends on the Kullback-Leibler divergence of the
approximating density to the posterior density.

Many existing approximating densities of the joint distribution of the log-
volatility are based on _local_ approximations, such as a second-order Taylor
expansion of the log target density around a point (e.g., the mode). As such,
these approximations are guaranteed to approximate the target density well
around the neighborhood of the point of expansion, but their accuracy
typically deteriorates rapidly away from the approximation point. In contrast,
we propose a _global_ approximation of the joint distribution of the log-
volatility that takes into account the entire support of the distribution. The
key idea is to set up a formal optimization problem to locate the ‘best’
density within a class of multivariate Gaussian distributions. More
specifically, we obtain the density within the family that is the closest to
the target posterior distribution, measured by the Kullback-Leibler divergence
between the two densities.

There are other Variational Bayesian methods that aim to ‘integrate out’ the
log-volatility using Monte Carlo draws from the joint conditional distribution
of the log-volatility. Examples include Tran, Nott, and Kohn (2017) and
Loaiza-Maya, Smith, Nott, and Danaher (2020), both of which consider
univariate stochastic volatility models. The key advantage of this approach is
that it provides a generally more accurate approximation of the posterior
distribution of model parameters. However, since this approach requires draws
of the log-volatility, it is more computationally intensive and is generally
not applicable to high-dimensional VARs.

To implement the proposed approach, we first reduce the difficult functional
optimization problem of finding the minimizer within a family of distributions
to a standard vector optimization problem by parameterizing the family of
Gaussian distributions. We then solve the associated optimization using the
Newton-Raphson method. Since the optimization problem is high-dimensional, we
carefully make use of fast band matrix routines to speed up computations.
Since the class of Gaussian distributions we construct includes some of the
existing Gaussian approximations, the optimal density located under this
proposed approach is guaranteed to be a better approximation — in the sense of
a smaller Kullback-Leibler divergence — than existing proposals.

We then demonstrate the superior approximation accuracy of the proposed
approach relative to existing methods in a Monte Carlo study. In particular,
the Monte Carlo results show that the mean squared errors under the proposed
global approximation can be more than an order of magnitude smaller than those
under existing local approximations.

The proposed methodology is illustrated using an application of a large VAR
with stochastic volatility to measure global bank network connectedness. More
specifically, we revisit the global bank network application in Demirer,
Diebold, Liu, and Yilmaz (2018) — they consider a 96-variable homoscedastic
VAR to measure global bank network connectedness. Since the connectedness
measures are functions of the error covariance matrix, how it is modeled might
be important for the analysis. We therefore extend their analysis by allowing
the error covariance matrix to vary over time via a stochastic volatility
process. Using data on bank returns volatility, we find qualitatively similar
results. In particular, our results show that North America and Europe are the
two largest net transmitters of future volatility uncertainty to the rest of
the world, whereas Asia is a large net receiver of future volatility
uncertainty from the rest of the world. We also find a substantial increase in
bank system-wide connectedness at the start of the Great Recession in late
2007.

The rest of the paper is organized as follows. We first present a
reparameterization of the reduced-form VAR with stochastic volatility in
Section 2, followed by some discussion of an adaptive Minnesota prior. Section
3 provides an overview of Variational Bayesian methods. We then introduce a
new global Gaussian approximation of the joint distribution of the log-
volatility in Section 4. Section 5 conducts a Monte Carlo study to provide
evidence of superior approximation accuracy of the proposed method relative to
existing approaches. It is followed by an application of using a VAR with
stochastic volatility to measure global bank network connectedness in Section
6. Lastly, Section 7 concludes and briefly discusses some future research
directions.

##  2 A VAR with Stochastic Volatility

In this section we describe a large VAR with stochastic volatility and then
outline an adaptive Minnesota prior. More specifically, we consider the
following reparameterization of the standard reduced-form VAR:

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

where
$\boldsymbol{\Sigma}_{t}=\text{diag}(\text{e}^{h_{1,t}},\ldots,\text{e}^{h_{n,t}})$
is a diagonal matrix and $\mathbf{B}_{0}$ is a lower triangular matrix with
ones on the main diagonal. Each log-volatility $h_{i,t},i=1,\ldots,n,$ evolves
as an independent random walk:

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

for $t=1,\ldots,T,$ where the initial condition $h_{i,0}$ is treated an
unknown parameter. It is straightforward to check that one can recover the
reduced-form intercepts and VAR coefficients by computing
$\widetilde{\mathbf{b}}=\mathbf{B}^{-1}_{0}\mathbf{b}$ and
$\widetilde{\mathbf{B}}_{j}=\mathbf{B}^{-1}_{0}\mathbf{B}_{j},j=1,\ldots,p$.
In addition, the implied reduced-form inverse covariance matrix, or precision
matrix, is
$\widetilde{\boldsymbol{\Sigma}}_{t}^{-1}=\mathbf{B}_{0}^{\prime}\boldsymbol{\Sigma}_{t}^{-1}\mathbf{B}_{0}$,
as considered in Cogley and Sargent (2005) and Carriero, Clark, and Marcellino
(2019).

Since the covariance matrix $\boldsymbol{\Sigma}_{t}$ in (1) is diagonal, we
can estimate this recursive system equation by equation without loss of
efficiency. Below we first rewrite (1) as $n$ separate univariate regressions.
For notational convenience, let $b_{i}$ denote the $i$-th element of
$\mathbf{b}$ and let $\mathbf{b}_{j,i}$ represent the $i$-th row of
$\mathbf{B}_{j}$. Then,
$\boldsymbol{\beta}_{i}=(b_{i},\mathbf{b}_{1,i},\ldots,\mathbf{b}_{p,i})^{\prime}$
is the intercept and VAR coefficients for the $i$-th equation. Furthermore,
let $\boldsymbol{\alpha}_{i}$ denote the free elements in the $i$-th row of
the impact matrix $\mathbf{B}_{0}$. Then, the $i$-th equation of the system in
(1) can be rewritten as:

$y_{i,t}=\widetilde{\mathbf{w}}_{i,t}\boldsymbol{\alpha}_{i}+\widetilde{\mathbf{x}}_{t}\boldsymbol{\beta}_{i}+\varepsilon_{i,t}^{y},\quad\varepsilon_{i,t}^{y}\sim\mathcal{N}(0,\text{e}^{h_{i,t}}),$

where $\widetilde{\mathbf{w}}_{i,t}=(-y_{1,t},\ldots,-y_{i-1,t})$ and
$\widetilde{\mathbf{x}}_{t}=(1,\mathbf{y}_{t-1}^{\prime},\ldots,\mathbf{y}_{t-p}^{\prime})$.
Note that this is a recursive system in which $y_{i,t}$ depends on the
contemporaneous variables $y_{1,t},\ldots,y_{i-1,t}$. Since the system is
recursive, the Jacobian of the change of variables from
$\boldsymbol{\varepsilon}_{t}^{y}$ to $\mathbf{y}_{t}$ has unit determinant.
Hence, the likelihood function has the usual Gaussian form.

Letting
$\mathbf{x}_{i,t}=(\widetilde{\mathbf{w}}_{i,t},\widetilde{\mathbf{x}}_{t})$,
we can further simplify the $i$-th equation as:

$y_{i,t}=\mathbf{x}_{i,t}\boldsymbol{\theta}_{i}+\varepsilon_{i,t}^{y},\quad\varepsilon_{i,t}^{y}\sim\mathcal{N}(0,\text{e}^{h_{i,t}}),$
(3)

where
$\boldsymbol{\theta}_{i}=(\boldsymbol{\alpha}_{i}^{\prime},\boldsymbol{\beta}_{i}^{\prime})^{\prime}$
is of dimension $k_{i}=np+i.$ We have therefore rewritten the VAR in (1) as a
system of $n$ separate univariate regressions. This representation facilitates
equation-by-equation estimation, which substantially speeds up the
computations.

To complete the model specification, we assume the following priors on the
parameters $\boldsymbol{\theta}_{i},h_{i,0}$ and
$\sigma_{h,i}^{2},i=1,\ldots,n$:

$\boldsymbol{\theta}_{i}\sim\mathcal{N}(\boldsymbol{\theta}_{0,i},\mathbf{V}_{\theta_{i}}),\quad
h_{i,0}\sim\mathcal{N}(0,V_{h_{i,0}}),\quad\sigma_{h,i}^{2}\sim\mathcal{IG}(\nu_{i},S_{i}),$

where $\mathcal{IG}(a,b)$ denotes the inverse-gamma distribution with mean
$b/(a-1)$. Since large VARs have a lot of parameters, a suitable shrinkage
prior on $\boldsymbol{\theta}_{i}$ is vital. Below we describe a version of
the Minnesota prior to elicit $\boldsymbol{\theta}_{0,i}$ and
$\mathbf{V}_{\theta_{i}}$.

We emphasize that even if other hierarchical shrinkage priors on
$\boldsymbol{\theta}_{i}$ are used, such as those in Bhattacharya, Pati,
Pillai, and Dunson (2015) and Griffin and Brown (2017), the proposed
variational Bayes method can be directly applied. Here we focus on the
Minnesota prior for two reasons. First, the Minnesota prior remains the most
popular shrinkage prior for large Bayesian VARs. Second, there is a growing
body of empirical evidence to suggest that it is more suitable for
macroeconomic data than other hierarchical shrinkage priors; see, for example,
Giannone, Lenza, and Primiceri (2017) and Cross, Hou, and Poon (2020).

For a general discussion of the Minnesota prior, we refer the readers to Koop
and Korobilis (2010), Karlsson (2013) or Chan (2020b). Below we outline how we
elicit $\boldsymbol{\theta}_{0,i}$ and $\mathbf{V}_{\theta_{i}}$. For growth
rates data, we set $\boldsymbol{\theta}_{0,i}=\mathbf{0}$ to shrink the VAR
coefficients to zero. For level data, $\boldsymbol{\theta}_{0,i}$ is also set
to be zero except for the coefficient associated with the first own lag, which
is set to be one. Next, for $\mathbf{V}_{\theta_{i}}$, we specify it to be
diagonal with the $k$-th diagonal element $V_{\theta_{i},k}$ set to be:

$V_{\theta_{i},k}=\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,\\\ \frac{s_{i}^{2}}{s_{j}^{2}},&\text{for
the $j$-th element of }\boldsymbol{\alpha}_{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,r=1,\ldots,n$.

Here the prior covariance matrix $\mathbf{V}_{\theta_{i}}$ depends on two key
hyperparameters: $\kappa_{1}$ and $\kappa_{2}$. The hyperparameter
$\kappa_{1}$ controls the overall shrinkage strength of the coefficients on
their own lags, while $\kappa_{2}$ controls those on lags of other variables.
In the application we select them by maximizing the variational lower bound of
the marginal likelihood on a two-dimensional grid.333The gold standard for
Bayesian model selection is the marginal likelihood. In our setting, however,
computing the marginal likelihood is computationally intensive, especially
over a two-dimensional grid. We instead use the variational lower bound, which
is readily available from the maximization, as a proxy. This is motivated by
papers such as Carriero, Clark, and Marcellino (2015) and Giannone, Lenza, and
Primiceri (2015), which show that one can substantially improve model fit and
forecast performance by selecting shrinkage hyperparameters in a data-based
fashion.

It is worth noting that the VAR in (1) is not order invariant for two reasons.
First, the VAR is written in structural form, and the prior on the structural-
form VAR coefficients induces a prior on the reduced-form parameters that
depend on the order of the variables. One can alleviate this issue by
eliciting prior means and variances on the reduced-form VAR coefficients, and
then derive the implied prior means and variances on the structural-form VAR
coefficients, along the lines suggested in Chan (2021). Second, the
multivariate stochastic volatility specification is constructed based on the
lower triangular impact matrix $\mathbf{B}_{0}$. As noted by Carriero, Clark,
and Marcellino (2019), since priors are independently elicited for
$\mathbf{B}_{0}$ and the stochastic volatility, the implied prior on the
covariance matrix
$\widetilde{\boldsymbol{\Sigma}}_{t}=\mathbf{B}_{0}^{-1}\boldsymbol{\Sigma}_{t}(\mathbf{B}_{0}^{-1})^{\prime}$
is not order invariant. This issue arises in all stochastic volatility models
that are based on a lower triangular parameterization, including the models in
Cogley and Sargent (2005) and Primiceri (2005). A few recent papers, such as
Bognanni (2018), Shin and Zhong (2020), Arias, Rubio-Ramirez, and Shin (2021)
and Chan, Koop, and Yu (2021), have considered order-invariant VARs with
multivariate stochastic volatility. Estimation of these order-invariant models
are more computationally intensive and applications typically involve small
and medium systems (e.g., 3-20 variables). It would therefore be useful to
develop Variational Bayesian methods for these models in the future.

##  3 Overview of Variational Bayes

The primary goal of Bayesian analysis is to characterize the posterior
distribution of the model parameters given the data, denoted as
$p(\boldsymbol{\theta}\,|\,\mathbf{y})$. Since this posterior distribution is
intractable for most econometric models, one often requires stochastic
simulation methods such as MCMC to characterize the posterior distribution.

In contrast, Variational Bayes is a collection of deterministic algorithms for
approximating the posterior distribution using a more tractable density. This
is done by first fixing a family of tractable densities. Then, we locate the
optimal density within this family by minimizing the Kullback-Leibler
divergence of the approximating density to the posterior density
$p(\boldsymbol{\theta}\,|\,\mathbf{y})$. Below we give a general overview of
this approach. For a more detailed discussion on Variational Bayesian methods,
see, e.g., Jordan, Ghahramani, Jaakkola, and Saul (1999), Chap. 10 of Bishop
(2006) and Ormerod and Wand (2010). Recent applications in econometrics
include Hajargasht and Woźniak (2018), Koop and Korobilis (2018), Gefang,
Koop, and Poon (2019) and Loaiza-Maya, Smith, Nott, and Danaher (2020).

Let $\mathcal{Q}$ denote a family of tractable densities within which to
locate the optimal approximating density. Recall that the Kullback-Leibler
divergence from a density $p_{1}$ to another density $p_{2}$ is defined as

$D_{KL}(p_{1}||p_{2})=\int
p_{1}(\mathbf{x})\log\frac{p_{1}(\mathbf{x})}{p_{2}(\mathbf{x})}\text{d}\mathbf{x}.$

Now, we locate the optimal approximating density $q^{*}(\boldsymbol{\theta})$
as the density in $\mathcal{Q}$ that minimizes the Kullback-Leibler divergence
to the posterior distribution $p(\boldsymbol{\theta}\,|\,\mathbf{y})$. More
precisely, $q^{*}(\boldsymbol{\theta})$ is the minimizer of the following
minimization problem:

$\min_{q\in\mathcal{Q}}D_{KL}(q||p(\boldsymbol{\theta}\,|\,\mathbf{y}))=\int
q(\boldsymbol{\theta})\log\frac{q(\boldsymbol{\theta})}{p(\boldsymbol{\theta}\,|\,\mathbf{y})}\text{d}\boldsymbol{\theta}.$

It turns out that minimizing the Kullback-Leibler divergence is equivalent to
maximizing the $q$-dependent lower bound on the marginal likelihood:

$\underline{p}(\mathbf{y};q)\equiv\exp\int
q(\boldsymbol{\theta})\log\frac{p(\mathbf{y},\boldsymbol{\theta})}{q(\boldsymbol{\theta})}\text{d}\boldsymbol{\theta}\leqslant
p(\mathbf{y}).$

Hence, this gives an alternative interpretation of the optimal density
$q^{*}(\boldsymbol{\theta})$ as the density in $\mathcal{Q}$ that has the
largest lower bound on the marginal likelihood $p(\mathbf{y})$. Moreover, the
lower bound is attained, i.e., $\underline{p}(\mathbf{y};q)=p(\mathbf{y})$ if
and only if $q(\boldsymbol{\theta})=p(\boldsymbol{\theta}\,|\,\mathbf{y})$.

In general both optimization problems are hard to solve as they involve a
typically high-dimensional integral. Nevertheless, one can substantially
simply the computations if the parameters $\boldsymbol{\theta}$ can be
naturally divided into $m$ blocks,
$\boldsymbol{\theta}_{1},\ldots,\boldsymbol{\theta}_{m}$, and the
approximating density $q$ is assumed to take the form

$q(\boldsymbol{\theta})=\prod_{i=1}^{m}q_{\boldsymbol{\theta}_{i}}(\boldsymbol{\theta}_{i}).$

This approach is known as mean field variational approximation. It amounts to
breaking up a high-dimensional optimization problem into $m$ lower-dimensional
ones. In particular, it can be shown that the optimal densities satisfy:

$q^{*}_{\boldsymbol{\theta}_{i}}(\boldsymbol{\theta}_{i})\propto\exp\left[\mathbb{E}_{-\boldsymbol{\theta}_{i}}\log
p(\mathbf{y},\boldsymbol{\theta})\right],\quad 1\leqslant i\leqslant m,$ (4)

where $\mathbb{E}_{-\boldsymbol{\theta}_{i}}$ denotes the expectation taken
with respect to the density $\prod_{j\neq
i}q_{\boldsymbol{\theta}_{j}}(\boldsymbol{\theta}_{j})$. This leads to an
iterative scheme that cycles through $i=1,\ldots,m$ via (4), until the
increase in the variational lower bound $\underline{p}(\mathbf{y};q)$ is
negligible. It can therefore be viewed as an instant of the coordinate ascent
method to find the maximizer in the function space.

For our VAR with stochastic volatility, the parameters are
$(\boldsymbol{\theta}_{i},h_{i,0},\sigma_{h,i}^{2},\mathbf{h}_{i}),i,\ldots,n$.
We approximate the posterior density
$p(\boldsymbol{\theta}_{i},h_{i,0},\sigma_{h,i}^{2},\mathbf{h}_{i}\,|\,\mathbf{y}_{i})$
using the approximating density of the form:

$q(\boldsymbol{\theta}_{i},h_{i,0},\sigma_{h,i}^{2},\mathbf{h}_{i})=q_{\boldsymbol{\theta}_{i}}(\boldsymbol{\theta}_{i})q_{h_{i,0}}(h_{i,0})q_{\sigma_{h,i}^{2}}(\sigma_{h,i}^{2})q_{\mathbf{h}_{i}}(\mathbf{h}_{i}),$

where the marginal densities $q_{\boldsymbol{\theta}_{i}},q_{h_{i,0}}$ and
$q_{\sigma_{h,i}^{2}}$ are unrestricted, whereas $q_{\mathbf{h}_{i}}$ is
assumed to be Gaussian. Since $q_{\mathbf{h}_{i}}$ is high-dimensional, it is
vital that the maximization step involving $q_{\mathbf{h}_{i}}$ is tractable,
and the Gaussian assumption provides a good trade-off between tractability and
flexibility. In the next section we propose a new global Gaussian
approximation $q_{\mathbf{h}_{i}}$ for $\mathbf{h}_{i}$. Detailed derivations
of the other densities and the variational lower bound are provided in
Appendix A.

##  4 A New Approximating Density of the Stochastic Volatility

In this section we introduce a _global_ Gaussian approximation of the joint
distribution of the log-volatility
$\mathbf{h}_{i}=(h_{i,1},\ldots,h_{i,T})^{\prime}$ in contrast to _local_
approximations that have been considered in the literature. For comparison, we
consider two local Gaussian approximations that have been recently used in
variational Bayes inference. The first Gaussian approximation is based on the
well-known approximation of the $\log$-$\chi^{2}_{1}$ distribution using the
$\mathcal{N}(-1.27,\pi^{2}/4)$ distribution. More specifically, let
$y_{i,t}^{*}=\log(y_{i,t}-\mathbf{x}_{i,t}\boldsymbol{\theta}_{i})^{2}$. Then,
one can rewrite (3) as

$y_{i,t}^{*}=h_{i,t}+\varepsilon_{i,t}^{y*},$

where $\varepsilon_{i,t}^{y*}$ follows the $\log$-$\chi^{2}_{1}$ distribution.
Given this transformation, Koop and Korobilis (2018) then replace the
$\log$-$\chi^{2}_{1}$ distribution with the $\mathcal{N}(-1.27,\pi^{2}/4)$
distribution. Consequently, the stochastic volatility model becomes
(approximately) a linear Gaussian state space model.444This approach of
approximating the stochastic volatility model can be traced back to Harvey,
Ruiz, and Shephard (1994), who suggest a quasi-maximum likelihood method to
estimate the linearized model based on the Kalman filter. One potential
problem of this approach is that the $\log$-$\chi^{2}_{1}$ distribution is
skewed and far from Gaussian, especially at the tails. In particular, the
$\log$-$\chi^{2}_{1}$ distribution has a heavier left tail but a much thinner
right tail compared to the Gaussian approximation. As a result, the stochastic
volatility estimates obtained using this approximation could be substantially
distorted.

In view of this problem, Gefang, Koop, and Poon (2019) suggest another
Gaussian distribution that is expected to provide a more accurate
approximation. More specifically, they approximate the optimal distribution of
$\mathbf{h}_{i}$ by a second-order Taylor approximation expanded around the
mode.555This type of local Gaussian approximation is first introduced to
estimate stochastic volatility models in the seminal papers by Durbin and
Koopman (1997) and Shephard and Pitt (1997). Specifically, they use a linear
Gaussian state space model to approximate the stochastic volatility model. The
distribution implied by the linear state space model is of course Gaussian,
which is then used as an importance sampling density (coupled with the Kalman
filter) to estimate the stochastic volatility model. Chan and Grant (2016) and
Chan and Eisenstat (2018) improve upon this approach by directly computing the
Gaussian approximating density using Newton-Raphson method based on fast band
matrix routines. They find that this Gaussian approximation works better in
their forecasting application than the one considered in Koop and Korobilis
(2018). Despite this improvement, Taylor expansion is a local approximation
that is guaranteed to work well only around the neighborhood of the point of
expansion (the mode of the optimal density in this case). Next, we introduce a
global Gaussian approximation that takes into account the entire support of
the distribution.

To set the stage, first note that the unrestricted optimal density of
$\mathbf{h}_{i}$ — i.e., not restricted to the class of Gaussian densities —
has the form:

$\widetilde{q}^{*}_{\mathbf{h}_{i}}(\mathbf{h}_{i})\propto\exp\left\\{\mathbb{E}_{-\mathbf{h}_{i}}\left[\log
p(\mathbf{h}_{i}\,|\,\mathbf{y}_{i},\boldsymbol{\theta}_{i},h_{i,0},\sigma_{h,i}^{2})\right]\right\\},$

where the expectation is taken with respect to the marginal density
$q_{-\mathbf{h}_{i}}(\boldsymbol{\theta}_{i},h_{i,0},\sigma_{h,i}^{2})=q_{\boldsymbol{\theta}_{i}}(\boldsymbol{\theta}_{i})q_{h_{i,0}}(h_{i,0})q_{\sigma_{h,i}^{2}}(\sigma_{h,i}^{2})$.
In particular, it can be shown that the log-density of
$\widetilde{q}^{*}_{\mathbf{h}_{i}}$ has the following explicit expression:

$\begin{split}\log\widetilde{q}^{*}_{\mathbf{h}_{i}}(\mathbf{h}_{i})&=\widetilde{c}_{\mathbf{h}_{i}}-\frac{1}{2}\sum_{t=1}^{T}h_{i,t}-\frac{1}{2}\sum_{t=1}^{T}\text{e}^{-h_{i,t}}\widehat{s}_{t}^{2}\\\
&\quad-\frac{1}{2}\mathbb{E}_{\sigma_{h,i}^{2}}\left[\frac{1}{\sigma_{h,i}^{2}}\right]\left(\sum_{t=2}^{T}(h_{i,t}-h_{i,t-1})^{2}+(h_{i,1}-\widehat{h}_{i,0})^{2}\right),\end{split}$
(5)

where $\widetilde{c}_{\mathbf{h}_{i}}$ and $\widehat{s}_{t}^{2}$ are constants
independent on $\mathbf{h}_{i}$ ($\widehat{s}_{t}^{2}$ depends on the
expectation and variance with respect to the density
$q_{\boldsymbol{\theta}_{i}}(\boldsymbol{\theta}_{i})$; its definition is
given in Appendix A). Unfortunately, $\widetilde{q}^{*}_{\mathbf{h}_{i}}$
given in (5) is a high-dimensional non-standard density that we cannot
directly use. To proceed, we approximate $\widetilde{q}^{*}_{\mathbf{h}_{i}}$
using a Gaussian distribution that is optimal in a well-defined sense.

The idea is to set up a formal optimization problem to locate the ‘best’
density within a parameterized class of Gaussian distributions. To that end,
consider the following family of Gaussian densities:

$\mathcal{G}=\left\\{f_{\mathcal{N}}(\cdot;\mathbf{m},\widehat{\mathbf{K}}_{\mathbf{h}_{i}}^{-1}):\mathbf{m}\in\mathbb{R}^{T}\right\\},$

where $f_{\mathcal{N}}(\cdot;\boldsymbol{\mu},\boldsymbol{\Sigma})$ is the
Gaussian density with mean vector $\boldsymbol{\mu}$ and covariance matrix
$\boldsymbol{\Sigma}$, and $\widehat{\mathbf{K}}_{\mathbf{h}_{i}}$ is the
negative Hessian of $\log\widetilde{q}^{*}_{\mathbf{h}_{i}}(\mathbf{h}_{i})$
evaluated at the mode of
$\log\widetilde{q}^{*}_{\mathbf{h}_{i}}(\mathbf{h}_{i})$.666One could expand
the class of Gaussian distributions by allowing the covariance matrix to vary
as well. By enlarging the class of distributions, the optimal density located
is expected to be better approximation to
$\widetilde{q}^{*}_{\mathbf{h}_{i}}(\mathbf{h}_{i})$. On the other hand, this
expanded class of distributions is much more challenging to handle due to the
complex nonlinear restrictions on the covariance matrix (i.e., symmetry and
positive-definiteness). We leave this possibility for future research. We then
locate the member in $\mathcal{G}$ that minimizes the Kullback-Leibler
divergence to $\widetilde{q}^{*}_{\mathbf{h}_{i}}$, say,
$f_{\mathcal{N}}(\cdot;\widehat{\mathbf{h}}_{i},\widehat{\mathbf{K}}_{\mathbf{h}_{i}}^{-1})$.
In other words, the optimal density for $\mathbf{h}_{i}$ we use, denoted as
$q^{*}_{\mathbf{h}_{i}}$, is then the
$\mathcal{N}(\widehat{\mathbf{h}}_{i},\widehat{\mathbf{K}}_{\mathbf{h}_{i}}^{-1})$
distribution. Since the class $\mathcal{G}$ includes the Gaussian
approximation proposed in Gefang, Koop, and Poon (2019), the optimal density
located under this proposed approach is guaranteed to be a better
approximation — in the sense of a smaller Kullback-Leibler divergence — than
the former approximation.

Now, to obtain the best Gaussian approximation within the class $\mathcal{G}$,
we consider the optimization problem:

$\min_{f\in\mathcal{G}}D_{KL}(f||\widetilde{q}^{*}_{\mathbf{h}_{i}})=\min_{\mathbf{m}\in\mathbb{R}^{T}}\mathbb{E}\left[\log\frac{f_{\mathcal{N}}(\mathbf{h}_{i};\mathbf{m},\widehat{\mathbf{K}}_{\mathbf{h}_{i}}^{-1})}{\widetilde{q}^{*}_{\mathbf{h}_{i}}(\mathbf{h}_{i})}\right],$
(6)

where the expectation is taken with respect to the density
$f_{\mathcal{N}}(\mathbf{h}_{i};\mathbf{m},\widehat{\mathbf{K}}_{\mathbf{h}_{i}}^{-1})$.
It turns out that the problem in (6) is a convex optimization problem with a
unique minimizer. In addition, we are able to derive analytical expressions
for the gradient and Hessian of the objective function, and therefore the
minimization problem can be quickly solved using the Newton-Raphson method.
Furthermore, the Hessian is a band matrix, and one can further speed up
computations by implementing fast band matrix routines.

Next, we provide some technical details in solving the optimization problem in
(6). For notational convenience, we write $f_{\mathbf{m}}(\cdot)\equiv
f_{\mathcal{N}}(\cdot;\mathbf{m},\widehat{\mathbf{K}}_{\mathbf{h}_{i}}^{-1})$.
Now, using equation (5), we can express
$\log\left[\frac{f_{\mathbf{m}}(\mathbf{h}_{i})}{\widetilde{q}^{*}_{\mathbf{h}_{i}}(\mathbf{h}_{i})}\right]$
as

$\displaystyle\log\left[\frac{f_{\mathbf{m}}(\mathbf{h}_{i})}{\widetilde{q}^{*}_{\mathbf{h}_{i}}(\mathbf{h}_{i})}\right]$
$\displaystyle=c_{1}-\frac{1}{2}(\mathbf{h}_{i}-\mathbf{m})^{\prime}\widehat{\mathbf{K}}_{\mathbf{h}_{i}}(\mathbf{h}_{i}-\mathbf{m})$
$\displaystyle\quad+\frac{1}{2}\left[\mathbf{1}_{T}^{\prime}\mathbf{h}_{i}+(\widehat{\mathbf{s}}^{2})^{\prime}\text{e}^{-\mathbf{h}_{i}}+\mathbb{E}_{\sigma_{h,i}^{2}}\left[\frac{1}{\sigma_{h,i}^{2}}\right](\mathbf{h}_{i}-\widehat{h}_{i,0}\mathbf{1}_{T})^{\prime}\mathbf{H}^{\prime}\mathbf{H}(\mathbf{h}_{i}-\widehat{h}_{i,0}\mathbf{1}_{T})\right],$

where $c_{1}$ is a constant independent of $\mathbf{h}_{i}$ and $\mathbf{m}$
and
$\widehat{\mathbf{s}}^{2}=(\widehat{s}_{1}^{2},\ldots,\widehat{s}_{T}^{2})^{\prime}$.
Then, taking expectation with respect to $f_{\mathbf{m}}$, we obtain:

$\mathbb{E}\log\left[\frac{f_{\mathbf{m}}(\mathbf{h}_{i})}{\widetilde{q}^{*}_{\mathbf{h}_{i}}(\mathbf{h}_{i})}\right]=c_{2}+\frac{1}{2}\left[\mathbf{1}_{T}^{\prime}\mathbf{m}+(\widehat{\mathbf{s}}^{2})^{\prime}\text{e}^{-\mathbf{m}+\frac{1}{2}\widehat{\mathbf{d}}_{i}}+\mathbb{E}_{\sigma_{h,i}^{2}}\left[\frac{1}{\sigma_{h,i}^{2}}\right](\mathbf{m}-\widehat{h}_{i,0}\mathbf{1}_{T})^{\prime}\mathbf{H}^{\prime}\mathbf{H}(\mathbf{m}-\widehat{h}_{i,0}\mathbf{1}_{T})\right],$

where $c_{2}$ is a constant independent of $\mathbf{m}$ and
$\widehat{\mathbf{d}}_{i}$ is a $T\times 1$ vector consisting of the diagonal
elements of $\widehat{\mathbf{K}}_{\mathbf{h}_{i}}^{-1}$. Given this
expression, it can be easily verify that
$\mathbb{E}\log\left[\frac{f_{\mathbf{m}}(\mathbf{h}_{i})}{\widetilde{q}^{*}_{\mathbf{h}_{i}}(\mathbf{h}_{i})}\right]$
is convex in $\mathbf{m}$. Hence, the minimization problem in (6) can be
solved readily. Furthermore, the gradient and the Hessian of the objective
function with respect to $\mathbf{m}$ can be computed easily:

grad
$\displaystyle=\mathbb{E}_{\sigma_{h,i}^{2}}\left[\frac{1}{\sigma_{h,i}^{2}}\right]\mathbf{H}^{\prime}\mathbf{H}(\mathbf{m}-\widehat{h}_{i,0}\mathbf{1}_{T})+\frac{1}{2}(\mathbf{1}_{T}-\widehat{\mathbf{s}}^{2}\odot\text{e}^{-\mathbf{m}+\frac{1}{2}\widehat{\mathbf{d}}_{i}}),$
Hess
$\displaystyle=\mathbb{E}_{\sigma_{h,i}^{2}}\left[\frac{1}{\sigma_{h,i}^{2}}\right]\mathbf{H}^{\prime}\mathbf{H}+\frac{1}{2}\text{diag}(\widehat{\mathbf{s}}^{2}\odot\text{e}^{-\mathbf{m}+\frac{1}{2}\widehat{\mathbf{d}}_{i}}),$

where $\odot$ denotes the component-wise product. Note that the Hessian is a
positive-definite matrix for all $\mathbf{m}\in\mathbb{R}^{T}$. Hence, Newton-
Raphson method can be used to quickly solve the minimization problem.
Furthermore, since the Hessian is also a band matrix, fast routines for band
matrices can be used to drastically speed up computations (see, e.g. Chan,
2017). Let $\widehat{\mathbf{h}}_{i}$ denote the unique minimizer. Finally, we
use
$f_{\mathcal{N}}(\cdot;\widehat{\mathbf{h}}_{i},\widehat{\mathbf{K}}_{\mathbf{h}_{i}}^{-1})$
as the optimal density $q^{*}_{\mathbf{h}_{i}}$.

##  5 Monte Carlo Experiments

In this section we conduct a Monte Carlo study to assess the accuracy of
approximating the posterior distribution of $\mathbf{h}$ using the proposed
global Gaussian approximation. We also document the runtimes of the
Variational Bayesian methods compared to MCMC for estimating VARs of different
dimensions.

First, we assess the accuracy of proposed variational approximation. As a
comparison, we include two local Gaussian approximations: 1) the Gaussian
distribution obtained by replacing the log-$\chi^{2}_{1}$ distribution in the
transformed observation equation with the $\mathcal{N}(-1.27,\pi^{2}/4)$
distribution; 2) a second-order Taylor approximation expanded around the mode
of the target posterior distribution of $\mathbf{h}$.

More specifically, we generate $R$ datasets from the following univariate
stochastic volatility model:

$\displaystyle z_{t}$ $\displaystyle=\text{e}^{{\frac{1}{2}h_{t}}}u_{t},$
$\displaystyle u_{t}$ $\displaystyle\sim\mathcal{N}(0,1),$ $\displaystyle
h_{t}$ $\displaystyle=h_{t-1}+v_{t},$ $\displaystyle v_{t}$
$\displaystyle\sim\mathcal{N}(0,0.1)$

for $t=1,\ldots,T$, and we set $h_{0}=0$. For each dataset
$\mathbf{z}^{(i)}=(z_{1}^{(i)},\ldots,z_{T}^{(i)})^{\prime},i=1,\ldots,R$ and
each Gaussian approximation with mean vector
$\widehat{\mathbf{h}}^{j}=(\widehat{h}_{1}^{j},\ldots,\widehat{h}_{T}^{j})^{\prime},j=1,\ldots,3$,
we compute the mean squared error (MSE) relative to the MCMC estimates,
$\text{MSE}_{i}(\widehat{\mathbf{h}}^{j})=\sum_{t=1}^{T}(\widehat{h}_{t}^{j}-\bar{h}_{t})^{2}/T$,
where $\bar{h}_{1},\ldots,\bar{h}_{T}$ are the posterior means obtained via
MCMC.

We first set the sample size to be $T=300$ and the number of Monte Carlo
replications to be $R=500$. The left panel in Figure 1 presents boxplots of
the MSEs for the three Gaussian approximations. Due to the differences in
scale, the right panel excludes the first approximation for better clarity. As
it is clear from the left panel, the first local approximation based on the
$\mathcal{N}(-1.27,\pi^{2}/4)$ approximation of the log-$\chi^{2}_{1}$
distribution (approx1) is typically an order of magnitude worse than the other
two Gaussian approximations. This approximation also performs poorly in
absolute terms for a substantial number of datasets.

Figure 1: Boxplots of the mean squared errors for the three Gaussian
approximations for $T=300$. Approx1 is the Gaussian distribution obtained by
the $\mathcal{N}(-1.27,\pi^{2}/4)$ approximation of the log-$\chi^{2}_{1}$
distribution and approx2 is the Gaussian distribution based on a second-order
Taylor approximation expanded around the mode of the target posterior
distribution. The central mark indicates the median, whereas the bottom and
top edges of the box indicate the 25-th and 75-th percentiles, respectively.
The whiskers extend to the minimum and the maximum.

Next, the right panel shows that the proposed global approximation is
substantially better than the second local approximation based on a second-
order Taylor expansion (approx2) — it provides another order of magnitude
reduction in MSE. For example, the median MSE of the proposed approximation is
only about 0.001, compared to about 0.012 for the second local approximation.
Moreover, for all datasets the proposed method provides a better approximation
— in terms of the lowest MSE — compared to the two alternatives.

We repeat the exercise with a longer sample of $T=2,000$ and the results are
reported in Figure 2. The main conclusion remains the same: the proposed
global approximation method is able to obtain a much more accurate
approximation of the log-volatility than existing local approximation methods.

Figure 2: Boxplots of the mean squared errors for the three Gaussian
approximations for $T=2,000$. Approx1 is the Gaussian distribution obtained by
the $\mathcal{N}(-1.27,\pi^{2}/4)$ approximation of the log-$\chi^{2}_{1}$
distribution and approx2 is the Gaussian distribution based on a second-order
Taylor approximation expanded around the mode of the target posterior
distribution. The central mark indicates the median, whereas the bottom and
top edges of the box indicate the 25-th and 75-th percentiles, respectively.
The whiskers extend to the minimum and the maximum.

After establishing that the proposed method is much more accurate than
existing approaches, we investigate its accuracy compared to MCMC. For this
purpose we compute the mean squared error against the true values of the log-
volatility:
$\text{MSE}(\widehat{\mathbf{h}})=\sum_{t=1}^{T}(\widehat{h}_{t}-h_{t})^{2}/T$,
where $\widehat{h}_{t}$ and $h_{t}$ are, respectively, estimates (obtained
from the proposed method or MCMC) and the true value. The results for $R=500$
Monte Carlo replications are reported in Figure 3. Each point in the scatter
plot corresponds to the MSEs of the proposed method and MCMC for one generated
dataset. While MCMC tends to be more accurate as expected, it is clear that
most of the points lie essentially on the diagonal line, indicating that the
approximation errors of the proposed method are fairly small.

Figure 3: Scatter plots of the mean squared errors for MCMC and the proposed
method for $T=300$.

Next, we document the runtimes of estimating VARs of different dimensions
using the proposed method versus MCMC. More specifically, Table 1 reports the
computation times to fit VARs of dimensions $n=25,50,100$ and sample sizes
$T=300,2,000$ using the two approaches (each value is the average of 100
datasets). For comparison we also include the runtimes of the two local
Gaussian approximations. The algorithms are implemented using
$\mathrm{M}\mathrm{{\scriptstyle ATLAB}}$ on a desktop with an Intel Core
i5-9600 @3.10 GHz processor and 16GB memory (when we implement the Variational
Bayesian methods, we do not use parallel computing for a fair comparison with
MCMC).

As it is evident from the table, the proposed approach is much faster than
conventional MCMC. For example, fitting a 100-variable VAR with $T=300$ using
the proposed approach takes only about 7 minutes, as opposed to about 70
minutes when MCMC is used. Among the three Variational Bayesian methods, the
local approximation based on the $\mathcal{N}(-1.27,\pi^{2}/4)$ distribution
(approx1) is noticeably faster than the other two approximations. But this is
at the cost of substantially larger approximation errors as presented earlier.
The proposed global approximation has similar runtimes as the local
approximation based on a second-order Taylor expansion (approx2), showing that
it offers better accuracy without sacrificing speed. In fact, in some
instances the proposed method is faster as it achieves convergence in fewer
iterations.

Table 1: The computation times (in minutes) to fit an $n$-variable VAR with a sample size $T$ using MCMC (to obtain 10,000 posterior draws) and the three Variational Bayesian methods: the Gaussian distribution obtained by the $\mathcal{N}(-1.27,\pi^{2}/4)$ approximation of the log-$\chi^{2}_{1}$ distribution (Approx1), the Gaussian distribution based on a second-order Taylor approximation expanded around the mode of the target posterior distribution (Approx2) and the proposed method. All VARs have $p=4$ lags. | $T=300$ | $T=2,000$
---|---|---
| $n=25$ | $n=50$ | $n=100$ | $n=25$ | $n=50$ | $n=100$
MCMC | 2.16 | 7.30 | 69.4 | 15.2 | 61.1 | 257.5
Approx1 | 0.01 | 0.01 | 0.16 | 0.05 | 0.11 | 0.24
Approx2 | 0.07 | 0.27 | 2.97 | 0.29 | 1.28 | 9.84
Proposed method | 0.06 | 0.22 | 6.99 | 0.21 | 0.93 | 6.78

##  6 Application: Bank Network Connectedness

In this section we illustrate the proposed methodology by using a large VAR
with stochastic volatility to measure global bank network connectedness with
the connectedness measures developed in the series of papers by Diebold and
Yilmaz (2009, 2014) and Demirer, Diebold, Liu, and Yilmaz (2018, DDLY
henceforth). In particular, we revisit the application in DDLY, who consider a
96-variable homoscedastic VAR regularized by the adaptive elastic net (Zou and
Zhang, 2009) to measure global bank network connectedness. Since the
connectedness measures are functions of the error covariance matrix, how it is
modeled is likely to be crucial. We therefore extend the analysis in DDLY by
allowing the error covariance matrix to vary over time via a stochastic
volatility process.777Korobilis and Yilmaz (2018) also use a time-varying VAR
to measure bank network connectedness. Their model is based on the discounted
Wishart process of Uhlig (1997) and West and Harrison (2006), which admits
efficient filtering and smoothing algorithms for estimation. However, the
discounted Wishart process seems to be too tightly parameterized for
macroeconomic data and it does not forecast well relative to standard
stochastic volatility models such as Cogley and Sargent (2005) and Primiceri
(2005). See, for example, Arias, Rubio-Ramirez, and Shin (2021) for a forecast
comparison exercise.

Another interesting aspect of the analysis is the impact of shrinkage methods
on the connectedness measures. While different shrinkage methods are expected
to have varying effects on the VAR coefficient estimates, their role on the
connectedness measures is less obvious. To regularize the large VAR, DDLY use
an adaptive elastic net — an average of Lasso and ridge penalties, where the
weights are inverses of the least squares estimates. The value of the overall
penalty weight is selected by 10-fold cross validation. In contrast, our
shrinkage method can be viewed as a combination of an adaptive and a
subjective ridge — the weights on individual VAR coefficients are subjectively
elicited according to the Minnesota prior, but the overall shrinkage parameter
is obtained by maximizing the marginal likelihood of the approximate model
(variational lower bound). It turns out that these two very different
shrinkage methods give quite similar results.

In what follows, we first define the connectedness measures as functions of
the VAR parameters. Then, we briefly discuss the dataset and present some
model comparison results. Finally, we report the connectedness measures
obtained from a VAR with stochastic volatility and an adaptive Minnesota
prior.

###  6.1 Connectedness Measures

In this section we define the connectedness measures that we use to
characterize bank network connectedness. These measures are based on variance
decompositions and are specifically designed to quantify how much individual
bank’s future uncertainty can be attributed to another specific bank or all
other banks as a whole.888The connectedness measures developed in Diebold and
Yilmaz (2009, 2014) and Demirer, Diebold, Liu, and Yilmaz (2018) are based on
a homoscedastic VAR. Here we extend these measures to a VAR with stochastic
volatility. They can be computed from the estimates of the VAR given in
(1)-(2).

More specifically, given the estimates (e.g., posterior means) of the
structural-form VAR in (1), we can recover the reduced-form estimates by
computing
$\widetilde{\mathbf{b}}=\mathbf{B}^{-1}_{0}\mathbf{b},\widetilde{\mathbf{B}}_{j}=\mathbf{B}^{-1}_{0}\mathbf{B}_{j},j=1,\ldots,p$
and
$\widetilde{\boldsymbol{\Sigma}}_{t}=\mathbf{B}_{0}^{-1}\boldsymbol{\Sigma}_{t}(\mathbf{B}_{0}^{-1})^{\prime}$.
Then, using these reduced-form estimates, we construct the corresponding
vector moving average matrices $\mathbf{A}_{h},h=0,1,2,\ldots$, with the
convention that $\mathbf{A}_{0}=\mathbf{I}_{n}$.

Next, following Diebold and Yilmaz (2014) we define the most granular
directional connectedness from one bank to another. More specifically, bank
$j$’s contribution to bank $i$’s $H$-step-ahead generalized forecast error
variance is defined as

$\theta_{ij,t}^{g}(H)=\frac{\widetilde{\sigma}_{jj,t}^{-1}\sum_{h=0}^{H-1}(\mathbf{e}_{i}^{\prime}\mathbf{A}_{h}\widetilde{\boldsymbol{\Sigma}}_{t}\mathbf{e}_{j})^{2}}{\sum_{h=0}^{H-1}(\mathbf{e}_{i}^{\prime}\mathbf{A}_{h}\widetilde{\boldsymbol{\Sigma}}_{t}\mathbf{A}_{h}^{\prime}\mathbf{e}_{i})^{2}},$

where $\widetilde{\sigma}_{jj,t}$ is the $j$-th diagonal elements of
$\widetilde{\boldsymbol{\Sigma}}_{t}$ and $\mathbf{e}_{i}$ is the selection
vector with one on the $i$-th position and and zeros otherwise. In contrast to
the static measure in DDLY, here the connectedness measure is time-varying.

Since these generalized forecast error variances might not sum to one, we
normalize them as follows:

$C_{i\leftarrow
j,t}^{H}=\frac{\theta_{ij,t}^{g}(H)}{\sum_{k=1}^{n}\theta_{ik,t}^{g}(H)}.$ (7)

Next, we aggregate these pairwise directional connectedness measures to form
total directional connectedness measures. The total directional connectedness
to bank $i$ from all other banks is:

$C_{i\leftarrow\bullet,t}^{H}=\frac{1}{n}\sum_{j=1,j\neq i}^{n}C_{i\leftarrow
j,t}^{H}.$ (8)

The total directional connectedness from bank $i$ to all other banks is
similarly defined as

$C_{\bullet\leftarrow i,t}^{H}=\frac{1}{n}\sum_{j=1,j\neq i}^{n}C_{j\leftarrow
i,t}^{H}.$ (9)

Finally, we can measure the total directional connectedness as

$C^{H}_{t}=\frac{1}{n}\sum_{i,j=1,i\neq j}^{n}C_{i\leftarrow j,t}^{H}.$ (10)

This measure is referred to as system-wide connectedness, as it aggregates the
total directional connectedness, both ‘to’ and ‘from’.

###  6.2 Data and Model Evidence

The dataset consists of daily stock prices of 96 banks from 29 developed and
emerging economies, and the sample period is from September 12, 2003 to
February 7, 2014. These 96 banks are those in the world’s top 150 by assets
that were publicly traded throughout the sample.999We thank Laura Liu for
providing us with the data and the associated R code. The dataset can also be
downloaded from the Journal of Applied Econometrics Data Archive. To measure
the connectedness in the global bank stock return volatility network, raw
daily stock prices (high, low, opening and closing prices) are used to compute
daily range-based realized volatility as proposed in Garman and Klass (1980).
This daily bank stock return volatility measure is then used as the dependent
variable. We refer the reader to DDLY for more details on the data.

We use the Minnesota prior described in Section 2, where the two shrinkage
priors $\kappa_{1}$ and $\kappa_{2}$ are selected by maximizing the
variational lower bound over a 2-dimensional grid. The optimal hyperparameter
values obtained are $\kappa_{1}=0.04$ and $\kappa_{2}=0.001$, again showing
much stronger shrinkage for coefficients on ‘other’ lags than on ‘own’ lags.
To see how well the VAR with stochastic volatility fits the data compared to a
standard homoscedastic VAR, we first obtain the variational lower bounds of
both models, and the results are reported in Table 2. For comparison, we also
include results based on the two local Gaussian approximations.

Table 2: Variational lower bounds of the standard homoscedastic VAR (VAR) and the VAR with stochastic volatility (VAR-SV). A larger value indicates a better model. VAR | VAR-SV
---|---
| approx1 | approx2 | proposed method
$-$286,050 |  $-$304,250 |  $-$283,520 |  $-$283,490

The results show that the VAR with stochastic volatility (approximated using
the proposed method) is strongly preferred by the data relative to its
homoscedastic counterpart. Since the variational lower bound of the marginal
likelihood has a built-in penalty for model complexity, the results indicate
that the increase in model-fit by allowing for time-varying variances
outweighs the cost of additional model complexity.

Among the three Variational Bayesian methods, the proposed method has the
largest variational lower bound, indicating that it gives the best
approximation. In particular, the variational lower bound associated with the
local approximation based on the $\mathcal{N}(-1.27,\pi^{2}/4)$ distribution
(approx1) is substantially smaller than the other two approximations,
suggesting that it provides a considerably worse model-fit. In addition, the
variational lower bound of the proposed method is about 30 in log scale larger
than that of the local approximation based on a second-order Taylor expansion
(approx2), highlighting that the proposed method offers much better
approximation accuracy. These findings are consistent with the Monte Carlo
results in Section 5 that show the proposed method delivers the lowest MSEs of
the log-volatility.

###  6.3 Connectedness Measures from a Large VAR with SV

Next, we report bank connectedness measures using a large VAR with stochastic
volatility. Since in our VAR the error covariance matrix is time-varying, the
connectedness measures defined in (7)–(10) are also time-varying. In contrast,
DDLY use rolling estimation with a 150-day window to characterize the global
banking network dynamically. In this section we compare the connectedness
measures under our VAR with stochastic volatility with their rolling-window
results.

We first report the bank network connectedness for the six-group aggregation
in Table 3. Since the network connectedness measures are time-varying under
the VAR with stochastic volatility, the values in the table are averages over
the whole sample 2003–2014.

The values of these connectedness measures are very similar to those under the
homoscedastic VAR in Table 4, even though now the model allows for stochastic
volatility. Consequently, the main message — that North America and Europe are
the two largest net transmitters of future volatility uncertainty and Asia is
a large net receiver of future volatility uncertainty from the rest of the
world — remains the same.

Table 3: Bank network connectedness for the six-group aggregation, 2003-2014, from the VAR-SV. | Africa | Asia | Europe | N. America | Oceania | S. America | From others
---|---|---|---|---|---|---|---
Africa | 0.00 | 7.56 | 24.79 | 21.81 | 2.18 | 2.33 | 58.67
Asia | 3.57 | 0.00 | 236.56 | 284.11 | 32.30 | 24.49 | 581.03
Europe | 5.16 | 64.37 | 0.00 | 678.02 | 29.32 | 33.26 | 810.13
N. America | 3.16 | 54.03 | 618.20 | 0.00 | 27.21 | 27.23 | 729.83
Oceania | 1.78 | 25.00 | 126.38 | 129.62 | 0.00 | 6.58 | 289.36
S. America | 1.12 | 13.35 | 51.39 | 52.32 | 2.94 | 0.00 | 121.12
To others | 14.80 | 164.30 | 1057.32 | 1165.88 | 93.95 | 93.88 | 2590.13

Next, we plot the dynamic system-wide connectedness in Figure 4. Our results
show an overall similar pattern compared to DDLY’s estimates from a
homoscedastic VAR obtained via a 150-day rolling window. In particular, our
dynamic estimates of the system-wide connectedness tend to increase from the
beginning of the sample to around 2008. The main difference is that our
measure peaks earlier — the first peak coincides with the start of the
liquidity crisis in August 2007, when the French bank BNP Paribas froze three
investment funds because of losses related to US subprime securities. This
caused the bond market to seize up, prompting the US Federal Reserve and the
European Central Bank to inject liquidity into the money markets to keep
interest rates down. Moreover, after this initial shock, volatility
connectedness remains high through the two waves of European Debt Crisis in
May 2010 and July-August 2011.

Figure 4: The dynamic system-wide connectedness measure from the VAR-SV.

To investigate the timing of the first peak, we follow DDLY to compute the
dynamic system-wide connectedness measure using a homoscedastic VAR with a
150-day rolling window. The results are reported in Appendix B. Despite the
very different shrinkage methods employed, our dynamic measure is remarkably
similar to that in DDLY. In particular, the dynamic system-wide measure
obtained using the rolling window increases substantially in 2007, but it does
not peak till around the collapse of Lehman Brothers in September 2008. To
better understand what drives these differences, we further estimate a VAR
with stochastic volatility using a 150-day rolling window (details are
reported in Appendix B). The rolling-window dynamic measure shows a similar
pattern: it increases drastically in 2007, but does not peak till September
2008. We thus conclude that the differences in the timing of the peak can be
attributed to the use of the rolling window. This could indicate that there
are structural breaks in the VAR coefficients.101010Another difference between
the rolling-window estimates and the estimates based on the full sample is
that the former exhibit more time variation. This is another indication that
there might be structural breaks in the VAR coefficients. Hence, developing
large time-varying parameter VARs, though computationally challenging, would
be a useful research direction.

Figure 5: Estimates of the dynamic system-wide connectedness measure from the
VAR-SV: computed using the variational approximation of the posterior means
(variational mean) and sample median of draws from the variational
approximation (sample median). The shaded areas represent the associated 68%
credit intervals (16- and 84-percentiles).

The dynamic system-wide connectedness measure presented in Figure 4 is
computed using the variational approximation of the posterior means of the VAR
coefficients and the log-volatility. Since the connectedness measure is a
nonlinear function of these model parameters, it is of interest to compare it
with the alternative approach of averaging draws from the variational
approximation. To that end, we sample 1,000 independent draws from the
variational approximation and compute the sample median of the connectedness
measure. An added advantage of this approach is that it also gives a measure
of parameter uncertainty. Figure 5 reports the sample median as well the the
associated 68% credit intervals (16- and 84-percentiles). While the estimates
based on the variational mean tend to be slightly smaller than those based on
the sample median, they track each other closely and have virtually the same
dynamic pattern.

Next, we decompose the dynamic system-wide connectedness into cross-country
and within-country components. More specifically, cross-country system-wide
connectedness is calculated as the sum of all pairwise connectedness across
banks located in different countries. Similarly, within-country system-wide
connectedness is the sum of pairwise connectedness across banks in the same
country. This decomposition allows us to explore the country origins of
volatility shocks and helps us better understand the dynamics of global bank
connectedness. The results of the decomposition are depicted in Figure 6.

Figure 6: The decomposition of the dynamic system-wide connectedness measure
from the VAR-SV into cross-country and within-country components.

Similar to the results reported in DDLY, we find that most variations in
system-wide connectedness are due to variations in cross-country system-wide
connectedness, whereas the within-country connectedness remains relatively
stable throughout the sample period. Our results show that cross-country
system-wide connectedness remains stable at around 50% from the beginning of
the sample till early 2017, and it then begins to fluctuate significantly.
Since the within-country connectedness remains relatively stable, cross-
country system-wide connectedness shows similar patterns as the total system-
wide connectedness depicted in Figure 4. In particular, cross-country system-
wide connectedness increases substantially around the liquidity crisis of
August 2007, and it remains high through the two waves of European Debt
Crisis.

##  7 Concluding Remarks and Future Research

We have developed a new variational approximation of the joint posterior
distribution of the log-volatility in the context of large VARs. In contrast
to existing approaches that are based on local approximations around a point
in the support of the distribution, the new method provides a global
approximation that takes into account the entire support. We have provided
evidence of superior approximation accuracy of the new proposal in a Monte
Carlo study compared to existing approximations.

Econometricians have only recently begun to use Variational Bayesian methods
as alternatives to MCMC for fitting high-dimensional models. In future
research, it would be useful to explore fitting other high-dimensional models,
such as large VARs with a factor stochastic volatility structure or with time-
varying VAR coefficients, using these methods. In addition, since the
variational lower bound can be obtained quite quickly, it would also be
interesting to explore using it to compare large stochastic volatility models
or shrinkage priors.

## Appendix A: Estimation Details

In this appendix we provide the technical details of the variational Bayes
approximation of the posterior distribution. Recall that the model can be
written as $n$ unrelated regressions for $i=1,\ldots,n$:

$y_{i,t}=\mathbf{x}_{i,t}\boldsymbol{\theta}_{i}+\varepsilon_{i,t}^{y},\quad\varepsilon_{i,t}^{y}\sim\mathcal{N}(0,\text{e}^{h_{i,t}}),$

where the log-volatility $h_{i,t}$ evolves as a random walk:

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

In addition, the priors on the parameters $\boldsymbol{\theta}_{i},h_{i,0}$
and $\sigma_{h,i}^{2}$ are also independent across equations. Specifically, we
assume for $i=1,\ldots,n$:

$\boldsymbol{\theta}_{i}\sim\mathcal{N}(\boldsymbol{\theta}_{0,i},\mathbf{V}_{\boldsymbol{\theta}_{i}}),\quad
h_{i,0}\sim\mathcal{N}(0,V_{h_{i,0}}),\quad\sigma_{h,i}^{2}\sim\mathcal{IG}(\nu_{i},S_{i}).$

Therefore, the joint distribution of the parameters and log-volatilities are
independent across equations. More specifically, the joint posterior
distribution of
$\boldsymbol{\theta}=(\boldsymbol{\theta}_{1}^{\prime},\ldots,\boldsymbol{\theta}_{n}^{\prime})^{\prime},$
$\mathbf{h}_{0}=(h_{1,0},\ldots,h_{n,0})^{\prime},$
$\boldsymbol{\sigma}_{h}^{2}=(\sigma_{h,1}^{2},\ldots,\sigma_{h,n}^{2})^{\prime}$
and
$\mathbf{h}=(\mathbf{h}_{1}^{\prime},\ldots,\mathbf{h}_{T}^{\prime})^{\prime}$
can be decomposed as:

$p(\boldsymbol{\theta},\mathbf{h}_{0},\boldsymbol{\sigma}_{h}^{2},\mathbf{h}\,|\,\mathbf{y})=\prod_{i=1}^{n}p(\boldsymbol{\theta}_{i},h_{i,0},\sigma_{h,i}^{2},\mathbf{h}_{i}\,|\,\mathbf{y}_{i}),$

where $\mathbf{y}_{i}=(y_{i,1},\ldots,y_{i,T})^{\prime}$. Hence, it suffices
to obtain a variational approximation of each of the $n$ components
$p(\boldsymbol{\theta}_{i},h_{i,0},\sigma_{h,i}^{2},\mathbf{h}_{i}\,|\,\mathbf{y}_{i}),i=1,\ldots,n$.

Now, we approximate
$p(\boldsymbol{\theta}_{i},h_{i,0},\sigma_{h,i}^{2},\mathbf{h}_{i}\,|\,\mathbf{y}_{i})$
using the product of densities:

$q(\boldsymbol{\theta}_{i},h_{i,0},\sigma_{h,i}^{2},\mathbf{h}_{i})=q_{\boldsymbol{\theta}_{i}}(\boldsymbol{\theta}_{i})q_{h_{i,0}}(h_{i,0})q_{\sigma_{h,i}^{2}}(\sigma_{h,i}^{2})q_{\mathbf{h}_{i}}(\mathbf{h}_{i}),$

where the marginal densities $q_{\boldsymbol{\theta}_{i}},q_{h_{i,0}}$ and
$q_{\sigma_{h,i}^{2}}$ are unrestricted, whereas $q_{\mathbf{h}_{i}}$ is
assumed to be Gaussian. Let
$q^{*}=q_{\boldsymbol{\theta}_{i}}^{*}q_{h_{i,0}}^{*}q_{\sigma_{h,i}^{2}}^{*}q_{\mathbf{h}_{i}}^{*}$
denote the optimal density. In what follows, we derive the explicit forms of
each of these optimal marginal densities and their associated parameters.

### The Optimal Density $q_{\boldsymbol{\theta}_{i}}^{*}$

The optimal density $q_{\boldsymbol{\theta}_{i}}^{*}$ has the form

$q_{\boldsymbol{\theta}_{i}}^{*}(\boldsymbol{\theta}_{i})\propto\exp\left\\{\mathbb{E}_{-\boldsymbol{\theta}_{i}}\left[\log
p(\boldsymbol{\theta}_{i}\,|\,\mathbf{y}_{i},\mathbf{h}_{i},h_{i,0},\sigma_{h,i}^{2})\right]\right\\},$

where the expectation is taken with respect to the marginal density
$q_{-\boldsymbol{\theta}_{i}}(h_{i,0},\sigma_{h,i}^{2},\mathbf{h}_{i})=q_{h_{i,0}}(h_{i,0})q_{\sigma_{h,i}^{2}}(\sigma_{h,i}^{2})q_{\mathbf{h}_{i}}(\mathbf{h}_{i})$.
To derive an explicit expression of $q_{\boldsymbol{\theta}_{i}}^{*}$, first
note that $\boldsymbol{\theta}_{i}$ is conditionally independent of
$(h_{i,0},\sigma_{h,i}^{2})$ given $(\mathbf{h}_{i},\mathbf{y}_{i})$. In
particular, the log-density is given by

$\log
p(\boldsymbol{\theta}_{i}\,|\,\mathbf{y}_{i},\mathbf{h}_{i})=c_{\boldsymbol{\theta}_{i}}-\frac{1}{2}\sum_{t=1}^{T}\text{e}^{-h_{i,t}}(y_{i,t}-\mathbf{x}_{i,t}\boldsymbol{\theta}_{i})^{2}-\frac{1}{2}(\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{0,i})^{\prime}\mathbf{V}_{\boldsymbol{\theta}_{i}}^{-1}(\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{0,i}),$

where $c_{\boldsymbol{\theta}_{i}}$ is a constant not dependent on
$\boldsymbol{\theta}_{i}$. Let
$\widehat{\mathbf{h}}_{i}=(\widehat{h}_{i,1},\ldots,\widehat{h}_{i,T})^{\prime}$
and $\widehat{\mathbf{K}}_{\mathbf{h}_{i}}$ denote respectively the mean
vector and precision matrix (i.e., inverse covariance matrix) of
$\mathbf{h}_{i}$ with respect to the density
$q_{\mathbf{h}_{i}}(\mathbf{h}_{i})$. Then, taking expectation of $\log
p(\boldsymbol{\theta}_{i}\,|\,\mathbf{y}_{i},\mathbf{h}_{i})$ gives

$\mathbb{E}_{-\boldsymbol{\theta}_{i}}\left[\log
p(\boldsymbol{\theta}_{i}\,|\,\mathbf{y}_{i},\mathbf{h}_{i})\right]=c_{\boldsymbol{\theta}_{i}}-\frac{1}{2}\sum_{t=1}^{T}\text{e}^{-\widehat{h}_{i,t}+\frac{1}{2}\widehat{d}_{i,t}}(y_{i,t}-\mathbf{x}_{i,t}\boldsymbol{\theta}_{i})^{2}-\frac{1}{2}(\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{0,i})^{\prime}\mathbf{V}_{\boldsymbol{\theta}_{i}}^{-1}(\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{0,i}),$

where $\widehat{d}_{i,t}$ is the $t$-th diagonal element of
$\widehat{\mathbf{K}}_{\mathbf{h}_{i}}^{-1}$. Since
$\mathbb{E}_{-\boldsymbol{\theta}_{i}}\left[\log
p(\boldsymbol{\theta}_{i}\,|\,\mathbf{y}_{i},\mathbf{h}_{i})\right]$ is a
(negative) quadratic form in $\boldsymbol{\theta}_{i}$, the optimal density
$q_{\boldsymbol{\theta}_{i}}^{*}$ is Gaussian. To derive the corresponding
mean vector and precision matrix, we rewrite
$\mathbb{E}_{-\boldsymbol{\theta}_{i}}\left[\log
p(\boldsymbol{\theta}_{i}\,|\,\mathbf{y}_{i},\mathbf{h}_{i})\right]$ in matrix
form as

$\mathbb{E}_{-\boldsymbol{\theta}_{i}}\left[\log
p(\boldsymbol{\theta}_{i}\,|\,\mathbf{y}_{i},\mathbf{h}_{i})\right]=c_{\boldsymbol{\theta}_{i}}-\frac{1}{2}(\mathbf{y}_{i}-\mathbf{X}_{i}\boldsymbol{\theta}_{i})^{\prime}\widehat{\mathbf{O}}_{\mathbf{h}_{i}}(\mathbf{y}_{i}-\mathbf{X}_{i}\boldsymbol{\theta}_{i})-\frac{1}{2}(\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{0,i})^{\prime}\mathbf{V}_{\boldsymbol{\theta}_{i}}^{-1}(\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{0,i}),$

where
$\widehat{\mathbf{O}}_{\mathbf{h}_{i}}=\text{diag}(\text{e}^{-\widehat{h}_{i,1}+\frac{1}{2}\widehat{d}_{i,1}},\ldots,\text{e}^{-\widehat{h}_{i,T}+\frac{1}{2}\widehat{d}_{i,T}})$
and
$\mathbf{X}_{i}=(\mathbf{x}_{i,1}^{\prime},\ldots,\mathbf{x}_{i,T}^{\prime})^{\prime}$.
Now, using standard linear regression results (see, .e.g, Chan, Koop, Poirier,
and Tobias, 2019, pp. 182-188), one can show that it is the
$\mathcal{N}(\widehat{\boldsymbol{\theta}}_{i},\widehat{\mathbf{K}}^{-1}_{\boldsymbol{\theta}_{i}})$
distribution, where

$\widehat{\mathbf{K}}_{\boldsymbol{\theta}_{i}}=\mathbf{V}_{\boldsymbol{\theta}_{i}}^{-1}+\mathbf{X}_{i}^{\prime}\widehat{\mathbf{O}}_{\mathbf{h}_{i}}\mathbf{X}_{i},\quad\widehat{\boldsymbol{\theta}}_{i}=\widehat{\mathbf{K}}_{\boldsymbol{\theta}_{i}}^{-1}(\mathbf{V}_{\boldsymbol{\theta}_{i}}^{-1}\boldsymbol{\theta}_{0,i}+\mathbf{X}_{i}^{\prime}\widehat{\mathbf{O}}_{\mathbf{h}_{i}}\mathbf{y}_{i}).$

### The Optimal Density $q^{*}_{h_{i,0}}$

Next, we derive the optimal density $q^{*}_{h_{i,0}}$, which takes the form

$q_{h_{i,0}}^{*}(h_{i,0})\propto\exp\left\\{\mathbb{E}_{-h_{i,0}}\left[\log
p(h_{i,0}\,|\,\mathbf{y}_{i},\boldsymbol{\theta}_{i},\mathbf{h}_{i},\sigma_{h,i}^{2})\right]\right\\},$

where the expectation is taken with respect to the marginal density
$q_{-h_{i,0}}(\boldsymbol{\theta}_{i},\sigma_{h,i}^{2},\mathbf{h}_{i})=q_{\boldsymbol{\theta}_{i}}(\boldsymbol{\theta}_{i})q_{\sigma_{h,i}^{2}}(\sigma_{h,i}^{2})q_{\mathbf{h}_{i}}(\mathbf{h}_{i})$.
First write

$\log
p(h_{i,0}\,|\,\mathbf{y}_{i},\boldsymbol{\theta}_{i},\mathbf{h}_{i},\sigma_{h,i}^{2})=\log
p(h_{i,0}\,|\,\mathbf{h}_{i},\sigma_{h,i}^{2})=c_{h_{i,0}}-\frac{1}{2\sigma_{h,i}^{2}}(h_{i,1}-h_{i,0})^{2}-\frac{1}{2V_{h_{i,0}}}h_{i,0}^{2},$

where $c_{h_{i,0}}$ is a constant not dependent on $h_{i,0}$. Then, taking
expectation with respect to the marginal density $q_{-h_{i,0}}$, we obtain

$\mathbb{E}_{-h_{i,0}}\left[\log
p(h_{i,0}\,|\,\mathbf{h}_{i},\sigma_{h,i}^{2})\right]=c_{h_{i,0}}-\frac{1}{2}\mathbb{E}_{\sigma_{h,i}^{2}}\left[\frac{1}{\sigma_{h,i}^{2}}\right]\left[(\widehat{h}_{i,1}-h_{i,0})^{2}+\widehat{d}_{i,1}\right]-\frac{1}{2V_{h_{i,0}}}h_{i,0}^{2},$

where $\widehat{d}_{i,1}$ is the first diagonal element of
$\widehat{\mathbf{K}}_{\mathbf{h}_{i}}^{-1}$ and the expectation
$\mathbb{E}_{\sigma_{h,i}^{2}}$ is taken with respect to the density
$q_{\sigma_{h,i}^{2}}(\sigma_{h,i}^{2})$ — this expectation can be computed
analytically as shown in the next subsection. Finally, using standard linear
regression results, one can show that $q^{*}_{h_{i,0}}$ is the
$\mathcal{N}(\widehat{h}_{i,0},\widehat{K}^{-1}_{h_{i,0}})$ distribution,
where

$\widehat{K}_{h_{i,0}}=V_{h_{i,0}}^{-1}+\mathbb{E}_{\sigma_{h,i}^{2}}\left[\frac{1}{\sigma_{h,i}^{2}}\right],\quad\widehat{h}_{i,0}=\widehat{K}_{h_{i,0}}^{-1}\mathbb{E}_{\sigma_{h,i}^{2}}\left[\frac{1}{\sigma_{h,i}^{2}}\right]\widehat{h}_{i,1}.$

### The Optimal Density $q^{*}_{\sigma_{h,i}^{2}}$

The kernel of the optimal density $q^{*}_{\sigma_{h,i}^{2}}$ is given by

$q^{*}_{\sigma_{h,i}^{2}}\propto\exp\left\\{\mathbb{E}_{-\sigma_{h,i}^{2}}\left[\log
p(\sigma_{h,i}^{2}\,|\,\mathbf{h}_{i},h_{i,0})\right]\right\\},$

where the expectation is taken with respect to the marginal density
$q_{-\sigma_{h,i}^{2}}(\boldsymbol{\theta}_{i},h_{i,0},\mathbf{h}_{i})=q_{\boldsymbol{\theta}_{i}}(\boldsymbol{\theta}_{i})q_{h_{i,0}}(h_{i,0})q_{\mathbf{h}_{i}}(\mathbf{h}_{i})$.
To derive an explicit expression for $q^{*}_{\sigma_{h,i}^{2}}$, first note
that

$\log
p(\sigma_{h,i}^{2}\,|\,\mathbf{h}_{i},h_{i,0})=c_{\sigma_{h,i}^{2}}-\frac{T}{2}\log\sigma_{h,i}^{2}-\frac{1}{2\sigma_{h,i}^{2}}(\mathbf{h}_{i}-h_{i,0}\mathbf{1}_{T})^{\prime}\mathbf{H}^{\prime}\mathbf{H}(\mathbf{h}_{i}-h_{i,0}\mathbf{1}_{T})-\nu_{i}\log\sigma_{h,i}^{2}-\frac{S_{i}}{\sigma_{h,i}^{2}},$

where $c_{\sigma_{h,i}^{2}}$ is a constant not dependent on
$\sigma_{h,i}^{2}$. Taking expectation with respect to the marginal density
$q_{-\sigma_{h,i}^{2}}$ gives

$\displaystyle\mathbb{E}_{-\sigma_{h,i}^{2}}\left[\log
p(\sigma_{h,i}^{2}\,|\,\mathbf{h}_{i},h_{i,0})\right]$
$\displaystyle=c_{\sigma_{h,i}^{2}}-\left(\nu_{i}+\frac{T}{2}\right)\log\sigma_{h,i}^{2}-\frac{S_{i}}{\sigma_{h,i}^{2}}$
$\displaystyle\quad-\frac{1}{2\sigma_{h,i}^{2}}\left[(\widehat{\mathbf{h}}_{i}-\widehat{h}_{i,0}\mathbf{1}_{T})^{\prime}\mathbf{H}^{\prime}\mathbf{H}(\widehat{\mathbf{h}}_{i}-\widehat{h}_{i,0}\mathbf{1}_{T})+\text{tr}(\mathbf{H}^{\prime}\mathbf{H}\widehat{\mathbf{K}}_{\mathbf{h}_{i}}^{-1})+\widehat{K}_{h_{i,0}}^{-1}\right].$

Exponentiating the term on the right-hand side, one recognizes that it is the
kernel of the $\mathcal{IG}(\widehat{\nu_{i}},\widehat{S}_{i})$ distribution,
where

$\widehat{\nu_{i}}=\nu_{i}+\frac{T}{2},\quad\widehat{S}_{i}=S_{i}+\frac{1}{2}\left[(\widehat{\mathbf{h}}_{i}-\widehat{h}_{i,0}\mathbf{1}_{T})^{\prime}\mathbf{H}^{\prime}\mathbf{H}(\widehat{\mathbf{h}}_{i}-\widehat{h}_{i,0}\mathbf{1}_{T})+\text{tr}(\mathbf{H}^{\prime}\mathbf{H}\widehat{\mathbf{K}}_{\mathbf{h}_{i}}^{-1})+\widehat{K}_{h_{i,0}}^{-1}\right].$

Since $q^{*}_{\sigma_{h,i}^{2}}$ is an inverse-gamma density, the expectation
of $1/\sigma_{h,i}^{2}$ can be obtained analytically as:

$\mathbb{E}_{\sigma_{h,i}^{2}}\left[\frac{1}{\sigma_{h,i}^{2}}\right]=\frac{\widehat{\nu_{i}}}{\widehat{S}_{i}}.$

### The Optimal Density $q^{*}_{\mathbf{h}_{i}}$

The unrestricted optimal density of $\mathbf{h}_{i}$ — i.e., not restricted to
the class of Gaussian densities — has the form:

$\widetilde{q}^{*}_{\mathbf{h}_{i}}(\mathbf{h}_{i})\propto\exp\left\\{\mathbb{E}_{-\mathbf{h}_{i}}\left[\log
p(\mathbf{h}_{i}\,|\,\mathbf{y}_{i},\boldsymbol{\theta}_{i},h_{i,0},\sigma_{h,i}^{2})\right]\right\\}.$

To derive an explicit expression for $\widetilde{q}^{*}_{\mathbf{h}_{i}}$,
note that

$\log
p(\mathbf{h}_{i}\,|\,\mathbf{y}_{i},\boldsymbol{\theta}_{i},h_{i,0},\sigma_{h,i}^{2})=c_{\mathbf{h}_{i}}-\frac{1}{2}\sum_{t=1}^{T}h_{i,t}-\frac{1}{2}\sum_{t=1}^{T}\text{e}^{-h_{i,t}}(y_{i,t}-\mathbf{x}_{i,t}\boldsymbol{\theta}_{i})^{2}-\frac{1}{2\sigma_{h,i}^{2}}\sum_{t=1}^{T}(h_{i,t}-h_{i,t-1})^{2},$

where $c_{\mathbf{h}_{i}}$ is a constant not dependent on $\mathbf{h}_{i}$.
Hence, taking expectation with respect to the marginal density
$q_{-\mathbf{h}_{i}}(\boldsymbol{\theta}_{i},h_{i,0},\sigma_{h,i}^{2})$ gives

$\displaystyle\mathbb{E}_{-\mathbf{h}_{i}}\left[\log
p(\mathbf{h}_{i}\,|\,\mathbf{y}_{i},\boldsymbol{\theta}_{i},h_{i,0},\sigma_{h,i}^{2})\right]$
$\displaystyle=c_{\mathbf{h}_{i}}-\frac{1}{2}\sum_{t=1}^{T}h_{i,t}-\frac{1}{2}\sum_{t=1}^{T}\text{e}^{-h_{i,t}}\widehat{s}_{t}^{2}-\frac{1}{2}\mathbb{E}_{\sigma_{h,i}^{2}}\left[\frac{1}{\sigma_{h,i}^{2}}\right]$
$\displaystyle\quad\times\left(\sum_{t=2}^{T}(h_{i,t}-h_{i,t-1})^{2}+(h_{i,1}-\widehat{h}_{i,0})^{2}+\widehat{K}_{h_{i,0}}^{-1}\right),$

where
$\widehat{s}_{t}^{2}=(y_{i,t}-\mathbf{x}_{i,t}\widehat{\boldsymbol{\theta}}_{i})^{2}+\text{tr}(\mathbf{x}_{i,t}^{\prime}\mathbf{x}_{i,t}\widehat{\mathbf{K}}^{-1}_{\boldsymbol{\theta}_{i}}).$
Therefore, the (log) kernel of $\widetilde{q}^{*}_{\mathbf{h}_{i}}$ has the
following explicit expression:

$\begin{split}\log\widetilde{q}^{*}_{\mathbf{h}_{i}}(\mathbf{h}_{i})=\widetilde{c}_{\mathbf{h}_{i}}-\frac{1}{2}\sum_{t=1}^{T}h_{i,t}-\frac{1}{2}\sum_{t=1}^{T}\text{e}^{-h_{i,t}}\widehat{s}_{t}^{2}-\frac{1}{2}\mathbb{E}_{\sigma_{h,i}^{2}}\left[\frac{1}{\sigma_{h,i}^{2}}\right]\left(\sum_{t=2}^{T}(h_{i,t}-h_{i,t-1})^{2}+(h_{i,1}-\widehat{h}_{i,0})^{2}\right),\end{split}$

where $\widetilde{c}_{\mathbf{h}_{i}}$ is a constant independent on
$\mathbf{h}_{i}$.

As discussed in the main text, we then locate the optimal Gaussian density
$q^{*}_{\mathbf{h}_{i}}$ within the following family parameterized by the mean
vector $\mathbf{m}$:

$\mathcal{G}=\left\\{f_{\mathcal{N}}(\cdot;\mathbf{m},\widehat{\mathbf{K}}_{\mathbf{h}_{i}}^{-1}):\mathbf{m}\in\mathbb{R}^{T}\right\\},$

where $\widehat{\mathbf{K}}_{\mathbf{h}_{i}}$ is the negative Hessian of
$\log\widetilde{q}^{*}_{\mathbf{h}_{i}}(\mathbf{h}_{i})$ evaluated at the mode
of $\log\widetilde{q}^{*}_{\mathbf{h}_{i}}(\mathbf{h}_{i})$. The optimal
density is obtained by solving the minimization problem defined in (6), which
for convenience we reproduce below:

$\min_{f_{\mathbf{m}}\in\mathcal{G}}D_{KL}(f_{\mathbf{m}}||\widetilde{q}^{*}_{\mathbf{h}_{i}})=\min_{\mathbf{m}\in\mathbb{R}^{T}}\mathbb{E}\log\left[\frac{f_{\mathbf{m}}(\mathbf{h}_{i})}{\widetilde{q}^{*}_{\mathbf{h}_{i}}(\mathbf{h}_{i})}\right],$

where the expectation is taken with respect to the density
$f_{\mathbf{m}}(\mathbf{h}_{i})$. This minimization problem can be quickly
solved using the Newton-Raphson method. We use
$f_{\mathcal{N}}(\cdot;\widehat{\mathbf{h}}_{i},\widehat{\mathbf{K}}_{\mathbf{h}_{i}}^{-1})$
as the optimal density $q^{*}_{\mathbf{h}_{i}}$, where
$\widehat{\mathbf{h}}_{i}$ is the unique minimizer.

### The Variational Lower Bound

Next, we derive the variational lower bound $\underline{p}(\mathbf{y}_{i};q)$.
To that end, we first compute the log ratio of the joint posterior density and
the variational approximation:

$\displaystyle\log$
$\displaystyle\left[\frac{p(\mathbf{y}_{i},\mathbf{h}_{i}\,|\,\boldsymbol{\theta}_{i},h_{i,0},\sigma_{h,i}^{2})p(\boldsymbol{\theta}_{i},h_{i,0},\sigma_{h,i}^{2})}{q(\boldsymbol{\theta}_{i},h_{i,0},\sigma_{h,i}^{2},\mathbf{h}_{i})}\right]=c_{i}-\frac{1}{2}\mathbf{1}_{T}^{\prime}\mathbf{h}_{i}-\frac{1}{2}(\mathbf{y}_{i}-\mathbf{X}_{i}\boldsymbol{\theta}_{i})^{\prime}\boldsymbol{\Omega}_{\mathbf{h}_{i}}^{-1}(\mathbf{y}_{i}-\mathbf{X}_{i}\boldsymbol{\theta}_{i})$
$\displaystyle-\frac{T}{2}\log\sigma_{h,i}^{2}-\frac{1}{2\sigma_{h,i}^{2}}(\mathbf{h}_{i}-h_{i,0}\mathbf{1}_{T})^{\prime}\mathbf{H}^{\prime}\mathbf{H}(\mathbf{h}_{i}-h_{i,0}\mathbf{1}_{T})-\frac{1}{2}(\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{0,i})^{\prime}\mathbf{V}_{\boldsymbol{\theta}_{i}}^{-1}(\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{0,i})$
$\displaystyle-\frac{1}{2V_{h_{i,0}}}h_{i,0}^{2}-(\nu_{i}+1)\log\sigma_{h,i}^{2}-\frac{S_{i}}{\sigma_{h,i}^{2}}+\frac{1}{2}(\mathbf{h}_{i}-\widehat{\mathbf{h}}_{i})^{\prime}\widehat{\mathbf{K}}_{\mathbf{h}_{i}}(\mathbf{h}_{i}-\widehat{\mathbf{h}}_{i})$
$\displaystyle+\frac{1}{2}(\boldsymbol{\theta}_{i}-\widehat{\boldsymbol{\theta}}_{i})^{\prime}\widehat{\mathbf{K}}_{\boldsymbol{\theta}_{i}}(\boldsymbol{\theta}_{i}-\widehat{\boldsymbol{\theta}}_{i})+\frac{\widehat{K}_{h_{i,0}}}{2}(h_{i,0}-\widehat{h}_{i,0})^{2}+(\widehat{\nu}_{i}+1)\log\sigma_{h,i}^{2}+\frac{\widehat{S}_{i}}{\sigma_{h,i}^{2}},$

where
$\boldsymbol{\Omega}_{\mathbf{h}_{i}}=\text{diag}(\text{e}^{h_{i,1}},\ldots,\text{e}^{h_{i,T}})$
and $c_{i}=-\frac{T}{2}\log(2\pi)-\frac{1}{2}\log
V_{h_{i,0}}-\frac{1}{2}\log|\mathbf{V}_{\boldsymbol{\theta}_{i}}|+\nu_{i}\log
S_{i}-\log\Gamma(\nu_{i})-\frac{1}{2}\log|\widehat{\mathbf{K}}_{\mathbf{h}_{i}}|-\frac{1}{2}\log|\widehat{\mathbf{K}}_{\boldsymbol{\theta}_{i}}|-\frac{1}{2}\log\widehat{K}_{h_{i,0}}-\widehat{\nu}_{i}\log\widehat{S}_{i}+\log\Gamma(\widehat{\nu}_{i})$.
Taking expectation of the above log ratio with respect to $q$, we obtain the
variational lower bound:

$\displaystyle\underline{p}(\mathbf{y}_{i};q)$
$\displaystyle=\mathbb{E}_{q}\left\\{\log\left[\frac{p(\mathbf{y}_{i},\mathbf{h}_{i}\,|\,\boldsymbol{\theta}_{i},h_{i,0},\sigma_{h,i}^{2})p(\boldsymbol{\theta}_{i},h_{i,0},\sigma_{h,i}^{2})}{q(\boldsymbol{\theta}_{i},h_{i,0},\sigma_{h,i}^{2},\mathbf{h}_{i})}\right]\right\\}$
$\displaystyle=c_{i}-\frac{1}{2}\mathbf{1}_{T}^{\prime}\widehat{\mathbf{h}}_{i}-\frac{1}{2}(\mathbf{y}_{i}-\mathbf{X}_{i}\widehat{\boldsymbol{\theta}}_{i})^{\prime}\widehat{\mathbf{O}}_{\mathbf{h}_{i}}(\mathbf{y}_{i}-\mathbf{X}_{i}\widehat{\boldsymbol{\theta}}_{i})-\frac{1}{2}\text{tr}(\mathbf{X}_{i}^{\prime}\widehat{\mathbf{O}}_{\mathbf{h}_{i}}\mathbf{X}_{i}\widehat{\mathbf{K}}_{\boldsymbol{\theta}_{i}}^{-1})-\frac{1}{2V_{h_{i,0}}}(\widehat{h}^{2}_{i,0}+\widehat{K}_{h_{i,0}}^{-1})$
$\displaystyle\quad-\frac{1}{2}(\widehat{\boldsymbol{\theta}}_{i}-\boldsymbol{\theta}_{i,0})^{\prime}\mathbf{V}_{\boldsymbol{\theta}_{i}}^{-1}(\widehat{\boldsymbol{\theta}}_{i}-\boldsymbol{\theta}_{i,0})-\frac{1}{2}\text{tr}(\mathbf{V}_{\boldsymbol{\theta}_{i}}^{-1}\widehat{\mathbf{K}}^{-1}_{\boldsymbol{\theta}_{i}})+\frac{1}{2}(T+k_{i}+1)$
$\displaystyle=c_{i}-\frac{1}{2}\mathbf{1}_{T}^{\prime}\widehat{\mathbf{h}}_{i}-\frac{1}{2}(\mathbf{y}_{i}-\mathbf{X}_{i}\widehat{\boldsymbol{\theta}}_{i})^{\prime}\widehat{\mathbf{O}}_{\mathbf{h}_{i}}(\mathbf{y}_{i}-\mathbf{X}_{i}\widehat{\boldsymbol{\theta}}_{i})-\frac{1}{2V_{h_{i,0}}}(\widehat{h}^{2}_{i,0}+\widehat{K}_{h_{i,0}}^{-1})$
$\displaystyle\quad-\frac{1}{2}(\widehat{\boldsymbol{\theta}}_{i}-\boldsymbol{\theta}_{i,0})^{\prime}\mathbf{V}_{\boldsymbol{\theta}_{i}}^{-1}(\widehat{\boldsymbol{\theta}}_{i}-\boldsymbol{\theta}_{i,0})+\frac{1}{2}(T+1),$

where
$\widehat{\mathbf{O}}_{\mathbf{h}_{i}}=\text{diag}(\text{e}^{-\widehat{h}_{i,1}+\frac{1}{2}d_{i,t}},\ldots,\text{e}^{-\widehat{h}_{i,T}+\frac{1}{2}d_{i,T}})$
and $k_{i}$ is the dimension of $\boldsymbol{\theta}_{i}$. Note that the
second equality in the above derivation follows from

$-\frac{1}{2}\text{tr}(\mathbf{X}_{i}^{\prime}\widehat{\mathbf{O}}_{\mathbf{h}_{i}}\mathbf{X}_{i}\widehat{\mathbf{K}}_{\boldsymbol{\theta}_{i}}^{-1})-\frac{1}{2}\text{tr}(\mathbf{V}_{\boldsymbol{\theta}_{i}}^{-1}\widehat{\mathbf{K}}^{-1}_{\boldsymbol{\theta}_{i}})=-\frac{1}{2}\text{tr}((\mathbf{V}_{\boldsymbol{\theta}_{i}}^{-1}+\mathbf{X}_{i}^{\prime}\widehat{\mathbf{O}}_{\mathbf{h}_{i}}\mathbf{X}_{i})\widehat{\mathbf{K}}_{\boldsymbol{\theta}_{i}}^{-1})=-\frac{1}{2}\text{tr}(\mathbf{I}_{k_{i}})=-\frac{1}{2}k_{i}.$

We summarize the variational Bayes iterative scheme below.

Algorithm 1 Iterative scheme for obtaining the parameters in the optimal
densities
$q^{*}(\boldsymbol{\theta}_{i},h_{i,0},\sigma_{h,i}^{2},\mathbf{h}_{i})$.

Initialize: $\widehat{\mathbf{K}}_{\boldsymbol{\theta}_{i}}$,
$\widehat{\boldsymbol{\theta}}_{i}$, $\widehat{\nu}_{i}$, $\widehat{S}_{i}$,
$\widehat{h}_{i,0}$. Set $\widehat{\nu}_{i}=\nu_{i}+\frac{T}{2}$.
Cycle:

$\displaystyle\widehat{\mathbf{K}}_{\mathbf{h}_{i}},\widehat{\mathbf{h}}_{i}$
$\displaystyle\xleftarrow[]{}\text{obtained by Newton-Raphson method}$
$\displaystyle\widehat{\mathbf{K}}_{\boldsymbol{\theta}_{i}}$
$\displaystyle\xleftarrow[]{}\mathbf{V}_{\boldsymbol{\theta}_{i}}^{-1}+\mathbf{X}_{i}^{\prime}\widehat{\mathbf{O}}_{\mathbf{h}_{i}}\mathbf{X}_{i}$
$\displaystyle\widehat{\boldsymbol{\theta}}_{i}$
$\displaystyle\xleftarrow[]{}\widehat{\mathbf{K}}_{\boldsymbol{\theta}_{i}}^{-1}(\mathbf{V}_{\boldsymbol{\theta}_{i}}^{-1}\boldsymbol{\theta}_{0,i}+\mathbf{X}_{i}^{\prime}\widehat{\mathbf{O}}_{\mathbf{h}_{i}}\mathbf{y}_{i})$
$\displaystyle\widehat{S}_{i}$
$\displaystyle\xleftarrow[]{}S_{i}+\frac{1}{2}\left[(\widehat{\mathbf{h}}_{i}-\widehat{h}_{i,0}\mathbf{1}_{T})^{\prime}\mathbf{H}^{\prime}\mathbf{H}(\widehat{\mathbf{h}}_{i}-\widehat{h}_{i,0}\mathbf{1}_{T})+\text{tr}(\mathbf{H}^{\prime}\mathbf{H}\widehat{\mathbf{K}}_{\mathbf{h}_{i}}^{-1})+\widehat{K}_{h_{i,0}}^{-1}\right]$
$\displaystyle\widehat{K}_{h_{i,0}}$
$\displaystyle\xleftarrow[]{}V_{h_{i,0}}^{-1}+\frac{\widehat{\nu_{i}}}{\widehat{S}_{i}}$
$\displaystyle\widehat{h}_{i,0}$
$\displaystyle\xleftarrow[]{}\widehat{K}_{h_{i,0}}^{-1}\frac{\widehat{\nu_{i}}}{\widehat{S}_{i}}\widehat{h}_{i,1}$

until the increase in $\underline{p}(\mathbf{y}_{i};q)$ is negligible.

## Appendix B: Additional Results

In this appendix we provide additional empirical results from the bank network
connectedness application.

We first report a range of bank connectedness measures using a large
homoscedastic VAR with the Minnesota prior described in Section 2. In our
implementation, the two shrinkage priors $\kappa_{1}$ and $\kappa_{2}$ are
selected by maximizing the variational lower bound over a 2-dimensional grid.
The optimal values obtained are $\kappa_{1}=0.04$ and $\kappa_{2}=0.001$,
indicating much stronger shrinkage toward zero for coefficients on ‘other’
lags than those on ‘own’ lags.

Following DDLY, instead of reporting individual bank connectedness measures,
we aggregate the network connectedness measures into six regions: Africa,
Asia, Europe, North America, Oceania and South America. The results are
reported in Table 4. These are the connectedness measures defined in (7) with
$H=10$, where each unit is a region rather than an individual bank.111111Note
that here we use a homoscedastic VAR, and consequently these connectedness
measures are time-invariant, i.e., $C_{i\leftarrow
j,1}^{H}=\cdots=C_{i\leftarrow j,T}^{H}.$ The row sums labeled ‘from others’
are the total directional connectedness from others defined in (8), and the
column sums labeled ‘to others’ are the total directional connectedness to
others defined in equation (9). Lastly, the lower right element is the system-
wide connectedness defined in (10).

Table 4: Bank network connectedness for the six-group aggregation, 2003-2014, from a homoscedastic VAR. | Africa | Asia | Europe | N. America | Oceania | S. America | From others
---|---|---|---|---|---|---|---
Africa | 0.00 | 7.54 | 22.45 | 22.11 | 2.12 | 2.39 | 56.62
Asia | 3.79 | 0.00 | 217.88 | 283.35 | 30.51 | 21.14 | 556.68
Europe | 4.93 | 67.35 | 0.00 | 734.41 | 32.60 | 33.35 | 872.64
N. America | 2.99 | 52.16 | 583.26 | 0.00 | 27.88 | 26.43 | 692.73
Oceania | 1.72 | 26.28 | 116.65 | 134.80 | 0.00 | 6.51 | 285.97
S. America | 1.22 | 12.07 | 50.02 | 54.00 | 2.83 | 0.00 | 120.14
To others | 14.66 | 165.40 | 990.26 | 1228.67 | 95.94 | 89.83 | 2584.77

Overall our results are qualitatively similar to those reported in DDLY,
despite their use of a different penalty (an adaptive elastic net). In
particular, our results also suggest that North America and Europe are the two
largest net transmitters of future volatility uncertainty (the difference
between ‘to others’ and ‘from others’) to the rest of the world. Moreover,
Asia has substantial total directional connectedness in both directions (i.e.,
total directional connectedness ‘to others’ and ‘from others’), and it is a
net receiver of future volatility uncertainty.

Next, Figure 7 plots the dynamic system-wide connectedness measure from the
96-variable homoscedastic VAR.

Figure 7: The dynamic system-wide connectedness measure from the homoscedastic
VAR using a 150-day rolling window.

The decomposition of the dynamic system-wide connectedness measure into cross-
country and within-country system-wide connectedness components is reported in
Figure 8.

Figure 8: The decomposition of the dynamic system-wide connectedness measure
from the homoscedastic VAR using a 150-day rolling window into cross-country
and within-country components.

Finally, we report the dynamic connectedness measures from the 96-variable VAR
with stochastic volatility using a 150-day rolling window in Figure 9.

Figure 9: The dynamic system-wide connectedness measure from the VAR with
stochastic volatility using a 150-day rolling window.

## References

  * (1)
  * Arias, Rubio-Ramirez, and Shin (2021) Arias, J. E., J. F. Rubio-Ramirez, and M. Shin (2021): “Macroeconomic forecasting and variable ordering in multivariate stochastic volatility models,” _Federal Reserve Bank of Philadelphia Working Papers_.
  * Banbura, Giannone, Modugno, and Reichlin (2013) Banbura, M., D. Giannone, M. Modugno, and L. Reichlin (2013): “Now-casting and the real-time data flow,” in _Handbook of Economic Forecasting_ , vol. 2, pp. 195–237. Elsevier.
  * Banbura, Giannone, and Reichlin (2010) Banbura, M., D. Giannone, and L. Reichlin (2010): “Large Bayesian vector auto regressions,” _Journal of Applied Econometrics_ , 25(1), 71–92.
  * Bhattacharya, Pati, Pillai, and Dunson (2015) Bhattacharya, A., D. Pati, N. S. Pillai, and D. B. Dunson (2015): “Dirichlet–Laplace priors for optimal shrinkage,” _Journal of the American Statistical Association_ , 110(512), 1479–1490.
  * Bishop (2006) Bishop, C. M. (2006): _Pattern Recognition and Machine Learning_. springer.
  * Bognanni (2018) Bognanni, M. (2018): “A class of time-varying parameter structural VARs for inference under exact or set identification,” _FRB of Cleveland Working Paper_.
  * Carriero, Clark, and Marcellino (2015) Carriero, A., T. E. Clark, and M. G. Marcellino (2015): “Bayesian VARs: Specification Choices and Forecast Accuracy,” _Journal of Applied Econometrics_ , 30(1), 46–73.
  * Carriero, Clark, and Marcellino (2016) (2016): “Common drifting volatility in large Bayesian VARs,” _Journal of Business and Economic Statistics_ , 34(3), 375–390.
  * Carriero, Clark, and Marcellino (2019) (2019): “Large Bayesian vector autoregressions with stochastic volatility and non-conjugate priors,” _Journal of Econometrics_ , 212, 137–154.
  * Carriero, Kapetanios, and Marcellino (2009) Carriero, A., G. Kapetanios, and M. Marcellino (2009): “Forecasting exchange rates with a large Bayesian VAR,” _International Journal of Forecasting_ , 25(2), 400–417.
  * Chan (2017) Chan, J. C. C. (2017): “The Stochastic Volatility in Mean Model with Time-Varying Parameters: An Application to Inflation Modeling,” _Journal of Business and Economic Statistics_ , 35(1), 17–28.
  * Chan (2020a) (2020a): “Large Bayesian VARs: A Flexible Kronecker Error Covariance Structure,” _Journal of Business and Economic Statistics_ , 38(1), 68–79.
  * Chan (2020b) (2020b): “Large Bayesian Vector Autoregressions,” in _Macroeconomic Forecasting in the Era of Big Data_ , ed. by P. Fuleky, pp. 95–125. Springer.
  * Chan (2021) (2021): “Asymmetric conjugate priors for large Bayesian VARs,” _Quantitative Economics_ , forthcoming.
  * Chan and Eisenstat (2018) Chan, J. C. C., and E. Eisenstat (2018): “Bayesian Model Comparison for Time-Varying Parameter VARs with Stochastic Volatility,” _Journal of Applied Econometrics_ , 33(4), 509–532.
  * Chan, Eisenstat, and Strachan (2020) Chan, J. C. C., E. Eisenstat, and R. W. Strachan (2020): “Reducing the State Space Dimension in a Large TVP-VAR,” _Journal of Econometrics_ , 218(1), 105–118.
  * Chan and Grant (2016) Chan, J. C. C., and A. L. Grant (2016): “On the Observed-Data Deviance Information Criterion for Volatility Modeling,” _Journal of Financial Econometrics_ , 14(4), 772–802.
  * Chan, Koop, Poirier, and Tobias (2019) Chan, J. C. C., G. Koop, D. J. Poirier, and J. L. Tobias (2019): _Bayesian Econometric Methods_. Cambridge University Press, 2 edn.
  * Chan, Koop, and Yu (2021) Chan, J. C. C., G. Koop, and X. Yu (2021): “Large order-invariant Bayesian VARs with stochastic volatility,” _Working Paper_.
  * Clark (2011) Clark, T. E. (2011): “Real-time density forecasts from Bayesian vector autoregressions with stochastic volatility,” _Journal of Business and Economic Statistics_ , 29(3), 327–341.
  * Cogley and Sargent (2005) Cogley, T., and T. J. Sargent (2005): “Drifts and volatilities: Monetary policies and outcomes in the post WWII US,” _Review of Economic Dynamics_ , 8(2), 262–302.
  * Cross, Hou, and Poon (2020) Cross, J., C. Hou, and A. Poon (2020): “Macroeconomic forecasting with large Bayesian VARs: Global-local priors and the illusion of sparsity,” _International Journal of Forecasting_ , 36(3), 899–915.
  * Cross and Poon (2016) Cross, J., and A. Poon (2016): “Forecasting structural change and fat-tailed events in Australian macroeconomic variables,” _Economic Modelling_ , 58, 34–51.
  * D’Agostino, Gambetti, and Giannone (2013) D’Agostino, A., L. Gambetti, and D. Giannone (2013): “Macroeconomic forecasting and structural change,” _Journal of Applied Econometrics_ , 28, 82–101.
  * Demirer, Diebold, Liu, and Yilmaz (2018) Demirer, M., F. X. Diebold, L. Liu, and K. Yilmaz (2018): “Estimating global bank network connectedness,” _Journal of Applied Econometrics_ , 33(1), 1–15.
  * Diebold and Yilmaz (2009) Diebold, F. X., and K. Yilmaz (2009): “Measuring financial asset return and volatility spillovers, with application to global equity markets,” _The Economic Journal_ , 119(534), 158–171.
  * Diebold and Yilmaz (2014) Diebold, F. X., and K. Yilmaz (2014): “On the network topology of variance decompositions: Measuring the connectedness of financial firms,” _Journal of Econometrics_ , 182(1), 119–134.
  * Durbin and Koopman (1997) Durbin, J., and S. J. Koopman (1997): “Monte Carlo maximum likelihood estimation for non-Gaussian state space models,” _Biometrika_ , 84, 669–684.
  * Ellahie and Ricco (2017) Ellahie, A., and G. Ricco (2017): “Government purchases reloaded: Informational insufficiency and heterogeneity in fiscal VARs,” _Journal of Monetary Economics_ , 90, 13–27.
  * Garman and Klass (1980) Garman, M. B., and M. J. Klass (1980): “On the estimation of security price volatilities from historical data,” _Journal of business_ , pp. 67–78.
  * Gefang, Koop, and Poon (2019) Gefang, D., G. Koop, and A. Poon (2019): “Variational Bayesian inference in large Vector Autoregressions with hierarchical shrinkage,” _CAMA Working Paper_.
  * Giannone, Lenza, and Primiceri (2017) Giannone, D., M. Lenza, and G. Primiceri (2017): “Economic predictions with big data: The illusion of sparsity,” _CEPR Discussion Paper No. DP12256_.
  * Giannone, Lenza, and Primiceri (2015) Giannone, D., M. Lenza, and G. E. Primiceri (2015): “Prior selection for vector autoregressions,” _Review of Economics and Statistics_ , 97(2), 436–451.
  * Griffin and Brown (2017) Griffin, J., and P. Brown (2017): “Hierarchical shrinkage priors for regression models,” _Bayesian Analysis_ , 12(1), 135–159.
  * Hajargasht and Woźniak (2018) Hajargasht, G., and T. Woźniak (2018): “Accurate computation of marginal data densities using variational Bayes,” _arXiv preprint arXiv:1805.10036_.
  * Harvey, Ruiz, and Shephard (1994) Harvey, A., E. Ruiz, and N. Shephard (1994): “Multivariate stochastic variance models,” _The Review of Economic Studies_ , 61(2), 247–264.
  * Jordan, Ghahramani, Jaakkola, and Saul (1999) Jordan, M. I., Z. Ghahramani, T. S. Jaakkola, and L. K. Saul (1999): “An introduction to variational methods for graphical models,” _Machine Learning_ , 37(2), 183–233.
  * Karlsson (2013) Karlsson, S. (2013): “Forecasting with Bayesian vector autoregressions,” in _Handbook of Economic Forecasting_ , ed. by G. Elliott, and A. Timmermann, vol. 2 of _Handbook of Economic Forecasting_ , pp. 791–897. Elsevier.
  * Kastner and Huber (2020) Kastner, G., and F. Huber (2020): “Sparse Bayesian vector autoregressions in huge dimensions,” _Journal of Forecasting_ , 37(7).
  * Koop (2013) Koop, G. (2013): “Forecasting with medium and large Bayesian VARs,” _Journal of Applied Econometrics_ , 28(2), 177–203.
  * Koop and Korobilis (2010) Koop, G., and D. Korobilis (2010): “Bayesian Multivariate Time Series Methods for Empirical Macroeconomics,” _Foundations and Trends in Econometrics_ , 3(4), 267–358.
  * Koop and Korobilis (2013) (2013): “Large time-varying parameter VARs,” _Journal of Econometrics_ , 177(2), 185–198.
  * Koop and Korobilis (2018) (2018): “Variational Bayes inference in high-dimensional time-varying parameter models,” _Available at SSRN 3246472_.
  * Korobilis and Yilmaz (2018) Korobilis, D., and K. Yilmaz (2018): “Measuring dynamic connectedness with large Bayesian VAR models,” _Available at SSRN 3099725_.
  * Loaiza-Maya, Smith, Nott, and Danaher (2020) Loaiza-Maya, R., M. Smith, D. Nott, and P. Danaher (2020): “Fast and Accurate Variational Inference for Models with Many Latent Variables,” _arXiv preprint arXiv:2005.07430_.
  * McCracken, Owyang, and Sekhposyan (2015) McCracken, M. W., M. Owyang, and T. Sekhposyan (2015): “Real-time forecasting with a large, mixed frequency, Bayesian VAR,” _FRB St. Louis Working Paper_ , 2015-30.
  * Morley and Wong (2020) Morley, J., and B. Wong (2020): “Estimating and accounting for the output gap with large Bayesian vector autoregressions,” _Journal of Applied Econometrics_ , 35(1), 1–18.
  * Ormerod and Wand (2010) Ormerod, J. T., and M. P. Wand (2010): “Explaining variational approximations,” _The American Statistician_ , 64(2), 140–153.
  * Primiceri (2005) Primiceri, G. E. (2005): “Time Varying Structural Vector Autoregressions and Monetary Policy,” _Review of Economic Studies_ , 72(3), 821–852.
  * Shephard and Pitt (1997) Shephard, N., and M. K. Pitt (1997): “Likelihood analysis of non-Gaussian measurement time series,” _Biometrika_ , 84, 653–667.
  * Shin and Zhong (2020) Shin, M., and M. Zhong (2020): “A new approach to identifying the real effects of uncertainty shocks,” _Journal of Business and Economic Statistics_ , 38(2), 367–379.
  * Tran, Nott, and Kohn (2017) Tran, M., D. J. Nott, and R. Kohn (2017): “Variational Bayes with intractable likelihood,” _Journal of Computational and Graphical Statistics_ , 26(4), 873–882.
  * Uhlig (1997) Uhlig, H. (1997): “Bayesian vector autoregressions with stochastic volatility,” _Econometrica_ , 65(1), 59–73.
  * West and Harrison (2006) West, M., and J. Harrison (2006): _Bayesian Forecasting and Dynamic Models_. Springer Science & Business Media.
  * You, Ormerod, and Mueller (2014) You, C., J. T. Ormerod, and S. Mueller (2014): “On variational Bayes estimation and variational information criteria for linear regression models,” _Australian & New Zealand Journal of Statistics_, 56(1), 73–87.
  * Zou and Zhang (2009) Zou, H., and H. Zhang (2009): “On the adaptive elastic-net with a diverging number of parameters,” _Annals of Statistics_ , 37(4), 1733\.