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

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

\onehalfspacing

\thispagestyle{empty}

Introduction

Since the influential work of \citet*{BGR10}, large Bayesian vector autoregressions (VARs) have been widely used to characterize the comovements of a large number of macroeconomic and financial variables.\footnote{Notable examples include \citet*{CKM09}, koop13, BGMR13, \citet*{CCM15}, \citet*{MOS15}, ER17 and MW20.} Since a vast empirical literature has demonstrated the importance of allowing for time-varying volatility in small systems,\footnote{See, for example, CS05, Primiceri05, clark11, DGG13, and CP16.} there is a lot of recent work that aims to develop stochastic volatility specifications for large VARs, including KK13, CCM16, CCM19, chan20, CES20 and KH20. 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 KK18 and GKP19, 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 {\em 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 {\em 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 TNK17 and LSND20, 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 DDLY18 --- 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 (ref), followed by some discussion of an adaptive Minnesota prior. Section (ref) 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 (ref). Section (ref) 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 (ref). Lastly, Section (ref) concludes and briefly discusses some future research directions.

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:

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:

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'\boldsymbol \Sigma_t^{-1}\mathbf{B}_0$, as considered in CS05 and CCM19.

Since the covariance matrix $\boldsymbol \Sigma_t$ in (ref) is diagonal, we can estimate this recursive system equation by equation without loss of efficiency. Below we first rewrite (ref) 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})'$ 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 (ref) 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}',\ldots, \mathbf{y}_{t-p}')$. 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:

where $\boldsymbol \theta_{i} = (\boldsymbol \alpha_{i}', \boldsymbol \beta_{i}')'$ is of dimension $k_i = np+i.$ We have therefore rewritten the VAR in (ref) 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 DLP15 and GB17, 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, \citet*{GLP17} and CHP20.

For a general discussion of the Minnesota prior, we refer the readers to KK10, karlsson13 or chan20b. 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\{

\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.\footnote{The 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 CCM15 and \citet*{GLP15}, 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 (ref) 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 chan21. Second, the multivariate stochastic volatility specification is constructed based on the lower triangular impact matrix $\mathbf{B}_0$. As noted by CCM19, 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})'$ is not order invariant. This issue arises in all stochastic volatility models that are based on a lower triangular parameterization, including the models in CS05 and Primiceri05. A few recent papers, such as Bognanni18, SZ20, ARRS21 and CKY21, 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.

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., Jordanetal99, Chap. 10 of Bishop06 and OW10. Recent applications in econometrics include HW18, KK18, GKP19 and LSND20.

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:

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 (ref), 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.

A New Approximating Density of the Stochastic Volatility

In this section we introduce a {\em global\/} Gaussian approximation of the joint distribution of the log-volatility $\mathbf{h}_i=(h_{i,1},\ldots, h_{i,T})'$ in contrast to {\em 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 (ref) 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, KK18 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.\footnote{This approach of approximating the stochastic volatility model can be traced back to \citet*{HRS94}, 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, GKP19 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.\footnote{This type of local Gaussian approximation is first introduced to estimate stochastic volatility models in the seminal papers by DK97 and SP97. 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. CG16b and CE18 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 KK18. 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:

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 (ref) 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)$.\footnote{One 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 GKP19, 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:

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 (ref) 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 (ref). 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 (ref), we can express $\log\left[\frac{f_{\mathbf{m}}(\mathbf{h}_i)}{\widetilde{q}^*_{\mathbf{h}_i}(\mathbf{h}_i)}\right]$ as

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)'$. 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'\mathbf{m} + (\widehat{\mathbf{s}}^2)'\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)'\mathbf{H}'\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 (ref) can be solved readily. Furthermore, the gradient and the Hessian of the objective function with respect to $\mathbf{m}$ can be computed easily:

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 chan17. 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}$.

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:

for $t=1,\ldots,T$, and we set $h_0=0$. For each dataset $\mathbf{z}^{(i)} = (z_1^{(i)},\ldots,z_T^{(i)})', i=1,\ldots, R$ and each Gaussian approximation with mean vector $\widehat{\mathbf{h}}^{j} = (\widehat{h}_{1}^{j},\ldots, \widehat{h}_{T}^{j})', 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 (ref) 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.

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 (ref). 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.

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 (ref). 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.

Next, we document the runtimes of estimating VARs of different dimensions using the proposed method versus MCMC. More specifically, Table (ref) 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.

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 DY09,DY14 and DDLY18. In particular, we revisit the application in DDLY, who consider a 96-variable homoscedastic VAR regularized by the adaptive elastic net ZZ09 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.\footnote{KY18 also use a time-varying VAR to measure bank network connectedness. Their model is based on the discounted Wishart process of Uhlig97 and WH06, 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 CS05 and Primiceri05. See, for example, ARRS21 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.

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.\footnote{The connectedness measures developed in DY09,DY14 and DDLY18 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 (ref)-(ref).

More specifically, given the estimates (e.g., posterior means) of the structural-form VAR in (ref), 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})'$. 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 DY14 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:

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:

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

Finally, we can measure the total directional connectedness as

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

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.\footnote{We 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 GK80. 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 (ref), 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 (ref). For comparison, we also include results based on the two local Gaussian approximations.

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 (ref) that show the proposed method delivers the lowest MSEs of the log-volatility.

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 (ref)--(ref) 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 (ref). 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 (ref), 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.

Next, we plot the dynamic system-wide connectedness in Figure (ref). 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.

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.\footnote{Another 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.

The dynamic system-wide connectedness measure presented in Figure (ref) 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 (ref) 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 (ref).

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 (ref). 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.

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.