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.
87,017 characters · 15 sections · 103 citation commands
Large Hybrid Time-Varying Parameter VARs
\onehalfspacing
\thispagestyle{empty}
Time-varying parameter vector autoregressions (TVP-VARs) developed by CS01, CS05 and Primiceri05 have become the workhorse models in empirical macroeconomics. These models are flexible and can capture many different forms of structural instabilities and the evolving nonlinear relationships between the dependent variables. Moreover, they often forecast substantially better than their homoskedastic or constant-coefficient counterparts, as shown in papers such as clark11, \citet*{DGG13}, KK13, CR15 and CP16. In empirical work, however, their applications are mostly limited to modeling small systems involving only a few variables because of the computational burden and over-parameterization concerns.
On the other hand, large VARs that use richer information have become increasingly popular due to their better forecast performance and more sensible impulse-response analysis, as demonstrated in the influential paper by \citet*{BGR10}. There is now a rapidly expanding literature that uses large VARs for forecasting and structural analysis. Prominent examples include \citet*{CKM09}, koop13, BGMR13, \citet*{CCM15}, ER17 and MW19. Since there is a large body of empirical evidence that demonstrates the importance of accommodating time-varying structures in small systems, there has been much interest in recent years to build TVP-VARs for large datasets. While there are a few proposals to build large constant-coefficient VARs with stochastic volatility CCM16, CCM19, KH18, chan20, chan21, the literature on large VARs with time-varying coefficients remains relatively scarce.
We propose a class of models we call hybrid TVP-VARs---VARs in which some equations have time-varying coefficients, whereas the coefficients are constant in others. More precisely, we develop an efficient Bayesian shrinkage and sparsification method that automatically decides, for each equation, (i) whether the VAR coefficients are constant or time-varying, and (ii) whether the parameters of the contemporaneous relations among variables are constant or time-varying. Given the importance of time-varying volatility, all equations feature stochastic volatility. Our framework nests many popular VARs as special cases, ranging from a constant-coefficient VAR with stochastic volatility on one end of the spectrum to the flexible but highly parameterized TVP-VARs of CS05 and Primiceri05 on the other end. More importantly, our framework also includes many hybrid TVP-VARs in between the extremes, allowing for a more nuanced modeling approach of the time-varying structures.
To formulate these large hybrid TVP-VARs, we use a reparameterization of the standard TVP-VAR in Primiceri05. Specifically, we rewrite the TVP-VAR in the structural form in which the time-varying error covariance matrices are diagonal. Hence, we can treat the structural TVP-VAR as a system of $n$ unrelated TVP regressions and estimate them one by one. This reduces the dimension of the problem and can substantially speed up computations. This approach is similar to the equation-by-equation estimation approach in \citet*{CCM19} that is designed for the reduced-form parameterization. But since under our parameterization there is no need to obtain the `orthogonalized' shocks at each iteration as in \citet*{CCM19}, the proposed approach is substantially faster. Moreover, under our parameterization the estimation can be parallelized to further speed up computations. This structural-form parameterization, however, raises the issue of variable ordering, that is, the assumed order of the variables might affect the model estimates compared to a standard reduced-form TVP-VAR. We investigate this issue empirically and find that the variability of the estimates from this structural-form parameterization is comparable to that of the TVP-VAR of Primiceri05.
Next, we adapt the non-centered parameterization of the state space model in FSW10 to our structural TVP-VAR representation. Further, for each equation we introduce two indicator variables, one determines whether the VAR coefficients are time-varying or constant, while the other controls whether the elements of the impact matrix are time-varying or not. Hence, each vector $\boldsymbol \gamma\in\{0,1\}^{2n}$, where $n$ is the number of endogenous variables, characterizes a hybrid TVP-VAR with a particular form of time variation. By treating these indicators as parameters to be estimated, we allow the data to determine the appropriate time-varying structures, in contrast to typical setups where time variation is assumed. The proposed approach therefore is not only flexible---it includes many state-of-the-art models routinely used in applied work as special cases---it also induces parsimony to ameliorate over-parameterization concerns. This data-driven hybrid TVP-VAR can also be interpreted as a Bayesian model average of $2^{2n}$ hybrid TVP-VARs with different forms of time variation, where the weights are determined by the posterior model probabilities $p(\boldsymbol \gamma\,|\,\mathbf{y})$. It follows that forecasts from such a model can be viewed as a forecast combination of a wide variety of hybrid TVP-VARs.
The estimation is done using Markov chain Monte Carlo (MCMC) methods. Hence, in contrast to earlier attempts to build large TVP-VARs, our approach is fully Bayesian and is exact---it simulates from the exact posterior distribution. There are, however, a few challenges in the estimation. First, the dimension of the model is large and there are thousands of latent state processes---time-varying coefficients and stochastic volatilities---to simulate. To overcome this challenge, in addition to using the equation-by-equation estimation approach described earlier, we also adopt the precision sampler of CJ09 to draw both the time-invariant and time-varying VAR coefficients, as well as the stochastic volatilities. In our high-dimensional setting the precision sampler substantially reduces the computational cost compared to conventional Kalman filter based smoothers. A second challenge in the estimation is that the indicators and the latent states enter the likelihood multiplicatively. Consequently, it is vital to sample them jointly; otherwise the Markov chain is likely to get stuck. We therefore develop algorithms to sample the indicators and the latent states jointly.
Using US datasets of different dimensions, we find evidence that the VAR coefficients and elements of the impact matrix in some, but not all, equations are time varying. In particular, in a formal Bayesian model comparison exercise, we show that there is overwhelming support for the (data-driven) hybrid TVP-VAR relative to a few standard benchmarks, including a constant-coefficient VAR with stochastic volatility and a full-fledged TVP-VAR in which all the VAR coefficients and error covariances are time varying. We further illustrate the usefulness of the hybrid TVP-VAR with a forecasting exercise that involves 20 US quarterly macroeconomic and financial variables. We show that the proposed model forecasts better than many benchmarks. These results suggest that using a data-driven approach to discover the time-varying structures---rather than imposing either constant coefficients or time-varying parameters---is empirically beneficial.
This paper contributes to the budding literature on developing large TVP-VARs. Earlier papers include KK13, KK18, who propose fast methods to approximate the posterior distributions of large TVP-VARs. BV18 and GH18 consider large VARs with only time-varying intercepts. \citet*{CES20} model the time-varying coefficients using a factor-like reduced-rank structure, whereas HKO19 develop a method that first shrinks the time-varying coefficients, followed by setting the small values to zero. As mentioned above, our estimation approach is exact and fully Bayesian, and the modeling framework is more flexible than many of those in earlier papers. There is also a growing literature on alternative, non-likelihood based approaches. Examples include GKP13 and Petrova19 that allow for the estimation of large TVP-VARs without imposing the Cholesky-type stochastic volatility, and hence they avoid the ordering issue. Nevertheless, one main advantage of the likelihood-based approach taken in this paper is that it is flexible and modular. In particular, it is straightforward to incorporate additional useful features into the proposed hybrid model, such as more sophisticated static and dynamic shrinkage priors for VARs pruser21,chan21 or more flexible error distributions to deal with outliers CCMM21,BH21.
The rest of the paper is organized as follows. We first introduce the proposed modeling framework in Section (ref). In particular, we discuss how we combine a reparameterization of the reduced-form TVP-VAR and the non-centered parameterization of the state space model to develop the hybrid TVP-VARs. We then describe the shrinkage priors and the posterior sampler in Section (ref). It is followed by a Monte Carlo study in Section (ref) that demonstrates that the proposed methodology works well and can select the correct time-varying or time-invariant structure. The empirical application is discussed in detail in Section (ref). Lastly, Section (ref) concludes and briefly discusses some future research directions.
We first introduce a class of models we call hybrid time-varying parameter VARs: VARs in which some equations have time-varying coefficients, whereas coefficients in other equations remain constant. To that end, let $\mathbf{y}_t= (y_{1,t},\ldots,y_{n,t})'$ be an $n \times 1$ vector of endogenous variables at time $t$. The TVP-VAR of Primiceri05 can be reparameterized in the following structural form:
where $\mathbf{b}_t$ is an $n\times 1 $ vector of time-varying intercepts, $\mathbf{B}_{1,t}, \ldots, \mathbf{B}_{p,t}$ are $n \times n$ VAR coefficient matrices, $\mathbf{A}_{t}$ is an $n \times n$ lower triangular matrix with ones on the diagonal and $\boldsymbol \Sigma_t = \operatorname{diag}(\exp(h_{1,t}), \ldots, \exp(h_{n,t}))$. The law of motion of the VAR coefficients and log-volatilites will be specified below. Since the system in (ref) is written in the structural form, the covariance matrix $\boldsymbol \Sigma_t$ is diagonal by construction. Consequently, we can estimate this recursive system equation by equation without loss of efficiency.
We note that \citet*{CCM19} pioneer a similar equation-by-equation estimation approach for a large reduced-form constant-coefficient VAR with stochastic volatility. The main advantage of the structural-form representation is that it allows us to rewrite the VAR as $n$ unrelated regressions, and it leads to a more efficient sampling scheme. The main drawback of this representation, however, is that the implied reduced-form estimates depend on how the variables are ordered in the system. We will investigate the extent to which these estimates depend on the ordering in Section (ref).
It is convenience to introduce some notations. Let $b_{i,t}$ denote the $i$-th element of $\mathbf{b}_t$ and let $\mathbf{B}_{i,j,t}$ represent the $i$-th row of $\mathbf{B}_{j,t}$. Then, $\boldsymbol \beta_{i,t} = (b_{i,t},\mathbf{B}_{i,1,t},\ldots,\mathbf{B}_{i,p,t})'$ is the intercept and VAR coefficients of the $i$-th equation and is of dimension $k_{\beta} \times 1$ with $k_{\beta} = np+1$. Moreover, let $\boldsymbol \alpha_{i,t} $ denote the free elements in the $i$-th row of the contemporaneous impact matrix $\mathbf{A}_t$ for $i=2,\ldots, n$. That is, $\boldsymbol \alpha_{i,t} = (A_{i1,t},\ldots, A_{i(i-1),t})'$ is of dimension $k_{\alpha_i}\times 1$ with $k_{\alpha_i} = i-1$. Then, the $i$-th equation of the system in (ref) can be rewritten as: \[ y_{i,t} = \widetilde{\mathbf{x}}_t \boldsymbol \beta_{i,t} + \widetilde{\mathbf{w}}_{i,t}\boldsymbol \alpha_{i,t} + \varepsilon_{i,t}^y, \quad \varepsilon_{i,t}^y \sim \mathcal{N}(0, \text{e}^{h_{i,t}}), \] where $\widetilde{\mathbf{x}}_t = (1, \mathbf{y}_{t-1}',\ldots, \mathbf{y}_{t-p}')$ and $\widetilde{\mathbf{w}}_{i,t} = (-y_{1,t},\ldots, -y_{i-1,t})$. Note that $y_{i,t}$ depends on the contemporaneous variables $y_{1,t},\ldots, y_{i-1,t}$. But since the system is triangular, when we perform the change of variables from $\boldsymbol \varepsilon_t^y$ to $\mathbf{y}_t$ to obtain the likelihood function, the density function remains Gaussian.
If we let $\mathbf{x}_{i,t} = (\widetilde{\mathbf{x}}_t, \widetilde{\mathbf{w}}_{i,t})$, we can further simplify the $i$-th equation as:
where $\boldsymbol \theta_{i,t} = (\boldsymbol \beta_{i,t}', \boldsymbol \alpha_{i,t}')'$ is of dimension $k_{\theta_i} = k_{\beta} + k_{\alpha_i} = np+i.$ Hence, we have rewritten the TVP-VAR in (ref) as $n$ unrelated regressions. Finally, the coefficients and log-volatilities are assumed to evolve as independent random walks:
where the initial conditions $\boldsymbol \beta_{i,0}, \boldsymbol \alpha_{i,0}$ and $h_{i,0}$ are treated as unknown parameters to be estimated. The system in (ref)--(ref) specifies a reparameterization of a standard TVP-VAR in which all equations have time-varying parameters and stochastic volatility.
Note that the innovations in (ref)-(ref) are assumed to be independent across equations. This assumption is partly motivated by the concern of proliferation of correlation parameters, especially when $n$ is large, if the correlations of the innovations are unrestricted. In addition, for $\boldsymbol \beta_{i,t}$ and $\boldsymbol \alpha_{i,t}$, this independence assumption is important for extending the setup later so that we can turn on and off the time variation in both equations. In contrast, it is feasible to allow the innovations to $h_{i,t}$ to be correlated across equations (with a slight increase of computational cost). In preliminary work we considered such an extension. While the estimation results suggest that the correlation parameters are sizable, this extension leads to only very modest forecast gains (see Appendix D for details). Therefore, in what follows we maintain the independence assumption in (ref)-(ref) as the baseline.
Next, we introduce a framework that allows the model to determine in a data-driven fashion whether the VAR coefficients and the contemporaneous relations among the endogenous variables in each equation are time varying or constant. For that purpose, we adapt the non-centered parameterization of FSW10 to our hybrid TVP-VARs. More specifically, for $i=1,\ldots, n, t=1,\ldots, T,$ we consider the following model:
where $\widetilde{\boldsymbol \beta}_{i,0} = \mathbf{0}$ and $ \widetilde{\boldsymbol \alpha}_{i,0} = \mathbf{0}$. Here $\gamma_i^{\beta}$ and $\gamma_i^{\alpha}$ are indicator variables that take values of either 0 or 1.
The model in (ref)-(ref) includes a wide variety of popular VAR specifications. For example, assuming that all indicators take the value of 1, the above model is just a reparameterization of the TVP-VAR in (ref)--(ref). To see that, define $\boldsymbol \beta_{i,t} = \boldsymbol \beta_{i,0} + \gamma_i^{\beta}\boldsymbol \Sigma_{\beta_i}^{\frac{1}{2}}\widetilde{\boldsymbol \beta}_{i,t}$ and $\boldsymbol \alpha_{i,t} = \boldsymbol \alpha_{i,0} + \gamma_i^{\alpha}\boldsymbol \Sigma_{\alpha_i}^{\frac{1}{2}}\widetilde{\boldsymbol \alpha}_{i,t}$. Then, when $\gamma_i^{\beta} = \gamma_i^{\alpha} = 1, i=1,\ldots, n$, it is clear that (ref) becomes (ref). In addition, we have
Hence, $\boldsymbol \beta_{i,t}$ and $\boldsymbol \alpha_{i,t}$ follow the same random walk processes as in (ref) and (ref), respectively. We have therefore shown that when $\gamma_i^{\beta} = \gamma_i^{\alpha} = 1, i=1,\ldots, n$, the proposed model reduces to a TVP-VAR with stochastic volatility.
For the intermediate case where $\gamma_i^{\beta} = 1$ and $ \gamma_i^{\alpha} = 0, i=1,\ldots, n$, the proposed model reduces to a structural-form reparameterization of the model in CS05, i.e., a TVP-VAR with stochastic volatility but the contemporaneous relations among the endogenous variables are restricted to be constant. In the extreme case where $\gamma_i^{\beta} = \gamma_i^{\alpha} = 0, i=1,\ldots, n$, the proposed model then becomes a constant-coefficient VAR with stochastic volatility---a reparameterization of the specification in CCM19. More generally, by allowing the indicators $\gamma_i^{\beta}$ and $\gamma_i^{\alpha}$ to take different values, we can have a VAR in which only some equations have time-varying parameters. Note that it is straightforward to include a few additional indicators to allow for more flexible forms of time variation. For example, one can replace $\gamma_i^{\beta}$ with two indicators, say, $\gamma_i^{\beta, \text{own}}$ and $\gamma_i^{\beta, \text{other}}$, which control the time variation in the elements of $\boldsymbol \beta_{i,t}$ that correspond to coefficients on own lags and lags of other variables, respectively. The posterior simulator in Section (ref) can be modified to handle this case, with a slight increase in computation time.
These indicators are not fixed but are estimated from the data. More precisely, we specify that each $\gamma^\beta_i$ follows an independent Bernoulli distribution with success probability $\mathbb P(\gamma^\beta_i = 1) = p^{\beta}_i, i=1,\ldots, n$. Similarly for $\gamma^{\alpha}_i$: $\mathbb P(\gamma^{\alpha}_i = 1) = p^{\alpha}_i$. These success probabilities $p^{\beta}_i$ and $p^{\alpha}_i, i=1,\ldots, n$, are in turn treated as parameters to be estimated. In contrast to typical setups where time variation in parameters is assumed CS01,CS05,Primiceri05, here the proposed model puts positive probabilities in simpler models in which the VAR coefficients and the contemporaneous relations among the variables are constant. The values of the indicators are determined by the data, and these time-varying features are turned on only when they are warranted. The proposed model therefore is not only flexible in the sense that it includes a wide variety of specifications popular in applied work as special cases, it also induces parsimony to combat over-parameterization concerns.
The proposed hybrid TVP-VAR can also be viewed as a Bayesian model average of a wide variety TVP-VARs with different forms of time variation. To see that, let $\boldsymbol \gamma = (\boldsymbol \gamma_1,\ldots, \boldsymbol \gamma_n)'$ denote the vector of indicator variables with $\boldsymbol \gamma_i = (\gamma^\beta_i, \gamma^\alpha_i)$. Note that each value of $\boldsymbol \gamma\in\{0,1\}^{2n}$ corresponds to a particular TVP-VAR in which the time variation of the $i$-th equation is characterized by $\boldsymbol \gamma_i$. For example, $\boldsymbol \gamma = \mathbf{0}$ corresponds to a constant-coefficient VAR with stochastic volatility. Then, the posterior distribution of any model parameters under the proposed model can be represented as the posterior average with respect to $p(\boldsymbol \gamma\,|\, \mathbf{y})$, i.e., the posterior model probabilities of the collection of $2^{2n}$ TVP-VARs with different forms of time variation, where $\mathbf{y}$ denotes the data. For example, the joint distribution of $\boldsymbol \beta$ and $\boldsymbol \alpha$, the time-varying VAR coefficients and free elements of the contemporaneous impact matrix, can be represented as \[ p(\boldsymbol \beta,\boldsymbol \alpha\,|\,\mathbf{y}) = \sum_{\mathbf{c}\in\{0,1\}^{2n}} p(\boldsymbol \beta,\boldsymbol \alpha \,|\,\mathbf{y},\boldsymbol \gamma=\mathbf{c}) p(\boldsymbol \gamma=\mathbf{c}\,|\,\mathbf{y}). \] For a small VAR with $n=3$ variables (and the additional assumption that $\gamma_i^{\beta} = \gamma_i^{\alpha}$), CE18b estimate all $2^3=8$ TVP-VARs and the corresponding posterior model probabilities. For larger $n$, this approach of computing $p(\boldsymbol \gamma=\mathbf{c}\,|\,\mathbf{y})$ and sampling from $p(\boldsymbol \beta,\boldsymbol \alpha \,|\,\mathbf{y},\boldsymbol \gamma=\mathbf{c})$ for all $2^{2n}$ possible models is clearly infeasible. In contrast, by including the model indicator $\boldsymbol \gamma$ in the estimation, we simultaneously explore the parameter space and the model space. This latter approach is convenient and computationally feasible for large systems.
It is also instructive to investigate how the value of the model indicator $\boldsymbol \gamma$ is determined by the data. To fix ideas, suppose we wish to compare two TVP-VARs, represented as $\boldsymbol \gamma=\mathbf{c}_1$ and $\boldsymbol \gamma=\mathbf{c}_2$. Let $p(\mathbf{y} \,|\,\boldsymbol \gamma=\mathbf{c}_j)$ denote the marginal likelihood under model $\boldsymbol \gamma=\mathbf{c}_j, j=1,2$, i.e.,
where $\boldsymbol \psi_j$ is the collection of model-specific time-invariant parameters and time-varying states (in our setting these parameters and states are common across models and $\boldsymbol \psi_1=\boldsymbol \psi_2$), $p(\mathbf{y} \,|\, \boldsymbol \psi_j, \boldsymbol \gamma=\mathbf{c}_j)$ is the (complete-data) likelihood and $p(\boldsymbol \psi_j\,|\, \boldsymbol \gamma=\mathbf{c}_j)$ is the prior density. Then, the posterior odds ratio in favor of model $\boldsymbol \gamma=\mathbf{c}_1$ against model $\boldsymbol \gamma=\mathbf{c}_2$ is given by: \[ \frac{p(\boldsymbol \gamma=\mathbf{c}_1 \,|\,\mathbf{y})}{p(\boldsymbol \gamma=\mathbf{c}_2\,|\,\mathbf{y})} = \frac{p(\boldsymbol \gamma=\mathbf{c}_1)}{p(\boldsymbol \gamma=\mathbf{c}_2)}\times \frac{p(\mathbf{y}\,|\, \boldsymbol \gamma=\mathbf{c}_1)}{p(\mathbf{y}\,|\, \boldsymbol \gamma=\mathbf{c}_2)}, \] where $p(\boldsymbol \gamma=\mathbf{c}_1)/p(\boldsymbol \gamma=\mathbf{c}_2)$ is the prior odds ratio. It follows that if both models are equally probable a priori, i.e., $p(\boldsymbol \gamma=\mathbf{c}_1) = p(\boldsymbol \gamma=\mathbf{c}_2)$, the posterior odds ratio between the two models is then equal to the ratio of the two marginal likelihoods, or the Bayes factor. More generally, under the assumption that each TVP-VAR has the same prior probability, the value of the model indicator $\boldsymbol \gamma$ is determined by the marginal likelihood $p(\mathbf{y}\,|\, \boldsymbol \gamma)$. That is, if the TVP-VAR represented by $\boldsymbol \gamma=\mathbf{c}$ forecasts the data better (as one-step-ahead density forecasts), the value $\mathbf{c}$ will have a higher weight.
In this section we first describe in detail the priors on the time-invariant parameters. We then outline the posterior simulator to estimate the model described in (ref)--(ref)
For notational convenience, stack $\mathbf{y}_{i} = (y_{i,1},\ldots, y_{i,T})'$, $\boldsymbol \beta_{i} = (\boldsymbol \beta_{i,1}',\ldots, \boldsymbol \beta_{i,T}')'$, $\boldsymbol \alpha_{i} = (\boldsymbol \alpha_{i,1}',\ldots, \boldsymbol \alpha_{i,T}')'$ and $\mathbf{h}_{i} = (h_{i,1},\ldots, h_{i,T})'$ over $t=1,\ldots, T$, and collect $\mathbf{y} = \{\mathbf{y}_{i}\}_{i=1}^n$, $\boldsymbol \beta = \{\boldsymbol \beta_{i}\}_{i=1}^n $, $\boldsymbol \alpha = \{\boldsymbol \alpha_{i}\}_{i=1}^n $ and $\mathbf{h}=\{\mathbf{h}_{i}\}_{i=1}^n$ over $i=1,\ldots,n$, and similarly define $\widetilde{\boldsymbol \beta}_i$ and $\widetilde{\boldsymbol \alpha}_i$. Furthermore, let $\boldsymbol \Sigma_{\theta_i} = \text{diag}(\boldsymbol \Sigma_{\beta_i},\boldsymbol \Sigma_{\alpha_i})$ and $\boldsymbol \gamma_i = (\gamma^\beta_i, \gamma^\alpha_i)$. In our model, the time-invariant parameters are $\boldsymbol \gamma = (\boldsymbol \gamma_1,\ldots, \boldsymbol \gamma_n)'$, $\boldsymbol \Sigma_\theta=\{\boldsymbol \Sigma_{\theta_i}\}_{i=1}^n$, $\boldsymbol \Sigma_h = \{\sigma_{i,h}^2\}_{i=1}^n$, $\boldsymbol \theta_0= (\boldsymbol \theta_{1,0}',\ldots, \boldsymbol \theta_{n,0}')'$, $\mathbf{h}_0 = (h_{1,0},\ldots, h_{n,0})'$, $\mathbf{p}^{\beta}=(p^{\beta}_1,\ldots,p^{\beta}_n)'$ and $\mathbf{p}^{\alpha}=(p^{\alpha}_1,\ldots,p^{\alpha}_n)'$. Below we give the details of the priors on these time-invariant parameters.
Since $\boldsymbol \theta_0= (\boldsymbol \beta_0',\boldsymbol \alpha_0')'$, the initial conditions of the VAR coefficients, is high-dimensional when $n$ is large, appropriate shrinkage is crucial. We assume a Minnesota-type prior on $\boldsymbol \theta_0$ along the lines in SZ98; see also DLS84, litterman86 and KK97. We refer the readers to KK10, DS12 and karlsson13 for a textbook discussion of the Minnesota prior. More specifically, consider $\boldsymbol \theta_0\sim\mathcal{N}(\mathbf{a}_{\boldsymbol \theta_0},\mathbf{V}_{\boldsymbol \theta_0})$, where the prior mean $\boldsymbol \theta_0$ is set to be zero when the variables are in growth rate to induce shrinkage and the prior covariance matrix $\mathbf{V}_{\boldsymbol \theta_0}$ is block-diagonal with $\mathbf{V}_{\boldsymbol \theta_0}=\text{diag}(\mathbf{V}_{\boldsymbol \theta_{1,0}},\ldots,\mathbf{V}_{\boldsymbol \theta_{n,0}})$---here $\mathbf{V}_{\boldsymbol \theta_{i,0}}$ is the prior covariance matrix for $\boldsymbol \theta_{i,0}, i=1,\ldots, n$. For each $\mathbf{V}_{\boldsymbol \theta_{i,0}},$ we in turn assume it to be diagonal with the $k$-th diagonal element $(V_{\boldsymbol \theta_{i,0}})_{kk}$ set to be: \[ (V_{\boldsymbol \theta_{i,0}})_{kk} = \left\{
\right. \] where $s_r^2$ denotes the sample variance of the residuals from regressing $y_{r,t}$ on $\mathbf{y}_{t-1},\ldots,\mathbf{y}_{t-4}$, $r=1,\ldots, n$. Here the prior covariance matrix $\mathbf{V}_{\boldsymbol \theta_0}$ depends on four hyperparameters---$\kappa_1,\kappa_2, \kappa_3$ and $\kappa_4$---that control the degree of shrinkage for different types of coefficients. For simplicity, we set $\kappa_3 = 1$ and $\kappa_4 = 100$. These values imply moderate shrinkage for the coefficients on the contemporaneous variables and no shrinkage for the intercepts.
The remaining two hyperparameters are $\kappa_1$ and $\kappa_2$, which control the overall shrinkage strength for coefficients on own lags and those on lags of other variables, respectively. Departing from SZ98, here we allow $\kappa_1$ and $\kappa_2$ to be different, as one might expect that coefficients on lags of other variables would be on average smaller than those on own lags. In fact, CCM15 and chan21 find empirical evidence in support of this so-called cross-variable shrinkage. In addition, we treat $\kappa_1$ and $\kappa_2$ as unknown parameters to be estimated rather than fixing them to some subjective values. This is motivated by a few recent papers, such as CCM15 and \citet*{GLP15}, which show that by selecting this type of overall shrinkage hyperparameters in a data-based fashion, one can substantially improve the forecast performance of the resulting VAR. In addition, this data-based Minnesota prior is also found to forecast better than many recently introduced adaptive shrinkage priors such as the normal-gamma prior, the Dirichlet-Laplace prior and the horseshoe prior. For example, this is demonstrated in a comprehensive forecasting exercise in CHP20.
We assume gamma priors for the hyperparameters $\kappa_1$ and $\kappa_2$: $\kappa_j\sim \mathcal{G}(c_{1,j},c_{2,j}), j=1,2$. We set $c_{1,1} = c_{1,2} = 1 $, $c_{2,1} = 1/0.04$ and $c_{2,2} = 1/0.04^2$. These values imply that the prior modes are at zero, which provides global shrinkage. The prior means of $\kappa_1$ and $\kappa_2$ are 0.04 and $0.04^2$ respectively, which are the fixed values used in \citet*{CCM15}. Next, following FSW10, the square roots of the diagonal elements of $\boldsymbol \Sigma_{\theta_i} = \operatorname{diag}(\sigma_{\theta_i, 1}^2 , \ldots, \sigma_{\theta_i, k_{\theta_i}}^2), i=1,\ldots, n,$ are independently distributed as mean 0 normal random variables: $\sigma_{\theta_i, j}\sim \mathcal{N}(0,S_{\theta_i, j}), i=1,\ldots, n, j = 1,\ldots, k_{\theta_i}$. We assume each $\sigma_{h,i}^2$ follow a conventional inverse-gamma priors: $\sigma_{h, i}^2\sim \mathcal{IG}(\nu_{h,i},S_{h, i}), i=1,\ldots, n$. The success probabilities $p^{\beta}_i$ and $p^{\alpha}_i$ are assumed to have beta distributions: $p^{\beta}_i \sim\mathcal{B}(a_{p^\beta},b_{p^\beta})$ and $p^{\alpha}_i\sim\mathcal{B}(a_{p^{\alpha}},b_{p^{\alpha}}), i=1,\ldots, n$. Finally, the elements of the initial condition $\mathbf{h}_0$ are assumed to be Gaussian: $h_{i,0}\sim\mathcal{N}(a_{h_{i,0}}, V_{h_{i,0}})$.
We now turn to the estimation of the model in (ref)--(ref) given the prior described in the previous section. There are a few challenges in the estimation. First, since $\boldsymbol \beta_i$ becomes degenerate when $\gamma^{\beta}_i=0$, making its sampling nonstandard (similarly for $\boldsymbol \alpha_i$). To sidestep this problem, we will use the parameterization in terms of $\widetilde{\boldsymbol \beta}_i$ and $\widetilde{\boldsymbol \alpha}_i$. Then, given the posterior draws of $\widetilde{\boldsymbol \beta}_i$, $\widetilde{\boldsymbol \alpha}_i$ and other parameters, we can recover the posterior draws of $\boldsymbol \beta_i$ and $\boldsymbol \alpha_i$ using the definitions $\boldsymbol \beta_{i,t} = \boldsymbol \beta_{i,0} + \gamma_i^{\beta}\boldsymbol \Sigma_{\beta_i}^{\frac{1}{2}}\widetilde{\boldsymbol \beta}_{i,t}$ and $\boldsymbol \alpha_{i,t} = \boldsymbol \alpha_{i,0} + \gamma_i^{\alpha}\boldsymbol \Sigma_{\alpha_i}^{\frac{1}{2}}\widetilde{\boldsymbol \alpha}_{i,t}$.
Second, since $\widetilde{\boldsymbol \beta}_i$ and the indicator $\gamma_i^\beta$ enter the likelihood in (ref) multiplicatively, it is vital to sample them jointly (similarly for $\widetilde{\boldsymbol \alpha}_i$ and $\gamma_i^{\alpha}$); otherwise the Markov chain might get stuck. To see this, consider a simpler sampling scheme in which we simulate $\widetilde{\boldsymbol \beta}_i$ given $\gamma_i^\beta$, followed by sampling $\gamma_i^\beta$ given $\widetilde{\boldsymbol \beta}_i$. Suppose $\gamma_i^\beta = 0$ in the last iteration. Given $\gamma_i^\beta = 0$, $\widetilde{\boldsymbol \beta}_i$ does not enter the likelihood and we simply sample it from its state equation. Since the sampled $\widetilde{\boldsymbol \beta}_i$ has no relation to the data, the implied time variation in the VAR coefficients would not match the data. Consequently, it is highly likely that the model would prefer no time variation, i.e., $\gamma_i^\beta=0$. Hence, it is unlikely for the Markov chain to move away from $\gamma_i^\beta = 0$ once it is there. It is therefore necessary to sample both $\widetilde{\boldsymbol \beta}_i$ and $\gamma_i^\beta$ in the same step. In addition, since the pair $(\widetilde{\boldsymbol \beta}_i, \gamma_i^\beta)$ and $(\widetilde{\boldsymbol \alpha}_i, \gamma_i^\alpha)$ enters the likelihood additively, we sample them jointly to further improve efficiency.
Next, define $\widetilde{\boldsymbol \theta}_{i} = (\widetilde{\boldsymbol \theta}_{i,1}',\ldots, \widetilde{\boldsymbol \theta}_{i,T}')'$ with $\widetilde{\boldsymbol \theta}_{i,t} = (\widetilde{\boldsymbol \beta}_{i,t}',\widetilde{\boldsymbol \alpha}_{i,t}')'$. Then, one can simulate from the joint posterior distribution using the following posterior sampler that sequentially samples from:
Step 2 to Step 7 mainly involve standard sampling techniques and we leave the details to Appendix A. Here we focus on the first step.
Step 1. We sample the four blocks of parameters $\widetilde{\boldsymbol \beta}_i, \gamma_i^\beta, \widetilde{\boldsymbol \alpha}_i$ and $\gamma_i^\alpha$ jointly to improve efficiency. This is done by first drawing the indicators $\boldsymbol \gamma_i = (\gamma^\beta_i, \gamma^\alpha_i)$ marginally of $\widetilde{\boldsymbol \theta}_{i,t} = (\widetilde{\boldsymbol \beta}_{i,t}', \widetilde{\boldsymbol \alpha}_{i,t}')'$---but conditional on other parameters---and then sample $\widetilde{\boldsymbol \beta}_i$ and $\widetilde{\boldsymbol \alpha}_i$ from their joint conditional distribution. The latter of these two steps is straightforward because given $\gamma^\beta_i$ and $\gamma^\alpha_i$, we have a linear Gaussian state space model in $\boldsymbol \theta_{i,t}$. Specifically, we stack the observation equation (ref) over $t=1,\ldots, T$:
where $\boldsymbol \Sigma_{\mathbf{h}_i} = \text{diag}(\text{e}^{h_{i,1}},\ldots, \text{e}^{h_{i,T}})$, \[ \mathbf{X}_i =
, \quad \mathbf{Z}_{\boldsymbol \gamma_i} =
. \] Here note that the matrix $\mathbf{Z}_{\boldsymbol \gamma_i}$ depends on the indicators $\boldsymbol \gamma_i=(\gamma_i^{\beta}, \gamma_i^{\alpha})$. Next, stack the state equations (ref)-(ref) over $t=1,\ldots, T$: \[ \mathbf{H}_{k_{\theta_i}}\widetilde{\boldsymbol \theta}_i = \boldsymbol \varepsilon_i^{\widetilde{\theta}}, \quad \boldsymbol \varepsilon_i^{\widetilde{\theta}} \sim \mathcal{N}(\mathbf{0},\mathbf{I}_{T k_{\theta_i}}), \] where $\mathbf{H}_{k_{\theta_i}}$ is the first difference matrix of dimension $k_{\theta_i}= k_{\beta} + k_{\alpha_i}$. Since $\mathbf{H}_{k_{\theta_i}}$ is a square matrix with unit determinant, it is invertible. It then follows that $\widetilde{\boldsymbol \theta}_i \sim \mathcal{N}(\mathbf{0},(\mathbf{H}_{k_{\theta_i}}'\mathbf{H}_{k_{\theta_i}})^{-1}).$ Finally, using standard linear regression results, we have
where
Since the precision matrix $\mathbf{K}_{\widetilde{\boldsymbol \theta}_i}$ is a band matrix, one can sample $(\widetilde{\boldsymbol \theta}_i \,|\, \mathbf{y}_i, \mathbf{h}_i, \boldsymbol \Sigma_{\theta_i}, \boldsymbol \theta_{i,0}, \boldsymbol \gamma_i)$ efficiently using the algorithm in CJ09. It is worth noting that one could include an additional step to accept or reject the draw $\widetilde{\boldsymbol \theta}_i$ to ensure stationarity by checking the roots of the characteristic polynomial associated with the implied reduced-form VAR coefficients along the lines suggested in CS05.
To sample $\boldsymbol \gamma_i=(\gamma_i^{\beta},\gamma_i^{\alpha})$ marginal of $\widetilde{\boldsymbol \theta}_i$, it suffices to compute the four probabilities that $\boldsymbol \gamma_i = (0,0), \boldsymbol \gamma_i=(0,1), \boldsymbol \gamma_i=(1,0),$ and $\boldsymbol \gamma_i=(1,1)$. To that end, note that \[ p(\boldsymbol \gamma_i \,|\, \mathbf{y}_i, \mathbf{h}_i, \boldsymbol \Sigma_{\theta_i}, \boldsymbol \theta_{i,0}) \propto \left[ \int_{\mathbb{R}^{Tk_{\theta_i}}} p(\mathbf{y}_i\,|\, \widetilde{\boldsymbol \theta}_i, \mathbf{h}_i, \boldsymbol \Sigma_{\theta_i}, \boldsymbol \theta_{i,0},\boldsymbol \gamma_i) p(\widetilde{\boldsymbol \theta}_i) \text{d} \widetilde{\boldsymbol \theta}_i \right] p(\boldsymbol \gamma_i), \] where both the conditional likelihood $p(\mathbf{y}_i\,|\, \widetilde{\boldsymbol \theta}_i, \mathbf{h}_i, \boldsymbol \Sigma_{\theta_i}, \boldsymbol \theta_{i,0},\boldsymbol \gamma_i)$ and the prior density $p(\widetilde{\boldsymbol \theta}_i)$ are Gaussian. It turns out that the above integral admits an analytical expression. In fact, using a similar derivation in CG16, one can show that
where $\widehat{\widetilde{\boldsymbol \theta}}_i$ and $\mathbf{K}_{\widetilde{\boldsymbol \theta}_i}$ are defined in (ref). Then, one can compute the relevant probabilities using the expression in (ref). For example, when $\boldsymbol \gamma_i = (0,0)$, $|\mathbf{K}_{\widetilde{\boldsymbol \theta}_i}|=1$ and $\widehat{\widetilde{\boldsymbol \theta}}_i = \mathbf{0}$. It follows that
Similarly, we have
where $\widehat{\widetilde{\boldsymbol \theta}}_i(1,1) $ and $\mathbf{K}_{\widetilde{\boldsymbol \theta}_i}(1,1)$ denote respectively $\widehat{\widetilde{\boldsymbol \theta}}_i$ and $\mathbf{K}_{\widetilde{\boldsymbol \theta}_i}$ evaluated at $\boldsymbol \gamma_i=(1,1)$. The probabilities that $\boldsymbol \gamma_i = (0,1)$ and $\boldsymbol \gamma_i = (1,0)$ can be computed similarly. A draw from this 4-point distribution is standard once we normalize the probabilities. The details of the remaining steps are provided in Appendix A.
In this section we first conduct a series of simulated experiments to assess how well the posterior sampler works in recovering the time-varying structure in the data generating process. We then document the runtimes of estimating the hybrid TVP-VARs of different dimensions to assess how well the posterior sampler scales to larger systems.
First, we generate 300 datasets from the hybrid VAR in (ref)--(ref) with $n=12$ variables and sample size $T=200$, $T=400$ or $T=800$. We set the vector of indicators $\boldsymbol \gamma$ by repeating the four combinations $(0,0), (0,1), (1,0), (1,1)$ three times --- that allows us to study the effect of different combinations of time-varying pattens as well as their positions in the system. We generate $\boldsymbol \beta_0$, the initial conditions of the VAR coefficients, stochastically as follows. The intercepts are drawn independently from the uniform distribution on the interval $(-10,10)$, i.e., $\mathcal{U}(-10, 10)$. For the VAR coefficients, the diagonal elements of the first VAR coefficient matrix are iid $\mathcal{U}(0,0.5)$ and the off-diagonal elements are from $\mathcal{U}(-0.2,0.2)$. All other elements of the $j$-th ($j > 1$) VAR coefficient matrices are iid $\mathcal{N}(0,0.1^2/j^2).$ Finally, the elements of $\boldsymbol \alpha_0$ are drawn independently from $\mathcal{U}(-0.5, 0.5)$.
If the coefficient $\theta_{ij,t}$ is time-varying (i.e., the associated indicator $\gamma_i^{\alpha}$ or $\gamma_i^{\beta}$ is 1), it is generated from the state equation (ref) or (ref) with $\sigma_{\theta_{i},j}^2 = 0.01^2$ if $\theta_{ij,t}$ is a VAR coefficient and $\sigma_{\theta_{i},j}^2 = 0.1^2$ if it is an intercept for $i=1,\ldots,n, j=1,\ldots, k_{\theta_i}$. Finally, for the log-volatility processes, we draw $h_{0,i}\sim \mathcal{U}(-2,2)$ and set $\sigma_h^2 = 0.1, i=1,\ldots,n$.
In the Monte Carlo study we use the priors described in Section (ref) with the following hyperparameters. The prior means of the initial conditions $\boldsymbol \theta_0$ and $\mathbf{h}_0$ are set to be zero $\mathbf{a}_{\boldsymbol \theta_0} = \mathbf{0}$ and $\mathbf{a}_{h}=\mathbf{0}$, and the prior covariance matrix of $\mathbf{h}_0$ is $\mathbf{V}_{h} = 10\times\mathbf{I}_n$. The hyperparameter of $\sigma_{\theta_i,j}$ is set so that the implied prior mean of $\sigma_{\theta_i,j}$ is $0.01^2$ if it is associated with a VAR coefficient and $0.1^2$ for an intercept. Finally, we set the hyperparameters of $\gamma_i^{\beta}$ and $\gamma_i^{\alpha}$ to be $a_{p^\beta} = b_{p^\beta} = a_{p^\alpha} = b_{p^\alpha} = 0.5$. These values imply that the prior modes are at 0 and 1, whereas the prior mean is 0.5.
Given a dataset and the priors described above, we estimate the hybrid VAR using the posterior sampler in Section (ref) and obtain the posterior mode of $\boldsymbol \gamma$. We repeat this procedure for all the datasets and compute the frequencies of $\gamma^\beta_i$ and $\gamma^\alpha_i$ being one, $i=1,\ldots, n$. The results are reported in Table (ref).
Overall, the posterior sampler works well and is able to recover the true time-varying structure in the simulated data on average. While it is harder to pin down the correct value of $\gamma^\alpha_i$ compared to $\gamma^\beta_i$, the frequencies of identifying the true value of $\gamma^\alpha_i$ are still reasonably good, even for a small sample of $T=200$. In addition, these results substantially improve when the sample size increases from $T=200$ to $T=800$. All in all, these Monte Carlo results confirm that the proposed hybrid model can recover salient patterns---such as time-varying conditional means and covariances---in the data.
To further investigate the effect of the beta prior on $\gamma_i^{\beta}$ and $\gamma_i^{\alpha}$, we repeat the Monte Carlo experiments but assume a uniform prior on the unit interval $(0,1)$, i.e., $\gamma_i^{\beta}, \gamma_i^{\alpha}\sim \mathcal{B}(1,1) = \mathcal{U}(0,1).$ Hence, both the prior means and modes are 0.5. The Monte Carlo results are similar to the baseline case and they are reported in Appendix D.
Next, we document the runtimes of estimating the hybrid TVP-VARs of different sizes to assess how well the posterior sampler scales to higher dimensions. More specifically, Table (ref) reports the runtimes (in minutes) to obtain 1,000 posterior draws from the hybrid models with $n= 10, 20, 30$ variables and $T=400, 800$ time periods. The posterior sampler is implemented in $\mathrm{M}\mathrm{{\scriptstyle ATLAB}}$ on a standard desktop with an Intel Core i7-7700 @3.60 GHz processor and 64 GB memory. As a comparison, we also include the corresponding runtimes of fitting the TVP-VAR of Primiceri05 using the algorithm in DP15. Note that the algorithm in DP15 samples all the time-varying VAR coefficients $\boldsymbol \beta$ in one block and it tends to be very computationally intensive for larger systems. One potential solution is to develop an equation-by-equation estimation procedure similar to that in CCCM21. Since the algorithm is designed for models with a constant contemporaneous impact matrix, extending it to handle the TVP-VAR of Primiceri05---which features a time-varying contemporaneous impact matrix---would be an interesting future research direction.
It is evident from the table that for typical applications with 15-30 variables, the proposed model can be estimated reasonably quickly. In addition, using the recursive representation that admits straightforward equation-by-equation estimation, fitting the proposed model is much faster than estimating the TVP-VAR of Primiceri05, even though the former is more flexible.
In this section we fit a large US macroeconomic dataset set to demonstrate the usefulness of the proposed model. After describing the dataset in Section (ref), we first investigate how different variable orderings affect the estimates from the proposed hybrid TVP-VAR relative to the TVP-VAR of Primiceri05 in Section (ref). We then present the full sample results in Section (ref). In particular, we conduct a formal Bayesian model comparison exercise to shed light on the time-varying patterns of the model parameters. We then consider a pseudo out-of-sample forecasting exercise in Section (ref). We show that the forecast performance of the proposed model compares favorably to a range of standard benchmarks.
The US dataset for our empirical application consists of 20 quarterly variables with a sample period from 1959Q1 to 2018Q4. It is sourced from the FRED-QD database at the Federal Reserve Bank of St. Louis as described in MN21. Our dataset contains a variety of standard macroeconomic and financial variables, such as Real GDP, industrial production, inflation rates, labor market variables, money supply and interest rates. They are transformed to stationarity, typically to annualized growth rates. The complete list of variables and how they are transformed is given in Appendix C.
We use the priors described in Section (ref). In particular, since the data are transformed to growth rates, we set the prior mean of $\boldsymbol \theta_0$ to be zero, i.e., $\mathbf{a}_{\boldsymbol \theta_0} = \mathbf{0}$. For the prior hyperparameters on $\kappa_1$ and $\kappa_2$, we set $c_{1,1} = c_{1,2} = 1 $, $c_{2,1} = 1/0.04$ and $c_{2,2} = 1/0.04^2$. These values imply that the prior means of $\kappa_1$ and $\kappa_2$ are respectively 0.04 and $0.04^2$. For the hyperparameters of the initial conditions $\mathbf{h}_0$, we set $\mathbf{a}_{h}=\mathbf{0}$ and $\mathbf{V}_{h} = 10\times\mathbf{I}_n$. Next, the hyperparameters of $\sigma_{h,i}^2$ are set so that the prior mean is $0.1$. Similarly, the hyperparameters of $\sigma_{\theta_i,j}$ are chosen so that the implied prior mean is $0.01^2$ if it is associated with a VAR coefficient and $0.1^2$ for an intercept. Finally, we set $a_{p^\theta} = b_{p^\theta} = a_{p^h} = b_{p^h} = 0.5$. These values imply prior modes at 0 and 1, whereas the prior mean is 0.5.
Since the proposed hybrid TVP-VAR is written in the recursive structural form---which is used as a computational device and not as an identification scheme---one naturally wonders how the assumed order of the variables affects the model estimates compared to standard TVP-VARs such as the model in Primiceri05. Conceptually, the proposed hybrid TVP-VAR is not order invariant due to two components. First, the VAR coefficients are in structural form, and their priors induce priors on the reduced-form parameters that depend on the order of the variables. Second, the multivariate stochastic volatility specification is constructed based on the lower triangular impact matrix $\mathbf{A}_t$. And since priors are independently elicited for $\mathbf{A}_t$ and the stochastic volatility, the implied prior on the covariance matrix is not order invariant.
Popular TVP-VARs such as CS05 and Primiceri05 share the second component but not the first (since they are formulated as reduced-form VARs); see, e.g., the discussion in Primiceri05 and CCM19. Hence, one might expect that estimates from the proposed hybrid TVP-VAR would be more sensitive to how the variables are ordered compared to those of Primiceri05. On the other hand, as discussed in Section (ref), the proposed hybrid TVP-VAR can be viewed as a Bayesian model average of a wide variety of VARs with many different forms of time variation. Since some of these VARs are more parsimonious and have restricted time variation---including the model of CCM19 that is less prone to the ordering issue---the resulting Bayesian model average estimates could in principle be more robust to different orderings compared to Primiceri05. Hence, whether the ordering issue is more severe in the proposed hybrid TVP-VAR relative to Primiceri05 is an empirical question, and we investigate this issue below.
It is worth noting that we choose the TVP-VAR of Primiceri05 to be our benchmark, even though it is not order invariant, because it is generally viewed as the state-of-the-art and it is widely used as a reduced-form VAR for both forecasting DGG13 and structural analysis using non-recursive identification schemes benati2008, BP13. There are a few recent papers that aim to develop order-invariant VARs with multivariate stochastic volatility, such as Bognanni18, SZ20, ARRS21 and CKY21. These models, however, are either designed for small TVP-VARs or they do not feature time-varying VAR coefficients.
Now, we investigate how different variable orderings impact the estimates from the proposed hybrid TVP-VAR relative to the TVP-VAR of Primiceri05. To that end, we use 6 variables---real GDP, PCE inflation, unemployment, Fed funds rate and industrial production and real average hourly earnings in manufacturing---and we consider all $6!=720$ possible orderings. For each ordering of these 6 variables, we fit the proposed model and obtain the fitted values (i.e., the time-varying conditional means). We then compute the mean squared errors (against the observed values) of the 6 variables. We repeat this exercise for the TVP-VAR of Primiceri05. Since there are 720 MSEs for each variable and each model, to better summarize the results we report the boxplots of the MSEs in Figure (ref). The middle line of each box denotes the median, while the lower and upper lines represent, respectively, the 25- and the 75-percentiles. The whiskers extend to the maximum and minimum.
Since the goal of this exercise is to compare the variability of the fitted means, we normalize the MSEs by the medians (so that the red line of each boxplot is one). Overall, the variability of the estimates from the proposed hybrid TVP-VAR is comparable to that of the model of Primiceri05. It is also interesting to note that for the majority of the variables the variability is relatively small. One exception is the unemployment rate---this is partly due to the very small base rate (e.g., the MSE of the unemployment is less than 1% of the MSE of the real GDP).
Next, we investigate the variability of the variance estimates due to different orderings. For that purpose, we obtain the fitted variance of each variable in each ordering, and compute the mean squared error against the `realized volatility' values (defined as $\text{RV}_{i,t} = (y_{i,t} - \widehat{y}_{i,t})^2,$ where $\widehat{y}_{i,t}$ is the fitted value from the regression of $y_{i,t}$ on an intercept and $y_{i,t-1}, \ldots, y_{i,t-4}$). The results are reported in Figure (ref). Again, the results show that the variability of the variance estimates from the proposed model is comparable to that of the TVP-VAR of Primiceri05, even though the former uses a recursive structural-form representation.
We also investigate the variability of point and density forecasts similar to the exercise in ARRS21. More specifically, for each variable ordering, we compute the root mean squared forecast error and the average of log predictive likelihoods (see Section 5.4 for more details) for each of the 6 variables from the proposed model. Consistent with the results in ARRS21, the variability of point forecasts is relatively small. In addition, the variability of the forecasts from the proposed model is similar to a full-fledged TVP-VAR where all the VAR coefficients and elements of the impact matrix are time varying. The details are reported in Appendix D.
In this section we report the full sample results of the hybrid TVP-VAR fitted using all $n=20$ variables, which are ordered as listed in Table (ref). Of particular interest are the posterior means of $\gamma_i^{\beta}$ and $\gamma_i^{\alpha}$, the indicators that control the time variation in the VAR coefficients and the elements of the impact matrix, respectively. These estimates are reported in Table (ref). The results clearly show that while many estimates are essentially 0, others are close to 1. In other words, while the time variation in many equations is essentially turned off, there is strong evidence for time-varying parameters in some equations. Since there is substantial heterogeneity in the time-variation pattern across equations, conventional approaches of either assuming time variation in all equations or restricting all parameters to be constant are unlikely to fit the time-variation pattern well.
In addition, the estimates of $\gamma_i^{\beta}$ and $\gamma_i^{\alpha}$ for the same equation are often of very different magnitudes, suggesting that time variation in one group of parameters does not necessarily imply time variation in the other. These results thus confirm the usefulness of having two separate indicators for each equation. Overall, while we find evidence for time variation in VAR coefficients and elements of the impact matrix in some equations, not all equations need both forms of time variation. These results therefore highlight the empirical relevance of the proposed hybrid TVP-VAR.
To assess the efficiency of the proposed posterior sampler, we compute the inefficiency factors of the posterior draws, which are reported in Appendix D. The values of the inefficiency factors are comparable to those of conventional TVP-VARs. The results thus show that the posterior sampler is efficient in terms of producing posterior draws that are not highly autocorrelated.
Next, we consider a Bayesian model comparison exercise to compare the proposed hybrid TVP-VAR to various VARs with different time-variation patterns. In general, to compare models using the Bayes factor, one would require the computation of the marginal likelihood given in (ref). Despite recent advances, obtaining the marginal likelihood for high-dimensional models with multiple latent states remains a nontrivial task. Fortunately, one much simpler approach is available when one wishes to compare nested models. More specifically, for nested models, the Bayes factor can be calculated using the Savage-Dickey density ratio VW95, which requires only the estimation of the unrestricted model. More importantly, no explicit computation of the marginal likelihood is needed. This approach has been used to compute the Bayes factor in many empirical applications, including KP99, DS09 and \citet*{KLS10}.
More specifically, suppose we wish to compare the proposed hybrid TVP-VAR against a TVP-VAR characterized by the vector of indicators $\boldsymbol \gamma = \mathbf{c}\in\{0,1\}^{2n}$ (e.g., $\mathbf{c} = \mathbf{0}$ represents a constant-coefficient VAR with stochastic volatility.) Since the latter is a restricted version of the former, the Bayes factor in favor of the unrestricted model can be obtained using the Savage-Dickey density ratio via \[ \text{BF}_{\text{u},\mathbf{c}} = \frac{p(\boldsymbol \gamma=\mathbf{c})}{p(\boldsymbol \gamma=\mathbf{c}\,|\, \mathbf{y})}, \] where the numerator and the denominator are, respectively, the marginal prior and posterior densities of $\boldsymbol \gamma$ evaluated at $\mathbf{c}$. Hence, to compute the relevant Bayes factor, one only needs to evaluate two densities at a point. Intuitively, if $\boldsymbol \gamma =\mathbf{c}$ is less likely under the posterior density relative to the prior density, i.e., $p(\boldsymbol \gamma=\mathbf{c}\,|\, \mathbf{y}) < p(\boldsymbol \gamma=\mathbf{c})$, then it is viewed as evidence against the restriction $\boldsymbol \gamma =\mathbf{c}$ and the unrestricted model is favored with $\log\text{BF}_{\text{u},\mathbf{c}} > 0$. Any two nested models can be compared similarly. We provide the technical details on evaluating the two densities at $\boldsymbol \gamma =\mathbf{c}$ in Appendix B.
We compare the proposed hybrid TVP-VAR to four VARs with stochastic volatility. For each of these VARs, $\gamma_1^{\beta}= \cdots = \gamma_n^{\beta}$ and $\gamma_1^{\alpha}=\cdots = \gamma_n^{\alpha}$, and they are denoted as HYB-$(\gamma_i^{\beta},\gamma_i^{\alpha})$. The first is a full-fledged TVP-VAR where all the VAR coefficients and elements of the impact matrix are time varying with $\gamma_i^{\beta} = \gamma_i^{\alpha} = 1, i=1,\ldots, n$; this is a structural-form version of the TVP-VAR in Primiceri05 and we denote this model as HYB-$(1,1)$. The second is a VAR with time-varying VAR coefficients but a constant impact matrix, i.e., $\gamma_i^{\beta} = 1, \gamma_i^{\alpha} = 0, i=1,\ldots, n$, which we denote as HYB-$(1,0)$; this is a structural-form parameterization of the TVP-VAR in CS05. The third is a constant-coefficient VAR with $\gamma_i^{\beta} = \gamma_i^{\alpha} = 0, i=1,\ldots, n$, which we denote as HYB-$(0,0)$; this is a structural-form variant of the VAR in CCM19. Lastly, we also consider a version in which the VAR coefficients are constant but the elements of the impact matrix are time varying, i.e., $\gamma_i^{\beta} = 0, \gamma_i^{\alpha} = 1, i=1,\ldots, n$, which we denote as HYB-$(0,1)$.
Table (ref) reports the log Bayes factors of the proposed hybrid TVP-VAR against these four VARs with stochastic volatility. It is clear that there is overwhelming support for the hybrid TVP-VAR relative to the alternatives. For instance, for $n=20$ the Bayes factor in favor of the hybrid TVP-VAR against the best alternative HYB-$(0,1)$ is $\text{e}^{401} \approx 1.42\times 10^{174}$, suggesting the posterior model probability of the hybrid TVP-VAR is practically 1. These model comparison results are consistent with the estimates of $\boldsymbol \gamma$ reported in Table (ref), which indicate substantial heterogeneity in the time-variation pattern across equations---hence, fixing these indicators to either 0 or 1 for all equations is likely to be too restrictive.
To investigate the support for the hybrid TVP-VAR across model dimensions, we also consider a small system ($n=3$) with only real GDP, PCE inflation and unemployment, as well as a medium system ($n=6$) with three additional variables: Fed funds rate, industrial production and real average hourly earnings in manufacturing. For both model dimensions, the hybrid TVP-VAR is strongly favored by the data compared to the alternatives.
It is interesting to note that the best model among the four benchmarks changes across model dimensions. More specifically, for the small system with $n=3$ variables, the best model is HYB-$(0,0)$, the constant-coefficient VAR with stochastic volatility, even though it is the most restrictive among the four VARs. This suggests that the additional flexibility in allowing time-varying coefficients does not sufficiently fit the data better to justify the added model complexity. This is inline with the results in CE18, who find that, for fitting a 3-variable US dataset, a constant-coefficient VAR with stochastic volatility performs well relative to various VARs with different forms of time variation according to the marginal likelihood.
While HYB-$(0,0)$ is the preferred model for the small system, when more variables are included HYB-$(0,1)$ is strongly favored. For instance, for $n=20$, allowing for time variation in the impact matrix increases the log marginal likelihood by 2,139---comparing HYB-$(0,0)$ and HYB-$(0,1)$. But further allowing for time variation in the VAR coefficients reduces the log marginal likelihood by 602---comparing HYB-$(0,1)$ and HYB-$(1,1)$. Taken together, these results suggest that it is generally useful to allow for time variation in the impact matrix. In contrast, adding time variation in the VAR coefficients should be done judiciously in large systems: while simply allowing all VAR coefficients to be time-varying can be detrimental, using the proposed data-based approach to add time variation equation-wise can substantially improve model-fit.
Next, we evaluate the forecast performance of the proposed hybrid TVP-VAR relative to a few standard benchmarks. In particular, we consider 1) a conventional homoscedastic and constant-coefficient VAR; 2) a constant-coefficient VAR with stochastic volatility and a constant impact matrix (by setting all $\gamma_i^{\beta}$ and $\gamma_i^{\alpha}$ to 0; denoted as HYB-$(0,0)$); and 3) a full-fledged TVP-VAR (by setting all the indicators to 1; denoted as HYB-$(1,1)$). All models use all $n=20$ variables. The sample period is from 1959Q1 to 2018Q4, and the evaluation period starts at 1985Q1 and runs till the end of the sample.
We perform a recursive forecasting exercise using an expanding window. More specifically, in each forecasting iteration $t$, we use only data up to time $t$, denoted as $\mathbf{y}_{1:t}$, to estimate the models. We then evaluate both point and density forecasts. We use the conditional expectation $\mathbb E(y_{i,t+m}\,|\, \mathbf{y}_{1:t})$ as the $m$-step-ahead point forecast for variable $i$ and the predictive density $p(y_{i,t+m} \,|\, \mathbf{y}_{1:t})$ as the corresponding density forecast.
The metric used to evaluate the point forecasts from model $M$ is the root mean squared forecast error (RMSFE) defined as \[ \text{RMSFE}_{i,m}^M = \sqrt{\frac{\sum_{t=t_0}^{T-m}( y_{i,t+m}^{\text{o}} - \mathbb E(y_{i,t+m}\,|\, \mathbf{y}_{1:t}))^2}{T-m-t_0+1}}, \] where $y_{i,t+m}^{\text{o}}$ is the actual observed value of $y_{i,t+m}$. For RMSFE, a smaller value indicates better forecast performance. To evaluate the density forecasts, the metric we use is the average of log predictive likelihoods (ALPL): \[ \text{ALPL}_{i,m}^M = \frac{1}{T-m-t_0+1}\sum_{t=t_0}^{T-m} \log p(y_{i,t+m}= y_{i,t+m}^{\text{o}}\,|\, \mathbf{y}_{1:t}), \] where $p(y_{i,t+m}= y_{i,t+m}^{\text{o}}\,|\, \mathbf{y}_{1:t})$ is the predictive likelihood. For this metric, a larger value indicates better forecast performance.
To compare the forecast performance of model $M$ against the benchmark $B$, we follow CCM15 to report the percentage gains in terms of RMSFE, defined as \[ 100 \times (1 - \text{RMSFE}_{i,m}^M/\text{RMSFE}_{i,m}^B), \] and the percentage gains in terms of ALPL: \[ 100 \times (\text{ALPL}_{i,m}^M - \text{ALPL}_{i,m}^B). \]
Figure (ref) reports the forecasting results of the hybrid TVP-VAR, where we use a conventional homoscedastic, constant-coefficient VAR as the benchmark. The top panel shows the percentage gains in RMSFE for all 20 variables, and the bottom panel presents the corresponding results in ALPL.
For both 1- and 4-quarter-ahead point forecasts, the hybrid TVP-VAR outperforms the benchmark for almost all variables (all but two for 1-quarter-ahead and one for 4-quarter ahead). For a few variables, such as the federal funds rate, real personal consumption expenditure and industrial production, the hybrid TVP-VAR outperforms the benchmark by more than 10% for 1-quarter-ahead forecasts (the differences in forecast performance are also statistically significant at 0.05 level according to the test of DM95). Overall, the median percentage gains in RMSFE for 1- and 4-quarter-ahead forecasts are, respectively, 5.1% and 6.0%.
For density forecasts, the hybrid TVP-VAR performs even better relative to the benchmark---it outperforms the benchmark for all variables in both forecast horizons. The median percentage gains in ALPL for 1- and 4-quarter-ahead forecasts are 14% and 13%, respectively. Moreover, for many variables the percentage gains are more than 20%. These results are consistent with numerous studies in the small VAR literature, such as clark11, \citet*{DGG13} and CR15, that show allowing for time-varying structures substantially improves forecast performance compared to VARs with constant parameters, especially for density forecasts.
Next, we compare the forecast performance of the hybrid TVP-VAR with that of HYB-$(0,0)$, the constant-coefficient VAR with stochastic volatility. The results are reported in Figure (ref). For both 1- and 4-quarter-ahead point forecasts, the hybrid TVP-VAR outperforms the HYB-$(0,0)$ for most variables. The median percentage gains in RMSFE are 1.7% and 4.1%, respectively. For density forecasts, the results are similar: the median percentage gains in ALPL for 1- and 4-quarter-ahead forecasts are, respectively, 0.8% and 3.1%. Overall, these results suggest that allowing for time variation in VAR coefficients---with appropriate shrinkage and sparsification---can further enhance the forecast performance of a VAR with stochastic volatility.
Finally, Figure (ref) compares the forecast performance of the the hybrid TVP-VAR with that of HYB-$(1,1)$, the full-fledged TVP-VAR where all the VAR coefficients and error variances are time varying. Again, for both point and density forecasts, the hybrid TVP-VAR performs better than the benchmark for most variables. In particular, the median percentage gains in RMSFE for 1- and 4-quarter-ahead forecasts are 1.0% and 1.4%, respectively; the median percentage gains in ALPL are 2.4% and 2.8%, respectively. These results suggest that imposing time variation in all equations is not necessary and would adversely impact the forecast performance.
Overall, these forecasting results show that the proposed hybrid TVP-VAR forecasts better than many state-of-the-art time-varying models. These forecasting results highlight the advantages of using a data-driven approach to discover the time-varying structures---rather than imposing either constant coefficients or time variation in parameters.
To better understand the sources of forecast gains, Figure (ref) compares the forecast performance of HYB-$(1,1)$, the full-fledged TVP-VAR, to HYB-$(0,0)$, the constant-coefficient VAR with stochastic volatility. The results are mixed: while HYB-$(1,1)$ does slightly better in terms of point forecasts for the majority of the variables, it performs worse in terms of density forecasts for many variables. These results suggest that allowing for time-varying VAR coefficients in all equations does not necessarily improve forecast performance. (The forecast performance of HYB-$(1,1)$ and HYB-$(1,0)$ are very similar, as reported in Appendix D). This finding is also consistent with the full-sample estimation results presented in Table (ref): while the data clearly favors time-varying VAR coefficients in a few equations, for the majority of the equations time variation is not needed.
Finally, one can also interpret the superior forecast performance of the hybrid TVP-VAR through the lens of Bayesian forecast combinations MZ93, AK08. As discussed in Section (ref), the proposed hybrid TVP-VAR can be viewed as a Bayesian model average of a wide variety TVP-VARs with different forms of time variation, where each component model is characterized by the vector of indicators $\boldsymbol \gamma$. Then, the forecasts from the hybrid TVP-VAR can be interpreted as a forecast combination weighted by the posterior model probabilities $p(\boldsymbol \gamma\,|\, \mathbf{y})$. Consistent with the large literature on forecast combinations, here we find that the forecast combination of the hybrid TVP-VAR performs better than many individual component models, including HYB-$(1,1)$ and HYB-$(0,0)$.
This paper has developed what we call hybrid TVP-VARs, i.e., VARs with time-varying parameters in some equations but not in others. Using US data, we found evidence that while VAR coefficients and error covariances in some equations are time varying, the data prefers constant coefficients in others. In a forecasting exercise that involves 20 macroeconomic and financial variables, we demonstrated the superior forecast performance of the proposed hybrid TVP-VARs compared to standard benchmarks.
In future work, it would be interesting to use hybrid TVP-VARs for structural analysis. Since they are formulated in the recursive structural-form, structural analysis using large VARs identified by recursive zero restrictions, such as the application in ER17, can directly use the proposed models.
In addition, developing an order-invariant version of these hybrid TVP-VARs would be an interesting and important extension. This would involve changing two components. First, one requires the use of the reduced-form VAR representation, but model indicators $\boldsymbol \gamma$ can be introduced similarly. An equation-by-equation estimation procedure can be developed along the lines in CCCM21, though computation would be more intensive. Second, one would need to replace the Cholesky stochastic volatility model with an order-invariant model suitable for large VARs, such as CCM16 or CKY21. The trade-off between the two stochastic volatility models is between computational speed and model flexibility: the former can be estimated quickly but the latter is more flexible. Finding a good modeling approach among all these choices---with an eye on the additional computational costs---would be an interesting research direction.