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.
86,581 characters · 14 sections · 38 citation commands
General Bayesian time-varying parameter VARs for predicting government bond yields
\thispagestyle{empty}
\doublespacing
Time-varying parameter vector autoregressive (TVP-VAR) models are commonly used in finance and macroeconomics to capture dynamic relations across variables, regime shifts and/or structural changes in economic processes primiceri2005,cogley2005drifts,dh2012. These models typically assume that the parameters evolve over time according to a simple stochastic process such as a random walk. While being rather flexible and parsimonious, this assumption does not allow to examine the extent to which covariates cause changes in the time-varying parameters (TVPs) over time. Moreover, wrongly assuming a random walk state equation could negatively impact predictive accuracy because it essentially implies a smoothness prior on the coefficients. This might be at odds with the rapid shifts we have observed in financial time series such as bond yields and thus could negatively impact predictive accuracy.
The literature has dealt with the issue of selecting the appropriate law of motion for time-varying parameters (TVPs) by estimating different models separately and then using model selection criteria to discriminate between competing specifications sims2006were, koop2009evolution, h2020. But, to the best of our knowledge, no attempt has been made to develop models that rely on a large set of competing laws of motion for the coefficients and decide which one describes the data best.
This paper proposes a flexible approach to TVP-VARs that efficiently integrates out uncertainty with respect to the state evolution equation. The approach assumes that the TVPs depend on a potentially large panel of covariates. These covariates are commonly labeled effect modifiers which can be partially latent and might feature their own state equations. In case they are observed, we obtain a model that is closely related to the varying coefficient model originally proposed in hastie1993varying. The main advantage of using observed, as opposed to latent effect modifiers, is that we can investigate the driving sources of parameter change. This feature is important if the researcher is interested in investigating why relations between variables in a VAR change over time, and to what extent these changes are explained by the observed effect modifiers.
Careful selection of these effect modifiers is crucial. On the one hand, deciding on the appropriate set of observed quantities is difficult, and a large set of candidates could arise. One key objective of this paper is to provide techniques to select promising subsets. On the other hand, appropriately selecting the latent effect modifiers allows us to capture situations where parts of the coefficients evolve smoothly whereas others move more abruptly. The latter behavior of the TVPs is often found for US macroeconomic data sims2006were, whereas the former is consistent with financial time series such as bond yields or stock returns dh2012, huber2019should.
Since large VARs often include both macroeconomic as well as financial quantities, a successful model should be able to accommodate both types of structural change or even rely on linear combinations of them. This is precisely what we aim to achieve in this paper. Our approach is capable of answering not only the question why coefficients change, but also to infer an appropriate law of motion using a broad set of latent quantities, each equipped with its own state evolution equation. These unobserved quantities range from latent factors that follow a random walk to Markov switching indicators that allow a subset of parameters to switch between a low number of regimes. Our model approach nests several alternatives proposed in the literature such as the TVP-VAR of primiceri2005 and cogley2005drifts or the reduced-rank model of chan2020reducing.
This large degree of flexibility, however, comes with two concerns. The first is that overfitting problems can easily arise. We overcome these by using Bayesian shrinkage priors. Our prior is a variant of the well-known Horseshoe prior carvalho2010horseshoe that allows us to shrink coefficients associated with irrelevant effect modifiers towards zero. The second concern relates to computation. Since inclusion of a large number of endogenous variables in a VAR quickly leads to a huge dimensional parameter space, we propose a computationally efficient Markov chain Monte Carlo (MCMC) algorithm. To circumvent mixing issues, we rely on different parameterizations of the model during MCMC sampling. The corresponding novel algorithm thus provides a second important contribution of the paper.
In our empirical work, we use the model approach to forecast the US term structure of interest rates. We investigate the empirical properties of our approach using two information sets in the underlying VAR. First, we include several interest rates at different maturities directly as endogenous variables. Second, we consider a three-factor Nelson-Siegel model for the term structure of interest rates as in diebold2006forecasting and diebold2008global. Adopting a long hold-out period that includes several recessionary episodes, our approach improves upon a wide range of competing models. While improvements for point forecasts are often muted, we find that our proposed model yields favorable density predictions. The predictive exercise is complemented by a comprehensive discussion of patterns in time-variation and their sources. Moreover, how our approach can be used to analyze low frequency relations between the observed quantities in the model over time.
The rest of the paper is structured as follows. Section (ref) introduces the econometric framework which includes the general form of the TVP-VAR, a flexible law of motion for the latent states as well as the effect modifiers which crucially impact the state dynamics. This section, moreover, introduces the Bayesian prior setup and techniques for posterior and predictive inference. Section (ref) applies the model approach to the term structure of US interest rates. It also serves to illustrate key model features and to highlight the predictive capabilities of the approach in an out-of-sample forecasting exercise. The last section summarizes and concludes the paper. Additional technical details and further empirical results are provided in the Appendix.
Let $\bm y_t$ denote an $M \times 1$ vector of macroeconomic and financial quantities at time $t=1, \dots, T$. We assume that $\bm y_t$ depends on its $P$ lags which we store in a $K=MP$-dimensional vector $\bm x_t = (\bm y'_{t-1}, \dots, \bm y'_{t-P})'$. Then the basic TVP-VAR can be written as a linear multivariate regression model:
where $\bm \beta_t$ represents a set of $k = MK$ dynamic regression coefficients and $\bm \epsilon_t \sim \mathcal{N}(\bm 0_M, \bm \Sigma_t)$ is a vector Gaussian shock process with time-varying $M \times M$-dimensional variance-covariance matrix $\bm \Sigma_t$. We assume that $\bm \Sigma_t$ can be decomposed as follows:
Here, $\bm Q_t$ denotes an $M \times M$ lower triangular matrix with unit diagonal with $v (= M (M-1))$ free elements denoted by $\bm q_t$. $\bm H_t = \text{diag}(e^{h_{1t}}, \dots, e^{h_{Mt}})$ is a diagonal matrix with $h_{jt}~(j=1, \dots, M)$ representing time-varying log-volatilities. These are assumed to evolve according to an AR(1) process:
$\mu_j$ is the unconditional mean, $\psi_j$ the persistence parameter and $\varsigma_j^2$ the variance of the log-volatility process for equation $j$.
In what follows, we rewrite the TVP-VAR using its non-centered parameterization sfs_wagner2010:
$\bm \beta$ denotes a $k$-dimensional vector of constant coefficients, and $\tilde{\bm \beta}_t = \bm \beta_t - \bm \beta$. This parameterization allows us to disentangle time-invariant (encoded by $\bm \beta$) from time-varying effects (encoded by $\tilde{\bm \beta}_t$) for the regressors. For the decomposed variance-covariance matrix, we have $\bm{Q}$, a lower triangular matrix with ones on the diagonal capturing the constant part of the covariances. $\tilde{\bm{Q}}_t = \bm{Q}_t - \bm{Q}$ is the corresponding lower triangular matrix with zero-diagonal elements containing the time-varying part. Their free elements are collected in the $v\times1$-vectors $\bm{q}$ and $\tilde{\bm{q}}_t$, respectively.
In this paper, the focus is on modeling the $N (=k+v)$-dimensional vector $\tilde{\bm{\gamma}}_t = (\tilde{\bm \beta}_t',\tilde{\bm{q}}_t')'$ and the constant part ${\bm{\gamma}} = ({\bm \beta}',{\bm{q}}')'$. The literature typically assumes that the transition distribution $p(\tilde{\bm \gamma}_t| \tilde{\bm \gamma}_{t-1})$ is given by:
This law of motion suggests that the expected value of $\tilde{\bm \gamma}_{t-1}$ equals $\tilde{\bm \gamma}_{t}$ and the amount of time-variation is determined by the $N \times N$-dimensional process innovation variance-covariance matrix $\bm V$. This matrix is often assumed to be diagonal. Notice that if selected elements in $\bm V$ are equal to zero, the corresponding regression coefficients are constant.
Estimation and inference is typically carried out via Bayesian methods. The recent literature proposes using shrinkage priors to allow for data-based selection of those coefficients which should be time-varying or constant. This already leads to substantial improvements in predictive accuracy but does not tackle the fundamental question whether the coefficients are better characterized by a random walk, a change-point process, or by mixtures of these.
As discussed in the previous sub-section, the typical assumption is that $\tilde{\bm \gamma}_t$ follows a random walk process. In addition, the shocks to the random walk state equation are often assumed to feature a positive error variance. We relax both assumptions to allow for more flexibility.
The random walk assumption is relaxed by assuming that the time-varying part stored in $\tilde{\bm{\gamma}}_t$ depends on a set of $R$ additional factors $\bm z_t$. These $\bm z_t$ are the effect modifiers mentioned in the introduction that can be observed or latent. The relationship between the TVPs and $\bm z_t$ is given by:
$\bm \Lambda$ denotes an $N \times R$ matrix of regression coefficients, and $\bm \eta_t \sim \mathcal{N}(\bm 0_N, \bm \Omega)$ is a Gaussian error term with diagonal error variance-covariance matrix $\bm \Omega = \text{diag}(\omega^2_1, \dots, \omega^2_N)$. If $R \ll N$, the coefficients feature a factor structure and co-move according to the effect modifiers in $\bm z_t$. The relationship between $\tilde{\bm{\gamma}}_t$ and $\bm z_t$ is determined by the factor loadings in $\bm \Lambda$. For instance, if the $j^{\text{th}}$ column of $\bm \Lambda, \bm \lambda_j$, is equal to zero, the corresponding $j^{\text{th}}$ factor in $\bm z_t$ does not enter the model and thus has no influence on $\tilde{\bm{\gamma}}_t$.
The specific selection of $\bm z_t$ is crucial for determining the dynamics of $\tilde{\bm{\gamma}}_t$. Appropriate choice of $\bm z_t$ yields a variety of important special cases that depend on the specific values of $\bm \Lambda$ and $\bm \Omega$ as well as on the composition of $\bm z_t$. In this sub-section, we briefly focus on special cases that arise independently of the choice of $\bm z_t$. The next sub-section deals with cases that arise if $\bm z_t$ is suitably chosen.
These two cases are the following. If $\bm \Lambda = \bm 0_{N \times R}$, with $\bm 0_{N \times R}$ being a $N \times R$ matrix of zeros, we obtain a random coefficients model that assumes that the regression coefficients follow a white noise process korobilis2019high, hhko2019. The second special case arises if both $\bm \Lambda = \bm 0_{N \times R}$ and $\bm \Omega = \bm 0_{N \times N}$. In this case, we obtain a standard constant parameter regression model.
Before we discuss the choice of $\bm z_t$, it is worth noting that if $\bm z_t$ is (partially) latent, the model in Eq. ((ref)) is not identified. Since our object of interest is $\bm \gamma_t$, this poses no greater issues. If we wish to structurally interpret elements in $\bm z_t$, standard identification strategies from the literature on dynamic factor models can be used geweke1996measuring,aguilar2000bayesian, stock2011dynamic.
The specific choice of $\bm z_t$ is crucial in determining how $\tilde{\bm{\gamma}}_t$ behaves over time. Hence, by suitably choosing the elements in $\bm z_t$, our model approach is related to the following specifications:
These examples show that our model, conditional on choosing a suitable set of effect modifiers, is capable of mimicking several prominent specifications in the literature. Since the question on the appropriate state evolution equation is essentially a model selection issue, we simply specify $\bm z_t$ to include most (with the exception of the $R=N$ setup) of the modifiers discussed above.
More precisely, we set $\bm z_t$ as follows:
Here we let $\bm r_t$ denote a set of $R_r$ observed factors and the dimension of $\bm z_t$ is thus $R = R_r + R_S + R_\tau$ with $R_S = M$. To allow for additional flexibility we assume that $\bm \tau_t = (\bm \tau_{1t}', \dots, \bm \tau_{Mt}')'$ with $\bm \tau_{jt}$ being equation-specific factors of dimension $R_{\tau j}$ (and thus $R_\tau = \sum_j R_{\tau j}$). Assuming that $R_{\tau i} = R_{\tau j} = \delta$ for all $i,j$, we have $R_\tau = \delta M$ latent random walk factors. Likewise, we estimate a separate Markov switching indicator $S_{jt}$ per equation (and thus $\bm S_t = (S_{1t}, \dots, S_{Mt})'$, with the corresponding transition probabilities matrix denoted by $\bm{P}_j$).
The corresponding loadings matrix $\bm \Lambda$ is structured such that the loadings in equation $j$ associated with the factors $\bm \tau_{it}$ and $S_{it}$ for $i \neq j$ equal zero. This assumption strikes a balance between assuming a large number of latent factors to achieve maximum flexibility (and thus risk overfitting) and using a rather parsimonious model (with the risk of being too simplistic). Recent contributions use similar assumptions on the state evolutions, koop2009evolution, maheu2018efficient. As opposed to these papers, our approach offers more flexibility since, if necessary, the presence of the idiosyncratic shocks to the TVPs allows for deviations if the factor structure does not represent the data well.
Our specification implies that, depending on the factor loadings $\bm \Lambda$, the evolution of $\tilde{\bm{\gamma}}_t$ might be a combination of a set of random walk factors, a Markov switching process and some observed quantities. To single out irrelevant elements in $\bm z_t$, one could simply set the corresponding columns in $\bm \Lambda$ equal to zero. In this paper, we achieve this through a Bayesian shrinkage prior. The next sub-section discusses our priors in more detail.
The discussion in Sub-section (ref) shows that our model approach nests a variety of competing models. To select the appropriate model variant and alleviate over-parameterization concerns, we opt for a Bayesian approach to introduce shrinkage. Here, we summarize the priors we impose on key parameters.
In light of the specific choice of $\bm z_t$, we introduce some additional notation to clarify details on our prior implementation. Let us assume that $\bm \Lambda$ is composed of the following matrices:
where $\bm \Lambda_r$ is an $N \times R_r$ matrix of loadings related to the observed quantities, $\bm \Lambda_{S}$ denotes an $N \times R_{S}$ matrix of loadings related to $\bm{S}_t$, and $\bm \Lambda_\tau$ represents an $N \times R_\tau$ factor loadings matrix associated with $\bm \tau_t$.
For imposing shrinkage we rely on variants of the horseshoe prior carvalho2010horseshoe. While in principle any global-local shrinkage prior may be used, we choose the horseshoe prior due to its excellent shrinkage properties and the lack of tuning parameters. In particular, we specify a column-wise horseshoe prior on the loadings matrix. Let $\bm{\Lambda}_j$ denote a sub-matrix of the free elements in $\bm{\Lambda}$ corresponding to the $j^{\text{th}}$ equation, and $\bm{\lambda}_{ji}$ mark the $i^{\text{th}}$ column of this matrix. $\lambda_{ji,\ell}$ refers to the $\ell^{\text{th}}$ element of this vector. The prior is given by:
Here, $\mathcal{C}^{+}(0, 1)$ denotes the half Cauchy distribution, and $\delta_{ji}$ is an equation- and column-specific global shrinkage factor, while $\kappa_{ji,\ell}$ is a local scaling parameter.
To further regularize our potentially huge-dimensional parameter space, we impose an equation-wise horseshoe prior on the constant part of the regression coefficients and covariances in $\bm{\gamma}_j$ corresponding to the $j^{\text{th}}$ equation. Let $\gamma_{ji}$ denote the $i^{\text{th}}$ element of the vector. The setup is similar to the one of the loadings matrix, and given by the hierarchical structure:
The hyperparameter $\zeta_j$ is an equation-specific global shrinkage factor, while the $\xi_{ji}$'s are local scalings.
We mentioned earlier that it is often assumed that shocks to the states feature positive error variances. To introduce shrinkage of $\omega_{ii}$ towards zero, we impose a horseshoe prior also on the square root of the innovation variances of the measurement errors in Eq. ((ref)). This prior is specified in an equation-specific manner. For equation $j$, let $\bm{\omega}_j$ denote a $v_j (=j-1 + k)$-dimensional vector which stores the diagonal elements in $\bm{\Omega}$ associated with the $j^{\text{th}}$ equation. This includes the process innovation variances on the $k$ regression coefficients and the $j-1$ covariance parameters in $\tilde{\bm Q}_t$. The square root of the $i^{\text{th}}$ element of $\bm \omega_j$, $\sqrt{\omega_{ji}}$, features the following prior hierarchy:
Choosing a Gaussian prior on the square root of the variance in the first level of the hierarchy implies a Gamma prior on $\omega_{ji}$, with $\omega_{ji}|\varpi_{ji},\vartheta_j\sim\mathcal{G}\left(1/2,\varpi_{ji}^{-2}\vartheta_j^{-2}/2\right)$, see also sfs_wagner2010. The hyperparameters $\vartheta_j$ and $\varpi_{ji}$ are again equation-specific global and local shrinkage parameters. Furthermore, we set $\bm{V}_\tau = \bm{I}_{R_\tau}$ and thus impose shrinkage through the factor loadings in $\bm{\Lambda}_\tau$ chan2020reducing.
On the parameters of the state equation of the log-volatility processes $\mu_j, \psi_j$ and $\varsigma_j^2$, we use the setup proposed in kastner2014ancillarity. That is, we assume a Gaussian prior on the unconditional mean, $\mu_j \sim \mathcal{N}(0, 10^2)$, a Beta prior on the transformed autoregressive parameter, $(\psi_j + 1)/2 \sim \mathcal{B}(5, 1.5)$, a Gamma prior on the state variance, $\varsigma_j^2 \sim \mathcal{G}(1/2, 1/2)$, while the prior on the initial state is the unconditional distribution for all $j=1,\dots, M$.
For the equation-specific transition probabilities $\bm{P}_j$ of Markov switching indicators, we assume that the $(i, i)^{\text{th}}$ element $p_{j,ii}$ arises from a Beta distribution given by:
and hence $p_{j, i\ell} = 1 - p_{j,i\ell}$ for $i \neq \ell$. In the empirical application, we define $e_{00} = e_{11} = 10$ and $e_{01} = e_{10} = 1$, in order to weakly push each $S_{jt}$ towards a single-state a priori.
To simulate from the full posterior distribution we develop an efficient MCMC algorithm. Since full-system estimation of the VAR quickly becomes computationally cumbersome, we rely on the equation-by-equation algorithm suggested in ccm2019.
Conditional on $\bm Q_t$, one can state the VAR as a system of (conditionally) independent equations. The first equation of this system is given by:
and the $j^{\text{th}}$ equation $ (j > 1)$:
${\bm{\beta}}_{j}$ and $\tilde{\bm{\beta}}_{jt}$ denote the $j^{\text{th}}$ subvectors of the constant and time-varying parts in $\bm \beta_t$ with $\bm \beta_t = (\bm \beta'_{1t}, \dots, \bm \beta'_{Mt})'$ and $\bm u_{jt} = (\varepsilon_{1t}, \dots, \varepsilon_{j-1, t})'$. The $(j-1)$-dimensional vectors $\bm q_j$ and $\tilde{\bm{q}}_{jt}$ store the constant and time-invariant part of the free elements in the $j^{\text{th}}$ row of $\bm Q_t$. This approach allows to estimate the different elements of $\bm \gamma_t$ that relate to the $M$ equations independently from each other conditional on the shocks to the preceding $j-1$ equations. This speeds up computation enormously.
Equation ((ref)) can be simplified to yield:
where $\bm m_{jt} = (\bm x'_t, \bm u_{jt}')'$, $\tilde{\bm \gamma}_{jt} = \bm \gamma_{jt} - \bm \gamma_{j}$, and $\bm \gamma_{jt}$ refers to the TVPs associated with the $j^{\text{th}}$ equation in $\bm \gamma_t$, and $\bm{\gamma}_j$ denotes the corresponding constant part. All the following steps will be carried out on an equation-by-equation basis and making use of the regression form in Eq. ((ref)). For notational simplicity, we assume that all elements in $\bm z_t$ are latent. In light of the discussion in Sub-section (ref), this implies that $\bm z_{jt} = (\bm \tau'_{jt},~ S_{jt})'$ and the extension to include observed factors is trivial.
Sampling $\bm z_{jt}$. Conditional on the remaining quantities of the model, we simulate the latent (random walk and Markov switching) components in $\bm z_{jt}$ by integrating out $\tilde{\bm{\gamma}}_{jt}$. This is achieved by rewriting Eq. ((ref)) as:
with $\tilde{y}_{jt} = y_{jt} - \bm m'_{jt} \bm \gamma_{j}$. Defining $\tilde{\bm m}'_{jt} = \bm m'_{jt} \bm \Lambda_j$ and $\hat{\varepsilon}_{jt} = \bm m'_{jt} \bm \eta_{jt} + \varepsilon_{jt}$ allows us to cast Eq. ((ref)) as a simple linear regression model:
This parameterization has the advantage that it does not depend on $\tilde{\bm \gamma}_{jt}$, and $\bm z_{jt}$ can thus be sampled marginally of $\tilde{\bm \gamma}_{jt}$. This improves mixing substantially since $\tilde{\bm \gamma}_{jt}$ and $\bm z_{jt}$ will often be highly correlated gerlach2000efficient, giordani2008efficient.
Depending on the precise law of motion for the elements in $\bm z_{jt}$, standard algorithms can now be used. In this paper, we use two different law of motions. For the latent random walk factors in $\bm \tau_{jt}$, we use the forward filtering backward sampling algorithm outlined in carter1994gibbs and fs1994. In case of the latent Markov switching factors in $S_{jt}$, we use the algorithm outlined in kim1999state. Both algorithms are well known and relevant details may be found in the original papers. Here, it suffices to note that in both cases, sampling the latent states is computationally easy since the state space is low dimensional with $R \ll N$. In this setting, sampling the factors equation-wise can be carried out in ${O}(R^3)$ steps, a substantial computational improvement relative to the ${O}(N^3)$ steps necessary to estimating an unrestricted TVP regression chan2020reducing.\footnote{It is worth mentioning that the ${O}(N^3)$ statement is true for the precision sampler and differs for forward-filtering backward-sampling algorithms.}
Sampling the state innovation variances. To obtain draws for the state innovation variances, reconsider Eq. ((ref)) and draw them conditional on the observed/latent states and the factor loadings using a generalized inverse Gaussian distribution. For further details and the moments of this distribution, see Appendix (ref).
Sampling $\bm \Lambda_j$ and $\bm \gamma_{j}$ jointly. Similarly to $\bm z_{jt}$, we sample the non-zero loadings in $\bm \Lambda$ and the time-invariant coefficients marginally of $\bm \gamma_{jt}$ by using equation-by-equation estimation. The observation equation for equation $j$ (conditional on $\bm z_{jt}$) can be written as a standard regression model:
where $\hat{\bm m}_{jt} = (\bm m_{jt}', (\bm z_{jt} \otimes \bm m_{jt})')'$ is an $R v_j^2$-dimensional vector of covariates, and $ \hat{\bm \gamma}_{j} = (\bm \gamma_{j}', \text{vec}(\bm \Lambda_j)')'$ denoting an $R v_j^2$-dimensional coefficient vector. The posterior of $\hat{\bm \gamma}_{j}$ is Gaussian with well known moments.
Sampling the stochastic volatilities. The latent log-volatility processes can again be sampled on an equation-by-equation basis. This step is implemented using the R-package stochvol.
Sampling the horseshoe prior hyperparameters. Our assumptions imply analogous horseshoe priors for the factor loadings matrix, the constant part of the coefficients and the square root of the state innovation variances. Posteriors are provided in Appendix (ref).
This section applies the model to predict the term structure of US interest rates. These time series are characterized by substantial non-linearities (e.g. during the period of the zero lower bound), feature substantial co-movement both in the level of the time series but also in the parameters describing their evolution. Our proposed model framework might thus be well suited to capture such features. We investigate this claim in a thorough forecasting exercise using several established benchmarks. After showing that our approach yields favorable forecasts, we discuss the driving forces behind parameter changes as well as discuss how key quantities that shape yield curve dynamics co-move over time at low frequencies.
Our aim is to predict monthly zero-coupon yields of US treasuries at different yearly maturities. The data is described in detail in gurkaynak2007us.\footnote{Available online at \href{https://www.federalreserve.gov/data/nominal-yield-curve.htm}{federalreserve.gov/data/nominal-yield-curve.htm}.} The target variables are $1,3,5,7,10 \text{ and } 15$ years maturities. All variables enter our model in first differences.
Estimation and forecasting is carried out recursively. Using data from 1973:01 to 1999:12, we produce one-month-ahead and three-months-ahead forecasts for 2000:01. After obtaining the predictive distributions we expand the sample and repeat this procedure until we reach 2019:12.\footnote{Note that for our set of financial indicators, data revisions and ragged edges arising from delays in the publication of the series do not matter. This is due to financial market data being available almost instantaneously, and the published quotes are not subject to revisions at later dates.} Point forecast performance is measured using Root Mean Squared Errors (RMSEs), while density forecasts are assessed in terms of Log Predictive Bayes Factors (LPBFs), averaged over the out-of-sample observations geweke2010comparing.
For analyzing the role of observed factors governing joint dynamics of the TVPs, we rely on three exogenous variables:
We use these three effect modifiers for simple reasons. First, there is strong evidence that yield curve dynamics differ across business cycle phases hevia2015estimating. Second, the RF interest rate serves as an early warning indicator which possesses predictive power for changes in the shape of the yield curve. Finally, the inclusion of the NFCI is motivated by the recent literature on forecasting tail risks in macroeconomic and financial time series adrian2019vulnerable, carriero2020capturing, adams2021forecasting.
The forecast exercise distinguishes between two model classes, with 15 distinct model specifications in each class. The first model class involves TVP-VARs that incorporate the six target variables as endogenous variables, that is $M=6$. This model class is labeled as VAR.
The second model class includes specifications based on the three-factor Nelson-Siegel (NS) model as in diebold2006forecasting that imposes a factor structure on the yields:
where $i_t(\theta)$ denotes the yield at maturity $\theta$ at time $t$, $L_t$ is a factor that controls the level, $\text{\textit{\v{S}}}_t$ determines the slope, and $\text{\textit{\c{C}}}_t$ represents the curvature factor of the yields. The parameter $\alpha$ governs the exponential decay rate. To maximize the loading on $\text{\textit{\c{C}}}_t$ we set $\alpha=0.7308\ (12\times0.0609)$.\footnote{See diebold2006forecasting for a discussion of this specific choice.} In what follows, we use the latent factors $L_t,\text{\textit{\v{S}}}_t \text{ and } \text{\textit{\c{C}}}_t$ as endogenous variables in the VAR specifications by defining $\bm{y}_t=(L_t, \text{\textit{\v{S}}}_t, \text{\textit{\c{C}}}_t)^\prime$, resulting in $M=3$. These latent factors are obtained by running OLS on a $t$-by-$t$ basis. This model class is subsequently labeled NS-VAR.
Model specifications are differentiated over a grid of effect modifier combinations. In particular, specifications within a model class differ in terms of three aspects (see Table (ref)): First, in terms of whether the three exogenous variables (collected in $\bm{r}_t$) are included in $\bm{z}_t$ or not ("x" marks inclusion, "--" indicates no observed factors); second, in terms of the number of latent random walk factors $R_{\tau j}$ included in $\bm{z}_t$ (which we assume to be equal across equations); and third, in terms of the presence of Markov switching indicators in $\bm{S}_t$, again with "x" marking their inclusion and "--" their absence. This setup implies that we have 15 time-varying parameter NS-VAR and 15 time-varying parameter VAR model specifications. For comparative purposes we also consider the two class-specific constant parameter model variants, labeled “Constant,” and a conventional independent random walk specification of the TVPs (i.e. we set the number of random walk factors equal to $K$ and exclude $\bm r_t$ and $\bm S_t$). Notice that we also have a specification which includes only $\bm S_t$ and another one which uses only observed factors. The latter one is closely related to a Markov switching model whereas the second one closely resembles a VAR with interaction terms.
Table (ref) shows the one-month and one-quarter-ahead out-of-sample forecasting results for US treasury yields at different maturities, using the 30 TVP model specifications and the two constant parameter model variants as described in the previous sub-section. Recall that $M=3$ (and $K=9$) in case of the NS-VAR and $M=6$ (and $K=18$) in case of the VAR model variants. $\bm{z}_t=(\bm{r}_t^\prime, \bm{S}_t^\prime, \bm{\tau}_t^\prime)$ where $\bm{r}_t$ denotes $R_r(=3)$-dimensional vector, $\bm{S}_t$ a $R_{S}$-dimensional vector, and $\bm{\tau}_t$ a $R_\tau$-dimensional vector. All specifications feature $P=3$ lags of the endogenous variables.
\afterpage{ \onehalfspacing
}
The performance of point forecasts is measured in terms of RMSEs and that of density forecasts in terms of LPBFs, relative to the constant parameter VAR model (shaded in yellow) that serves as benchmark. RMSEs are presented in ratios, and LPBFs in differences are given below the RMSEs in parentheses. RMSEs below one indicate superior performance, relative to the benchmark, as LPBF figures greater than zero do. The best performing model specification by column is given in bold, highlighting the specification with the smallest RMSE ratio and the largest positive LPBF difference, respectively.
The vast range of competing specifications, loss functions used to evaluate forecasts and maturities makes it hard to identify a single best performing model. We first provide a general overview on model performance and then zoom into differences in predictive accuracy for point and density forecasts.
At a very general level, Table (ref) suggests a pronounced degree of heterogeneity in forecast accuracy across models and for both the NS-VAR and VAR specifications. While differences at some maturities are substantial, they are muted or non-existent for others. It is also worth mentioning that the TVP variants of the NS-VAR and VAR models outperform the constant parameter specifications in most cases (apart from one-quarter ahead point forecasts of treasuries with a maturity of five years). This general observation holds true irrespective of whether only point forecasts or the full predictive distribution are considered. These accuracy premia point towards the necessity of addressing structural breaks in the dynamic evolution of the yield curve.
Comparing the NS-VAR and VAR specifications indicates that the latter usually perform better for density forecasts, while the former are often superior in terms of point forecasts. The better point forecasting performance of the NS-VAR suggests that the three factors contain relevant information for the first moment of the predictive distribution. When we also consider higher order moments, this story changes. The better density performance of the VAR models is most likely driven by two sources. The first is that, as opposed to the conditional mean, the strong implicit assumptions of the NS-VAR on the prediction error variance-covariance matrix seems overly restrictive. Allowing for richer dynamics in the covariances by explicitly modeling the shocks to a panel of yields thus yields better density forecasts. The second fact is that the small-scale NS-VARs might feature a too tight predictive variance since relevant information is ignored. And this could harm density forecasting accuracy during turbulent times.
Next, we zoom into specific model classes. Within these, differences in performance along sub-divisions (provided by the inclusion/exclusion of the exogenous variables or the Markov switching processes and the number of latent random walk factors) are often negligible. This finding indicates that our flexible approach to probabilistically selecting the most adequate state evolution via Bayesian shrinkage is successful, and thus not susceptible to overfitting concerns.
We now turn to considering model-specific forecasting accuracy with the aim to find the best performing models for point and density forecasts and for the two different forecast horizons. The overall winner for point forecasts at the one-month horizon, on average, is the flexible NS-VAR model specification featuring the exogenous variables, one latent factor per equation and the Markov switching processes. While this specification yields RMSEs that are $23$ percent lower than those of the benchmark specification, it must be acknowledged that most competing specifications exhibit values that are similar in magnitude. The same is true for the best-performing model specification at the one-quarter-ahead horizon, the one-factor NS-VAR without exogenous variables and Markov switching processes, albeit at much smaller margins versus the benchmark. Here, improvements are about nine percent in terms of relative RMSEs. Assessing predictive performance for individual maturities allows to identify which segment of the yield curve drives the overall results. While gains for shorter maturities in the case of one-month-ahead forecasts are muted, we observe large gains at the long end of the yield curve. Relative to the best performing specification, improvements are about $40$ percent in RMSE terms for the best performing specification. The same is true for one-quarter-ahead forecasts, however, at smaller margins of about $20$ percent.
We proceed with our findings for density forecasts. As mentioned above, the VAR specifications overall exhibit more favorable relative LPBFs compared to the NS-VAR. This is easily observable by noting that most bold values (in parentheses) are located in the lower panel of Table (ref). In terms of average performance at the one-month horizon, we find that the TVP-VAR specification with one factor but no exogenous variables or Markov switching performs best, closely followed by the most flexible specification with one unobserved factor including Markov switching. This again serves as an example that even though it is extremely flexible, our approach to shrinking the parameter space avoids overfitting and does not harm predictive performance. In fact, for one-quarter-ahead forecasts, the TVP-VAR with exogeneous variables, one latent factor and Markov switching shows the largest gains in terms of density forecasts. Again, it must be acknowledged that margins within this model class are rather small. As it is the case for point forecasts, these gains are mostly obtained in terms of the long end of the yield curve, while gains in forecast accuracy at shorter horizons are negligible.
Summing up, while improvements relative to established models are often small, our proposed framework is competitive for most maturities at both the one-month and one-quarter ahead horizons. Main differences arise from the considered model class, with the NS-VAR exhibiting promising results in terms of point forecasts, and superior density forecasts for the VAR. Moreover, we detect the largest improvements at the long end of the yield curve. The additional flexibility of our proposed approach does not lead to severe overfitting since we use shrinkage priors to regularize several parts of the parameter space. And it almost never harms forecast accuracy while improving performance in some cases.
\doublespacing
The previous sub-section established that the proposed approach yields more favorable predictions than conventional TVP-VARs. Including observed and latent effect modifiers allows for investigating the sources of time-variation in coefficients, and thus the driving factors of improvements in predictive accuracy. We carry out a detailed analysis of these determinants in this sub-section.
Because of its favorable forecasting properties, we choose the NS-VAR with one latent factor per equation to investigate the driving forces of parameter variation over the full estimation sample. Recall that the choice of $\bm z_t$ for this specification translates into a single equation-specific latent factor ($\tau_{jt}$), an equation-specific Markov switching indicator $S_{jt}$ and the three observed early warning indicators.
To illustrate the observed and latent factors, we transform each column vector in $\bm z_t$ ex-post such that it is bounded between zero and one. This allows for comparisons even though some blocks of the respective matrices are not econometrically identified.\footnote{In particular, this is the case for the product $\bm \Lambda_\tau \bm \tau_t$. Note that we do not face this issue for $\bm \Lambda_r \bm r_t$ and $\bm \Lambda_S \bm S_t$ because both are either observed or already bounded between zero and one, and their scale and sign are identified.} We therefore exploit the fact that introducing any invertible $R \times R$-dimensional matrix does not alter the likelihood of the model, since $\bm \Lambda \bm U^{-1} \bm U \bm z_t = \bm \Lambda \bm z_t$. Define $\bm U$ as diagonal matrix such that the maximum of each $\bm z_j$ (for $j = 1, \dots, R$) corresponds to one, the minimum is zero, $\tilde{\bm \Lambda} = \bm \Lambda \bm U^{-1}$ and $\tilde{\bm z}_t = \bm U \bm z_t$. This simple linear transformation allows for assessing the relative movement of the indicators in $\tilde{\bm z}_t$ without affecting overall dynamics.
Figure (ref) displays the evolution of the normalized indicators in $\tilde{\bm z}_t$ over time. First, we focus on the features of the observed effect modifiers depicted in the upper panel (a). The normalized NFCI index peaks during the oil crisis in the mid $1970$s, while exhibiting a stable evolution at a quite low level during the Great Moderation (the period from $1990$ until to the onset of the global financial crisis in $2007$). Spikes in $r_{\text{NFCI}, t}$ tend to coincide with recessionary episodes, indicated by the binary recession indicator $r_{\text{REC},t}$. The risk-free interest rate $r_{\text{RF},t}$ can be related to the monetary policy stance. Early in the sample we observe substantial increases, peaking during the Volcker disinflation in the early $1980$s. Subsequently, large and abrupt decreases are notable during recessions, while an overall decreasing trend is observable. Before turning to the latent indicators in $\bm S_t$ and $\bm \tau_t$, note that the respective elements $S_{jt}$ and $\tau_{jt}$ are equation-specific with $j\in\{L, \text{\textit{\v{S}}},\text{\textit{\c{C}}}\}$ referring to the level, slope and curvature of the yield curve. The middle panel (b) indicates the posterior median of the three Markov switching factors collected in $\bm S_t$. The lower panel (c) shows the posterior medians of the three (transformed) gradually changing latent random walk factors in $\bm \tau_t$.
Several features of the latent indicators are worth highlighting. Each of the latent quantities exhibits distinct dynamics and thus carries information in addition to the observed indicators. Examining the posterior means, which are essentially the unconditional posterior probabilities that a given Markov indicator equals one, we observe that both $S_{L,t}$ and $S_{\text{\textit{\v{S}}},t}$ evolve comparatively gradually over time. Before the Volcker disinflation, both tend to increase with posterior medians above $0.5$, indicating that regime $1$ is more likely. After $1985$, we observe a major shift towards regime $0$. Interestingly, for $S_{\text{\textit{\v{S}}},t}$ this transition appears immediately after $1985$, while we detect a notable delay in $S_{L,t}$.
During the Great Moderation both indicators tend to remain associated with regime $0$. The indicator associated with the middle segment of the yield curve, $S_{\text{\textit{\c{C}}},t}$, by contrast, transitions between regimes at a higher frequency. Particularly during the Volcker disinflation we observe mixed patterns and no clear or steady tendency towards a single regime. This changes between $1990$ and $1997$, where the posterior mean of $S_{\text{\textit{\c{C}}},t}$ is consistently above $0.5$, albeit with several high-frequency movements. In the aftermath of the global financial crisis, with short-term rates approaching the zero lower bond, $S_{\text{\textit{\c{C}}},t}$ switches abruptly into regime $0$. Conditional on the respective loadings in $\bm \Lambda_S$ being non-zero, this feature would directly relate to the observed narrowing spread of the yield curve and accompanying structural breaks in coefficients of the equation related to the curvature of the yield curve.
We observe several interesting features of the unobserved factors in $\bm{\tau}_t$. While $\tau_{L,t}$ is noisy and indicates substantial high-frequency movements, $\tau_{\text{\textit{\v{S}}},t}$ and $\tau_{\text{\textit{\c{C}}},t}$ are much smoother. The factor governing coefficients in the level-equation of the yield curve peaks early in the sample, followed by a decline between 1980 and 1990. After a brief increase and stabilization between 1990 and 2000, we see gradual declines until the global financial crisis starting in 2007. Since then, the factor shows upward trending movement, with several high-frequency troughs. By contrast, the unobserved factor related to $\mathcal{S}_t$ exhibits approximately linear trending behavior from the beginning of the sample until the early 2000s, where it plateaued. After 2010, a gradual but moderate decrease is visible. $\tau_{\text{\textit{\c{C}}},t}$ is comparable to $\tau_{L,t}$, albeit with several differences. While several peaks coincide, we also find adverse movements, for instance in the brief early 1980s recession and after 2000. Interestingly, high-frequency movements are muted when compared to $\tau_{L,t}$.
The preceding discussion of $\tilde{\bm z}_t$ must be considered in light of the rescaled loadings in $\tilde{\bm \Lambda} = \bm \Lambda \bm U^{-1}$. $\tilde{\bm \Lambda}$ translates the law of motion captured in $\tilde{\bm z}_t$ to the coefficients in $\bm \beta_t$ by acting either as amplifier or attenuator. Figure (ref) shows the posterior mean of the rescaled loadings $\bm \Lambda \bm U^{-1}$ and allows to assess which elements in $\tilde{\bm z}_t$ determine the time variation in the TVPs. We differentiate the loadings along two dimensions. Panel (a) shows the effect modifiers related to the VAR coefficients, while panel (b) depicts the block of $\tilde{\bm \Lambda}$ related to the covariances (stored in $\bm q_t$).
Assessing the loadings in $\bm \Lambda_r$ for the $L_t$-equation reveals that a major part of the coefficients loads strongly positive on $r_{\text{NFCI},t}$. In the case of $r_{\text{RF},t}$ and $r_{\text{REC},t}$ the patterns are more mixed. For $r_{\text{NFCI},t}$ we find a maximum loading of around $0.25$ for $\text{\textit{\v{S}}}_{t-1}$, $L_{t-2}$, $\text{\textit{\c{C}}}_{t-2}$ and $L_{t-3}$, while for $r_{\text{RF},t}$ and $r_{\text{REC},t}$ some coefficients load moderately positive (e.g., the first own lag $L_{t-1}$) and others moderately negative (e.g., loadings related to lags of the curvature factor). The loadings $\bm \Lambda_S$ and $\bm \Lambda_\tau$ related to estimated latent factors exhibit only modest relevance for defining the law of motion in coefficients related to $L_t$. In particular, the amplifiers for $\bm{\tau}_t$ are shrunk heavily towards zero, implying that low frequency movements in the respective coefficients are either already captured by other measures, or irrelevant for coefficients in the $L_t$-equation. The same is true for the indicators in $\bm{S}_t$ and $\bm{\tau}_t$ for the case of the $\text{\textit{\v{S}}}_t$-equation.
Turning to observed factors, the corresponding factor loadings are positive for lower-order lags while strongly negative for higher-order lags. By contrast, we find mostly negative loadings for observed factors for first-lags in the $\text{\textit{\c{C}}}_t$-equation. Loadings related to higher-order lags are mixed, and no clear patterns are observable. Interestingly, the curvature factor is the only equation where we detect substantial loadings on the latent factors, with most loadings showing a positive sign.
While unobserved factors appear to be less important for coefficients in the conditional mean of the model, they play an important role for the covariances, as indicated in panel (b). The covariance between $\text{\textit{\v{S}}}_t$ and $\text{\textit{\c{C}}}_t$, $w(\text{\textit{\v{S}}}_t,\text{\textit{\c{C}}}_t)$, shows pronounced loadings on all other effect modifiers than NFCI, where we observe a modest negative loading. A mixed pattern emerges for the $L_t$ and $\text{\textit{\c{C}}}_t$ covariance, $w(L_t,\text{\textit{\c{C}}}_t)$, with some positive and negative measures. The contemporaneous relationship $w(L_t,\text{\textit{\v{S}}}_t)$ between $L_t$ and $\text{\textit{\v{S}}}_t$, marks a particularly interesting case, with modestly negative loadings on NFCI and RF, while REC and the unobserved factor loadings are close to zero.
Summarizing, we find that observed factors often load strongly on the coefficients of all equations. While the Markov switching indicator appears to be important particularly for the $L_t$- and $\text{\textit{\c{C}}}_t$-equations, loadings are muted for the coefficients of $\text{\textit{\v{S}}}_t$. Another interesting aspect is that conditional on observed effect modifiers, the gradually evolving coefficients captured by $\bm{\tau}_t$ are mostly irrelevant for the $L_t$- and $\text{\textit{\v{S}}}_t$-equations, different to strong positive loadings in the case of the $\text{\textit{\c{C}}}_t$-equation and the covariances. Additional results showing the actual regression coefficients are provided in Appendix (ref).
To assess what our framework implies on the relations between the factors that determine the yield curve, we compute the low-frequency relationship between the $L_t$, $\text{\textit{\v{S}}}_t$ and $\text{\textit{\c{C}}}_t$. We choose this long-run correlation measure for two reasons. First, it allows to illustrate movements in time-varying coefficients (i.e., transmission channels) and changes in the error variances in a single indicator over time. Second, this measure compresses information of all coefficients in a structurally meaningful way since it isolates long-run trends and correlations from short-run fluctuations sargent2011two, kliem2016low. Additional results for the reduced form coefficients are provided in Appendix B.
To construct the measure, we transform the TVP-VAR($P$) in Eq. ((ref)) to its state-space TVP-VAR($1$) form. In what follows, the observation equation is given by $\bm y_t = \bm J \bm Y_t$, while the state equation is defined as $\bm Y_t = \bm B_t \bm Y_{t-1} + \bm \omega_t$ with $\bm \omega_t \sim \mathcal{N}(\bm 0, \bm \Omega_t)$. Here, $\bm J$ maps a $K$-dimensional vector $\bm Y_t = (\bm y_t', \dots, \bm y_{t-P+1}')'$ to $\bm y_t$, $\bm B_t$ collects the elements in $\bm \beta_t$ in the upper $M \times K$ block and defines identities otherwise. Similar, the $K \times K$-dimensional variance-covariance matrix $\bm \Omega_t$ collects elements in $\bm \Sigma_t$ in the upper-left $M \times M$ block and is zero otherwise. We follow sargent2011two and first calculate the spectral density $\bm \Phi_t(0)$ of $\bm y_t$ at a zero frequency which coincides with the unconditional variance-covariance matrix of $\bm y_t$:
Next, we transform the covariances of $\bm \Phi_t$ into a correlation measure for each period $t$ and each variable combination $i,j~(i \neq j \text{ and } i,j = 1, \dots, M)$:
The measure $\phi_{ij,t}$ describes the long-run relations between variable $i$ and $j$ at each point in time, and is displayed in Figure (ref) for the model specification TVP-NS-VAR with $\bm{z}_t=(\bm{r}_t^\prime, \bm{S}_t^\prime, \bm{\tau}_t^\prime)^\prime$, $R_{\tau j}=1$ and $P=3$. Note that the variables enter our model in differences. Hence, Figure (ref) depicts the low frequency relations of changes in the level, slope and curvature of the yield curve.
We observe several interesting periods characterized by structural breaks. First, the relationship between the level and slope of the yield curve was close to zero until the Volcker disinflation of the 1980s. After this period, we can identify an abrupt decrease to significantly negative values. The long-run correlation stays negative until the end of the sample, with minor low-frequency movements.
Second, for most of the sample the long-run coefficient between level and curvature of the yield curve is insignificant. Substantial structural breaks are detectable in the early 1980s. Again, this coincides with a shift in the US monetary policy regime. Between the two recessions in the early 1980s and in the recovery period afterwards, we note a strong positive relationship between the level and curvature of the yield curve. Afterwards, during the Great Moderation, there are mostly insignificant values. This ends during the Great Recession, where the relationship is estimated to be significantly negative. This finding may be linked to short-term interest rates approaching the zero lower bound rapidly. The trend reverses late in the sample, with the Federal Reserve conducting several subsequent rate hikes starting at the end of 2015. Turning to the relationship between the slope and curvature of the yield curve, we again find an insignificant relationship for most of the sample. Large breaks are observable in the period between the two 1980s recessions, but different to the relationship between the level and curvature of the yield curve, this is not visible in the subsequent recovery period. Interestingly, we estimate a significantly positive relationship in the period of the aforementioned rate hikes by the Federal Reserve starting in 2015.
This paper proposes methods for automatically selecting adequate state equations in TVP-VAR models in a data-driven fashion. The TVPs are assumed to depend on a set of observed and unobserved covariates, also known as effect modifiers. As unobserved covariates, we consider a set of low dimensional latent factors that follow a random walk, alongside Markov switching indicators that allow for abrupt structural breaks. Our model nests several alternatives commonly used in the literature on modeling macroeconomic and financial time series. To choose between state equations, we use a hierarchical Bayesian global-local shrinkage prior on the most flexible specification.
We apply our econometric framework to US yield curve data. Carrying out a thorough predictive exercise, we show that our techniques produce favorable point and density forecasts vis-\`{a}-vis a set of established benchmark models (which are nested variants of our proposed modeling approach). The performance is specific to the information set used in the underlying TVP-VAR and appears to be more pronounced for density forecasts. This exercise illustrates that our approach produces very competitive forecasts without increasing the risk of overfitting, while providing a framework to trace the sources of time-variation -- a key advantage compared to conventional TVP-VARs. This predictive exercise is complemented by a full-sample analysis of structural breaks in the relationship between the level, slope and curvature of the US yield curve. We detect several interesting patterns in abrupt and gradual time-variation patterns in long-run cross-variable relations. These changes appear to be specific to the monetary regime and the state of the business cycle.\\
Acknowledgments: The authors gratefully acknowledge financial support by the Jubil\"aums-fonds of the Oesterreichische Nationalbank (OeNB, project 18127) and by the Austrian Science Fund (FWF, project ZK 35).
{\setstretch{2} \addcontentsline{toc}{section}{References} }