Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
118,442 characters · 16 sections · 52 citation commands
Approximate Bayesian inference and forecasting in huge-dimensional multi-country VARs
\thispagestyle{empty}
\footnotetext{ We thank the editor Jesús Fernández-Villaverde and two anonymous referees for their valuable suggestions and constructive comments. We also thank Maximilian B\"ock, Todd Clark, Niko Hauzenberger, Ed Knotek, James Mitchell, Anna Stelzer, Saeed Zaman, and participants of the IAAE 2021 Annual Conference in Rotterdam. Florian Huber and Michael Pfarrhofer gratefully acknowledge financial support from the Austrian Science Fund (FWF, grant no. ZK 35) and the Jubiläumsfonds of the Oesterreichische Nationalbank (grant no. 18304). Please address correspondence to: Florian Huber. Department of Economics, University of Salzburg. Address: M\"{o}nchsberg 2a, 5020 Salzburg, Austria. Email: [email removed].}
\doublespacing
There is much evidence that working with multi-country time series models improves macroeconomic forecasting and structural analysis psw2004,pss2009,cc2016. This is to be expected in the modern globalized economy where countries are linked together through trade and financial flows and events in one country can spill over into others. However, the relevant data sets can be enormous. In this paper, we work with a $38$ country data set that contains $487$ variables. If all of these are treated as endogenous variables in an unrestricted multi-country vector autoregression (VAR), the number of equations in the VAR will be huge, as will the number of right hand side variables in each equation. The resulting model will be over-parameterized.
Much of the existing literature deals with this problem by imposing restrictions on the model or compressing the data (e.g., by using factor methods). For instance, the popular class of global VARs psw2004,pss2009,cfh2016,dfh2016,huber2016 assumes that information from all other countries impacts a country solely through a single weighted average of other country information. The weights in the average typically are based on bilateral trade flows. By contrast, the literature on panel VARs (PVARs) mainly deals with over-parameterization issues through constraints on the parameters describing the dynamic and static relations across countries canova2009estimating.
We propose an unrestricted PVAR specification where any variable can affect any or all other variables either contemporaneously or with a lag. This implies that the influence of foreign variables is neither a priori restricted (e.g., using trade weights) nor based on the output of a dimension-reduction procedure such as principal components. In other words, we want to let the data decide the exact nature and extent of linkages between countries. This feature of our approach turns out to be a substantial improvement over multi-country models commonly used in the literature.
Bayesian shrinkage methods are increasingly used for overcoming the over-parameterization problems which arise when working with unrestricted PVARs (or large VARs in general). Influential early contributions such as canova2009estimating used simple methods for choosing the prior (e.g., subjective elicitation or training sample methods). More recently, Bayesians have been working with global-local shrinkage or variable selection priors kk2016, korobilis2016prior, bai2022. These priors are commonly used in Big Data problems where models involve a large number of parameters. They automatically sort through all the parameters and decide which ones to shrink to zero and which ones to estimate freely. Our goal in this paper is to use a global-local shrinkage prior in a large unrestricted PVAR and let it decide which cross-country linkages are important and which can be ignored. In our empirical work, we mainly use the Horseshoe prior of horseshoe although the econometric methods developed in this paper will work with any hierarchical shrinkage prior that takes a conditionally Gaussian form including the LASSO and Dirichlet-Laplace priors of lasso and bhattacharya2015dirichlet, respectively.
Bayesian estimation and forecasting in VARs using global-local shrinkage priors typically requires the use of computationally demanding Markov Chain Monte Carlo (MCMC) methods. These algorithms are impractical in the very large PVARs that arise when working with many countries and many variables. For this reason, the existing Bayesian literature which uses unrestricted PVARs has focused on relatively small models. For instance, kk2016 use a PVAR involving four variables for each of seven European countries which is much smaller than the one considered in this paper.
To overcome the computational hurdle, we develop an Integrated Rotated Gaussian Approximation (IRGA) for the PVAR. IRGAs were recently developed by irga as a machine learning tool to speed up computation in high-dimensional models. These methods build on the intuition that some parameters are more important than other parameters. The other parameters, which in our case are mainly coefficients associated with other countries' lagged endogenous variables and covariance terms, are estimated using efficient approximations. The more important parameters are then estimated conditional on these approximations using precise MCMC techniques. We adapt these methods for use with PVARs. The resulting IRGA-based algorithm leads to vast improvements in speed of computation and has appealing approximation properties which we illustrate through simulations.
In our empirical work, we estimate a huge model of the world economy that contains 487 endogenous variables. In a forecasting exercise, we compare the performance of our unrestricted PVAR with a Horseshoe prior, estimated using IRGA methods (PVAR-IRGA), to a range of alternatives including a GVAR, factor-augmented VARs with the factors constructed using other country variables, and single country Bayesian VARs. In terms of computation, our key finding is that the computational improvements are large enough to enable Bayesian forecasting and structural analysis to be done even in huge PVARs. In terms of empirical results, we find our large approximate model forecasts well and often outperforms competing models. The forecast improvements are particularly strong for short-run density forecasts of stock market returns and longer-run inflation and output predictions. In these cases, we find PVAR-IRGA, with few exceptions, to forecast substantially better than all the alternatives.
We then proceed by analyzing the properties of the forecasts of the PVAR-IRGA via recent techniques used for analyzing social networks holland1983stochastic, karrer2011stochastic, vziberna2014blockmodeling, pati2015optimal. These stochastic block models build on the correlation matrix and sort correlations into clusters which facilitate interpretation. This analysis provides novel insights on the properties of the forecasts which are consistent with actual developments during the global financial crisis and the euro area sovereign debt crisis.
Our forecasting exercise is complemented by additional empirical results on the degree of cross-country spillovers. Considering different variants of the dy2009 spillover index shows that our model detects sizable international relations across output, prices, long-term interest rates and stock markets which sharply increase throughout the hold-out period. These increases are particularly pronounced during times of economic turmoil. Our findings hence provide novel insights on truly global connectivity since comparable large scale analyses are not feasible with standard econometric techniques.
The remainder of the paper is structured as follows. In Section (ref) we define the PVAR likelihood and the prior we use to carry out Bayesian inference and prediction. Section (ref) discusses computation and develops our IRGA methods which allow for fast computation. It also includes a theoretical discussion on the approximation properties of IRGA. In Section (ref) we carry out a simulation exercise to complement the theoretical discussion while in Section (ref) we present results for our forecasting exercise which compares our PVAR-IRGA to a range of other approaches. This section also includes information on the extent of cross-country spillovers. The final section summarizes and concludes the paper. The appendix provides additional technical details and further empirical results such as robustness checks.
This section develops the basic PVAR and briefly discusses the main specification issues commonly faced by researchers. We then discuss how Bayesian techniques can be used to deal with over-parameterization concerns and thus, in an automatic fashion, solve several of these specification issues.
Our goal is to model dynamic and static relations in an international panel of macroeconomic and financial time series which are stored in an $n$-dimensional vector $\bm y_t = (\bm y'_{1t}, \dots, \bm y'_{Nt})'$ for $t=1,\hdots,T$. This vector is composed of $N$ country-specific sub-vectors $\bm y_{it}$ which are $M \times 1$ dimensional.\footnote{Note that $M$ may differ across countries but is used here to simplify notation. Our approach naturally allows for different covariates across equations and countries. In the empirical work, variable coverage differs across countries.} Assuming that each $\bm y_{it}$ depends on the lagged values of $\bm y_t$, we obtain a PVAR given by:
where $\bm \Gamma_{ij}$ are $M \times M$ coefficient matrices associated with the lagged endogenous variables of country $i$. Lags of variables from countries other than $i$ are denoted by $\bm z_{it} = (\bm y'_{-i, t-1}, \dots, \bm y'_{-i, t-p})'$ with $\bm y_{-i, t}=(\bm y'_{1t}, \dots, \bm y'_{i-1, t}, \bm y'_{i+1, t}, \dots, \bm y'_{Nt})'$. The coefficient matrix on other country lags, $\bm \Xi_i$, is an $M \times K_{other}$ matrix where $K_{other} = (N-1)Mp$. Note that $\bm \Xi_i$ will contain an enormous number of parameters unless $N$ and/or $M$ are small. The matrix $\bm \Xi_i$ encodes the dynamic relations across countries (which are commonly referred to as dynamic interdependencies in the literature) while the $M \times k(=Mp)$ matrix $\bm \Gamma_i = (\bm \Gamma_{i1}, \dots, \bm \Gamma_{ip})$ captures domestic dynamics.
The usual VAR representation in terms of $\bm y_t$ is obtained by stacking all country-specific models and reshuffling the columns of $\bm \Gamma_i$ and $\bm \Xi_i$ appropriately:
The coefficient matrices $\tilde{\bm \Gamma}_j$ are of dimension $n \times n$ and the errors $\bm \epsilon_t = (\bm \epsilon'_{1t}, \dots, \bm \epsilon'_{Nt})'$ are i.i.d. Gaussian with $\bm \Sigma$ being an $n \times n$-dimensional variance-covariance matrix. The off-diagonal elements of this matrix determine both contemporaneous relations across variables within a country and instantaneous dependencies across countries. The latter relations are typically referred to as static interdependencies in the literature. Notice that unrestricted estimation of all contemporaneous relations across countries implies estimating $n (n-1) /2$ covariances. For large $n$, this adds to the already huge number of parameters in the $\bm \Gamma_{i}$'s and $\bm \Xi_i$'s.
In the literature on PVARs, estimation is often facilitated by introducing restrictions canova2009estimating, cc2016 on the coefficients in ((ref)) and ((ref)). For instance, the so-called cross-sectional homogeneity restriction arises if $\bm \Gamma_{i} = \bm \Gamma_{s}$ for $i \neq s$. This implies that domestic dynamics across countries are identical -- a rather restrictive assumption if the panel of countries includes, e.g., developed and developing economies. Another restriction often introduced is $\bm \Xi_i = \bm 0$ for some (or even all) $i$. This rules out dynamic relations across some countries but substantially reduces the number of free parameters. GVARs are also not restriction free. The assumption here is that cross-country linkages can be approximated by cross-country weighted averages and hence restrict $\bm \Xi_i$.\footnote{These weights have to be specified exogenously and are often based on measures of economic connectivity such as bilateral trade flows. For an overview, see Feldkircher2016a.} Another restriction sometimes considered assumes that shocks across countries are uncorrelated, i.e., $Cov(\bm \epsilon_{it}, \bm \epsilon_{st}) = \bm 0$. This implies introducing zero restrictions on the relevant elements in $\bm \Sigma$. All these restrictions, however, have serious implications for forecasting and structural inference and potentially introduce mis-specification if chosen wrongly. These considerations inspire us to use Bayesian variable selection methods via a global-local shrinkage prior so as to choose the appropriate restrictions in a data based manner.
If left unrestricted, estimation of the PVAR using traditional Bayesian MCMC methods is computationally cumbersome. For large data sets such as the one used in this paper, the computational burden becomes impractical. Hence, a goal of this paper is to speed up computation. In a first step, we greatly simplify computation by transforming the PVAR to allow for equation-by-equation estimation. This can be achieved by taking a Cholesky-type decomposition of $\bm \Sigma = \bm U \bm H \bm U'$. Here, we let $\bm U$ denote a lower triangular matrix with unit diagonal and $\bm H$ is a diagonal matrix with main diagonal $\bm \sigma^2 = (\bm \sigma_1^{2'}, \dots, \bm \sigma_N^{2'})'$. The $M$-dimensional vector $\bm \sigma_i^{2} = (\sigma_{\varepsilon, i1}^2, \dots, \sigma_{\varepsilon, iM}^2)'$ stores the idiosyncratic variances $\sigma_{\varepsilon, ij}^2$ associated with the shock in country $i$ and equation $j$. We can use this decomposition to recover the structural form of ((ref)):
where $\bm W = (\bm I - \bm U^{-1})$ encodes the contemporaneous relations across the shocks in the system and the matrices $\bm A_j~(j=1,\dots, p)$ denote structural coefficients. Within a given country, we can easily obtain a representation similar to ((ref)) by reshuffling the explanatory variables:
with $\bm W_i$ denoting the $M$ rows of $\bm W$ associated with the $i^{th}$ country. This allows us to rewrite the $j^{th} (>1)$ equation in the country-specific model $i (>1)$ as follows:\footnote{For $j=1$ and $i=1$, the equation simplifies and only $\bm x_{it}$ and $\bm z_{it}$ appear as regressors.}
where $y_{ij, t}$ denotes the $j^{th}$ element of $\bm y_{it}$, $\bm x_{it} = (\bm y'_{it-1}, \dots, \bm y'_{it-p})'$ while $\bm A_{ij, \bullet}$ and $\bm B_{ij, \bullet}$ denote the $j^{th}$ rows of $\bm A_i$ and $\bm B_i$, respectively. $\bm u_{ij, \bullet} = (w_{i1}, \dots, w_{i, j-1}, \bm u'_{i1}, \dots, \bm u'_{i i-1})'$ are the covariance parameters associated with the relevant row in $\bm W_i$.
Note that the errors are now independent across equations (i.e., both within and across countries) with error variance given by $\sigma^2_{\varepsilon, ij}$. This independence allows for estimating one equation at a time which greatly speeds up computation. Cholesky-based formulations such as this have been used in many recent papers, including ccm2019, kkp2019, hko2020 and carriero2021corrigendum. As opposed to carriero2021corrigendum, our approach includes the contemporaneous values of the endogenous variables and is thus not order invariant. In the appendix we show that the results are robust to different orderings of the countries in $\bm y_t$, implying only negligible empirical differences.
Equation ((ref)) is a simple regression model which regresses $y_{ij,t}$ on the lags of $\bm y_{it}$, the lags of the other countries' endogenous variables in $\bm z_{it}$, the contemporaneous values of the preceding $j-1$ variables domestic variables in country $i$ as well as the contemporaneous values of the preceding $i-1$ countries. Our approach builds on the notion that $\bm x_{it}$ is more important in explaining $y_{ij,t}$ than all other quantities in ((ref)).
To simplify notation, let $\tilde{\bm z}_{ij,t} = (\bm z'_{it}, y_{i1,t}, \dots, y_{ij-1,t}, \bm y'_{1t}, \dots, \bm y'_{i-1 t})'$ denote a $K_{ij} (= K_{other} + j-1 + (i-1)M)$ vector which stores the international quantities (both lagged and contemporaneously) as well as the time $t$ values of the endogenous variables up to equation $i$. Moreover, we stack the corresponding regression coefficients in a $K_{ij}$ vector $ \tilde{\bm B}_{ij, \bullet} = (\bm B'_{ij, \bullet}, \bm u'_{ij, \bullet})'$. Notice that $K_{ij}$ is much larger than $k$ which implies that $\tilde{\bm B}_{ij, \bullet}$ is difficult to estimate for large $M$, $N$ and $p$. We can rewrite ((ref)) in full-data form by stacking the $T$ observations into vectors to obtain the PVAR equations which define the likelihood function of our model:
with $\bm x_i$ being a $T \times k$ matrix where $k=Mp$ with $t^{th}$ row given by $\bm x'_{it}$. $\tilde{\bm z}_{ij}$ is a $T \times K_{ij}$ matrix with $t^{th}$ row $\tilde{\bm z}'_{ijt}$. This equation is a regression model which discriminates between a high-dimensional set of predictors related to covariances and international quantities in $\tilde{\bm z}_{ij}$ and a low-dimensional set of domestic quantities in $\bm x_i$.
The methods developed in this paper apply for any prior which has a hierarchical Gaussian form and thus leads to a full conditional posterior which is Gaussian. This is due to the fact that IRGA methods exploit the property that Gaussian distributions are invariant to rotations. A popular class of priors which has this form is the class of Gaussian global-local shrinkage priors. These can be represented as scale mixtures of Gaussians.\footnote{triplegamma provides a taxonomy of a range of priors in this class and discusses their properties.}
At a general level, consider the $j^{th}$ coefficient in a model, $\phi_j$. A global-local shrinkage prior can be written as:
where $f$ and $g$ are mixing densities and many different choices for them have been proposed. In a global-local shrinkage prior, $\lambda$ controls global shrinkage (common to all parameters). Having global shrinkage has often been found useful in Bayesian VARs (e.g., the Minnesota prior has a global shrinkage parameter) to reduce over-fitting concerns.\footnote{In Appendix (ref), we show how to implement a hierarchical Minnesota-type prior within this framework.} $\psi_j$ does local shrinkage (specific to the $j^{th}$ parameter). That is, if $\psi_j$ is estimated to be close to zero then $\phi_j$ is shrunk to be close to zero.
Suitably chosen mixing densities $f$ and $g$ result in a wide range of popular shrinkage priors such as the LASSO lasso, Normal-Gamma griffin2010inference, Dirichlet-Laplace bhattacharya2015dirichlet or the Horseshoe horseshoe. Due to its empirical success and ease of implementation, we use the Horseshoe prior. It takes the form:
whereby $\mathcal{C}^+$ denotes the half-Cauchy distribution. MS2016 show that the Horseshoe can be equivalently stated in terms of inverse Gamma distributions using suitable auxiliary variables. Specifically,
with $\mathcal{G}^{-1}$ denoting the inverse Gamma distribution and $\psi_j$, $\xi$ being auxiliary parameters. These auxiliary parameters are merely used to simplify posterior inference.
Up to this point, we have assumed that all parameters of the model are forced to zero through a single global shrinkage parameter $\lambda$. However, in the PVAR model we have different sets of parameters, countries, and variable types (equations within a given country). Specifying a single shrinkage parameter would imply that global shrinkage is symmetric across these different dimensions, a rather restrictive assumption. In the PVAR, we will assume that the $\lambda$ does not shrink all parameters towards zero but is specified to differ across countries, equations and types of parameters. This implies that for each coefficient vector $\bm A_{ij, \bullet}$ and $\tilde{\bm B}_{ij, \bullet}$ we will replace $\lambda$ with $\lambda_{A,ij}$ and $\lambda_{B,ij}$, respectively. In other words, each equation will have its own global shrinkage parameters and there will be two of them: one for own country coefficients and one for other country and contemporaneous coefficients.
Since $\tilde{\bm B}_{ij, \bullet}$ includes the covariance parameters as well, our discussion suggests that we use a Horseshoe also on the off-diagonal elements of $\bm U$. This implies that our prior allows for detecting whether static interdependencies across countries and variables are present and, if not, introduces shrinkage. We choose independent weakly informative inverse Gamma prior on the variances collected on the main diagonal of $\bm{H}$. In particular, $\sigma_{\epsilon,ij}^2\sim\mathcal{G}^{-1}(a_\sigma,b_\sigma)$ with $a_\sigma=b_\sigma=0.01$ for all countries and equations.
In this section we provide a framework that is capable of estimating huge PVAR models at reasonable computational cost. The next sub-section introduces IRGA to the PVAR case and provides some information on how posterior simulation can be carried out. Since this approach relies on approximating certain regions of the parameter space we then discuss the theoretical properties of the approximation.
As written in ((ref)), the PVAR simply involves $MN$ regression models. Bayesian MCMC methods for posterior and predictive inference in the regression model using a Horseshoe prior are well-established MS2016. In theory, we could simply use such methods with our PVAR. However, the problem is that MCMC methods are simply too slow for dealing with the high-dimensional parameter spaces that arise with PVARs. The main source of this high-dimensionality is that $\tilde{\bm B}_{ij, \bullet}$ (the coefficients on contemporaneous and other country variables in the equation for the $j^{th}$ variable for country $i$) potentially contains tens of thousands of coefficients. The matrix of own country coefficients, $\bm A_{ij, \bullet}$, is much smaller.
One empirical regularity often found in multi-country data sets is that own country effects are usually more important than other country effects. Hence, the literature sometimes sets $\bm B_i$ equal to a zero matrix to rule out dynamic interdependencies. This, however, could translate into a mis-specified model. In this paper, we allow the data to speak about the degree of sparsity in $\tilde{\bm B}_{ij, \bullet}$. A second empirical regularity is that the time series in $\tilde{\bm z}_{ij}$ often display substantial co-movement kose2003international. A potential solution would be to extract a low number of principal components from $\tilde{\bm z}_{ij}$. Another solution, which we adopt here, builds on the notion that if elements in $\tilde{\bm z}_{ij}$ are very similar to each other, it might pay off to not include all of them and effectively control for collinearity. This is also consistent with $\tilde{\bm B}_{ij, \bullet}$ being very sparse.
By contrast, $\bm A_{ij, \bullet}$ is likely non-sparse. This consideration motivates the way we implement IRGA with the PVAR. The general idea of IRGA is to use MCMC methods on important (low-dimensional) parameters but compute a fast approximation of the posterior for other (high-dimensional) parameters of less importance. In the PVAR, we consider $\bm A_{ij, \bullet}$ as the important parameters and $\tilde{\bm B}_{ij, \bullet}$ the less important ones.\footnote{It is worth noting that other choices of important and less important coefficients are possible. For instance, if it is felt likely that the U.S. is a dominant unit, then coefficients on U.S. variables could always be treated as important even for countries other than the U.S. The key consideration is that the number of parameters in $\bm A_{ij, \bullet}$ should not be too large so as to prevent practical use of MCMC methods on it.}
We now provide details of how we implement IRGA in the PVAR. Let $\bm Q_i$ be the $T \times T$ rotation matrix obtained from the QR-decomposition of $\bm x_i$ and partition it as $\bm Q_i = (\bm Q_{i1}, \bm Q_{i2})$ with $\bm Q_{i1}$ being $T \times k$ and $\bm Q_{i2}$ being $T \times (T-k)$.
Multiplying the equation for the $j^{th}$ variable in country $i$ by $\bm Q_i$ and exploiting the rotation invariance of the Gaussian distribution yields an equivalent representation of ((ref)):
The second equation follows since $\bm Q_{i2}' {\bm x}_{i} = \bm 0$. Note that $\bm A_{ij, \bullet}$ does not appear in it. This gives rise to a computational strategy which estimates $\tilde{\bm B}_{ij, \bullet}$ independently of $\bm A_{ij, \bullet}$. IRGA involves calculating the posteriors based on the two likelihood functions defined by ((ref)) and ((ref)). An approximate posterior for the (high-dimensional) $\tilde{\bm B}_{ij, \bullet}$ and $\sigma^2_{\varepsilon, ij}$ is obtained using ((ref)). Conditional on this approximate posterior $\hat{p}(\tilde{\bm B}_{ij, \bullet}|\bm Q'_{i2} \bm y_{ij})$, the posterior of $\bm A_{ij, \bullet}$ is obtained using MCMC based on ((ref)).
Any Gaussian approximation can be used for $\hat{p}(\tilde{\bm B}_{ij, \bullet}|\bm Q'_{i2} \bm y_{ij})$. We use vector approximate message passing (VAMP). Our choice of VAMP is driven by its scalability in huge dimensions and the fact that recent papers in machine learning and econometrics korobilis2019high have shown that it works extremely well for forecasting purposes. The specific implementation of the VAMP algorithm is the one proposed in rangan2019vector and details are given in Appendix (ref).\footnote{We have also experimented with variational Bayes to approximate $\hat{p}(\tilde{\bm B}_{ij, \bullet}|\bm Q'_{i2} \bm y_{ij})$. However, computation of the expected lower bound on evidence in these dimensions becomes computationally and numerically cumbersome. In addition, in our experiments the VAMP-based algorithm led to more precise predictions of the model.} The result is a Gaussian approximation: $\mathcal{N}(\overline{\bm B}_{ij, \bullet}, \overline{\bm V}_{ij, \bullet})$.
Rewriting ((ref)) and plugging in the approximate moments of $\hat{p}(\tilde{\bm B}_{ij, \bullet}|\bm Q'_{i2} \bm y)$ yields:
This gives us a Gaussian likelihood for a regression model which can be combined with any (conditionally) Gaussian prior on $\bm A_{ij, \bullet}$ leading to a textbook form for the posterior of $\bm A_{ij, \bullet}$ which can be estimated using MCMC methods. Additional details on the full conditional posterior distributions and how we approximate the error variances as well as the hyperparameters of the prior are given in Appendix (ref).
In this sub-section we briefly discuss why a Gaussian approximation is reasonable and how this approximation impacts the posterior distribution of $\bm A_{ij, \bullet}$. Intuitively, the accuracy of the approximate posterior of the domestic coefficients $\hat{p}(\bm A_{ij, \bullet}|\bm y_{ij})$ depends on the goodness of the approximation to ${p}(\tilde{\bm B}_{ij, \bullet}|\bm Q'_{i2} \bm y_{ij})$.
To investigate this relationship more formally, let $\text{KL}(p(a) || \hat{p}(a))$ denote the Kullback-Leibler (KL) divergence between an exact and an approximating distribution. irga show that the expected (with respect to the conditional distribution of $\bm y_{ij}$ given $\bm Q'_{i2}\bm y_{ij})$ KL divergence between $p(\bm A_{ij, \bullet}|\bm y_{ij})$ and $\hat{p}(\bm A_{ij, \bullet}|\bm y_{ij})$ is bounded from above by the approximation error to ${p}(\tilde{\bm B}_{ij, \bullet}|\bm Q'_{i2} \bm y_{ij})$.
Formally, irga establish a link between the approximation quality of $\hat{p}(\tilde{\bm B}_{ij, \bullet}|\bm Q'_{i2} \bm y_{ij})$ and how this impacts the approximate full conditional posterior $\hat{p}(\bm A_{ij, \bullet}|\bm y_{ij})$:
with $\mathcal{N}_{\sigma_{ij}} = \mathcal{N}(0, \sigma^2_{\varepsilon, ij} \bm I_k)$. This equation has three main implications. First, if the Gaussian approximation $\hat{p}(\tilde{\bm B}_{ij, \bullet}|\bm Q'_{i2} \bm y_{ij})$ is close to the exact full conditional posterior ${p}(\tilde{\bm B}_{ij, \bullet}|\bm Q'_{i2} \bm y_{ij})$ the corresponding conditional distribution $\hat{p}(\bm A_{ij, \bullet}|\bm y_{ij})$ will be close to $p(\bm A_{ij, \bullet}|\bm y_{ij})$. Second, the prior on $\bm A_{ij, \bullet}$ does not impact the error bound. Third, it does not depend on concentration properties around the true value of $\tilde{\bm B}_{ij,\bullet}$ which makes the result relevant for settings with $T$ being small.
In the next step, we justify our Gaussian approximation. Under certain mild assumptions on ${p}(\tilde{\bm B}_{ij, \bullet}|\bm Q'_{i2} \bm y_{ij})$ and $\tilde{\bm z}_{ij}$ and if $k \ll K_{ij}$, a multivariate central limit theorem implies that $\bm Q'_{i1} \tilde{\bm z}_{ij} \tilde{\bm B}_{ij, \bullet}$ is close to a Gaussian distribution diaconis1984asymptotics even if elements in $\tilde{\bm B}_{ij, \bullet}$ are non-Gaussian. This motivates a Gaussian approximation to ${p}(\tilde{\bm B}_{ij, \bullet}|\bm Q'_{i2} \bm y_{ij})$. irga, in Theorem 2, show that the expected KL divergence (with respect to $\bm Q'_{i2} \tilde{\bm z}_{ij}$) between the actual full conditional and the Gaussian approximation is bounded by two constants $\varpi_1$ and $\varpi_2$:\footnote{The constants $\varpi_1$ and $\varpi_2$ take a complicated form. To avoid introducing additional notation we summarize their main properties here. Precise definitions can be found in irga.}
$\varpi_1$ depends on the concentration properties of $p(\tilde{\bm B}_{ij, \bullet}|\bm Q'_{i2} \bm y_{ij})$ around its mean and $\varpi_2$ on the average correlation between the elements in $\tilde{\bm B}_{ij, \bullet}$.
If the posterior covariance of $\tilde{\bm B}_{ij, \bullet}$ is small relative to $\sigma^2_{\varepsilon, ij}$ and when $k \ll K_{ij}$, $\varpi_1$ approaches zero. The second constant measures approximation errors between the first two moments of ${p}(\tilde{\bm B}_{ij, \bullet}|\bm Q'_{i2} \bm y_{ij})$ and the approximating density. This quantity depends on the posterior covariance of the approximating density to the posterior of $\tilde{\bm B}_{ij, \bullet}$ and thus the error bound can be small even if the approximate variance-covariance matrix differs sharply from the true covariance of the posterior of $\tilde{\bm B}_{ij, \bullet}$. This finding also has important implications for our estimates of $\bm A_{ij, \bullet}$ since ((ref)) and ((ref)) can be combined to arrive at an upper bound for the approximation error to the posterior distribution of $\bm A_{ij, \bullet}$.
The theoretical discussion in the previous section builds on certain assumptions about the rows of $\tilde{\bm z}_{ij}$ and the correlation properties of the posterior of ${p}(\tilde{\bm B}_{ij, \bullet}|\bm Q'_{i2} \bm y_{ij})$. In this section, we will use synthetic data to illustrate how our approach performs under a realistic data generating process (DGP).
Our DGP assumes that $N=10, K=2, n=20$ and $T=500$ and features a single lag. It is given by:
where $ \bm H = \sigma^2_\varepsilon \times \bm I_{20}$ and $u_{ij} \sim \mathcal{N}(0, \sigma^2_\varepsilon/10) \text{ for } i=2, \dots, n; j=1, \dots, i-1$. The matrix $\hat{ \bm \Gamma}_1$ is obtained as follows. The blocks referring to the domestic coefficients are centred around the same mean vector and we add Gaussian shocks:
for all countries $i$. We simulate the elements in $\bm \Xi_{i}$ from $\mathcal{N}(0, \sigma_\beta^2)$. To obtain sparsity in $\bm \Xi_i$ we randomly zero out elements such that we have around $60$ percent zeroes in $\bm \Xi_i$. The initial value of $\bm y_t$ is sampled from a multivariate zero-mean Gaussian with variance $0.01$. To analyze how the goodness of our approximation changes with different measurement errors and cross-country heterogeneity, we consider $\sigma_\varepsilon, \sigma_\beta \in \{0.01, 0.025, 0.05\}$.
The main question is whether our IRGA-based approach yields estimates of the domestic coefficients which are close to the ones obtained from exact methods. Since exact methods become prohibitively slow in truly large data sets we make this comparison operational by simulating from a moderately large DGP and estimate a PVAR model with a Horseshoe prior and without using IRGA (i.e., all coefficients are estimated using MCMC). For both models we include $p=2$ lags of $\bm y_t$.
To avoid mixing up approximation errors arising from using IRGA to any errors in estimates that come from our equation-by-equation estimation approach, both the exact and approximate models are estimated based on the same (correct) ordering and in what follows we compare the coefficients $\bm A_{i1}$ and $\bm B_i$ to their true (implied) values.
Table (ref) shows averages of relative mean absolute error (MAE) ratios between the IRGA-PVAR and the PVAR estimated through MCMC across $100$ replications of the DGP. The white rows refer to MAE ratios for the domestic coefficients whereas the gray shaded rows denote MAE ratios for all regression coefficients. For both models we use the posterior median as our point estimator.
When we consider the results for the domestic coefficients (i.e., the white rows) we find a great deal of relative MAEs close to $1$ (except for cases in which both $\sigma_\beta$ and $\sigma_\varepsilon$ are very small). In principle, when we fix a given row and consider increasing values of $\sigma_\beta$ we find that the relative MAEs approach unity. This is consistent with the theoretical predictions in the previous section and shows that, at least when the posterior mean is considered, approximation accuracy of our IRGA-based approach increases with the ratio $\sigma_\beta / \sigma_\varepsilon$. When we fix a given column, as long as $\sigma_\varepsilon$ is not too small, we find no discernible differences across values of $\sigma_\beta$.
Turning to the relative MAEs across all coefficients reveals that our approximate method also yields estimates of $\bm B_i$ which are competitive to the ones obtained using MCMC. In fact, for $\sigma_\beta=0.05$ we find that the mean estimates essentially equal to the ones obtained from the MCMC-based method.
To illustrate our approach using a single draw from the DGP for $\sigma_\beta = \sigma_\varepsilon=0.05$, (ref) shows the marginal posterior distributions of the domestic coefficients for our IRGA-based PVAR (in solid black) and the MCMC-based estimates (in dashed black). This figure shows that in most cases, posterior distributions are similar. Especially when we focus on the mean/median we observe only small differences across coefficients (with some few cases suggesting a larger disagreement between MCMC and approximate estimates). When we focus on the higher moments of the marginal distributions we find similar variances, tail behavior and skewness properties. This small discussion has shown that, at least when synthetic data is considered, our approach yields reasonable estimates.
In this section we develop a huge-dimensional model of the world economy. The model is used to forecast output (measured by industrial production), inflation, long-term interest rates and stock prices for a large panel of countries. We moreover analyze the properties of the forecasts using novel stochastic block models as well as Diebold-Yilmaz (DY) spillover indices.
We have collected macroeconomic and financial data from the OECD's short-term indicator data base. The data are monthly and span the period from 2001m2 to 2019m12. In principle, we have a panel of 18 series for 38 OECD countries but not all variables are available with the same country coverage (see Table (ref)). In total we have $487$ variables in our PVAR.
The series fall into three categories: macroeconomic, financial and leading indicators. For macroeconomic variables, we consider measures of economic activity (industrial production growth, the output gap and the unemployment rate), export and import growth, consumer and producer price inflation as well as changes in earnings in the manufacturing sector.
Financial data cover short- and long-term interest rates (overnight, 3-months money market rates and long-term government bond yields), changes in stock prices and broad money growth. For euro area countries, we include the 3-months euribor as a measure of short-term interest rates.
In addition to this rather standard macro-financial data set, we gathered data on confidence/sentiment and short-term leading indicators. These comprise the OECD's leading indicator, which is a constructed measure to provide early signals of turning points in business cycles,\footnote{This measure shows fluctuations of an economy's activity around its long-term trend (i.e., potential output). The main difference to the output gap measure we employ is that the leading indicator is amplitude adjusted. For more information, see \href{http://www.oecd.org/sdd/leading-indicators/oecdcompositeleadingindicatorsreferenceturningpointsandcomponentseries.htm}{oecd.org/sdd/leading-indicators/oecdcompositeleadingindicatorsreferenceturningpointsandcomponentseries.htm}.} measures of manufacturers' and consumers' confidence, as well as changes in passenger car registration and newly permitted dwellings. Confidence measures potentially contain additional information to predict economic activity Batchelor1998, Ludvigson2004 as do car registrations and new dwellings for consumer expenditures. Changes or growth rates refer to either year-on-year or month-on-month growth rates. For a detailed overview, see Table (ref).
We carry out a forecasting exercise comparing an unrestricted PVAR with a Horseshoe prior to a range of alternatives. In what follows, we focus on predictions for four target variables: consumer price inflation (labeled Infl. in the tables and figures), industrial production (Ind. prod.), stock returns (Equities) and long-term interest rates (LT-IR).
Our forecasting design is recursive, implying that we use the data from 2001m2 to 2006m12 as an initial estimation period. We use data through 2006m12 to produce one-month up to twelve-months-ahead forecast distributions. This initial estimation period is then extended by data for 2007m1 and forecast densities for 2007m2 (up to 2008m1) are constructed. We repeat this procedure until we reach the end of the hold-out period.
To compare point forecasts across models we use MAEs. Since this disregards higher-order moments of the predictive distribution we also use Log Predictive Likelihoods (LPLs) to compare the density forecast performance of these alternatives.
The set of competing models is chosen not only to reflect a variety of popular approaches, but also to be computationally practical. In particular, alongside the proposed PVAR-IRGA approach, we consider the following models. First, to assess the role of allowing for cross-country spillovers, we estimate single-country Bayesian VARs (BVARs) that rule out both static and dynamic interdependencies. These are estimated country-by-country. Second, to compare our approach to other models incorporating international information, we consider two types of specifications. The first is a factor-augmented VAR (FAVAR-10) model, which augments the single-country BVARs with $10$ factors extracted from the non-domestic country variables. This procedure serves to obtain a lower dimensional representation of the international information set. As a second option to include international information, we use a Bayesian GVAR.
To ensure that differences in forecast performance are not driven by the respective priors on the VAR coefficients, we estimate all models with Horseshoe priors and set the number of lags equal to two. As a robustness check we also repeat the forecasting exercise replacing the Horseshoe prior with a conventional Minnesota prior for each of the PVAR-IRGA, BVAR, FAVAR-10 and the GVAR. Further specification details and results are given in Appendices (ref), (ref) and (ref).
Other than PVAR-IRGA, all of the models are estimated using exact MCMC methods. We stress that approaches which involve using MCMC methods for an unrestricted $487$ dimensional VAR would simply be computationally impractical. Computation times (average estimation time per model over all periods in the hold-out) for doing the pseudo real time forecasting exercise are provided in Table (ref). It can be seen that PVAR-IRGA is substantially faster than any of the competing approaches, even though all of the latter are much more parsimonious models. Notice that estimating the PVAR-IRGA is even faster than estimating a set of $N$ country-specific VARs of the same size. This is because MCMC sampling of the domestic coefficients based on the IRGA posterior is faster since it only relies on a part of the likelihood function to form the conditional posterior distributions.
The results of our forecasting exercise are summarized in Table (ref). It presents absolute values of MAEs and LPLs for our PVAR-IRGA approach (rows shaded in red) for two forecast horizons. Results for the other approaches are benchmarked relative to these. To be precise, for MAEs we take ratios relative to PVAR-IRGA (with numbers exceeding unity implying that the PVAR improves upon the competitors), for LPLs we take differences relative to PVAR-IRGA (with values below zero suggesting that the PVAR is outperforming the respective model). {The numbers in the table are GDP-weighted averages of country-specific LPLs computed by using GDP in 2015 U.S. dollars averaged over the period 2002 to 2019. To provide a rough gauge of model performance across countries, the numbers in parentheses represent the percentage of countries in which a given model performs best in absolute terms.}
The most important finding is that PVAR-IRGA works -- it produces sensible forecasts quickly. Bayesian estimation of huge dimensional PVARs has been made possible through the use of IRGA methods. {The other main finding is that (with some exceptions) PVAR-IRGA works well and is highly competitive with competing approaches. These improvements are limited for point forecasts but sometimes very pronounced for LPLs which measure density forecast performance. For short-term forecasts of equity returns in particular, PVAR-IRGA is producing strong improvements in LPLs while its performance is slightly weaker for the remaining variables under consideration.}
{When we focus on higher-order forecasts the relative performance of PVAR-IRGA improves. While we observe that predictive accuracy is deteriorating for equities at the twelve-months-ahead horizon, forecasts of industrial production and inflation improve considerably. For the latter two variables, the PVAR-IRGA is the single best performing model.}
BVAR is the only alternative that does not allow for any cross-country spillovers. On average, it forecasts fairly well for both forecast horizons. Improvements, however, are more more pronounced at the one-month-ahead horizon. For twelve-months-ahead, we find that taking cross-country linkages into account helps forecast accuracy. This suggests that for short-run forecasts, cross-country spillovers are not that strong (or at least do not significantly help in predicting output, inflation and long-term rates).
The fact that the single country BVARs improve upon the GVAR indicates that having a smaller-sized sparse model seems to be more important than taking into account cross-country linkages for improving forecasts. The good performance of our PVAR-IRGA suggests that the Horseshoe prior successfully strikes a balance between exploiting cross-country information and sparsity. Since the model is unrestricted it also allows the data to speak about the precise form of cross-country linkages.
The statements in the preceding paragraphs are based on an examination of LPLs. Analyzing MAEs reveals similar patterns, but to a weaker extent. This indicates that the benefits of unrestricted modeling of the high-dimensional PVAR offers some benefits in terms of point forecast performance, but the benefits are larger for density forecasts.
{The discussion has focused on averages across countries (and time periods). However, it could be that a model (such as our PVAR-IRGA) yields lower average LPLs but still provides the best performance for individual countries in our sample. Considering the percentage of wins for each model corroborates the findings based on LPLs and MAEs. But it is worth emphasizing that even though the PVAR-IRGA is sometimes outperformed by simpler competitors such as BVAR, for most variables we still find a sizable fraction of wins across countries for our proposed model. In the case of one-month-ahead inflation density forecasts, this share is about 34 percent whereas it is around 47 percent for longer-run forecasts of equity returns (in which BVAR outperforms the PVAR-IRGA if we consider LPLs). These sizable shares suggest that it might be worthwhile to carefully analyze country-specific results.}
In the preceding sub-section, we compared the average (over countries and time) forecast performance of PVAR-IRGA to various alternatives. In this sub-section, we look behind the average performance to investigate forecast performance at the country level and see how it changes over time. We do so through heatmaps of cumulative LPLs for PVAR-IRGA for the individual countries. {To assess when it pays off to allow for cross-country linkages (and because it is the strongest competitor to our approach) we benchmark our results to the LPLs of the BVAR.}
Figures (ref) and (ref) contain these heatmaps for the two forecast horizons. The figures are grouped into four categories: Advanced European, Emerging European, Advanced Other and Emerging Other and individual countries are labeled using ISO country codes. Intensifying shades of blue (red) indicate stronger support for PVAR-IRGA (BVAR).
{We first focus on one-month-ahead forecasts. Starting with inflation predictions, we observe that the stronger performance of the BVAR is mainly driven by a weak performance of the PVAR-IRGA in emerging economies (most notably Brazil and Turkey, with some exceptions, such as India) while it performs best for developed economies such as the U.S., France and Ireland. This brief discussion shows why it is important to also consider results at the country level. If the researcher is interested in short-term forecasting of U.S. inflation and has to decide on one of the models we consider, focusing on overall LPLs masks the particularly strong U.S.-specific forecasting performance.}
{When we focus on industrial production we find a somewhat different pattern, with PVAR-IRGA being outperformed by the BVAR for several countries in Advanced Europe (except for Austria, Finland, Norway and the Netherlands) and some gains in several countries located in Emerging Europe (Slovakia, Russia and the Baltics). For some countries (such as Sweden and Italy) the PVAR-IRGA performs well prior to the global financial crisis. The rapid decline in output in the final half of 2008, however, led to a deterioration in forecast performance. This is because PVAR-IRGA yields predictive distributions which are sometimes too tight and thus capturing outliers becomes increasingly difficult.}
{Turning to the results for equity returns reveals a great number of blue-colored cells. In principle, our model works well for most economies (with a slightly weaker performance for, e.g., Portugal, Norway, Slovenia and Japan). Interestingly, we also find some heterogeneity with respect to model performance over time. In the case of the U.S., for instance, our model only improves upon the BVAR from 2012 onward. We conjecture that the slightly weaker performance prior and during the financial crisis is, again, driven by too tight predictive intervals. But these tight intervals then help in predicting returns after the financial crisis, a period characterized by steady increases in U.S. stock markets.}
{The PVAR-IRGA displays the weakest performance for long-term interest rates. For some few countries (e.g., Portugal, Greece, Latvia and Israel) the PVAR performs well. In general, the weak performance for long-term rates is driven by the fact that these display a downward trend during the hold-out period for most countries. PVAR-IRGA captures this downward trend rather well but the predictive variance is considerably smaller than the one of the single-country BVAR. Hence, under the predictive distribution of the PVAR-IRGA, even relatively small changes in long rates have strong effects on LPLs. The countries which depart from this general pattern (such as Portugal and Greece) feature large spikes in long-term interest rates. The PVAR captures this well and quickly adjusts the predictive variance. Since it takes slightly longer for the BVAR to adjust we conjecture that this quick increase is mostly driven by the large information set.}
{In the previous sub-section we have shown that on average, the PVAR-IRGA produces the most precise density forecasts for inflation and industrial production when the forecast horizon is increased while the performance for equity returns deteriorates. When we focus on twelve-month-ahead predictions (see (ref)) we find that the strong overall performance for inflation is mostly driven by excellent forecasts in major developed economies such as the U.S. or Japan. A similar pattern is found for industrial production. Again, we find that the PVAR-IRGA produces precise density forecasts for most developed economies located in Europe as well as the United States and Japan. However, it is also worth stressing that the PVAR also produces accurate forecasts for developing economies such as Turkey as well as several countries located in Central Eastern Europe. This strong performance is driven by more precise point forecasts but also by the fact that the predictive distributions for multi-step-ahead forecasts seem to be heavy tailed and thus make observing outliers more probable.}
{For equities and long-term interest rates we find that the PVAR performs slightly weaker than the BVAR. Especially for equities, this is driven by a particularly bad performance in the U.S. For long-term interest rates we again find that the PVAR is competitive when used to forecast long-term rates in Portugal and Israel but it appears to be outperformed in countries such as Denmark.}
{In summary, PVAR-IRGA yields precise equity return predictions for short-term forecasts and shows good performance when used to forecast inflation in major economies such as the U.S. For higher-order forecasts, the results somewhat reverse and the PVAR-IRGA works particularly well when it is used to predict industrial production and inflation. In general, and this is consistent with the findings based on average LPLs, we find that more predictive evidence in favor of cross-country spillovers increases with the forecast horizons. We will provide additional evidence on the increasing importance of cross-country spillovers in the following sub-sections.}
In the previous section we have shown that our PVAR yields highly competitive forecasts and often improves upon other single- or multi-country models. In this section, our aim is to quantitatively analyze the country-specific forecasts to investigate the role of cross-country correlations in the point predictions. For brevity, we focus on one-month-ahead forecasts for industrial production and long-term interest rates.
Since the corresponding correlation matrices are high-dimensional and thus difficult to interpret, we use techniques from network analysis vziberna2014blockmodeling to search for clusters in the correlation matrices of the posterior median and the posterior standard deviation of the forecasts.\footnote{These are implemented through the R package blockmodeling.} Intuitively speaking, we reorder the rows and columns such that correlations between countries are grouped into $R=8$ distinct blocks.\footnote{The choice of 8 blocks is arbitrary. For our application, it yields a good balance between a too granular and a too coarse approach.} The relations (i.e., positive correlations) within a given block are maximized whereas the relations of countries within a block to economies outside of a block are much less important (or even negative).
Since this algorithm needs a correlation matrix of forecasts, we compute the initial correlation matrix based on the first 12 observations and then expand this window until we reach the end of the hold-out period which is used to compute forecast distributions. This yields a sequence of correlation matrices which we then analyze using a stochastic block model.
In what follows, we focus on one-month-ahead forecasts for three distinct periods in our hold-out sample. First, we examine correlation structures among our forecasts for 2009m1, the onset of the global financial crisis. The second period we consider is 2012m6, the month prior to Mario Draghi's famous “whatever it takes” speech, which marks the height of the euro area sovereign debt crisis. The final period is 2019m12, the end of our sample. Using all available information is a natural choice to investigate how forecasts are related. Figure (ref) shows the correlation matrices (multiplied by 10) for industrial production forecasts sorted using a stochastic block model.
The figure shows the $R=8$ clusters on the main diagonal of the matrix. The first two blocks are each defined by a single country (Norway and Portugal, respectively). Forecasts between the blocks are (modestly) negatively correlated, as indicated by the red shading in the off-diagonal elements between the blocks. Considering the remaining countries, we see that forecasts for Norway are not only negatively correlated with those for Portugal, but also with the rest of the sample. As an oil-based economy, Norway was considerably less affected by the global financial crisis than the rest of the sample. Looking at the remaining economies reveals two large clusters. One of them contains forecasts for Greece, Italy, Spain and the Baltics -- countries that showed massive contractions in output during the global financial crisis. The Irish economy, which also significantly contracted during the global financial crisis, appears in another cluster. Another large cluster contains countries that were comparably less affected by the crisis. Importantly though and with the exception of Norway, all clusters are positively correlated indicated by the grey to black colour shading in Figure (ref). This reflects the global nature of the financial crisis.
Considering the period of the euro area sovereign debt crisis, depicted in the upper right panel of Figure (ref) reveals a very similar picture: Norway stands out, and forecasts of the remaining countries are positively correlated. One big cluster emerges which covers mostly European economies and Russia -- the latter which shares strong trade ties with the European Union. Other European countries that were not affected by the sovereign debt crisis, such as the Baltics, are allocated into different blocks and Ireland again appears in a separate cluster. Turning now to one-month ahead forecasts for the end of the hold-out sample, we still find isolated countries such as Norway and Ireland, two European clusters and one large, international cluster. The latter one shows only a modest within-correlation, indicating that at the end of our sample, with the pandemic yet to fully unfold, forecasts are less strongly correlated as in periods of severe downturns.
Figure (ref) shows the same analysis for long-term interest rate forecasts. During the period of the global financial crisis, long-term rates spiked in countries like Greece, but have been downward-trending in countries considered as safe havens: the United States, Canada and Japan. PVAR-IRGA forecasts are consistent with these historical observations putting the aforementioned country groups into separate blocks. Two further clusters emerge: one solely consisting of European economies, while the other one contains both Advanced European and non-European economies. At the height of the sovereign debt crisis, long-term interest rates shot up for most economies but to a different extent. This is mirrored in the first four correlation blocks depicted on the main diagonal of the correlation matrix. Forecasts for crisis-stricken economies, such as Ireland, Italy and Spain are strongly correlated and appear in a separate block. Interestingly, Greece and Portugal appear in one block, but forecasts between these clusters are positively correlated. Both blocks of crisis stricken countries are clearly separated from the rest of the sample. Last and looking at the end of sample period, we see a very different picture. With the exception of South Africa and to some extent Greece, global long-term interest forecasts are very homogeneous and positively correlated
Summing up, examining correlation structures of one-month ahead forecasts revealed insights as to which extent the model is capable of mirroring correlation structures present in the data. Correlations are strongest during periods of simultaneous contractions such as witnessed during the global financial crisis. Looking at the episode of the sovereign debt crisis, which unfolded on a more regional basis, also reveals clusters but they differ from those identified during the global financial crisis. This highlights the overall flexibility of the model. For some countries forecasts are always separated from the rest of the sample (e.g., Norway as an oil-exporting economy). This implies that the PVAR-IRGA can take cross-country links into account when they are important but at the same time does not enforce them on the whole set of countries -- a flexibility which is of ample importance when dealing with large, heterogeneous cross-sections.
Another way of looking at cross-country correlations is to use the Diebold-Yilmaz (DY) spillover index dy2009. This index distinguishes to which extent forecast error variance can be explained by its own history as opposed to effects through all other variables in the system. The latter effects are dubbed "spillovers" and serve as a measure of overall connectivity. For instance, in this paper, we calculate the index based on a generalized forecast error variance decomposition (GFEVD) which avoids order dependence with respect to the elements in $\bm y_t$. We report the total share of effects from non-domestic variables (spillovers from one variable to another) in the GFEVD. The calculation of the underlying GFEVD is recursive (in the same manner as was done for the forecasting exercise), and we show results for a forecast horizon of $12$ and $24$ months. Our focus on higher-order spillovers is motivated by evidence reported in the literature on GVARs Feldkircher2016a which shows that cross-country spillovers (measured through FEVDs) become sizable only after several quarters.
{Figure (ref) shows the DY index over the hold-out period for both horizons. Most importantly, our results point at a sizable degree of cross-country connectedness. At the end of our sample period, the DY indices amount to about 45 at the twelve-months-ahead forecast horizon and to 68 percent at the 24-months-ahead horizon. Investigating the DY index over time, reveals a steady increase of the index until 2016 after which the indicator levels out.\footnote{Note that consistent with the forecasting exercise, we use an expanding window to calculate the DY index which automatically introduces a certain degree of persistence. Appendix (ref) provides additional empirical results for estimates using a rolling window of observations.} A particularly pronounced increase of the index can be observed during the period of the euro area sovereign debt crisis (between 2010 and 2012). In general, economic variables tend to co-move more strongly during turbulent times Pham2021 -- a pattern that we also observe with our data. It is worth stressing that our estimates are surrounded by considerable posterior uncertainty, which tends to attribute more posterior mass above than below the posterior median. This behavior is even more pronounced during the period of the sovereign debt crisis.}
The discussion up to this point focused on an overall measure of cross-country connectivity considering all focus variables jointly. To investigate whether connectivity plays a larger role for certain variables, we display the DY index for each of the focus variables separately in Figure (ref).
The figure reveals some interesting variable-specific differences. For example, the evolution of the DY index for industrial production, inflation and long-term interest rates is very similar to the behavior of the overall index displayed in Figure (ref). The degree of connectedness, is however, comparably smaller (35 percent). It is worth stressing that the posterior distribution is again strongly tilted towards higher levels of connectivity. A high degree of cross-country dependence between inflation rates is in line with Borio2007, Ciccarelli2010, km2018 who stress the importance of a global component in determining domestic inflation rates.
Spillovers between equity prices show a distinct dynamic. The degree of connectivity in financial markets is generally higher compared to that of the remaining variables. The associated DY index is about 25 percent at the beginning of the hold-out period and rises sharply during the global financial crisis. This finding is in line with demirer2018estimating, who demonstrate a strong increase in equity connectedness between banks during periods of financial stress. At the end of the sample period, the DY index amounts to about 75 percent. Notably, there is also little difference between the two forecast horizons.
In Figure (ref), we repeat the calculation and measure the role of cross-country effects for each country. In general, most of the countries considered display sustained increases in their respective DY indices peaking at about 75 percent. In some emerging economies, such as China and Indonesia, the index increases sharply right after 2008. Some notable exceptions are Turkey and Mexico which experienced a decrease in connectivity after the global financial crisis surrounding the Fed's tapering statement in June 2013. In the case of Greece, we observe a gradual decline in connectivity during the euro area sovereign debt crisis. This finding reflects a strong domestic component which complies with the fact that Greece was the epicenter of the crisis at the time.
{Summing up and jointly considered with the forecasting results, the cross-variable DY indices paint a similar and consistent picture. Benefits from estimating a large multi-country system tend to increase with the forecast horizon. This finding is mirrored by increasing levels of connectivity if higher-order forecasts are considered. Our results indicate that the PVAR is capable of, in light of heavy shrinkage introduced through the Horseshoe prior, capturing cross-country relations flexibly. As opposed to other models (such as the FAVAR or the GVAR) the PVAR-IRGA introduces no particular restrictions on the coefficients. This feature is crucial to appropriately control for cross-country heterogeneity in light of a diverse set of countries such as the one we have in our data set.}
Appendix (ref) provides additional results for the spillover indices computed using rolling windows of varying length. Using a rolling instead of an expanding window implies that past observations do not impact the estimates at some point. Thus, parameters are quicker to adjust to new information (at the cost of discarding past information). Results for both the overall and the variable-specific DY indices tell a similar story to the ones based on an expanding window. The main difference is that the estimates feature more movements in periods of economic turmoil such as the global financial crisis and the sovereign debt crisis in Europe.
Multi-country VARs have the potential to be enormous and simply working with unrestricted versions of them leads to over-parameterization and computational problems. The existing literature typically deals with these problems by imposing restrictions or reducing the dimension of the data. But the former strategy risks mis-specification and the latter risks losing information. Accordingly, in this paper, we have developed Bayesian methods for working with unrestricted VARs and rely on the Horseshoe prior to gain in parsimony by imposing shrinkage in a data-based fashion. Existing Bayesian work with VARs with such a prior has been done using MCMC methods. These are too computationally demanding to be used with PVARs with hundreds of dependent variables. In this paper, we have used IRGA methods to overcome this computational hurdle. We show that these allow for practical inference even in PVARs of huge dimension. Our macroeconomic empirical application demonstrates the benefits of being able to work with such large PVARs.
\addcontentsline{toc}{section}{References}