# General Bayesian time-varying parameter VARs for predicting government bond
yields

Manfred M. Fischer*1 *Corresponding author: Manfred M. Fischer, Vienna
University of Economics and Business,
Welthandelsplatz 1, A-1020 Vienna, Austria. E-mail<EMAIL_ADDRESS>
1 Vienna University of Economics and Business, 2 University of Salzburg  Niko
Hauzenberger1, 2
Florian Huber2 and Michael Pfarrhofer2

Time-varying parameter (TVP) regressions commonly assume that time-variation
in the coefficients is determined by a simple stochastic process such as a
random walk. While such models are capable of capturing a wide range of
dynamic patterns, the true nature of time variation might stem from other
sources, or arise from different laws of motion. In this paper, we propose a
flexible TVP VAR that assumes the TVPs to depend on a panel of partially
latent covariates. The latent part of these covariates differ in their state
dynamics and thus capture smoothly evolving or abruptly changing coefficients.
To determine which of these covariates are important, and thus to decide on
the appropriate state evolution, we introduce Bayesian shrinkage priors to
perform model selection. As an empirical application, we forecast the US term
structure of interest rates and show that our approach performs well relative
to a set of competing models. We then show how the model can be used to
explain structural breaks in coefficients related to the US yield curve.

JEL: C11, C30, E37, E43
Keywords: Bayesian shrinkage, interest rate forecasting, latent effect
modifiers, MCMC sampling, time-varying parameter regression

##  1 Introduction

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 (see
Primiceri, 2005; Cogley and Sargent, 2005; Dangl and Halling, 2012). 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 (see, e.g., Sims and Zha, 2006; Koop et al., 2009;
Hauzenberger, 2020). 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 Hastie and Tibshirani (1993). 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 (see Sims and Zha, 2006), whereas the former is consistent
with financial time series such as bond yields or stock returns (see Dangl and
Halling, 2012; Huber et al., 2019).

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 Primiceri (2005) and Cogley and Sargent
(2005) or the reduced-rank model of Chan et al. (2020).

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 (see Carvalho et al., 2010) 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 Diebold and Li (2006) and Diebold et al.
(2008). 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 2 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 3 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.

##  2 Flexible Bayesian Inference in TVP-VARs

###  2.1 The TVP-VAR

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}^{\prime}_{t-1},\dots,\bm{y}^{\prime}_{t-P})^{\prime}$.
Then the basic TVP-VAR can be written as a linear multivariate regression
model:

$\bm{y}_{t}=(\bm{I}_{M}\otimes\bm{x}^{\prime}_{t})\bm{\beta}_{t}+\bm{\epsilon}_{t},$

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:

$\bm{\Sigma}_{t}=\bm{Q}_{t}\bm{H}_{t}\bm{Q}^{\prime}_{t}.$

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:

$h_{jt}=\mu_{j}+\psi_{j}(h_{jt-1}-\mu_{j})+\nu_{jt},\quad\nu_{jt}\sim\mathcal{N}(0,\varsigma^{2}_{j}).$

$\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 (Frühwirth-Schnatter and Wagner, 2010):

$\bm{y}_{t}=(\bm{I}_{M}\otimes\bm{x}^{\prime}_{t})(\bm{\beta}+\tilde{\bm{\beta}}_{t})+({\bm{Q}}+\tilde{\bm{Q}}_{t})\bm{\varepsilon}_{t},\quad\bm{\varepsilon}_{t}\sim\mathcal{N}(\bm{0}_{M},\bm{H}_{t}).$

$\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\times 1$-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}^{\prime},\tilde{\bm{q}}_{t}^{\prime})^{\prime}$
and the constant part
${\bm{\gamma}}=({\bm{\beta}}^{\prime},{\bm{q}}^{\prime})^{\prime}$. The
literature typically assumes that the transition distribution
$p(\tilde{\bm{\gamma}}_{t}|\tilde{\bm{\gamma}}_{t-1})$ is given by:

$\tilde{\bm{\gamma}}_{t}|\tilde{\bm{\gamma}}_{t-1}\sim\mathcal{N}(\tilde{\bm{\gamma}}_{t-1},\bm{V})\text{
with }\tilde{\bm{\gamma}}_{0}=\bm{0}_{N}.$

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.

###  2.2 A General Specification for the TVPs

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:

$\tilde{\bm{\gamma}}_{t}=\bm{\Lambda}\bm{z}_{t}+\bm{\eta}_{t}.$ (1)

$\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 (for some recent papers that follow this approach, see
Korobilis, 2019; Hauzenberger et al., 2019). 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. (1) 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 (see, e.g., Geweke and Zhou, 1996; Aguilar and West, 2000;
Stock and Watson, 2011).

###  2.3 Possible Choices for the Effect Modifiers

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:

  * •

Chan et al. (2020): We assume that $\bm{z}_{t}$ consists exclusively of a
sequence of $R=R_{\tau}$ latent factors $\bm{\tau}_{t}$, which follow a
multivariate random walk:

$\bm{\tau}_{t}=\bm{\tau}_{t-1}+\bm{\nu}_{t},\quad\bm{\nu}_{t}\sim\mathcal{N}(0,\bm{V}_{\tau}).$

$\bm{V}_{\tau}=\text{diag}(v_{1}^{2},\dots,v^{2}_{R_{\tau}})$ denotes a
diagonal variance-covariance matrix with $v_{j}^{2}$ being process innovation
variances that determine the smoothness of the elements in $\bm{\tau}_{t}$.
Note that setting $v_{j}^{2}$ close to zero effectively implies that
$\tau_{jt}$, the $j^{\text{th}}$ element in $\bm{\tau}_{t}$, is constant. This
model implies a factor structure in $\tilde{\bm{\gamma}}_{t}$ if $R_{\tau}\ll
N$.

  * •

Primiceri (2005): If $R=R_{\tau}=N$, the elements in $\bm{z}_{t}$ are random
walks, and $\bm{\Lambda}=\bm{I}_{N}$, we obtain a standard time-varying
parameter model. Assuming that the covariances are constant we obtain the
model put forth in Cogley and Sargent (2005).

  * •

Sims and Zha (2006): A Markov switching model can be obtained by setting
$\bm{z}_{t}=S_{t}$, with $S_{t}\in\\{0,1\\}$ denoting a binary indicator with
transition probabilities given by:

$p(S_{t}=i|S_{t-1}=j)=p_{ij}\quad\text{for }i,j=0,1,$

with $p_{ij}$ denoting the $(i,j)^{\text{th}}$ element of a $2\times
2$-dimensional transition probability matrix $\bm{P}$. Inclusion of this
random quantity allows to capture structural breaks in
$\tilde{\bm{\gamma}}_{t}$ that are common to all coefficients.

  * •

Caggiano et al. (2017): Assuming that $\bm{z}_{t}$ is exclusively composed of
observed quantities we obtain a regression model with interaction terms.

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:

$\bm{z}_{t}=(\bm{r}^{\prime}_{t},\bm{S}^{\prime}_{t},\bm{\tau}^{\prime}_{t})^{\prime}.$

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}^{\prime},\dots,\bm{\tau}_{Mt}^{\prime})^{\prime}$
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})^{\prime}$,
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, (see Koop et al., 2009; Maheu and Song, 2018). 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.

###  2.4 The Prior Setup

The discussion in Sub-section 2.1 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:

$\bm{\Lambda}=[\bm{\Lambda}_{r}\quad\bm{\Lambda}_{S}\quad\bm{\Lambda}_{\tau}],$

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 (Carvalho et
al., 2010). 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:

$\lambda_{ji,\ell}|\kappa_{ji,\ell},\delta_{ji}\sim\mathcal{N}(0,\kappa_{ji,\ell}^{2}\delta_{ji}^{2}),\quad\kappa_{ji,\ell}\sim\mathcal{C}^{+}(0,1),\quad\delta_{ji}\sim\mathcal{C}^{+}(0,1).$

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:

$\gamma_{ji}|\xi_{ji},\zeta_{j}\sim\mathcal{N}(0,\xi_{ji}^{2}\zeta_{j}^{2}),\quad\xi_{ji}\sim\mathcal{C}^{+}(0,1),\quad\zeta_{j}\sim\mathcal{C}^{+}(0,1).$

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. (1). 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:

$\sqrt{\omega_{ji}}|\varpi_{ji},\vartheta_{j}\sim\mathcal{N}(0,\varpi_{ji}^{2}\vartheta_{j}^{2}),\quad\varpi_{ji}\sim\mathcal{C}^{+}(0,1),\quad\vartheta_{j}\sim\mathcal{C}^{+}(0,1).$

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 Frühwirth-Schnatter and Wagner (2010). 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}$
(see Chan et al., 2020).

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
Kastner and Frühwirth-Schnatter (2014). 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:

$p_{j,ii}\sim\mathcal{B}(e_{i0},e_{i1}),\quad\text{for }i=0,1\text{,
}j=1,\dots,M,$

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.

###  2.5 Full Conditional Posterior Simulation

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 Carriero et al. (2019).

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:

$y_{1t}=\bm{x}^{\prime}_{t}(\bm{\beta}_{1}+\tilde{\bm{\beta}}_{1t})+\varepsilon_{1t},$

and the $j^{\text{th}}$ equation $(j>1)$:

$\displaystyle
y_{jt}=\bm{x}^{\prime}_{t}(\bm{\beta}_{j}+\tilde{\bm{\beta}}_{jt})+\bm{u}^{\prime}_{jt}(\bm{q}_{j}+\tilde{\bm{q}}_{jt})+\varepsilon_{jt}.$
(2)

${\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}^{\prime}_{1t},\dots,\bm{\beta}^{\prime}_{Mt})^{\prime}$
and $\bm{u}_{jt}=(\varepsilon_{1t},\dots,\varepsilon_{j-1,t})^{\prime}$. 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 (2) can be simplified to yield:

$y_{jt}=\bm{m}^{\prime}_{jt}(\bm{\gamma}_{j}+\tilde{\bm{\gamma}}_{jt})+\varepsilon_{jt},$
(3)

where $\bm{m}_{jt}=(\bm{x}^{\prime}_{t},\bm{u}_{jt}^{\prime})^{\prime}$,
$\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. (3). For
notational simplicity, we assume that all elements in $\bm{z}_{t}$ are latent.
In light of the discussion in Sub-section 2.3, this implies that
$\bm{z}_{jt}=(\bm{\tau}^{\prime}_{jt},~{}S_{jt})^{\prime}$ 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. (3) as:

$\tilde{y}_{jt}=\bm{m}^{\prime}_{jt}\bm{\Lambda}_{j}\bm{z}_{jt}+\bm{m}^{\prime}_{jt}\bm{\eta}_{jt}+\varepsilon_{jt},$
(4)

with $\tilde{y}_{jt}=y_{jt}-\bm{m}^{\prime}_{jt}\bm{\gamma}_{j}$. Defining
$\tilde{\bm{m}}^{\prime}_{jt}=\bm{m}^{\prime}_{jt}\bm{\Lambda}_{j}$ and
$\hat{\varepsilon}_{jt}=\bm{m}^{\prime}_{jt}\bm{\eta}_{jt}+\varepsilon_{jt}$
allows us to cast Eq. (4) as a simple linear regression model:

$\tilde{y}_{jt}=\tilde{\bm{m}}^{\prime}_{jt}\bm{z}_{jt}+\hat{\varepsilon}_{jt},\quad\hat{\varepsilon}_{jt}\sim\mathcal{N}(0,\bm{m}^{\prime}_{jt}\text{diag}(\bm{\omega}_{j})\bm{m}_{jt}+e^{h_{jt}}).$
(5)

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
(for a detailed discussion of this issue, see Gerlach et al., 2000; Giordani
and Kohn, 2008).

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 Carter and Kohn
(1994) and Frühwirth-Schnatter (1994). In case of the latent Markov switching
factors in $S_{jt}$, we use the algorithm outlined in Kim and Nelson (1999).
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 (see also the discussion in Chan et al., 2020).111It 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. (1) 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 A.

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:

$y_{jt}=\hat{\bm{m}}_{jt}^{\prime}\hat{\bm{\gamma}}_{j}+\hat{\varepsilon}_{jt},$
(6)

where
$\hat{\bm{m}}_{jt}=(\bm{m}_{jt}^{\prime},(\bm{z}_{jt}\otimes\bm{m}_{jt})^{\prime})^{\prime}$
is an $Rv_{j}^{2}$-dimensional vector of covariates, and
$\hat{\bm{\gamma}}_{j}=(\bm{\gamma}_{j}^{\prime},\text{vec}(\bm{\Lambda}_{j})^{\prime})^{\prime}$
denoting an $Rv_{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 A.

##  3 Forecasting US Government Bond Yields

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.

###  3.1 Data and Design of the Forecasting Exercise

Our aim is to predict monthly zero-coupon yields of US treasuries at different
yearly maturities. The data is described in detail in Gürkaynak et al.
(2007).222Available online at 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.333Note 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 (Geweke and Amisano, 2010).

For analyzing the role of observed factors governing joint dynamics of the
TVPs, we rely on three exogenous variables:

  1. (a)

a binary recession indicator (labeled REC) for the US, dated on a monthly
basis by the Business Cycle Dating Committee of the National Bureau of
Economic Research,

  2. (b)

the National Financial Conditions Index (NFCI), maintained by the Federal
Reserve Bank of Chicago, downloaded from the FRED database of the Federal
Reserve Bank of St. Louis (available online at fred.stlouisfed.org), and

  3. (c)

the Risk-Free (RF) interest rate from the Fama-French Portfolios and Factors
database, provided on the web page of Kenneth R. French as another important
early warning indicator.444Available online at
mba.tuck.dartmouth.edu/pages/faculty/ken.french. Thus,
$\bm{r}_{t}=({r}_{\text{NGCI},t-1},{r}_{\text{REC},t-1},{r}_{\text{RF},t-1})^{\prime}$
and $R_{r}=3$. Exogenous variables enter the model as first order lags.
Higher-order forecasts involving exogenous variables are based on random walk
predictions of these quantities.

We use these three effect modifiers for simple reasons. First, there is strong
evidence that yield curve dynamics differ across business cycle phases (see,
e.g., Hevia et al., 2015). 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 (see, e.g., Adrian et al., 2019; Carriero et al., 2020; Adams et al.,
2021).

###  3.2 Competing Model Specifications

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 Diebold and Li (2006) that imposes a factor
structure on the yields:

$i_{t}(\theta)=L_{t}+\left(\frac{1-\exp(-\theta\alpha)}{\theta\alpha}\right)\text{{\v{S}}}_{t}+\left(\frac{1-\exp(-\theta\alpha)}{\theta\alpha}-\exp(-\theta\alpha)\right)\text{{\c{C}}}_{t}.$

where $i_{t}(\theta)$ denotes the yield at maturity $\theta$ at time $t$,
$L_{t}$ is a factor that controls the level, $\text{{\v{S}}}_{t}$ determines
the slope, and $\text{{\c{C}}}_{t}$ represents the curvature factor of the
yields. The parameter $\alpha$ governs the exponential decay rate. To maximize
the loading on $\text{{\c{C}}}_{t}$ we set $\alpha=0.7308\ (12\times
0.0609)$.555See Diebold and Li (2006) for a discussion of this specific
choice. In what follows, we use the latent factors
$L_{t},\text{{\v{S}}}_{t}\text{ and }\text{{\c{C}}}_{t}$ as endogenous
variables in the VAR specifications by defining
$\bm{y}_{t}=(L_{t},\text{{\v{S}}}_{t},\text{{\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 LABEL:tab:LPS): 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.

###  3.3 Results

Table LABEL:tab:LPS 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.

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 LABEL:tab:LPS 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 LABEL:tab:LPS. 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.

###  3.4 Determinants of Time-Variation in the Coefficients

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.666In 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 1: Evolution of the normalized effect modifiers
$\tilde{\bm{z}}_{t}=\bm{U}\bm{z}_{t}$ over time.

(a) Observed factors

(b) Equation-specific latent Markov switching factors

(c) Equation-specific latent random walk factors

Notes: Results are based on a TVP-NS-VAR model specification
$\bm{z}_{t}=(\bm{r}^{\prime}_{t},\bm{S}_{t}^{\prime},\bm{\tau}_{t}^{\prime})^{\prime}$
and $\delta=R_{\tau}/M=1$ and $P=3$. Panel (a) shows the normalized observed
factors $r_{\text{NFCI},t},r_{\text{REC},t}$ and $r_{\text{RF},t}$ (collected
in $\bm{r}_{t}$), panel (b) the posterior mean of the latent Markov switching
factors, $S_{L,t},S_{\text{{\v{S}}},t}$ and $S_{\text{{\c{C}}},t}$ (collected
in $\bm{S}_{t}$) and panel (c) the posterior median of the latent random walk
factors, $\tau_{L},\tau_{\text{{\v{S}}}}$ and $\tau_{\text{{\c{c}}}}$
(collected in $\bm{\tau}_{t}$). Note that $S_{jt}$ and $\tau_{jt}$ are
equation-specific latent quantities with $M(=3)$ endogenous variables,
$j\in\\{L,\text{{\v{S}}},\text{{\c{C}}}\\}$. $S_{L,t}$ and $\tau_{L,t}$
correspond to the first equation, $S_{\text{{\v{S}}},t}$ and
$\tau_{\text{{\v{S}}},t}$ to the second, while $S_{\text{{\c{C}}},t}$ and
$\tau_{\text{{\c{C}}},t}$ to the third. Results are based on the TVP-NS-VAR
model variant with $\delta=R_{\tau}/M=1$ and using 15,000 MCMC draws. The gray
shaded vertical bars represent recessions dated by the NBER Business Cycle
Dating Committee. Sample period: 1973:01 to 2019:12. Vertical axis: normalized
values. Front axis: months.

Figure 1 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{{\v{S}}},\text{{\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{{\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{{\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{{\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{{\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{{\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{{\v{S}}},t}$ and $\tau_{\text{{\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{{\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}$.

Figure 2: Heat maps for rescaled loadings in $\tilde{\bm{\Lambda}}=\bm{\Lambda
U}^{-1}$.

(a) Coefficients

${L}_{t}$-equation

$\text{{\v{S}}}_{t}$-equation

$\text{{\c{C}}}_{t}$-equation

Space keeper

(b) Covariances

Notes: $\tilde{\bm{\Lambda}}$ translates the law of motion captured in
$\tilde{\bm{z}}_{t}$ to (a) the VAR coefficients in $\bm{\beta}_{t}$, and (b)
the covariances stored in $\bm{q}_{t}$, based on a TVP-NS-VAR model
specification with
$\bm{z}_{t}=(\bm{r}^{\prime}_{t},\bm{S}_{t}^{\prime},\bm{\tau}_{t}^{\prime})^{\prime}$,
$\delta=R_{\tau}/M=1$ and $P=3$. $\iota_{t}$ denotes the time-varying
intercept, and $L_{t-p},\text{{\v{S}}}_{t-p},\text{{\c{C}}}_{t-p}$ for
$p=1,2,3$ the lagged variables $L_{t},\text{{\v{S}}}_{t}$ and
$\text{{\c{C}}}_{t}$, respectively. $\bm{\lambda}_{r}$ relates to elements of
$\bm{\tilde{\Lambda}}$ associated with $r_{\text{NFCI},t},r_{\text{REC},t}$
and $r_{\text{RF},t}$ and the respective equations, $\bm{\lambda}_{S}$ denotes
loadings related to a single Markov switching factor collected in
$\bm{S}_{t}$, and $\bm{\lambda}_{\tau}$ loadings corresponding to latent
random walk factors in $\bm{\tau}_{t}$. Note that $S_{jt}$ and $\tau_{jt}$ for
$j\in\\{L,\text{{\v{S}}},\text{{\c{C}}}\\}$ are equation-specific quantities,
while $\bm{r}_{t}$ stays fixed across equations.
$w(\text{{\v{S}}}_{t},\text{{\c{C}}}_{t})$ denotes the contemporaneous
relation between the $\text{{\v{S}}}_{t}$\- and
$\text{{\c{C}}}_{t}$-equations, $w(L_{t},\text{{\c{C}}}_{t})$ and
$w(L_{t},\text{{\v{S}}}_{t})$ are defined analogously. Sample period: 1973:01
to 2019:12.

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 2 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{{\v{S}}}_{t-1}$, $L_{t-2}$, $\text{{\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{{\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{{\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{{\v{S}}}_{t}$ and $\text{{\c{C}}}_{t}$,
$w(\text{{\v{S}}}_{t},\text{{\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{{\c{C}}}_{t}$ covariance,
$w(L_{t},\text{{\c{C}}}_{t})$, with some positive and negative measures. The
contemporaneous relationship $w(L_{t},\text{{\v{S}}}_{t})$ between $L_{t}$ and
$\text{{\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{{\c{C}}}_{t}$-equations, loadings are muted for the coefficients of
$\text{{\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{{\v{S}}}_{t}$-equations, different to strong positive loadings in the
case of the $\text{{\c{C}}}_{t}$-equation and the covariances. Additional
results showing the actual regression coefficients are provided in Appendix B.

###  3.5 Low Frequency Relations between the Nelson-Siegel Factors

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{{\v{S}}}_{t}$ and $\text{{\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 (Sargent and Surico, 2011; Kliem et al., 2016). Additional
results for the reduced form coefficients are provided in Appendix B.

To construct the measure, we transform the TVP-VAR($P$) in Eq. (1) 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}^{\prime},\dots,\bm{y}_{t-P+1}^{\prime})^{\prime}$ 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 Sargent and Surico (2011) 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}$:

$\bm{\Phi}_{t}(0)=\bm{J}\left(\bm{I}-\bm{B}_{t}\right)\bm{\Omega}_{t}(\bm{I}-\bm{B}^{\prime}_{t})^{-1}\bm{J}^{\prime},\quad\text{for
}t=1,\dots,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)$:

$\phi_{ij,t}=\frac{\Phi_{ij,t}(0)}{\Phi_{jj,t}(0)}.$

Figure 3: Posterior median and the 68 percent credible set for pairwise long-
run correlations.

(a) $\phi_{L\text{{\v{S}}},t}$

(b) $\phi_{L\text{{\c{C}}},t}$

(c) $\phi_{\text{{\v{S}}}\text{{\c{C}}},t}$

Notes: Panel (a) $\phi_{L\text{{\v{S}}},t}$, (b) $\phi_{L\text{{\c{C}}},t}$
and (c) $\phi_{\text{{\v{S}}}\text{{\c{C}}},t}$, based on a TVP-NS-VAR with
$\bm{z}_{t}=(\bm{r}^{\prime}_{t},\bm{S}_{t}^{\prime},\bm{\tau}_{t}^{\prime})^{\prime}$,
$\delta=R\tau/M=1$ and $P=3$. The colored solid lines denote the posterior
medians, the black dashed line the zero line, and the colored shaded band the
68 percent posterior coverage interval, while the gray vertical bars represent
recessions dated by the NBER Business Cycle Dating Committee. $\phi_{ij,t}$ is
the long-run correlation for the variables $i$ and $j$ at time $t$, for
$i,j\in\\{L,\text{{\v{S}}},\text{{\c{C}}}\\}$. Results are based on 15,000
MCMC draws. Sample period: 1973:01 to 2019:12. Vertical axis: correlation
measurements. Front axis: months.

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 3 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 3 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.

##  4 Conclusions

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-à-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äums-fonds of the Oesterreichische Nationalbank (OeNB, project 18127) and
by the Austrian Science Fund (FWF, project ZK 35).

## References

  * Adams et al. (2021) Patrick A Adams, Tobias Adrian, Nina Boyarchenko, and Domenico Giannone.  Forecasting macroeconomic risks.  _International Journal of Forecasting_ , in press, published online 6 Feb, 2021.
  * Adrian et al. (2019) Tobias Adrian, Nina Boyarchenko, and Domenico Giannone.  Vulnerable growth.  _American Economic Review_ , 109(4):1263–89, 2019\.
  * Aguilar and West (2000) Omar Aguilar and Mike West.  Bayesian dynamic factor models and portfolio allocation.  _Journal of Business & Economic Statistics_, 18(3):338–357, 2000.
  * Caggiano et al. (2017) Giovanni Caggiano, Efrem Castelnuovo, and Giovanni Pellegrino.  Estimating the real effects of uncertainty shocks at the zero lower bound.  _European Economic Review_ , 100:257–272, 2017.
  * Carriero et al. (2019) Andrea Carriero, Todd Clark, and Massimiliano Marcellino.  Large vector autoregressions with stochastic volatility and flexible priors.  _Journal of Econometrics_ , 212(1):137–154, 2019\.
  * Carriero et al. (2020) Andrea Carriero, Todd E Clark, and Massimiliano Giuseppe Marcellino.  Capturing macroeconomic tail risks with Bayesian vector autoregressions.  Working Paper 202002R, Federal Reserve Bank of Cleveland, 2020.
  * Carter and Kohn (1994) Chris K Carter and Robert Kohn.  On Gibbs sampling for state space models.  _Biometrika_ , 81(3):541–553, 1994.
  * Carvalho et al. (2010) Carlos M Carvalho, Nicholas G Polson, and James G Scott.  The horseshoe estimator for sparse signals.  _Biometrika_ , 97(2):465–480, 2010.
  * Chan et al. (2020) Joshua CC Chan, Eric Eisenstat, and Rodney Strachan.  Reducing the state space dimension in a large TVP-VAR.  _Journal of Econometrics_ , 218(1):105–118, 2020\.
  * Cogley and Sargent (2005) Timothy Cogley and Thomas J. Sargent.  Drifts and volatilities: Monetary policies and outcomes in the post WWII US.  _Review of Economic Dynamics_ , 8(2):262 – 302, 2005.
  * Dangl and Halling (2012) Thomas Dangl and Michael Halling.  Predictive regressions with time-varying coefficients.  _Journal of Financial Economics_ , 106:157–181, 2012.
  * Diebold and Li (2006) Francis X Diebold and Canlin Li.  Forecasting the term structure of government bond yields.  _Journal of Econometrics_ , 130(2):337–364, 2006\.
  * Diebold et al. (2008) Francis X Diebold, Canlin Li, and Vivian Z Yue.  Global yield curve dynamics and interactions: A dynamic Nelson–Siegel approach.  _Journal of Econometrics_ , 146(2):351–363, 2008\.
  * Frühwirth-Schnatter (1994) Sylvia Frühwirth-Schnatter.  Data augmentation and dynamic linear models.  _Journal of Time Series Analysis_ , 15(2):183–202, 1994.
  * Frühwirth-Schnatter and Wagner (2010) Sylvia Frühwirth-Schnatter and Helga Wagner.  Stochastic model specification search for Gaussian and partial non-Gaussian state space models.  _Journal of Econometrics_ , 154(1):85–100, 2010\.
  * Gerlach et al. (2000) Richard Gerlach, Chris Carter, and Robert Kohn.  Efficient Bayesian inference for dynamic mixture models.  _Journal of the American Statistical Association_ , 95(451):819–828, 2000.
  * Geweke and Amisano (2010) John Geweke and Gianni Amisano.  Comparing and evaluating Bayesian predictive distributions of asset returns.  _International Journal of Forecasting_ , 26(2):216–230, 2010.
  * Geweke and Zhou (1996) John Geweke and Guofu Zhou.  Measuring the pricing error of the arbitrage pricing theory.  _Review of Financial Studies_ , 9(2):557–587, 1996.
  * Giordani and Kohn (2008) Paolo Giordani and Robert Kohn.  Efficient Bayesian inference for multiple change-point and mixture innovation models.  _Journal of Business & Economic Statistics_, 26(1):66–77, 2008.
  * Gürkaynak et al. (2007) Refet S Gürkaynak, Brian Sack, and Jonathan H Wright.  The US treasury yield curve: 1961 to the present.  _Journal of Monetary Economics_ , 54(8):2291–2304, 2007.
  * Hastie and Tibshirani (1993) Trevor Hastie and Robert Tibshirani.  Varying-coefficient models.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 55(4):757–779, 1993.
  * Hauzenberger (2020) Niko Hauzenberger.  Flexible mixture priors for time-varying parameter models.  _arXiv preprint_ , 2006.10088, 2020.
  * Hauzenberger et al. (2019) Niko Hauzenberger, Florian Huber, Gary Koop, and Luca Onorante.  Fast and flexible Bayesian inference in time-varying parameter regression models.  _arXiv preprint_ , 1910.10779, 2019.
  * Hevia et al. (2015) Constantino Hevia, Martin Gonzalez-Rozada, Martin Sola, and Fabio Spagnolo.  Estimating and forecasting the yield curve using a Markov switching dynamic Nelson and Siegel model.  _Journal of Applied Econometrics_ , 30(6):987–1009, 2015.
  * Huber et al. (2019) Florian Huber, Gregor Kastner, and Martin Feldkircher.  Should I stay or should I go? A latent threshold approach to large-scale mixture innovation models.  _Journal of Applied Econometrics_ , 34(5):621–640, 2019.
  * Kastner and Frühwirth-Schnatter (2014) Gregor Kastner and Sylvia Frühwirth-Schnatter.  Ancillarity-sufficiency interweaving strategy (ASIS) for boosting MCMC estimation of stochastic volatility models.  _Computational Statistics & Data Analysis_, 76:408–423, 2014.
  * Kim and Nelson (1999) Chang-Jin Kim and Charles R. Nelson.  _State-space models with regime switching: Classical and Gibbs-sampling approaches with applications_.  MIT Press, Cambridge and London, 1999.
  * Kliem et al. (2016) Martin Kliem, Alexander Kriwoluzky, and Samad Sarferaz.  On the low-frequency relationship between public deficits and inflation.  _Journal of Applied Econometrics_ , 31(3):566–583, 2016.
  * Koop et al. (2009) Gary Koop, Roberto Leon-Gonzalez, and Rodney W Strachan.  On the evolution of the monetary policy transmission mechanism.  _Journal of Economic Dynamics and Control_ , 33(4):997–1017, 2009.
  * Korobilis (2019) Dimitris Korobilis.  High-dimensional macroeconomic forecasting using message passing algorithms.  _Journal of Business & Economic Statistics_, in press, published online 22 Nov, 2019.
  * Maheu and Song (2018) John M Maheu and Yong Song.  An efficient Bayesian approach to multiple structural change in multivariate time series.  _Journal of Applied Econometrics_ , 33(2):251–270, 2018.
  * Makalic and Schmidt (2015) Enes Makalic and Daniel F Schmidt.  A simple sampler for the horseshoe estimator.  _IEEE Signal Processing Letters_ , 23(1):179–182, 2015.
  * Primiceri (2005) Giorgio Primiceri.  Time varying structural autoregressions and monetary policy.  _Review of Economic Studies_ , 72(3):821–852, 2005.
  * Sargent and Surico (2011) Thomas J Sargent and Paolo Surico.  Two illustrations of the quantity theory of money: Breakdowns and revivals.  _American Economic Review_ , 101(1):109–28, 2011\.
  * Sims and Zha (2006) Christopher A Sims and Tao Zha.  Were there regime switches in US monetary policy?  _American Economic Review_ , 96(1):54–81, 2006\.
  * Stock and Watson (2011) James H Stock and Mark Watson.  Dynamic factor models.  In Michael P. Clements and David F. Hendry, editors, _The Oxford Handbook of Economic Forecasting_ , pages 35–60. Oxford University Press, 2011\.

Appendices

##  Appendix A Technical appendix

###  A.1 Sampling the state innovation variances

For sampling the state innovation variances based on Eq. (1), we let
${\eta}_{ji,t}$ denote the shock to the $i^{\text{th}}$ coefficient in
$\tilde{\bm{\gamma}}_{t}$ with respect to the $j^{\text{th}}$ equation. The
posterior of the state innovation variances is a generalized inverse Gaussian
(GIG) distribution:777The generalized inverse Gaussian distribution is
specified such that its density function is proportional to
$f(x)=x^{\lambda-1}\exp(-(\chi/x+\psi x)/2)$ for a random variable
$x\sim\mathcal{GIG}(\lambda,\chi,\psi)$.

$\omega_{ji}|\varpi_{ji},\vartheta_{j}\sim\mathcal{GIG}\left(\frac{1-T}{2},\sum_{t=1}^{T}\eta_{ji,t}^{2},\varpi_{ji}\vartheta_{j}\right).$

###  A.2 Posterior for the horseshoe prior

Our specification of the horseshoe prior on $\bm{\Lambda}_{j}$,
$\bm{\gamma}_{j}$ and the square root of $\bm{\omega}_{j}$ in Section 2.4 for
a generic parameter $b_{i}$ for $i=1,\ldots,K$ is:

$b_{i}|c_{i},d\sim\mathcal{N}(0,c_{i}^{2}d^{2}),\quad
c_{i}\sim\mathcal{C}^{+}(0,1),\quad d\sim\mathcal{C}^{+}(0,1).$

We rely on this prior in its auxiliary representation as in Makalic and
Schmidt (2015) for efficient sampling of the local ($c_{i}$) and global ($d$)
shrinkage parameters:

$c_{i}^{2}|e_{i}\sim\mathcal{G}^{-1}(1/2,1/e_{i}),\quad
d^{2}|f\sim\mathcal{G}^{-1}(1/2,1/f),\quad
e_{i}\sim\mathcal{G}^{-1}(1/2,1),\quad f\sim\mathcal{G}^{-1}(1/2,1).$

Here, $\mathcal{G}^{-1}$ denotes the inverse Gamma distribution. This setup
yields the following conditional posterior distributions:

$\displaystyle c_{i}^{2}|b_{i},d,e_{i}$
$\displaystyle\sim\mathcal{G}^{-1}\left(1,\frac{1}{e_{i}}+\frac{b_{i}^{2}}{2d^{2}}\right),\quad
d^{2}|b_{i},c_{i},f\sim\mathcal{G}^{-1}\left(\frac{K+1}{2},\frac{1}{f}+\sum_{i=1}^{K}\frac{b_{i}^{2}}{2c_{i}^{2}}\right),$
$\displaystyle e_{i}|c_{i}$
$\displaystyle\sim\mathcal{G}^{-1}\left(1,1+c_{i}^{-2}\right),\quad
f|d\sim\mathcal{G}^{-1}\left(1,1+d^{-2}\right).$

##  Appendix B Further empirical results

This Appendix contains additional results for the reduced form coefficients.
While Sub-section 3.4 provides posterior estimates of the lower dimensional
effect modifiers, Figs. B.1–B.3 display the coefficients obtained by
multiplying the factor loadings with the observed/latent factors based on the
relationship established in Eq. (1).

Figure B.1: Posterior median of the coefficients associated with the own lags
of $L_{t}$, $\text{{\v{S}}}_{t}$ and $\text{{\c{C}}}_{t}$.

(a) Coefficients of first own lag

(b) Coefficients of second own lag

(c) Coefficients of third own lag

Notes: Panels (a), (b) and (c) show the dynamic evolution of the coefficients
related to the variables’ own lags $p\in\\{1,2,3\\}$ of the respective
equation for $L_{t}$, $\text{{\v{S}}}_{t}$ and $\text{{\c{C}}}_{t}$. Results
are based on the TVP-NS-VAR model variant with $\delta=R_{\tau}/M=1$ and using
15,000 MCMC draws. The black dashed line denotes the zero line, while the gray
shaded vertical bars represent recessions dated by the NBER Business Cycle
Dating Committee. Sample period 1973:01 to 2019:12. Vertical axis: posterior
median estimate. Front axis: months. Figure B.2: Posterior median of the
coefficients associated with cross-variable lags of $L_{t}$,
$\text{{\v{S}}}_{t}$ and $\text{{\c{C}}}_{t}$.

(a) Coefficients of first other lags

(b) Coefficients of second other lags

(c) Coefficients of third other lags

Notes: Panels (a), (b) and (c) show the dynamic evolution of the coefficients
related to cross-variable lags $p\in\\{1,2,3\\}$ of the respective equation
for $L_{t}$, $\text{{\v{S}}}_{t}$ and $\text{{\c{C}}}_{t}$. Results are based
on the TVP-NS-VAR model variant with $\delta=R_{\tau}/M=1$ and using 15,000
MCMC draws. The black dashed line denotes the zero line, while the gray shaded
vertical bars represent recessions dated by the NBER Business Cycle Dating
Committee. Sample period 1973:01 to 2019:12. Vertical axis: posterior median
estimate. Front axis: months.

Figure B.3: Posterior median of the contemporaneous relationships between
$L_{t}$, $\text{{\v{S}}}_{t}$ and $\text{{\c{C}}}_{t}$.

Notes: $w(\text{{\v{S}}}_{t},\text{{\c{C}}}_{t})$ denotes the contemporaneous
relation between the $\text{{\v{S}}}_{t}$\- and
$\text{{\c{C}}}_{t}$-equations, $w(L_{t},\text{{\c{C}}}_{t})$ and
$w(L_{t},\text{{\v{S}}}_{t})$ are defined analogously. Results are based on
the TVP-NS-VAR model variant with $\delta=R_{\tau}/M=1$ and using 15,000 MCMC
draws. The black dashed line denotes the zero line, while the gray shaded
vertical bars represent recessions dated by the NBER Business Cycle Dating
Committee. Sample period 1973:01 to 2019:12. Vertical axis: posterior median
estimate. Front axis: months.