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.
112,594 characters · 21 sections · 59 citation commands
Gaussian Process Vector Autoregressions and Macroeconomic Uncertainty
\affil[a]{ University of Salzburg} \affil[b]{ Bocconi University, IGIER, CEPR, Baffi-Carefin and BIDSA} \thispagestyle{empty}\linespread{1.5}
\doparttoc \faketableofcontents \part
Economic relations can change over time for a variety of reasons, such as technological progress, institutional changes, major policy interventions, but also wars, terrorist attacks, stock market crashes and pandemics. Standard econometric models, such as linear single and multivariate regressions assume instead stability of the parameters characterizing the conditional first and second moments of the dependent variables. When stability is formally tested, it is often rejected stock1996evidence. This has led to the development of a variety of methods to handle structural change in econometric models.
Parameter evolution is assumed to be either observable (i.e., driven by the behavior of observable economic variables) or unobservable, and either discrete and abrupt or continuous and smooth. Examples include threshold and smooth transition models tong1990non, terasvirta1994specification, Markov switching models hamilton1989new, and double stochastic models nyblom1989testing. In all these models, a specific and fully parametrized type of parameter evolution is assumed, and then linear or non-linear filters are used for estimation in a classical context or Markov chain Monte Carlo (MCMC) methods in a Bayesian framework. Examples of economic applications of all these methods include koop2013large, accm2017changing, aastveit2017economic, caggiano2017estimating, alessandri2019financial, and caggiano2021uncertainty.\footnote{A special mention is due to primiceri2005time who popularized the use of time-varying parameters and stochastic volatility in macroeconometrics.}
Assuming a specific type of parameter evolution increases estimation efficiency but can lead to mis-specification. A more flexible alternative allows for a smooth evolution of parameters without specifying the form of parameter time variation. In a classical context, the evolution can be either deterministic robinson1991time, chen2012testing, or stochastic giraitis2014inference, giraitis2018inference, kapetanios2019large. Kernel estimators are the main tool used in this literature. Alternative approaches, which can also capture non-linear relationships between the target and explanatory variables, are, e.g., based on using functional-coefficient regressions cai2000functional, kowal2017bayesian, regression trees chipman2010bart, huber2020nowcasting, coulombe2020macroeconomy, neural networks hornik1989multilayer,gu2021autoencoder, coulombe2022neural {or infinite mixtures hirano2002semiparametric, bassetti2014beta, kalli2018bayesian, billio2019bayesian, jin2022infinite. Most of these approaches, however, are fairly different from the VAR models that are the workhorse of modern time series econometrics, making interpretation of the estimation results and computation of quantities such as impulse response functions difficult. In addition, they typically focus on extending the specification of the conditional mean, while assuming a constant conditional variance, which can be restrictive for macroeconomic and financial data. Finally, some of the methods do not scale well into high dimensions and are thus not particularly suited for large datasets nowadays used in macroeconomics.}
In this paper, we propose a new model that belongs to the non-parametric class and is capable of capturing, in a very flexible way, not only parameter evolution but also general non-linear relationships. Our model can be applied in a large data context while keeping the flexibility and ease of use of VARs{, and allowing for time-varying conditional variances}. Specifically, we combine the statistical literature on Gaussian process (GP) regressions crawford2019variable, with that on VARs to obtain a GP-VAR model. Borrowing ideas from the literature on Bayesian Minnesota-type VARs, the model assumes, for each endogenous variable, a different non-linear relationship with its own lags and with the lags of all the other variables (and possibly of additional exogenous regressors). Gaussian processes are used to model non-linearities in a flexible but efficient way. They can be viewed as a non-parametric alternative to the adaptive Minnesota-type shrinkage proposed in chan2021ijof. {They are also similar to neural networks, in the sense that they are universal approximators based on infinite mixtures of Gaussian distributions.\footnote{{In fact, specific choices of the kernels underlying Gaussian processes can produce a variety of neural network models, see novak2018bayesian for details. }}}
We develop an efficient (Bayesian) estimation procedure, based on the structural form of the GP-VAR, which has the additional benefit that its complexity is linear in the number of endogenous variables and does not depend on the number of lags. Hence, estimation can be parallelized and, in addition, the conjugate structure of the model we use allows for pre-computing various matrix multiplications and kernel operations, which further speeds up computation. As a result, estimation is feasible also for very large models.
As in all Bayesian procedures, an assumption on the distribution of the errors is required. As is common in the Bayesian VAR literature, we assume that the errors are Gaussian. Yet, we also permit the variance of the errors to change over time, adopting a stochastic volatility (SV) specification. Already in linear BVARs, the use of SV permits to have time variation in the conditional distribution of the variables, so much so that BVAR-SV are empirically a good alternative to quantile regressions carriero2022specification. Moreover, multi-step predictive densities, used to produce forecasts and generalized impulse responses, are non-Gaussian. The use of SV in the GP-VAR adds flexibility, and prevents overfitting in the sense of avoiding that large realizations of the shocks are interpreted as changes in the conditional mean.
We illustrate the GP-VAR with synthetic data generated from a highly non-linear multivariate data generating process (DGP) that have both Gaussian and non-Gaussian shocks. This DGP assumes that some equations feature parameters that exhibit structural breaks while others depend non-linearily on the lags of the endogenous variables. To assess whether the GP-VAR is also capable of recovering linear relations, one equation is a standard linear regression model. In all these cases, our approach works reasonably well in terms of detecting the nature of non-linearities.
Our GP-based model has a vast range of applicability, for both reduced form and structural analysis. This mirrors the possible applications of standard VARs but allows for much more general dynamic relationships across the variables. Besides the evaluation with synthetic data, we consider a forecasting application and a more structural economic analysis. In the forecasting application, we compare the performance of the GP-VAR with other linear and non-linear competitors that allow for parameter change, and with the BVAR with SV. Predicting US output, inflation and interest rates, we show that the GP-VAR improves upon all competing models, with gains that are particularly pronounced at the four-quarter-ahead horizon.
As an example of a structural economic analysis, and to gather new insights on a topic that has recently attracted considerable attention bloom2014fluctuations, we use the GP-VAR to investigate the effects of exogenous uncertainty shocks on US macroeconomic and financial time series. Comparing the responses of the GP-VAR with the ones of a standard linear BVAR reveals that our model produces sensible responses for real activity and stock markets. Differences between the responses relate to the shape and magnitudes, with the GP-VAR producing stronger reactions of uncertainty, real GDP growth, and stock markets returns while yielding similar reactions of employment growth. Considering models of differing sizes shows that the impulse responses do not differ markedly across model sizes.
Our proposed framework naturally allows for analyzing potential asymmetries in transmission channels. The responses to a positive uncertainty shock (higher unexpected uncertainty) are typically much stronger than the ones to a negative shock. Interestingly, the shape of the IRFs also differ, with positive shocks leading to responses that peak later. In addition, our findings suggest that the relationship between real activity and uncertainty becomes proportionally slightly smaller for large shocks, while financial markets react relatively more strongly to larger increases in uncertainty. Finally, our framework also allows us to investigate whether transmission mechanisms have changed over time. Doing so reveals that the effects of uncertainty have been smaller in the great inflation period ($1970$Q$1$ to $1984$Q$4$), more pronounced through the great moderation ($1985$Q$1$ to $2006$Q$4$) before again turning more muted during the post great moderation phase ($2007$Q$1$ to $2019$Q$4$).
The paper is structured as follows. Section (ref) provides an introduction to Gaussian process regression. In Section (ref) we develop the GP-VAR model. This section also provides necessary details on the prior setup and posterior computation. We then analyze model performance using synthetic data in Section (ref). Sections (ref) and (ref) include our empirical work. In the former section we briefly discuss the dataset and provide some in- and out-of-sample model evidence. In the latter, we focus on the macroeconomic implications of an uncertainty shock. The final section briefly summarizes and concludes the paper. The Online Appendix contains additional empirical results and technical details on the specification, estimation and use of the GP-VAR model.
In this section we briefly discuss Gaussian process (GP) regressions with a focus on time series data.\footnote{For a textbook treatment, see williams2006gaussian.} GP regressions are a non-parametric technique to establish a flexible relationship between a scalar time series $y_t$ and a set of $K$ predictors $\bm x_t$ in period $t$. The key advantage of this approach is that it does not rely on parametric assumptions on the precise functional relationship between $y_t$ and $\bm x_t$.
In general, a non-parametric regression is given by:
with $f$ being some unknown regression function $f: \mathbb{R}^K \to \mathbb R$ and $\varepsilon_t$ denoting an independent Gaussian shock with zero mean and constant variance $\sigma^2$. We relax this assumption in Sub-section (ref) to allow for heteroskedastic shocks. An assumption on the error distribution is needed in a Bayesian context and Gaussianity is the most common one, though different distributions can be easily accommodated by exploiting a scale-location mixture of Gaussians representation escobar1995bayesian.
In standard regression models, the function $f$ is assumed to be linear with $f(\bm x_t) = \bm \beta' \bm x_t$ where $\bm \beta$ is a $K \times 1$ vector of linear coefficients. If mean relations are non-linear, this assumption might be too restrictive. To gain more flexibility one can embed the covariates in $\bm x_t$ into a higher dimensional space such as the space of powers $\bm x_t \to \psi(\bm x_t) =(\bm x'_t, (\bm x^2_t)', \dots, (\bm x^R_t)')'$, with $\bm x^2_t = (\bm x_t \odot \bm x_t)$ and higher orders defined recursively. Conditional on choosing a sufficiently large integer $R$, this would provide substantial flexibility to approximate any smooth function $f$. However, adequately selecting $R$ is key and the mapping, moreover, is ad-hoc in the sense that there exist infinitely many non-linear mappings $\psi$.
Standard Bayesian methods place a prior on the coefficients associated with the covariates (and possible non-linear transformations thereof) and thus control for uncertainty with respect to these basis functions but at the cost of remaining within a class of functions (such as linear, polynomial or trigonometric functions). By contrast, in GP regressions we treat the function $f$ as an unknown quantity and let the data decide about the appropriate form (and degree) of non-linearities.
The key inferential goal in GP regression is to infer the function $f$ from the data under relatively mild assumptions. This is achieved by specifying a prior on $f(\bm x_t)$. A typical assumption is to assume that $f(\bm x_t)$ follows a Gaussian process prior:
with $\mu(\bm x_t)= \mathbb{E}[f(\bm x_t)]$ being the mean function and
denoting a kernel (or covariance) function that determines the relationship between $f(\bm x_t)$ and $f(\bm x_\tau)$ for periods $t$ and $\tau$. The kernel is typically parameterized by a low dimensional vector of hyperparameters $\bm \vartheta$ and controls the behavior of the function $f$. This kernel needs to be positive semidefinite and symmetric.
In what follows, we will set the function $\mu(\bm x_t) = 0$ for all $t$. This is without loss of generality, since any explicit basis function for $\mu(\bm x_t)$ can be used to model the mean process. If the focus is on modeling stationary data, $\mu(\bm x_t) = 0$ implies that a priori the process is centered around a white noise process. In case one would like to model persistent or non-stationary data it would be straightforward to implement a prior that forces the system towards a set of random walk processes. This can be achieved by setting $\mu(y_{t-1}) = \rho y_{t-1}$, where $\rho$ denotes a persistence parameter with prior mean $\mathbb{E}[\rho] = 1$. Alternatively, one could specify the prior on $f$ to imply persistence in $y_t$. This possibility is discussed in much more detail in Section (ref) of the Online Appendix.
A common choice in GP regressions is the Gaussian (or squared exponential) kernel function:
with $\xi$ denoting a scaling parameter and $\kappa$ the (inverse) length scale and thus $\bm \vartheta = (\xi, \kappa)'$. Larger values of $\kappa$ lead to a GP which displays more high frequency variation whereas lower values imply a slowly varying mean function. The parameter $\xi$ controls the prior variance of the function $f$. To see this, note that if $\bm x_t = \bm x_\tau$, we obtain $\text{Var}[f(\bm x_t)] =\xi$.
This specification is quite flexible and fulfills several convenient conditions. For instance, williams2006gaussian show that the use of the Gaussian kernel implies that $f(\bm x_t)$ is mean square continuous and differentiable. Moreover, this kernel function represents a positive semidefinite and symmetric covariance function. Furthermore, Mercer's theorem mercer1909xvi, under this kernel, states that the GP regression can be written in terms of an infinite number of basis functions. These basis functions are Gaussians with different means and variances. This suggests a connection to the literature on Bayesian non-parametrics escobar1995bayesian, neal2000markov,kalli2018bayesian, fruhwirth2019here that relies on infinite mixtures of Gaussians to estimate unknown densities. The link between GPs and infinite mixtures of Gaussians shows that the Gaussian assumption on $\varepsilon_t$ is not too restrictive as the model allows to recover non-Gaussian features in the data.
The GP prior represents an infinite dimensional prior over the space of functions. This implies that the estimation problem is infinite dimensional as well. However, since we sample data in a discrete manner, the GP prior becomes a multivariate Gaussian prior on $\bm f = (f(\bm x_1), \dots, f(\bm x_T))'$:
with $\bm 0_T$ being a $T \times 1$ vector of zeros, $K_{\bm \vartheta}(\bm X, \bm X)$ a $T \times T$ kernel matrix with typical element $k_{\bm \vartheta}(\bm x_t, \bm x_\tau)$ and $\bm X = (\bm x_1, \dots, \bm x_T)'$. This implies that, in terms of full data matrices, the GP regression is given by:
where $\bm I_T$ denotes a $T \times T$ identity matrix.
Assuming for the moment that $\sigma^2$ is known, the posterior of $\bm f$ follows a multivariate Gaussian distribution:
with variance-covariance matrix $\overline{\bm V}_{\bm f}$ and posterior mean vector $\overline{\bm f}$:
The mean function $\overline{\bm f}$ can be interpreted as a weighted average of the values of the endogenous variable:
where $\bm \alpha = (\alpha_1, \dots, \alpha_T)' = (K_{\bm \vartheta}(\bm X, \bm X) + \sigma^2 \bm I_T)^{-1}\bm y$. This (finite dimensional) representation shows how one moves from an infinite dimensional problem to a finite dimensional one.
The expression for the variance-covariance matrix $\overline{\bm V}_{\bm f}$ also has an intuitive interpretation. The first term is the prior variance (i.e., the kernel matrix). The second term measures how much of the variance is expressed through the covariates in $\bm X$ and thus the posterior covariance indicates how much the model learns from $\bm X$.
The predictive distribution of $f(\bm x_{T+h})$ can be easily derived by exploiting basic properties of the multivariate Gaussian:
whereby
Similar to the posterior mean $\overline{\bm f}$, the predictive mean $\overline{f}_{T+h}$ is a weighted average of the values of the endogenous variables $\bm y$ with the weights depending on the relationship between $\bm X$ and a realization of the vector of covariates $\bm x_{T+h}$ related to the $h$-step-ahead horizon. The predictive variance $\overline{V}_{T+h}$, again, depends on a term that is purely driven by the prior evaluated at $\bm x_{T+h}$ minus a term that measures the informational content in the covariates.
Before proceeding to the discussion on how to set the kernel it is worth noting that what we have discussed above is often labeled the function-space view of the GP. This is because the prior is elicited directly on $f$. Another way of analyzing GPs is based on the weight-space view. Under the weight-space view one can rewrite the GP regression as a standard regression model as follows:
with $\bm W_{\bm \vartheta}$ denoting the lower Cholesky factor of $K_{\bm \vartheta}(\bm X, \bm X) = \bm W_{\bm \vartheta} \bm W'_{\bm \vartheta}$ and $\bm \eta$ is a Gaussian shock vector with zero mean and unit variance. This is a standard regression model with $T$ regressors, a coefficient vector $\bm \eta$ and a Gaussian prior on $\bm \eta$. Standard textbook formulas for the Bayesian linear regression model koop2003bayesian can be used to carry out posterior inference.
This also shows that if we set $K_{\bm \vartheta}(\bm X, \bm X) = \bm X \bm V_{\bm \vartheta} \bm X'$, we obtain a linear regression model that features a Gaussian prior with zero mean and a typical prior variance-covariance matrix $\bm V_{\bm \vartheta}$. This kernel implies many more parameters than the parsimony inducing Gaussian kernel. Hence, the resulting fit and forecasts can be expected to have posterior distributions with larger variances than those associated with the Gaussian kernel, though bias would be lower if the true model is linear and features stable parameters.
In the previous sub-section the quantities for the posterior of $\bm f$ and the predictive density for future values of $f(\bm x_{T+h})$ suggest that the kernel and its hyperparameters play an important role. In this sub-section, we discuss this issue in more detail.
One of the key advantages of GPs is that by constructing suitable kernels, one can determine the space of possible functions. This gives rise to substantial flexibility and allows for capturing a large range of competing models within a single econometric model. For instance, williams2006gaussian discuss how kernels can be constructed to mimic the behavior of neural networks, regression splines, polynomial and linear regressions. Tree-based techniques such as Bayesian additive regression trees chipman2010bart can be cast in this framework by exploiting the ANOVA-representation of the model and then the weight-space view of the GP. In principle, and we will build on this feature later, summing over the corresponding kernels gives rise to another kernel and suitable weights could be constructed to select, in a data-driven way, which model summarizes the data best.
As stated in the previous sub-section, our focus will be on the Gaussian kernel due to its excellent empirical properties and analytical tractability. The two hyperparameters $\kappa$ and $\xi$ control the curvature and the marginal variance of the function, respectively. We illustrate the effect of $\kappa$ on the prior and posterior of $\bm f$ in Figures (ref) and (ref) by means of two simple univariate examples. The first example models quarterly US inflation (in year-on-year terms) and sets $x_t=t$ for periods ranging from $2005$Q$1$ to $2015$Q$4$. The second example models US GDP growth as a function of the first lag of a macroeconomic uncertainty measure for the same sub-sample.\footnote{Throughout the paper, we use the macroeconomic uncertainty measure of jurado2015measuring provided (and regularly updated) on the web page of Sydney C. Ludvigson (available online via \href{https://www.sydneyludvigson.com/macro-and-financial-uncertainty-indexes}{sydneyludvigson.com/macro-and-financial-uncertainty-indexes}). Detailed information on this index and the econometric techniques used to obtain this measure can be found in jurado2015measuring.} The figures then show (for both the prior and posterior) the value of the function $f(x_t)$ on the y-axis and $x_t$ on the x-axis. Both figures display in the left (right) panel the $5^{th}$ and $95^{th}$ prior (posterior) percentiles (with the area, the $90\%$ credible set, in between shaded in light red) as well as three random draws from the prior (dashed red lines) in the left panel, and the posterior median (solid red lines) in the right panel.
(ref) reveals that if $y_t$ is a (possibly non-linear) function of time and the inverse length scale parameter is set small, the model generates functions that track the trend in inflation rather well. This is similar to the unobserved components model of stock2007has, which features a persistent stochastic trend in inflation. Once we increase $\kappa$ we observe that the draws from the prior display more high frequency variation with shorter cycles between peaks and troughs. Once this prior is combined with the data, the estimated mean functions display much more curvature and fit the actual data increasingly well.
Once $\kappa$ is set too high, the functions arising from the prior vary substantially and are likely to capture also very small deviations of inflation from its trend. This translates into a close-to-perfect fit of the posterior mean of the functions and gives rise to serious overfitting concerns.
To see how a GP regression captures a possibly non-linear relationship between $y_t$ and $x_t$, (ref) shows the functional relationship between output growth and lagged macroeconomic uncertainty. If $\kappa = 0.01$, the regression relationship is almost linear and suggests that high levels of (lagged) macroeconomic uncertainty are accompanied by negative output growth rates.
When we set $\kappa=0.1$ we observe much more curvature (both in the prior and the posterior) in the relationship, indicating that if uncertainty is between $0$ and around $1.7$, GDP growth is between 2.3 and 2.5 percent. However, once a certain threshold in the first lag of uncertainty is reached, the relationship becomes strongly negative until it becomes essentially flat for very high levels of uncertainty. A similar finding, but slightly more pronounced, arises if we set $\kappa = 4$. In this case GDP growth does not change much as long as lagged uncertainty is between 0 and 1.7 and then the relationship becomes, again, strongly negative.
As is clear from these stylized examples, the role of the kernel and its hyperparameters crucially impacts the posterior estimates of the function $f$. Setting $\kappa$ too small leads to a model which might miss important (higher frequency) information whereas a $\kappa$ set too large translates into an overfitting model which might yield a very strong in-sample fit but poor ouf-of-sample predictions. Setting $\kappa$ is thus of crucial importance and in all our empirical work we will infer it through Bayesian updating.
Another key question is whether the estimated function converges to the true underlying function. The literature deals with this question using several assumptions on the error distributions (mostly setting $\sigma^2 = 0$) or how the GP regression behaves if the underlying function $f$ differs in terms of smoothness from the GP prior controlled by the kernel stone1982optimal, van2008rates, yang2017frequentist, teckentrup2020convergence. stone1982optimal, by focusing on iid data, shows that the optimal rate of estimation of a $\zeta-$smooth function is $T^{- \zeta/(2\zeta + K)}$ and thus decreases in $K$ while it increases in the smoothness of the true function. Building on this finding, teckentrup2020convergence analyzes the contraction properties of a Gaussian process regression under a general Mat\'{e}rn kernel function and provides error bounds that also depend on the relationship of the smoothness of the true and estimated functions. If these agree, one can achieve a convergence rate of $T^{-\zeta/K}$.
After having provided the necessary foundations on Gaussian process regression, we will now focus on developing a model that is suitable for macroeconomic analysis.
In this section, we first develop the GP-VAR in Sub-section (ref). Next, Sub-sections (ref) to (ref) are devoted to the development of efficient MCMC schemes to carry out posterior and structural inference. Finally, Sub-section (ref) details how to compute forecasts and (generalized) impulse response functions for the GP-VAR.
In the following discussion, let $\bm y_t = (y_{1t}, \dots, y_{Mt})'$ denote an $M \times 1$ vector of macroeconomic and financial variables.\footnote{We assume that the elements in $\bm y_t$ are demeaned. In our empirical application we include a constant term with an uninformative prior.} Moreover, $\bm x_t = (\bm x'_{1t}, \dots, \bm x'_{Mt})'$ denotes an $Mp \times 1$ vector with $\bm x_{jt} = (y_{jt-1}, \dots, y_{jt-p})'$ storing the “own" lags of the $j^{th}$ endogenous variable and $\bm z_t =(\bm z'_{1t}, \dots, \bm z'_{Mt})'$ an $(M-1) M p \times 1$ vector of “other" lags. Hence, $\bm z_{jt}=(\bm y'_{-j t-1}, \dots, \bm y'_{-j t-p})'$, where $\bm y_{-j t}$ denotes the vector $\bm y_t$ with the $j^{th}$ element excluded.
We discriminate between own and other lags of $\bm y_t$ because we assume that lags of other endogenous variables impact a given endogenous variable differently from its own lags. The literature on Bayesian VARs banbura2010large, koop2013forecasting, huber2019adaptive, chan2021ijof has captured this through shrinkage priors that treat coefficients on own and other lags differently. We wish to capture this equation-specific asymmetry by specifying our GP-VAR to depend on two latent processes: one driven by $\bm x_t$ and one by $\bm z_t$. The structural form of the resulting GP-VAR is then given by:
with $F(\bm x_t) = (f_1 (\bm x_{1t}), \dots ,f_M (\bm x_{Mt}) )'$ and $G(\bm z_t) = (g_1 (\bm z_{1t}), \dots ,g_M (\bm z_{Mt}) )'$ and $f_j$ and $g_j$ being equation-specific functions. The function $f_j$ controls how $y_{jt}$ depends on its own lags while $g_j$ encodes the relationship between $y_{jt}$ and the lags of the other endogenous variables. The functions $f_j$ and $g_j$, and hence $F$ and $G$, differ because we construct different kernels with distinct hyperparameters.\footnote{It is worth stressing that one could also think of our decomposition in terms of a new function with a kernel that is given by the sum of the kernels of the functions $f_j$ and $g_j$.} The matrix $\bm Q$ is an $M \times M$ lower triangular matrix with zeros along its main diagonal. This matrix defines the contemporaneous relations across the elements in $\bm y_t$.
Finally, $\bm \varepsilon_t$ is an $M \times 1$ vector of Gaussian shocks with zero mean and an $M \times M$ time-varying variance-covariance matrix $\bm H_t = \text{diag}(\omega_{1t},\dots,\omega_{Mt})$. We will assume that $\omega_{jt}$ follows a flexible stochastic volatility (SV) model:
with the logarithm of $h_{jt} = \log \omega_{jt}$ being assumed to evolve according to a stationary AR($1$) state equation. We let $\rho_{hj}$ denote the persistence parameter, $\sigma^2_{hj}$ the error variance, and $h_{j0}$ the initial state of the log-volatility process.
Allowing for time variation in the shock variances provides additional flexibility and enables us to capture non-Gaussian features in the shocks (not only, but also due to the fact that $h_{jt}$ enters the model non-linearly).\footnote{One could also introduce additional scaling factors that arise from inverse Gamma distributions to obtain a model with t-distributed shocks.} In principle, we could also allow for unknown functional relations between the contemporaneous terms of the preceding $j-1$ equations and the response of equation $j$. However, this would lead to a complicated non-linear covariance structure. Since we are interested in carrying out structural identification based on zero impact restrictions we opt for choosing this simpler approach which implies multivariate Gaussian reduced form shocks, but with a time-varying covariance matrix. Given that the literature on GPs typically assumes the shocks to be Gaussian and homoskedastic, this is already a substantial increase in flexibility.\footnote{A rare exception is jylanki2011robust, who propose a GP regression with heavy tailed errors and mainly focus on fast and robust approximate inference of a posterior that is analytically intractable due to a $t$-distributed likelihood.}
The model in (ref) assumes that the shocks in $\bm \varepsilon_t$ are, conditional on $\bm Q \bm y_t$, orthogonal and hence estimation can be carried out equation-by-equation. We will exploit this representation for simplicity and computational tractability. The $j^{th}$ equation, in terms of full-data matrices, is given by:
with $\bm Y_j =(y_{j1}, \dots, y_{jT})', \bm f_j = (f_j(\bm x_{j1}), \dots, f_j(\bm x_{jT}))', \bm g_j = (g_j(\bm z_{j1}), \dots, g_j(\bm z_{jT}))'$, $\bm \Omega_j = \text{diag}(\omega_{j1}, \dots, \omega_{jT}), \bm \epsilon_j = (\varepsilon_{j1}, \dots, \varepsilon_{jT})'$ and $q_{jk}$ denoting the $(j,k)^{th}$ element of $\bm Q$. We will use this form to carry out inference about the unknown functions $f_j$ and $g_j$ as well as the remaining parameters and latent states of the model.
Notice that our estimation strategy is not invariant with respect to reordering the elements in $\bm y_t$, a common problem if this orthogonalization strategy is used. In Sub-section (ref) of the Online Appendix, we show that different orderings have only a small impact on the estimated impulse responses.
In this sub-section, our focus will be on the priors on $f_j$ and $g_j$. The priors on the remaining, linear quantities are standard and thus not discussed in depth. We use a Horseshoe prior carvalho2010horseshoe on the free elements in $\bm Q$, a Beta prior on the (transformed) persistence parameter $(\rho_{h j}+1)/2 \sim \mathcal{B}(25, 5)$, and an inverse Gamma prior on the state innovation variances $\sigma^2_{h j}$. This prior is specified to have mean $0.1$ and variance $0.01$.
For equation-specific functions $f_j$ and $g_j$, we specify two GPs with one conditional on $\bm X_j = (\bm x_{j1}, \dots, \bm x_{jT})'$ and one conditional on $\bm Z_j = (\bm z_{j1}, \dots, \bm z_{jT})'$:
We let $\sqrt{\bm \Omega_j} = \text{diag}(\sqrt{\omega_{j1}}, \dots, \sqrt{\omega_{jT}})$ while $K_{\bm \vartheta_{j1}}(\bm X_j, \bm X_j)$ and $K_{\bm \vartheta_{j2}}(\bm Z_j, \bm Z_j)$ denote two suitable kernels with typical elements given by:
For $j=1, \dots, M$, $\bm \vartheta_{j1}$ and $\bm \vartheta_{j2}$ are equation and kernel-specific hyperparameters and the matrices $\bm D_{\bm X_j}, \bm D_{\bm Z_j}$ are diagonal matrices with typical $i^{th}$ element $\hat{\sigma}^2_{\bm X_j i}, \hat{\sigma}^2_{\bm Z_j i}$. These are set equal to the empirical variances of the $i^{th}$ column of $\bm X_j$ and $\bm Z_j$, respectively. Inclusion of the diagonal scaling matrices $\bm D_{\bm X_j}$ and $\bm D_{\bm Z_j}$ serves to control for differences in the scaling of the explanatory variables. Notice that since the hyperparameters are allowed to differ, we essentially treat own and other lags asymmetrically through different functional approximations $f_j$ and $g_j$.
The kernel is scaled with the error variances in $\bm \Omega_j$. A typical diagonal element of the corresponding re-scaled kernel is given by $\omega_{jt} \times k_{\bm \vartheta_{j1}}(\bm x_{jt}, \bm x_{jt}) = \omega_{jt} \xi_{j1}$ and $\omega_{jt} \times k_{\bm \vartheta_{j2}}(\bm z_{jt}, \bm z_{jt}) = \omega_{jt} \xi_{j2}$. Typical off-diagonal elements are given by $\sqrt{\omega_{jt}} \sqrt{\omega_{j \tau}} \times k_{\bm \vartheta_{j1}}(\bm x_{jt}, \bm x_{j\tau})$ and $\sqrt{\omega_{jt}} \sqrt{\omega_{j \tau}} \times k_{\bm \vartheta_{j2}}(\bm z_{jt}, \bm z_{j\tau})$. The interaction between the kernel and the error variances gives rise to convenient statistical and computational properties.
First, note that if $\omega_{jt}$ is large, the corresponding prior on the unknown functions is more spread out. In macroeconomic data, $\omega_{jt}$ is typically large in crisis periods when the $\bm x_{jt}$ and $\bm z_{jt}$ are far away from their previous values. Since the diagonal elements of the kernels are effectively determined by $\bm \xi_j = (\xi_{j1}, \xi_{j2})'$ the presence of $\omega_{jt}$ allows for larger values in the marginal prior variance and thus makes large shifts in the unknown functions more likely. Second, the interaction between $\omega_{jt}$ and $\omega_{jt-1}$ implies that the covariances are scaled down if $\omega_{jt} \gg \omega_{jt-1}$, suggesting that the informational content decreases if increases in uncertainty are substantial (i.e., $\Delta \omega_{jt}$ is large). If $\omega_{jt} \approx \omega_{j\tau}$ and both are large, the corresponding covariance will be scaled upwards. This implies that our model learns from previous crisis episodes as well. Third, as we will show in Sub-section (ref), interacting the kernel with the error variances leads to a conjugate Gaussian process structure which implies that we can factor out the error volatilities and do not need to update several quantities during MCMC sampling. This speeds up computation enormously and allows for estimating large models.
Before discussing how we select the hyperparameters, it is worth highlighting a possible identification problem of our model. In our baseline specification we center $\bm f_j$ and $\bm g_j$ around zero a priori. If we introduce an additional intercept term (or a simpler mean function) no identification issues arise. However, if we believe that $\bm f_j$ and $\bm g_j$ are centered on non-zero values, we can not separately identify them. In our empirical work, we normalize the grand mean of $\bm g_j$ to be equal to zero.\footnote{Notice that this only concerns the posterior distribution since, under the prior, this condition is automatically fulfilled.}
It is worth stressing, however, that if interest is on predictions or impulse responses, this does not cause any additional issues since the conditional mean function (which is the sum over $\bm f_j$ and $\bm g_j$) is identified. Exploiting basic properties of the Gaussian distribution one can easily show that the sum of $\bm f_j$ and $\bm g_j$ in (ref) gives rise to a new latent process $\bm m_j$ which is, again, Gaussian:
Hence, our model can be also viewed as a standard Gaussian process that combines information in $\bm X_j$ and $\bm Z_j$ by summing over two different kernels. This immediately implies that if we are interested in sampling from the posterior predictive distribution of $\bm y_{t+h}$ (and related functions such as impulse responses) it is sufficient to estimate $\bm m_j$.
So far, we always conditioned on the hyperparameters that determine the shape of the Gaussian kernel. A simple way of specifying $\bm \vartheta_{j1}$ and $\bm \vartheta_{j2}$ is the median heuristic approach stipulated in chaudhuri2017mean. This choice works well in a wide range of applications featuring many covariates crawford2019variable. The median heuristic fixes $\xi_{j1} = \xi_{j2} = 1$ and defines the inverse of the bandwidth parameter as:
for $j = 1, \dots, M$. This simple approach has the convenient property that it automatically selects a bandwidth which is consistent with the time series behavior of the elements in $\bm y_t$. To illustrate this, suppose that $y_{jt}$ is a highly persistent process (e.g., inflation or short-term interest rates). In this case, for $\tau=t-1$, the Euclidean distance $\lVert \bm x_{jt} - \bm x_{j\tau} \rVert$ will be quite small and, hence, the mean function $\bm f_j$ smoothly adjusts. If $y_{jt}$ is less persistent and displays large fluctuations (e.g., stock market or exchange rate returns), the Euclidean distance $\lVert \bm x_{jt} - \bm x_{j\tau} \rVert$ will be large and, thus, $\bm f_j$ allows for capturing this behavior. The dispersion in $\bm z_{jt}$ might have important implications for $y_{jt}$ if the aim is to model a trend in $y_{jt}$ that depends on other covariates. This could arise in a situation where the prior on $\bm f_j$ is set very tight (i.e., the posterior of $\bm f_j$ will be centered on zero) and information not coming from $\bm x_{jt}$ would then determine the behavior of $y_{jt}$. This discussion highlights how the median heuristic allows for flexibly discriminating between signal and noise and thus acts as a non-linear filter which purges the time series from high frequency variation, if necessary.
Given that we work with potentially large panels of time series, it is questionable that the median heuristic works equally well for all elements in $\bm y_t$. As a solution, we propose to use the median heuristic to set up a discrete grid for both $\xi_{j1}$ ($\xi_{j2}$) and $\kappa_{j1}$ ($\kappa_{j2}$). For each element in this grid we specify a hyperprior. We use Gamma priors on all elements. For $j = 1, \dots, M$, that is
for the linear shrinkage hyperparameters and
for the bandwidth parameters. Here, $c_{\xi1}, c_{\xi2}, c_{\kappa 1}$ and $c_{\kappa 2}$ are scalars that define the tightness of the hyperprior. In the empirical application, we set $c_{\xi1} = c_{\xi2} = c_\xi$ and $c_{\kappa 1} = c_{\kappa 2} = c_\kappa$. These parameters strongly influence the shape of the conditional mean and are crucial modeling choices and we set them through cross-validation. In our empirical application, we find that small values of $c_\kappa$ work reasonably well, yielding an informative prior that forces $\kappa_{j1}$ and $\kappa_{j2}$ towards zero.
Based on this set of priors we can derive the conditional posterior distribution. Since $\xi_{j1}$ ($\xi_{j2}$) and $\kappa_{j1}$ ($\kappa_{j2}$) are placed on a grid, we can pre-compute several quantities related to the kernel (such as inverses and Cholesky factors) while at the same time infer them from the data with sufficient accuracy, which is crucial for precise inference. In what follows, we center these grids around the median heuristic and additionally take into account the considerations of the informative Gamma priors:
Here, the intervals indicate the minimum (maximum) value supported for each hyperparameter. Within this two dimensional range, we define a discrete grid of around $1000$ combinations with equally sized increments along each dimension.\footnote{Implicitly, this two dimensional grid results in a prior view in which any hyperparameter combination not included in the grid has zero support.} The corresponding posterior is discrete and we can use inverse transform sampling to carry out posterior inference. Further details are provided in Sub-section (ref).
Posterior inference for the GP-VAR is carried out using a novel yet conceptually simple MCMC algorithm which cycles between several steps. In this section we will focus on how to sample from the posterior of $\bm f_j$, $p(\bm f_j | \bullet)$, with $\bullet$ denoting conditioning on everything else, and $\bm \vartheta_{j1}$. Sampling from $p(\bm g_j|\bullet)$ and $p(\bm \vartheta_{j2}|\bullet)$ works analogously with some adjustments. These relate to the fact that we introduce a linear restriction that $(\bm \iota' \bm \iota)^{-1} \bm \iota' \bm g_j = 0$, with $\bm \iota$ denoting a $T \times 1$ vector of ones. The corresponding conditional posterior distribution is a hyperplane truncated Gaussian where efficient sampling algorithms are available cong2017fast. Further details can be found in Sub-section (ref) of the Online Appendix. It is worth stressing that we sample $\bm f_j$ and $\bm g_j$ separately. This increases the computational burden slightly but allows us to consider both latent processes separately from each other. In case our focus is purely on prediction or impulse response analysis, one can also simulate the process $\bm m_j$ defined in Sub-section (ref) without any additional restriction. Both procedures yield exactly the same results.
{Generalizing the results in Section (ref), it can be shown that} the posterior of $\bm f_j$ is Gaussian for all $j$:
with posterior moments given by:
where for $j=1$, the term $\sum_{k=1}^{j-1} q_{jk} \bm Y_k$ is excluded.
In principle, computing the inverse and the Cholesky factor of $\overline{\bm V}_{\bm f_j}^{-1}$ constitutes the main bottleneck when it comes to sampling from $p(\bm f_j|\bullet)$. This is especially so if $T$ is large. But, in common macroeconomic applications which use quarterly US data, $T$ is moderate and thus computation is feasible. In our case, even if interest centers on using monthly data or even higher frequencies, we can exploit the convenient fact that, conditional on the hyperparameters $\xi_{j1}$ and $\kappa_{j1}$,
as well as its Cholesky factor $\bm B_{\bm f_j}$ can be pre-computed. In addition, notice that
These practical properties (due to the conjugate structure) substantially speed up computation in terms of sampling from $p(\bm f_j|\bullet)$.
These results are conditional on the hyperparameters. As outlined in the previous sub-section, we will estimate them by defining a discrete two dimensional grid of $1000$ combinations. For each hyperparameter combination on this grid, we compute the corresponding kernel $K_{\bm \vartheta_{j1}}(\bm X_j, \bm X_j)$ as well as all relevant quantities (i.e., $\bm C_{\bm f_j}$). Based on these values we jointly evaluate the conditional posterior ordinate by applying Bayes theorem. The exact form of the conditional likelihood is given by:
Note that the shape of $K_{\bm \vartheta_{j1}}(\bm X_j, \bm X_j)$ depends on the hyperparameters $\bm \vartheta_{j1} = (\xi_{j1}, \kappa_{j1})'$, which we want to update. For each pair of values $\bm \vartheta_{j1}^{(s)}$ on our two dimensional grid (with $s$ denoting a specific combination), we compute the corresponding kernel $K_{\bm \vartheta_{j1} = \bm \vartheta_{j1}^{(s)}}(\bm X_j, \bm X_j)$ as well as $\text{det}\left(K_{\bm \vartheta_{j1} = \bm \vartheta_{j1}^{(s)}}(\bm X_j, \bm X_j)\right)$ and $\left(K_{\bm \vartheta_{j1} = \bm \vartheta_{j1}^{(s)}}(\bm X_j, \bm X_j)\right)^{-1}$ prior to MCMC sampling. Hence, within our sampler evaluating the likelihood is straightforward and computationally efficient. All that remains is to multiply the likelihood with the prior. The corresponding posterior ordinates for each $\bm \vartheta_{j1}^{(s)}$ are used to compute probabilities to perform inverse transform sampling to sample from $p(\bm \vartheta_{j1}|\bullet)$.
Conditional on $\bm f_j$ and $\bm g_j$, the remaining parameters (i.e., the free elements in $\bm Q$, the log-volatilities and the associated coefficients in the corresponding state equations) can be sampled through (mostly) standard steps. One modification relates to how we sample the volatilities in $\bm \Omega_j$. The main difference stems from the fact that the volatilities in $\bm \Omega_j$ also show up in the prior on $\bm f_j$ and $\bm g_j$. To circumvent this issue we integrate out the latent processes $\bm f_j$ and $\bm g_j$. This calls for a minor adjustment of the original sampler by integrating out the latent processes $\bm f_j$ and $\bm g_j$ first and then sampling the log-volatilities using an independent Metropolis Hastings update similar to the one proposed in chan2017stochastic. We provide additional details and the full posterior simulator in Section (ref) of the Online Appendix.
In non-parametric models such as the GP regression described in Section (ref), the effect of the covariates on $\bm y_t$ are typically analyzed through so-called partial dependence plots friedman2001greedy. Our large dimensional setting and the fact that we have a VAR-type structure in the conditional mean, imply that partial dependence plots are difficult to compute and visualize since they would require integration over a large number of covariates. Moreover, VARs are dynamic models and partial dependence plots are difficult to employ in dynamic settings. To capture non-linear model dynamics and possible relations across variables, impulse responses are used to investigate the effects of structural shocks on $\bm y_t$. In this paper, we will follow this route as well. Because the model is highly non-linear, we need to resort to generalized impulse responses (GIRFs) originally proposed in koop1996impulse. {As GIRFs are based on (multi-step-ahead) forecasts, and since forecasting with the GP-VAR is of interest by itself, we start with a discussion on forecast computation.}
We begin by computing the predictive distribution of the one-step-ahead forecasts $p({\bm y}_{t+1}|\mathcal{I}_t)$, with $\mathcal{I}_t$ denoting all available information up to time $t$. The one-step-ahead predictive distribution is obtained by simulating from the predictive distribution of $\bm m_{t+1}$, $p(\bm m_{t+1}|\mathcal{I}_t)$, and sampling from the marginal distribution of the shocks $\bm \varepsilon_{t+1} \sim \mathcal{N}(\bm 0, \bm H_{t+1})$. The draw from the one-step-ahead predictive density is used to set up $\bm x_{t+2}$ and $\bm z_{t+2}$. Based on these, we can compute the corresponding kernels and obtain a draw $\bm m_{t+2}$ from the density $p(\bm m_{t+2}|\mathcal{I}_t)$. Again, a draw from $\bm y_{t+2} \sim p(\bm y_{t+2}|\mathcal{I}_t)$ is obtained by sampling from the marginal shock distribution $\bm \varepsilon_{t+2} \sim \mathcal{N}(\bm 0, \bm H_{t+2})$ and adding this draw to $\bm m_{t+2}$. Higher order forecasts are obtained analogously. The resulting predictive distribution of $\bm y_{t+h}$ will be highly non-Gaussian and might feature heavy tails and/or asymmetries. This forms the baseline.
The GIRFs are computed as follows. To analyze the effects of a structural disturbance (such as an uncertainty shock) we assume that the uncertainty indicator is (without loss of generality) in the $j^{th}$ position in $\bm y_t$. A corresponding shock of size $\varsigma$ in time $t$ to the uncertainty indicator shifts all elements in $\bm y_t$ by $\varsigma ~ \bm q_j$, i.e., the $j^{th}$ column of $(\bm I - \bm Q)^{-1}$ scaled by a scalar that reflects the shock size $\varsigma$. The other shocks are sampled, again, from their marginal distributions, i.e., for all $i \neq j$ we have that $\varepsilon_{it} \sim \mathcal{N}(0, \omega_{it})$. Based on this we draw from the predictive distribution conditional on the uncertainty shock $\hat{\bm y}_{t+1} \sim p(\bm y_{t+1}|\mathcal{I}_t, \varepsilon_{jt}=\varsigma)$ and use this draw to compute $\hat{\bm x}_{t+2}$ and $\hat{\bm z}_{t+2}$. For higher order conditional forecasts we proceed as in the case of the unconditional forecast distribution by simulating from the marginal distribution of the structural shocks which are added to the conditional mean forecasts $\hat{\bm m}_{t+h}$.
The corresponding dynamic responses are then obtained by subtracting the mean of the unconditional predictive distribution from the conditional (on the uncertainty shock) predictive density. This yields:
Notice that $\bm \delta_{ht}$ is state-dependent and, due to the non-linear nature of the conditional mean function, allows for asymmetries in how $\bm y_t$ reacts to shocks.\footnote{The full details on computing generalized impulse response functions can be found in Section (ref) of the Online Appendix.} This gives rise to two inferential opportunities. First, one can assess how a given shock has impacted the economy in a given point in time. This allows us to investigate whether transmission mechanisms depend on the underlying state of the economy. Second, the non-linear mean function directly implies that shock transmission can be asymmetric, so that positive shocks might feed through the economy differently than negative shocks, and non-proportional, so that larger shocks can have proportionally different effects than smaller shocks. In all our empirical work we will exploit both dimensions and focus on asymmetries in the sign, and non-proportionality in the size, of the shock as well as explicitly consider state dependencies by computing $\bm \delta_{ht}$ over time. Finally, we also integrate out uncertainty with respect to the state of the economy by averaging over all values of $t$.
In this section, we illustrate the computational merits of our approach and evaluate whether it successfully recovers different features of a highly non-linear DGP.
To illustrate our methods, we simulate $T=200$ observations from a highly non-linear small-scale VAR with $M=3$ equations. The three equations differ in terms of whether they are linear or non-linear in the parameters but also with respect to the distribution of the shocks. Non-linearities are captured in two ways. First, we assume a break point and second we assume non-linear relations between the response variable and the lags of the other variables. In all these equations, we assume that the functions $f_j$ and $g_j$ differ to assess whether our approach is capable of discriminating between the two. The precise form of our DGP is given by:
with
where $\bm H_t = \text{diag}(\omega_{1t},\omega_{2t},\omega_{3t})$, $\phi_{11,1} = 0.8$, $\phi_{22,1} = 0.65$, $\phi_{ij,k} \sim \mathcal{N}\left(0, \left(\frac{0.3}{k}\right)^2\right)$ for $i \neq j$ and $k = 1, \dots p$, and $\mathcal{I}(\bullet)$ denotes the indicator function that equals one if its argument is true and zero otherwise. The free elements in $\bm Q$, are also simulated from a Gaussian distribution with $q_{jk} \sim \mathcal{N}(0,0.1^2)$. Moreover, we introduce an SV specification with heavy tails for the structural error variances $\omega_{jt} = \lambda_{jt} \tilde{\omega}_{jt}$ with each $\tilde{h}_{jt} = \log \tilde{\omega}_{jt}$ following an independent random walk law of motion: $\tilde{h}_{jt} = \tilde{h}_{jt-1} + \sigma_{\tilde{h} j} u_{\tilde{h} t}$, with $u_{\tilde{h} t} \sim \mathcal{N}(0, 1)$. For each equation, we set the initial state $\tilde{\omega}_{j0} = \exp \tilde{h}_{j0} = 0.01$ and the state innovation variance $\sigma_{\tilde{h} j} = 0.01$. We consider (conditionally) $t_3$-distributed errors with three degrees of freedom in the first equation by simulating $\lambda_{1t} \sim \mathcal{G}^{-1}(3/2, 3/2)$ and (conditionally) Gaussian-distributed errors for the second and third equation by setting $\lambda_{2t} = \lambda_{3t} = 1$ for all $t$.
Figure (ref) shows results obtained from simulating a single realization from the DGP. Horizontal panels refer to each equation of the DGP while vertical panels show the different components. The red shaded areas represent the $90\%$ posterior credible set of {$\bm f_j, \bm g_j$, $\bm m_j (= \bm f_j + \bm g_j$}) and $\bm y_j$, while the solid black line denotes the actual outcome of these quantities. The figure suggests that our approach is capable of detecting different functional relations between $\bm y_t$, $\bm x_t$ and $\bm z_t$. Across all three equations, we find that the estimated conditional mean function $\bm m_j$ tracks the actual value rather well. It is also worth stressing that if the DGP features non-Gaussian shocks, the model recovers the true mean function particularly well (see the first row of Figure (ref)). Zooming into the estimates for the different latent components reveals that most of this strong fit is driven by accurate estimates of $\bm f_j$. Considering the estimates for $\bm g_j$ suggests that, if the actual relationship between $\bm z_t$ and $\bm y_t$ is non-existent, our approach accurately detects this behavior. However, there are some cases where $\bm f_j$ soaks up variation in $\bm g_j$. This, however, does not impact the mean estimate $\bm m_j$.
Since this discussion has been based on a single draw from the DGP, one might ask whether the strong performance of the GP-VAR is due to a particularly favorable realization from the DGP. To briefly investigate whether this is the case we show, in the first row of Table (ref), the average correlations between the true and posterior median of the functions for 100 draws from the DGP. Numerical standard errors (across these $100$ repetitions) are shown in the second row. The first row indicates that correlations are high, reaching $0.95$ for the first, $0.93$ for the second and $0.75$ for the third equation. The fact that mean correlations slightly decrease are mainly driven by the fact that equations two and three feature substantial high frequency movements which, in our framework, is mostly picked up by the stochastic volatility component.
One key advantage of the GP-VAR is that we can remain agnostic on the form of non-linearities the conditional mean function might take. This is confirmed for synthetic data where, irrespective of the non-linear nature of the model, the estimated sum of $\bm f_j$ and $\bm g_j$ (i.e., $\bm m_j$) closely tracks the dynamics of the actual outcome. This finding holds both for the linear case (i.e., the first equation){, even with fat-tailed errors,} and for highly non-linear situations (i.e., the second and third equations).
We have stressed that our approach is computationally efficient and scalable to large datasets. To investigate this claim more carefully, Figure (ref) shows the time required to generate 1000 draws from the joint posterior for a given equation across different values of $K$ and for $T=200$. We show the computation times for our GP-VAR and a VAR with SV. Since our approach is embarrassingly parallel the actual times for generating a draw from the joint posterior of the full system are approximately $M$ times the runtimes reported in the figure.\footnote{Since we need to augment the $j^{th}$ equation with the contemporaneous values of the preceding $j-1$ equations, this statement is only approximately valid.}
The most striking take away from the figure is that the computation time of the GP-VAR does not depend on $K$. This implies that increasing the number of lags and/or endogenous variables does not impact estimation times considerably. By contrast, the time necessary to generate a draw from the joint posterior of the VAR rises rapidly in $K$. This shows that our approach scales well in high dimensions and, in fact, is much faster than competing approaches to non-linear VAR models such as TVP-VARs or regime-switching VARs.
To provide a rough gauge on actual estimation times for practitioners, MCMC estimation of the GP-VAR with eight endogenous variables and five lags takes around 30.5 minutes on a standard desktop computer (for 10,000 MCMC draws). These are gross estimation times and thus include pre-computation of matrices used during MCMC simulation and are not based on parallel computation of the individual equations (which is possible due to the structural form). Estimating larger models (such as the 64 variable GP-VAR) takes around four hours.
In this section, we apply our GP-VAR to US macroeconomic data. Sub-section (ref) provides details on the data and Sub-section (ref) discusses whether the GP-VAR fits the data well, shows some of its key in-sample features, and investigates its forecast performance.
We use the quarterly version of the dataset proposed in mccracken2016fred and consider time series that range from $1960$Q$1$ to $2019$Q$4$. We exclude the years of the Covid-19 pandemic to make the comparison with a linear VAR similar to that used in jurado2015measuring (JLN) sensible. In the following, we consider four different specifications that differ in the number of endogenous variables used. These specifications are:
All variables are transformed to be approximately stationary and we include five lags of the endogenous variables. The precise variables included (and transformations applied to each variable) are shown in Table (ref) in the Online Appendix. We consider these different model sizes for several reasons. First, we would like to assess how adding additional information impacts the forecasting performance and the responses of key variables to an uncertainty shock. Second, we are interested in the relationship between non-linearities and the size of the model.
In this section, we start by providing some predictive evidence of our GP-VAR and investigate how the role of the parameters associated with the kernel impact predictive accuracy. To this end, we employ a recursive forecasting design. Our initial training period goes from $1960$Q$1$ to $1999$Q$4$. After computing the one- and four-step-ahead predictive distributions, we add an additional observation, recompute all models and simulate from the corresponding predictive densities. This procedure is repeated until the end of the sample ($2019$Q$4$) is reached.
We consider three ways of specifying the equation-specific ($j = 1, \dots, M)$ and kernel-specific $(k = \{1, 2\})$ hyperparameters. The first one is a semi-automatic approach that is based on putting the (inverse) length scale and the linear scaling parameter on a two dimensional grid $\kappa_{jk} \in [0.1 \bar{\kappa}_{jk}, 2 \bar{\kappa}_{jk}]$ and $\xi_{jk} \in [0.04, 4]$. Notice that the median heuristic is used to determine the lower and upper bound of the grid. The second approach (labeled “semi-automatic w/o linear scaling" in Table (ref)) puts $\kappa_{jk} \in [0.1\bar{\kappa}_{jk}, 2 \bar{\kappa}_{jk}]$ on a grid and sets $\xi_{jk} = 1$. Finally, we also consider a specification that does not rely on the median heuristic (labeled “naive" in the table). The grid is $\kappa_{jk} \in [0.1,2]$ and $\xi_{jk}$ is again set equal to $1$.
As competing models, we include the standard Minnesota BVAR with SV, a TVP-VAR-SV similar to the one used in primiceri2005time and the BART-VAR with SV proposed in huber2022inference. All models are estimated for different model sizes and benchmarked to a small-scale Minnesota BVAR with SV {for the three variables we focus on: real GDP growth (RGDP), CPI inflation (CPI), and the Fed funds rate (FFR).} For the three focus variables, we compute log predictive Bayes factors (LPBFs) relative to the small-scale BVAR with SV, so that positive numbers indicate that a given model works better than the benchmark while negative values suggest a weaker forecasting performance. The LPBFs do not only take into account how well a given model predicts the realization of a given variable but also factor in forecasting performance for higher order features of the predictive distribution geweke2010comparing. To investigate whether controlling for heteroskedasticity pays off in the GP-VAR, we also consider homoskedastic variants of the GP-VAR in the upper part of the table.
(ref) shows the results of our forecasting exercise across different model sizes. The columns “Joint" show the joint LPBF for the three focus variables and thus provide a comparable (across model sizes) metric of overall predictive accuracy. Considering joint LBPFs reveals that GP-VAR SV with $M=16$ and a semi-automatic approach for hyperparameter elicitation improves upon all competing models for both forecast horizons. The smallest model ($M=8$) also yields competitive predictions. Once we further increase the size of the dataset, predictive accuracy slightly deteriorates. Interestingly, and consistent with findings in clark2021tail, we find that the gains in predictive accuracy increase when we focus on higher forecast horizons. Comparing the models with SV to their homoskedastic counterparts paints a very consistent picture. Models which do not control for time variation in the error variances perform consistently worse than the models that have SV in the error terms.
To drill deeper into which variables drive the overall forecasting performance, we now focus on the marginal LPBFs for the three focus variables. Starting with one-quarter-ahead predictions of GDP growth, we observe that the GP-VARs with SV are beaten by the BVAR-SV with $M=32$. When we turn off SV, predictive accuracy sometimes increases by small margins. This result, however, changes if we focus on higher order forecasts. For one-year-ahead predictions of GDP growth, the single best performing model is the largest ($M=64$) homoskedastic GP-VAR with the semi-automatic approach that fixes the linear scaling parameter to one. Strikingly, at that horizon and for this specific variable, using a homoskedastic specification is almost uniformly better than the corresponding SV setup. For inflation and the Fed funds rate, the smaller-sized GP-VARs again yield the best density forecasting performance, outperforming all competing models. Finally, focusing on interest rate forecasts shows that GP-VARs with SV do well (for all model sizes) but once we turn off SV forecasts become highly imprecise. This is driven by the fact that during the zero lower bound, the conditional variance of the interest rate equation approaches zero and a model which assumes homoskedasticity fails to take that into account.
After having established that the different GP-VARs with SV do well when used to forecast US macroeconomic quantities, we focus on some in-sample features for the GP-VAR with eight endogenous variables. To get an impression on how the linear shrinkage parameter evolves over time, Figure (ref) plots the product of the error variances times the linear shrinkage parameter that determines the kernel of $\bm f_j$ and $\bm g_j$. Two interesting features emerge. First, less shrinkage (larger parameter values) is applied to the own lags than to the lags of the other variables. This holds for most variables under scrutiny (except for CPI inflation and S&P 500 returns). Second, less shrinkage is typically applied during recessionary times. In particular, we introduce little shrinkage during the recessions in the early 1980s and the financial crisis. Notice, however, that there are also some exceptions from this pattern (such as stock market returns, hours worked or hours employed). This finding is, again, in line with results from the forecasting literature showing that more information is particularly useful during problematic times koop2013forecasting.
Finally, we investigate differences in $\kappa_{j1}$ and $\kappa_{j2}$ across equations ($j = 1, \dots M$) and variable types. Boxplots that show the posterior distribution of the inverse of the length scale parameters for own and other lags are in (ref). Recall that large values of $\kappa_{jk}~(k=1,2)$ imply more variation in the latent processes whereas values of $\kappa_{jk}$ close to zero imply less variation in $\bm f_j$ and $\bm g_j$. A general pattern is that for all variables the hyperparameters associated with the kernel on own lags are considerably larger than the ones for the kernel related to the other lags. This indicates that the own lags of a given endogenous variable require more flexibility (i.e., functions that allow for much more variation) whereas the effect of other lags appears to be more linear. For the majority of variables (except for hours worked, CPI inflation and the Fed funds rate), the posterior distribution of the hyperparameter looks similar. For the three exceptions, $\kappa_{j1}$ is much smaller and more precisely estimated.
We now analyze the effects of uncertainty shocks using our GP-VAR and focus on assessing how macroeconomic uncertainty feeds through the economy. In Sub-section (ref), as a benchmark exercise, we compare our impulse responses to those obtained from a model similar to that used by JLN. In Sub-section (ref) we leverage the non-linear nature of the GP-VAR and analyze how the effects of uncertainty shocks change according to the sign or size of the shocks, and over time.
We benchmark the IRFs of our GP-VAR-8 to the ones of a BVAR with SV that is closely related to the original JLN specification.\footnote{While they use a classical homoskedastic VAR estimated on monthly data, we work with a quarterly BVAR with SV.} In what follows, our focus will be on the variables discussed in JLN: year-on-year growth rates of output (measured through real GDP) and employment. We also show the responses of the uncertainty indicator and the quarter-on-quarter returns of the S&P 500 and include the responses of the other variables in Section (ref) of the Online Appendix.
In (ref) we report the (average over time in the case of the GP-VAR) posterior quantiles ($16^{th}$, $50^{th}$ and $84^{th}$) of the responses to a macroeconomic uncertainty shock in the JLN model (in gray) and in the corresponding GP-VAR (in orange) with eight endogenous variables and uncertainty ordered second after stock market returns. Uncertainty responses to its own shock differ slightly between the GP-VAR and the linear model. These differences relate to responses within the first five quarters after the shock hit the system. The BVAR yields uncertainty reactions that peak after one quarter, declining steadily afterwards. As opposed to this swift reaction in uncertainty, the GP-VAR generates endogenous uncertainty reactions which peak after five quarters, declining sharply afterwards. After around eight quarters, both IRFs (almost) coincide.
This uncertainty reaction has direct implications on how the other variables in the model react. Real GDP growth reacts in an hump-shaped manner under both models. However, driven by the somewhat later peak in uncertainty, the GP-VAR produces much stronger output growth reactions that peak slightly later (after around five quarters). When we focus on employment growth the IRFs differ less. In principle, both models suggest a peak decline of around 0.8 percentage points, with the GP-VAR generating a somewhat slower response, reaching its trough after about six to seven quarters. But in principle, responses between the linear and non-parametric model tell a similar story. Finally, financial market reactions measured through the S&P500 suggest a much stronger decline in stock prices under the GP-VAR. Interestingly, the shape of the IRFs suggests that the linear model generates the strongest reaction after around two quarters. In the GP-VAR, we find that stock markets react faster and stronger to uncertainty shocks, with substantial reactions within the first year after the shock hit the system.
To conclude, in (ref) we report the GIRFs to the uncertainty shock for the same four variables displayed in (ref) but obtained from GP-VARs of different dimensions (with 8, 16, 32, and 64 variables). Differences across model sizes are small (or non-existent) for most variables. Small differences arise for employment growth, with the magnitude of the responses increasing with the model size. Stock market reactions also differ slightly across datasets, with no clear-cut pattern. Since the GIRFs are very similar across model size and given its excellent forecasting properties, we will focus on asymmetries generated by the GP-VAR-8 model in the following sections. Results for the larger models are provided in the Online Appendix.
The non-linear and non-parametric nature of our models allows for asymmetries in the impulse response functions. This implies that shocks propagate non-linearily through the model, giving rise to differences in the GIRFs both over time but also for different shock magnitudes or signs.
The first aspect we consider relates to whether positive and negative uncertainty shocks trigger different responses of the economy. In (ref) we report the responses to negative and positive uncertainty shocks from the GP-VAR-8, averaged over time. The figure thus shows GIRFs to a positive (in orange), negative (in blue) and a negative shock multiplied by -1 (in gray, to ease comparison).
From the figure we observe some differences. These differences mostly relate to peak reactions as well as short-run (i.e. within two years) responses. In general, we find that positive shocks (higher uncertainty) trigger a stronger reaction of uncertainty, which in turn translates into more pronounced reactions of real activity and stock market quantities.
More specifically, considering the endogenous reaction of the uncertainty indicator shows that responses to a positive uncertainty shock peak after around a year and quickly die out afterwards. However, if uncertainty unexpectedly declines, the peak happens on impact and is much smaller as opposed to an adverse uncertainty shock.
Turning to real GDP and employment growth, we find that positive shocks trigger stronger reactions for both variables. Interestingly, the timing of the peak responses is similar for negative and positive shocks but reactions appear much more pronounced for the latter. Stock market reactions also differ markedly across positive and negative shocks. For positive shocks we, again, find that the peak effect happens after one year and that it is more pronounced as compared to the negative shock. Overall, the picture that emerges from the GP-VAR is that higher unexpected uncertainty has stronger effects on the economy than lower uncertainty, a feature that is a priori ruled out in linear VARs.
Our GP-VAR also allows for analyzing how shocks of different sizes impact the economy. As opposed to a standard VAR which assumes that shocks enter linearly (and thus responses to shocks of different sizes are exactly proportional to each other) our GP-VAR is more flexible and allows for investigating whether shocks of different magnitudes trigger different dynamics in the GIRFs.
In (ref) we consider two shock sizes: a one standard deviation and a two standard deviation shock. To permit straightforward comparison of the shapes of the responses to differently sized shocks, we also add impulses to a two standard deviation shock which are then re-scaled to match the impact of the one standard deviation shock (the gray shaded area in the figures).
This figure gives rise to at least two observations. First, when we compare the shape of the responses to a one standard deviation to the ones of a two standard deviation shock we find differences in the timing (and more generally in the shape) of the IRFs. A stronger shock triggers a faster peak reaction of GDP and employment growth. Stock market reactions display a somewhat different shape. After a sharp immediate reaction (for both shock sizes) the peak effect happens to be on impact if the size of the shock is large whereas it turns out to materialize after one year if the shock size is smaller.
Second, in terms of the magnitudes we find that a two standard deviation shock triggers peak responses with magnitudes that are less than twice the magnitudes to a one standard deviation shock. This is particularly visible for employment and output reactions. For stock market responses, the impact reactions are (almost) proportional to each other.
After showing that the economy reacts asymmetrically with respect to the sign and size of the uncertainty shock, this section asks whether the effect of uncertainty shocks changes over time (see Section 2.5 of castelnuovo2022jes, or mumtaz2018changing for some evidence using TVP-VARs).
We start by considering impulse responses averaged over certain sub-periods in (ref). The classification into sub-periods is mostly taken from d2012century, and it is such that the main events in each sub-period include, respectively, the great inflation ($1970$Q$1$ to $1984$Q$4$), the great moderation ($1985$Q$1$ to $2006$Q$4$), and the post great moderation period ($2007$Q$1$ to $2019$Q$4$). To also get a rough feeling about whether asymmetries between positive and negative shocks have changed over time, all figures include the IRFs to positive (in orange) and negative (in blue) shocks.
The main feature emerging from the figure is the different behavior of the response of uncertainty across sub-samples. In the final two sub-samples, uncertainty responses increase up to four quarters after the shock, with peak effects being strongest in the great moderation period and becoming slightly weaker in the final sub-sample. Moreover, sign effects of uncertainty responses increase appreciably in the last two sub-samples.
These differences in the responses of uncertainty trigger differences in the IRFs of the other quantities which relate not only to the magnitudes but also to the shapes of the responses. We find that GDP growth, employment and the S&P 500 display the strongest reactions in the great moderation regime, becoming slightly weaker in the post great moderation period. The weaker reaction of real activity over time corroborates findings in mumtaz2018changing who also report smaller responses of real activity to uncertainty shocks. As opposed to their findings, we observe that stock market reactions do not change much in magnitude but the shape differs (in accordance with the different shape in the uncertainty reaction described above). We, moreover, observe that asymmetries in terms of the sign of the shocks have decreased over time for GDP and employment growth. Only for stock market reactions these sign asymmetries have increased, with benign uncertainty shocks yielding a much weaker positive reaction of stock markets during the post great moderation regime.
Finally, to conclude this section we raise the issue that considering GIRFs averaged over sub-samples possibly still masks important differences over time within sub-samples. To shed light on whether IRFs change within regimes, (ref) displays yearly averages of posterior medians of the IRFs over time during each of the three periods. Yellow IRFs refer to the beginning of the respective sub-sample and red ones denote IRFs computed towards the end of the sub-sample. This figure suggests substantial heterogeneity in responses during the great inflation period. Especially towards the end of this sample, reactions of GDP growth and employment point towards a substantial real activity overshoot. During the great moderation, the intra-period variation of the IRFs becomes much smaller, yielding patterns more consistent with the common wisdom in the uncertainty literature: real activity and stock markets decline in response to increases in economic uncertainty. In the years from 2007 to 2019, we find that IRFs differ especially in the beginning of the sample (from 2007 to 2009). For the remaining years, there is much less variation in responses and these appear to be similar to those observed in the 1985 to 2006 period.
Overall, we can conclude that the effects of uncertainty change both during sub-samples defined by economic considerations and sometimes also within each sub-sample. This kind of time variation is a priori ruled out in linear VAR models, which can therefore lead to biased estimates of the effects of uncertainty.
In this paper, we have developed a flexible multivariate model that uses Gaussian processes to model the unknown relationship between a panel of macroeconomic time series and their lagged values. Our GP-VAR is a very flexible model which remains agnostic on the precise relations between the endogenous variables and the predictors. This model can be viewed as a very flexible and general extension of the linear VAR commonly used in empirical macroeconomics. We also control for changes in the error variances by introducing a stochastic volatility specification. While a more flexible conditional mean can reduce the need of a time-varying conditional variance, empirically we find heteroskedasticity to be relevant also for GP-VARs.
We develop efficient MCMC estimation algorithms for the GP-VAR, which are scalable to high dimensions, so much so that for large models estimation is even faster than for the corresponding BVAR-SV. Scaling the covariance of the Gaussian process by the latent volatility factors is particularly helpful to achieve computational gains, as it permits to pre-compute several quantities before MCMC sampling. This speeds up computation enormously.
To illustrate the practical working of the GP-VAR, we first test it on simulated data from different linear and non-linear models, finding that it is capable of reproducing a variety of non-linear patterns (but also a linear behavior). Then, we show in a forecasting exercise that our model yields favorable density forecasts of US output, inflation and short-term interest rates with respect to both linear and other non-parametric and time-varying specifications.
In the main part of our empirical work we re-assess the effects of uncertainty shocks by replicating and extending the analysis carried out by jurado2015measuring based on linear VARs with the GP-VAR. Overall, our empirical results suggest that the measurement of uncertainty and its effects with a simple linear VAR can lead to several incorrect conclusions. Not only the effects of uncertainty can be over-stated, but they can also be treated as stable over time, symmetric for positive and negative shocks, and proportional to the shock size. Instead the GP-VAR model, which is preferred to the linear VAR in terms of fit and forecasting performance, returns time variation in the responses, asymmetry and non-proportionality. Hence, the empirical features we uncover should be also replicated by theoretical models about uncertainty and its effects, which instead at the moment typically assume stability and symmetry bloom2014fluctuations.
{\setstretch{0.85} \addcontentsline{toc}{section}{References} }
\setcounter{page}{0} \thispagestyle{empty} \setcounter{footnote}{0}