Dynamic Shrinkage Priors for Large Time-Varying Parameter Regressions Using Scalable Markov Chain Monte Carlo Methods

2.1 A TVP Regression

The static representation of a TVP regression model involving a T-dimensional dependent variable, y , and a T × K-dimensional matrix of predictors, X is:

where α is a K-dimensional vector of time-invariant coefficients, β _t is a K-dimensional vector of time-varying coefficients and L  = diag(σ ₁, …, σ _T ) with σ _t denoting time-varying error volatilities. The TVP part of this model arises through the Wβ term. W is a T × k(=TK) matrix given by:

with x _t denoting a K-dimensional sub-vector of X . Equation (1) is simply a regression which leads to the terminology static representation. But it is a regression with an enormous number of explanatory variables.

Note that (2) implies that the TVPs are mean zero and uncorrelated over time. However, extensions to other forms can be trivially done through a re-definition of W . For instance, if we are interested in random walk-type behavior in the TVPs, we can set

This specification implies that β can be interpreted as the changes in the parameters and multiplication with W yields the cumulative sum over β . In our empirical exercise, we consider both of these specifications for W and refer to the former as the flexible (FLEX) and the latter as the random walk (RW) specification.

The existing literature using Bayesian shrinkage techniques typically uses MCMC methods. Exact MCMC sampling, however, quickly becomes computationally cumbersome since k is extremely large even for moderate values of T and K.

Various solutions to this have been proposed in the literature. The standard solution is simply not to work with the static representation, but instead make some parametric assumption about how the TVPs evolve (e.g. assume they follow random walks or a Markov switching process). Unless K is extremely large, exact MCMC methods are feasible. However, with macroeconomic data it is common to find strong evidence of changes in the conditional variance of a series, but much less evidence in favor of change in the conditional mean of a series, (see, e.g. Clark 2011). When K is large, it is plausible to assume that only some of the predictors have time-varying coefficients and, even for these, coefficient change may only rarely happen. Common conventional approaches are not suited for data sets which exhibit such sparsity in the TVPs. If changes in the conditional mean of the parameters happen only rarely then a random walk assumption, which assumes change is continually happening, is not appropriate. If changes in the conditional mean only occur for a small sub-set of the K variables (or occur at different times for different variables), then a Markov switching model which assumes all coefficients change at the same time is not appropriate. These considerations motivate our use of the static representation and the development of a dynamic shrinkage prior suited for the case of TVP sparsity.

The literature has proposed a few ways of overcoming the computational hurdle that arises if the static representation is used. Korobilis (2021) uses message passing techniques to estimate large TVP regressions and shows that these large models outperform a range of competing models. Similarly, Huber, Koop, and Pfarrhofer (2020) approximate the TVPs using message passing techniques based on a rotated model representation and sample from the full conditional posterior of α using MCMC methods. Both approaches have the drawback that the quality of the approximation inherent in the use of message passing techniques might be questionable. In another recent paper, Hauzenberger et al. (2022) propose using the singular value decomposition of W in combination with a conjugate shrinkage prior on β to ensure computational efficiency. However, this method has the potential drawback that conjugate priors might be too restrictive for discriminating signals and noise in high dimensional models.

In this paper, we develop another approach which should work particularly well when β is extremely sparse. This is the scalable MCMC method, based on posterior perturbations, of Johndrow, Orenstein, and Bhattacharya (2017).

2.2 Extension to a TVP-VAR

Before discussing the scalable MCMC algorithm, we note that methods developed for the TVP regression can also be used for the TVP-VAR if it is written in equation-by-equation form (see, for instance, Carriero, Clark, and Marcellino 2019; Huber Koop, and Onorante 2021). In particular, we can use the following structural representation of the TVP-VAR:

with y _t being an M-dimensional vector of endogenous variables, c _t denoting an M-dimensional vector of intercepts, A _pt , for p = 1, …, P, denoting an M × M-dimensional time-varying coefficient matrix that may be stacked in a matrix A _t  = ( A _1t , …, A _Pt ). Furthermore, ϵ _t is an M-dimensional vector of errors and Σ t = diag σ 1 t 2 , … , σ M t 2 refers to its diagonal time-varying covariance matrix. Finally, A _0t defines contemporaneous relationships between the elements of y _t and is lower-triangular with zeros on the diagonal.

The ith (i = 2, …, M) equation of y _t can be written as a standard TVP regression model:

Here, x _it is a K _i (= MP + i)-dimensional vector of covariates with x i t = 1 , { y j t } j = 1 i − 1 , y t − 1 ′ , … , y t − P ′ ′ , γ i t = ( α i + β i t ) = c i t , { a i j , 0 t } j = 1 i − 1 , A i • , t ′ denotes a K _i -dimensional vector of time-varying coefficients, with c _it referring to the ith element in c _t , a _ij,0t denoting the (i,j)th element of A _0t and A _i•,t referring to the ith row of A _t . For i = 1, x 1 t = 1 , y t − 1 ′ , … , y t − p ′ ′ and γ _1t  = (c _1t , A _1•,t )′. Thus, the TVP-VAR can be written as a set of M independent TVP regressions which can be estimated separately using the MCMC methods described in the following section. An additional computational advantage arises in that the M equations can be estimated in parallel using multiple CPUs.

Depending on the particular choice of W , this model nests a variety of commonly used specifications in the literature. For instance, if W implies a random walk behavior of the latent states we arrive at a TVP-VAR closely related to the one proposed in Primiceri (2005). As we will show below, the main difference is that we have a more flexible state equation by allowing for heteroskedasticity in the shocks to the states through dynamic shrinkage priors. Another model that is closely related to ours is the one proposed in Cogley, Primiceri, and Sargent (2010). This model assumes that the variances of the state innovations evolve according to independent stochastic volatility models.

4.1 Dynamic Global-Local Shrinkage Priors

For the time-invariant coefficients, α , we use a horseshoe shrinkage prior (Carvalho, Polson, and Scott 2010). Since the properties of this prior are familiar and posterior simulation methods for this prior are standard, we do not discuss it further here. See the appendix for additional details.

The important contribution of the present paper lies in the development of a dynamic extension of the horseshoe prior for β . We modify methods outlined in Kowal, Matteson, and Ruppert (2019) to design a prior which reflects our beliefs about what kinds of parameter change are commonly found in macroeconomic applications. In particular, we want to allow for a high degree of sparsity in the TVPs. That is, we want a prior that allows for the possibility that parameter change is rare and may occur for only some coefficients in the regression. There may be periods of instability when parameters change and times of stability when they do not. A dynamic global-local shrinkage prior which has these properties is:

where β _t  = (β _1t , …, β _Kt )′ denotes the coefficients at time t, τ denotes a global shrinkage parameter that pushes all elements in β towards zero, λ _t is a time-specific shrinkage factor that pushes all elements in β _t towards zero and ϕ _jt is a coefficient and time-specific shrinkage term that follows a half-Cauchy distribution.

Thus, the prior covariance matrix of β _t is given by:

which implies that λ _t acts as a common factor that aims to detect periods characterized by substantive amounts of time variation.

The main innovation of this paper lies in our treatment of this common factor. Before we discuss the precise specifications for λ _t , it is worth summarizing the key innovation of this prior. As opposed to the dynamic horseshoe of Kowal, Matteson, and Ruppert (2019), we only introduce persistence in the common shrinkage factor λ _t . The key point to note here is that, as opposed to assuming a dynamic law of motion for the coefficient-specific prior scaling parameters, we borrow strength from the cross-sectional dimension and by doing this we substantially reduce the computational burden necessary.

For the global shrinkage parameter we consider four different laws of motion. The first and second of these involve setting g _t  = log(τλ _t ) and assuming it follows an AR(1) process:

with μ = log  τ. We consider two possible distributions for ν _t . In the first of these it follows a four parameter Z-distribution, Z ( 1 / 2,1 / 2,0,0 ) , leading to a variant of the dynamic horseshoe prior proposed in Kowal, Matteson, and Ruppert (2019) (henceforth labeled dHS svol-Z). The second of these follows a Gaussian distribution, leading to a standard stochastic volatility model for this prior variance (labeled dHS svol-N). This model resembles the one stipulated in Cogley, Primiceri, and Sargent (2010) but with a single dynamic volatility process. Both of these processes imply a gradual evolution of g _t and thus a smooth transition from times of rapid parameter change to times of less parameter change.

The third and fourth specifications allow for more abrupt change between times of stability and times of instability. They assume that λ _t is a regime switching process with:

Here, d _t denotes an indicator that either follows a Markov switching model (labeled dHS MS) or a mixture specification (labeled dHS Mix) and κ ₀, κ ₁ denote prior variances with the property that κ ₁ ≫ κ ₀. For the Markov switching model, we assume that d _t is driven by a (2 × 2)-dimensional transition probability matrix P with transition probabilities from state i to j denoted by p _ij (with p i i ∼ B ( a i , M S , b i , M S ) , for i = 0, 1, following a Beta distribution a priori). The mixture model assumes that p ( d t = 1 ) = p ̲ , with p ̲ ∼ B ( a Mix , b Mix ) . In the empirical application we specify κ ₁ = 100/K, κ ₀ = 0.01/K, a _Mix = a _1,MS  = b _0,MS  = 3 and b _Mix = a _0,MS  = b _1,MS  = 30.

We also include a fifth specification by setting λ _t  = 1 for all t. We refer to this setup as the static horseshoe prior (abbreviated as sHS). For these last three specifications (i.e. the ones that do not assume λ _t to evolve according to an AR(1) process), we use a half-Cauchy prior on τ ∼ C + ( 0,1 ) .

4.2 Markov Chain Monte Carlo (MCMC) Algorithm

For all of these models, Bayesian estimation and prediction can be done using MCMC methods. In this sub-section we mainly focus on how to sample λ _t under the assumption that it evolves according to an AR(1) process. For this step we propose a simple and accurate approximation that renders the corresponding hierarchical model linear and conditionally Gaussian. We only briefly discuss the remaining steps since most of them are standard in the literature.

For the time varying regression coefficients, the scalable algorithms (with or without sparsification) of the preceding section, based on Johndrow, Orenstein, and Bhattacharya (2017), can be used. The only modification is that we construct D ₀ as follows:

with λ _t depending on the specific law of motion adopted. Most of the prior hyperparameters introduced in this section have posterior conditionals of standard forms. These are given in the appendix.

Sampling λ _t for the specifications that assume it to be binary is also straightforward and can be carried out using standard algorithms. To sample from the posterior of λ _t under the assumption that it evolves according to an AR(1) process, the algorithm proposed in Jacquier, Polson, and Rossi (1995) can be used. However, since this algorithm simulates the λ _t ’s one at a time mixing is often an issue. A second option would be to view the prior (after squaring each element of β _t and taking logs) as the observation equation of a dynamic factor model. This strategy, however, would be computationally challenging for moderate to large values of K. As a solution, we propose a new algorithm that is straightforward to implement and, if K is large, has good properties.

Let β ̂ t be a K-dimensional vector of normalized TVPs with typical element β ̂ j t = β j t / ϕ j t τ 1 / 2 . Using (5) and squaring yields:

with ν t = v t ′ v t for v t ∼ N ( 0 K , I K ) . Notice that ν _t follows a χ ² distribution with K degrees of freedom, denoted by χ K 2 . This implies that sampling algorithms that rely on the Gaussian mixture approximation proposed in Kim, Shephard, and Chib (1998) cannot be used. Instead we approximate the χ K 2 using a well-known limit theorem that implies, as K → ∞,

This approximation works if K is large. In our case, K is often large. For instance, in the largest TVP-VAR model we consider, K is around 100. Since we estimate the TVP-VAR one equation at a time, values of this order of magnitude hold in each equation and the approximation is likely to be good. But if one were to do full system estimation of the TVP-VAR, there are on the order of MK VAR coefficients at each point in time and the approximation would be even better.

Substituting the Gaussian approximation into (7) and taking logs yields:

Finally, under the assumption that ( 2 K q t + K ) > 0 and by using a Taylor series expansion,^[3] we approximate log v ̂ t with a N log ( K ) − 1 / K , 2 / K to render (8) conditionally Gaussian. This implies that any of the standard algorithms proposed in the literature on Gaussian linear state space models can be used. In this paper, we simulate log  λ _t using the precision sampler outlined, for example, in Chan and Jeliazkov (2009) and McCausland, Miller, and Pelletier (2011).

The accuracy of this approximation for different values of K is illustrated in Figure 1. From this figure it is clearly visible that, if K is greater than 5, our approximation works extremely well. In these cases, there is hardly any difference visible between the log ⁡ χ K 2 and the single-component Gaussian distribution. For K = 1 (the most extreme case) and K = 5, some differences arise which mainly relate to the left tail of the distribution. However, already for K = 5 these differences are so small that we do not expect them to have any serious consequences on our estimates of λ _t , even for small values of K.

Figure 1:

This figure shows the approximation error of a single-component Gaussian used to approximate a log ⁡ χ K 2 distribution. It illustrates the approximation error resulting from approximating the error distribution (which is log ⁡ χ K 2 ) with a single-component Gaussian with mean log(K) − 1/K and variance 2/K. For different values of K, the blue shaded areas show the exact error distribution, while the red shaded areas indicate the approximate error distribution.

5.1 How Does our Algorithm Compare to Exact MCMC?

We start by showing that using our approximate (sparsified) algorithm yields estimates that are close to the exact ones in terms of precision. This is achieved by considering five different data generating processes (DGPs). These are all based on Equation (1) but make different assumptions about the density and nature of parameter change. Dense DGPs are characterized by having time-variation in a large number of parameters (with sparse DGPs being the opposite of dense). The nature of parameter change can be gradual (e.g. characterized by constant evolution of the parameters) or abrupt. For each of the five DGPs, we simulate a time series of length T = 250 and with K = 50.

The different DGPs assume that the states evolve as follows:

1. dense gradual: β t ∼ N ( β t − 1 , 1 100 × I K ) ,

2. dense mixed: β t ∼ N β t − 1 , d t + ( 1 − d t t ) 100 × I K with Prob(d _t  = 1) = 0.1,

3. medium-dense gradual: β t ∼ N ( β t − 1 , d t 100 × I K ) with Prob(d _t  = 1) = 0.3,

4. sparse abrupt: β t ∼ N ( β t − 1 , I K ) with Prob(d _t  = 1) = 0.02,

5. no TVPs: β _t  =  0 _K×1 for all t.

We investigate the accuracy of our scalable approximate MCMC methods relative to the exact MCMC algorithm of Bhattacharya, Chakraborty, and Mallick (2016) (i.e. it is the version of our algorithm which imposes δ _j  = 1 for all j). Table 1 shows the ratio of mean absolute errors (MAEs), computed using the posterior mean of { β t } t = 1 T and the true parameters, for the approximate relative to the exact approach for the five priors averaged over the five DGPs. With one exception, MAE ratios are essentially one indicating that the approximate and exact algorithms are producing almost identical results. The one exception is for the DGP which does not have any TVPs. For this case, the approximate algorithm is substantially better than the exact one. This is because our approximate algorithm uses SAVS which (correctly for this DGP) can set the TVPs to be precisely zero. In this case, draws from the posterior will coincide with draws from the prior that induce heavy shrinkage. Hence, compared to the exact model, the likelihood does not influence the prior and more shrinkage can be achieved.

Thus, Table 1 shows that, where there is substantial time variation in parameters, the approximation inherent in our scalable MCMC algorithm is an excellent one, yielding results that are virtually identical to the slower exact algorithm. The table also shows the usefulness of SAVS in cases of very sparse DGPs.

5.2 How Big are the Computational Gains of our Algorithm?

Our second artificial data experiment is designed to investigate the computational gains of our algorithm relative to exact MCMC for various choices of K, T, degrees of sparsity and data configurations. Since we are only interested in computation time we just generate one artificial data set for each of two different ways of specifying W . The random numbers refered to below are drawn from the standard Gaussian distribution.

For K = 1, …, 400 and T ∈ {100, 200} we randomly draw a y and an X . The W is drawn in two ways which correspond to the flexible and random walk specifications of equations (2) and (3), respectively.

In terms of sparsity, we consider four scenarios based on how we choose W ̃ S :

1. 100 % dense: W ̃ S = W . This is the exact algorithm.

2. 50 % dense: W ̃ S contains 50 % of the columns of W (i.e. s = 0.5k).

3. 10 % dense: W ̃ S contains 10 % of the columns of W (i.e. s = 0.1k).

4. 1 % dense: W ̃ S contains 1 % of the columns of W (i.e. s = 0.01k).

Figure 2 depicts the computational advantages of our approximate MCMC algorithm relative to the exact algorithm of Bhattacharya, Chakraborty, and Mallick (2016). It shows the time necessary to obtain a draw of β . It can be seen that when the TVPs are highly correlated over time as with the random walk specification, then our scalable algorithm has substantial computational advantages relative to the exact algorithm particularly for large K and in sparse data sets. When the TVPs are uncorrelated the computational advantages of our approach relative to the exact algorithm are smaller, but still appreciable.^[4]

Figure 2:

This figure shows the time necessary to obtain a draw for K time-varying coefficients. It depicts the estimation time in seconds required to obtain a draw for K time-varying coefficients for different degrees of overall sparsity (i.e. 1 %, 10 %, 50 %, and 100 % dense). The dots refer to the empirical run times for which we fit a nonlinear trend (indicated by the solid lines). The red colored dots and red solid lines indicate run times of the exact algorithm (100 % dense, with s = k).

6.1 Data Overview and Specification Issues

We illustrate our methods using a monthly data set of 30 government bond yields in the euro area (EA). As opposed to forecasting standard US macroeconomic time series such as output, inflation and unemployment rates, forecasting EA government bond yields is challenging due to, at least, three reasons. The first is that the researcher has to decide on the segment of yield curve she is interested in or use techniques that allow for analyzing the full term structure of government bond yields. Following the latter approach leads to overfitting issues whereas the former approach might suffer from omitted variable bias. The second challenge is that these time series are often subject to outliers as well as sharp shifts in the conditional variance. The final reason is that the time series we consider are rather short and in such circumstances TVP-VARs risk overfitting if the estimates of the TVPs are not regularized sufficiently. We expect that the techniques proposed in this paper are capable of handling both issues well.

We use monthly yield curve data obtained from Eurostat. This dataset includes the yield to maturity of a (hypothetical) zero coupon bond on AAA-rated government bonds of eurozone countries for 30 different maturities. These maturities range from one-year to 30-years and span the period from 2005:01 to 2019:12.

If we wish to model all 30 yields jointly we have to estimate a TVP-VAR with M = 30 equations, a challenging statistical and computational task which we will take on in the next sub-section. Since the parameter space of such a model is vast and difficult to interpret, in this sub-section where we present some in-sample results, we will use a small-scale example. This model is based on the Nelson-Siegel three factor model (see, e.g. Diebold, Rudebusch, and Aruoba 2006; Nelson and Siegel 1987) and assumes that the yield on a security with maturity t , labeled r t ( t ) , features a factor structure:

Here, L _t , S _t and C _t refer to the level, slope and curvature factor, respectively, while η t ( t ) denotes maturity-specific measurement errors which are independent across maturities and feature variance σ η 2 ( t ) . ζ denotes a parameter that controls the shape of the factor loadings. Following Diebold, Rudebusch, and Aruoba (2006), we set ζ = 0.7308 (12 × 0.0609). Since the loading of the level factor is one for all maturities and does not feature a discount factor, it defines the behavior at the long end of the yield curve. Moreover, the slope factor mainly shapes the short end of the yield curve and the curvature factor defines the middle part of the curve. The latent yield curve factors are obtained by running OLS on a t-by-t basis. These estimates are then consequently used as our endogenous variables by setting y _t  = (L _t , S _t , C _t )′ and estimating the TVP-VAR defined in (4). We use the flexible specification for W in (2) and the approximate algorithm to estimate the model. In addition, we set the lag length to two. After obtaining forecasts for y _t , we use (9) to map the factors back to the observed yields. It is worth noting that (9) constitutes an observation equation which links the observed yields to the latent Nelson-Siegel factors. To compute predictive densities, we also take the corresponding measurement errors into account by estimating the measurement error variance independently for each observed series.

6.2 In-Sample Results

To provide some information on the amount of time variation, Figures 3 and 4 depicts heatmaps of the posterior inclusion probability (PIPs) for a Nelson-Siegel model with panels a)–d) referring to the four different dynamic priors for λ _t . These PIPs are the posterior means of the elements of δ .

Figure 3:

This figure shows heatmaps of posterior inclusion probability (PIPs) for time-variation in structural TVP-VAR coefficients with a gradually changing common shrinkage factor. Gray shaded areas indicate coefficients which do not appear in the model due to the lower triangularity of A _0t .

Figure 4:

This figure shows heatmaps of posterior inclusion probability (PIPs) for time-variation in structural TVP-VAR coefficients with a regime-switching common shrinkage factor. Gray shaded areas indicate coefficients which do not appear in the model due to the lower triangularity of A _0t .

The main impression provided by Figures 3 and 4 is that there is little evidence of strong time-variation in the parameters when using this data set. However, there does seem to be some in the sense that there are many variables and time periods where the PIPs are appreciably above zero. That is, even though the figures contain a lot of white (PIPs essentially zero) and just a handful of deep reds (PIPs above one half), there is a great deal of pink of various shades (e.g. PIPs 20 %–30 %). This is consistent with time-variation being small, episodic and only occurring in some coefficients.

Results for our four different dynamic horseshoe priors are slightly different indicating the dynamic prior choice can have an impact on results. A clear pattern emerges only for the dynamic horseshoe prior with a mixture specification. It is finding that small amounts of time-variation occur only for the coefficients on the curvature factor. If the mixture part of the prior is replaced by a Markov switching specification, we tend to find short-lived periods where a small amount of time-variation occurs for all of the coefficients in an equation. But, interestingly, dHS MS finds that different equations have time-variation occuring at different periods of time. Evidence for TVPs is the least when we use stochastic volatility specifications in the dynamic horseshoe priors. For these priors, tiny amounts of time variation (i.e. tiny PIPs) are spread much more widely throughout the sample and across variables.

6.3 Forecast Exercise

The dataset covers the entire yield curve and includes yields from one-year to thirty-year bonds in one-year steps. We choose {1y, 3y, 5y, 7y, 10y, 15y, 30y} maturities as our target variables that we wish to forecast and consider one-month and one-quarter ahead as forecast horizons. We use a range of competing models that differ in terms of how they model time-variation in coefficients and the number of endogenous variables they have. All models feature stochastic volatility in the measurement errors and have two lags. We also offer comparison between the two MCMC algorithms: exact and approximate.

In terms of VAR dimension, we have large TVP-VARs and VARs with all 30 maturities (M = 30) as well as the three factor Nelson-Siegel model described in the previous sub-section (M = 3).

In terms of time variation specified through the likelihood function (i.e. through the definition of W ), we consider the flexible (FLEX) and random walk (RW) specifications defined in (2) and (3). In terms of time variation specified through the prior, we consider the five global-local shrinkage priors (four dynamic and one static) given in Sub-Section 4.1. In addition, we consider as a competitor the conventional TVP-VAR setup of Primiceri (2005). We estimate the TVP-VAR only for the Nelson-Siegel model since the original prior overfits in higher dimensions.^[5]

We also have VAR models where coefficients are constant over time. For these we do two versions, one with a Minnesota prior (MIN) and the other a horseshoe prior (HS). These models are estimated by setting β  = 0 and then using the sampling steps for α detailed in the appendix. For the Minnesota prior, we use a non-conjugate version that allows for asymmetric shrinkage patterns and integrate out the corresponding hyperparameters within MCMC.

To evaluate one-month and one-quarter ahead forecasts, we use a recursive prediction design and split the sample into an initial estimation period that ranges from 2005:01 to 2008:12 and a forecast evaluation period from 2009:01 to 2019:12. We use Root Mean Squared Forecast Errors (RMSEs) as the measure of performance for our point forecasts and Continuous Ranked Probability Scores (CRPSs, Gneiting and Raftery 2007) as the measure of performance of our density forecasts. Both are presented in ratio form relative to the benchmark model which is the large VAR with Minnesota prior. Values less than one indicate an approach is beating the benchmark.

We present our forecasting results in two tables. Table 2 shows the one-month ahead forecast performance of the different models while Table B.1 in the appendix shows the one-quarter ahead forecasting results. Our focus on one-step ahead forecasts is predicated by the fact that the density forecast measures based on proper scoring rules (such as CRPSs) can be viewed as a training sample marginal likelihood and thus enables model comparison (see Gneiting and Raftery 2007).

Overall, the evidence in Table 2 (and Table B.1) is mixed, with no single approach being dominant. In principle, one robust pattern is that models with TVPs tend to produce more accurate forecasts than the large VAR with stochastic volatility benchmark. These gains range from being rather small (particularly at the short-end of the yield curve) to appreciable (when the focus is on the long-end of the yield curve). This is consistent with recent findings in Fischer et al. (2023) who document that flexible models work well for this particular dataset when longer maturities are considered.

If we compare results for the large TVP-VARs to results for the smaller TVP-VARs based on the Nelson-Siegel factors reveals that both specifications produce forecasts of similar quality. When forecasting one-month ahead and focusing on the CRPS as a measure of forecast performance, the best average forecast performance is produced by one of the large TVP-VARs. But when we focus on point forecasting performance, one of the NS models emerges as the best forecasting model. This finding indicates that using more information in an unrestricted manner seems to exert benign effect on higher order moments of the predictive density whereas for the first moment the effect is negligible (or even negative). Interestingly, this finding only holds for one-month ahead predictive densities. When we focus on one-quarter ahead forecasts (see Table B.1 in the appendix), this result is reversed with CRPSs indicating one of the NS models is forecasting best and RMSEs indicating one of the large TVP-VARs is forecasting best.

The comparison of the different choices for W also yields a mixed pattern of results. At the short end of the yield curve the RW specification tends to forecast better, but at the longer end the FLEX specification does better. It is interesting to note, however, that the good performance for RW occurs with a large TVP-VAR whereas for the FLEX specification it occurs for a Nelson-Siegel version of the model.

In terms of which of our dynamic horseshoe priors forecasts best, it does seem to be the priors which assume λ _t to exhibit rapid change between values forecast better than the gradual change of the stochastic volatility specifications. That is, the Markov switching or mixture versions of the prior, dHS MS and dHS Mix, tend to forecast better than dHS svol-Z or dHS svol-N. Although there are several exceptions to this pattern. At this point it is also worth highlighting that the original Primiceri (2005) model is outperformed by our shrinkage specifications in all segments of the yield curve. This suggests that using proper shrinkage priors on the state innovation variances and allowing for dynamic shrinkage pays off.

Thus, overall (and with several exceptions) we have a story where, in this data set, time variation in the regression coefficients is present and there are gains to be made from capturing them. This can be seen by noting that the constant parameter VARs with stochastic volatility are never the best performing specifications across the different maturities and also for both time horizons we consider. As we have shown in the previous sub-section, this time variation is episodic (rather than gradually evolving) and only occurs occasionally and for some of the coefficients. However, ignoring this time variation and using constant parameter models leads to a deterioration in forecasts in almost all situations.

In terms of computation, our scalable algorithm does seem to work well. If we compare results from the exact MCMC algorithm to our approximate (non-sparsified) algorithm, it can be seen that using the computationally-faster approximation is not leading to a deterioration in forecast performance. In fact, there are some cases where the approximate forecasts are better than their exact counterparts.