Inducing Sparsity and Shrinkage in Time-Varying Parameter Models

Inducing Sparsity and Shrinkage in Time-Varying Parameter Models

Introduction

\doublespace Time-varying parameter (TVP) regressions and Vector Autoregressions (VARs) have enjoyed great popularity among econometricians in recent years as a way of modelling the parameter change that occurs in many macroeconomic and financial time series variables. These are state space models which have been found to work well in forecasting giannone2013 and been successfully used for structural economic analysis in a changing environment cogley2005drifts, primiceri2005. They are flexible and capable of modelling almost any nonlinear relationship between explanatory and dependent variables. However, this flexibility comes at a cost: TVP models can be over-parameterized and suffer from the curse of dimensionality, particularly when the number of potential explanatory variables is large. This can lead to very good in-sample fit, but poor out-of-sample forecast performance.

There is a large and growing literature that proposes various methods for overcoming these over-parameterization concerns using Bayesian methods fs_wagner, bkk, kg2014, kowal, uribe_lopes, rockova, koopkorobilisvb, bitto_fs, huberetalJAE, eisenstat_TVPVAR. These papers propose different approaches to obtain more precise inference. Much of this literature uses hierarchical global-local shrinkage priors. A key property of these priors is that they ensure shrinkage in the sense that they pull coefficients towards zero. However, they do not impose them to be exactly zero and, thus, estimation uncertainty remains. In contrast to shrinkage approaches, selection approaches seek to choose a single sparse specification. That is, they select a particular set of explanatory variables and, by doing so, impose coefficients on non-selected explanatory variables to be zero.\footnote{In the Bayesian literature, there are some global-local priors, such as the spike and slab prior, which do select variables, but these are less popular since Markov Chain Monte Carlo (MCMC) algorithms tend to mix poorly.}

Which is better: shrinkage or sparsity? The answer to this question depends on the empirical application. In macroeconomics, there is evidence that shrinkage and sparsity can both play a role. For instance, in the case of constant coefficient regressions and VARs, there is debate among Bayesian econometricians as to whether models are sparse (in which case sparsification methods are appropriate) or dense (in which case shrinkage is appropriate). A recent paper, giannone2017economic, considers a range of data sets in macroeconomics, microeconomics and finance and finds evidence mostly in favor of dense models, a finding reinforced by chp. But there are exceptions to this pattern where sparse models do better. But, instead of choosing one of sparsity or shrinkage, why not do both? This is exactly what recent papers such as hahncarvalho2015dss propose. That is, first shrinkage is done using a Bayesian global-local shrinkage prior and then sparsification is done on the resulting estimates. Such an approach could add the benefits of sparsity (i.e. the reduction in estimation error that is important for improving forecasts) along with the benefits of shrinkage which are so useful with dense data sets. {Recent contributions in finance provide evidence that this works well if interest centers on non-linear modeling of expected returns of companies fisher2018monotonic or constructing optimal portfolios puelz2019portfolio.}

One consideration that arises in some approaches is computation. Bayesian inference with hierarchical shrinkage priors requires computationally-burdensome MCMC methods. Adding a second sparsification step can greatly increase the burden if this step uses cross-validation methods for choosing key tuning parameters. However, in a recent contribution, bhattacharya2018signal propose a simple algorithm, the signal adaptive variable selector (SAVS), for doing the sparsification step. This involves no tuning parameters and is computationally trivial. bhattacharya2018signal provides a theoretical justification for SAVS and shows it to have good empirical performance in simulated and real data contexts.

The papers cited in the preceding paragraph all relate to constant coefficient regression or VAR models rather than the TVP state space models which are the focus of this paper. We develop Bayesian methods for inference and forecasting in TVP regressions and TVP-VARs which both shrink and sparsify. The shrinkage step can be done using any of the hierarchical shrinkage priors that have been used with TVP regressions. In this paper, we use the Dirichlet-Laplace prior bhattacharya2015dirichlet, a fully hierarchical variant of the stochastic search variable selection prior george1993variable, ishwaran2005spike, george2008bayesian, the Horseshoe carvalho2010horseshoe, the Bayesian Lasso of park2008bayesian and the Normal-Gamma prior of griffin2010inference. The sparsity step is done using the SAVS method of bhattacharya2018signal.

Another extension we make in this paper relative to the shrink-then-sparsify methods of hahncarvalho2015dss and bhattacharya2018signal is that we allow for uncertainty in the sparsified estimates. That is, hahncarvalho2015dss and bhattacharya2018signal take the posterior mean from the shrinkage step and use only this in the sparsification step. We sparsify every MCMC draw in the shrinkage step, thus allowing for parameter uncertainty. This feature is crucial if interest centers on computing non-linear functions of the parameters (such as higher order predictive distributions) and allows for uncertainty quantification with respect to the chosen model. Our methods are illustrated with simulated and real data and we find them to improve estimation accuracy and forecast performance.

The remainder of this paper is organized as follows: Section (ref) discusses various global-local shrinkage priors in the context of the regression model with constant coefficients. It describes how the sparsification strategy of bhattacharya2018signal works in the regression model. Section (ref) extends these methods to TVP regressions and TVP-VARs. Section (ref) investigates the performance of our methods relative to non-sparsified alternatives using simulated data from a range of sparse and dense TVP regressions. Section (ref) carries out a forecasting exercise using TVP-VARs. A comparison of forecasts which are both shrunk and sparsified to those which are only shrunk shows the benefits of doing both. Section (ref) concludes the paper and a technical appendix provides further details on the specific prior setup and the posterior simulation algorithms.

Shrinkage and Sparsification in Regression Models

In this section we describe the shrinkage and sparsification methods for regression which we build on in this paper. In the next section, we will show how they can be adapted for dynamic regressions and multiple equation models such as VARs. Consider the regression model:

for $t=1,\ldots,T$ where $y_{t}$ is a scalar dependent variable, $\bm{X}_{t} = (X_{1t}, \dots, X_{Kt})'$ is a $K\times 1$ vector that stacks the explanatory variables $X_{jt}~(j=1,\dots,K)$ and $\bm \beta$ is a $K$-dimensional vector of regression coefficients. The errors are assumed to be independent and follow a zero mean Gaussian distribution with variance $\sigma^2_\varepsilon$.

When $K$ is large relative to $T$, Bayesians increasingly use hierarchical priors so as to induce shrinkage. Global-local shrinkage priors are particularly popular polson2010shrink. These contain shrinkage that is both global (i.e. common to all parameters) and local (i.e. specific to each parameter). We consider priors which can be represented as scale mixtures of Gaussians. In particular, for the $j^{th}$ regression coefficient we assume:

Global shrinkage is controlled by $\lambda$ while $\phi_j$ handles the shrinkage of coefficient $j$. $f$ and $g$ are mixing densities and many different choices have been proposed for them. In this paper, we consider the Horseshoe (HS) prior of carvalho2010horseshoe, the Bayesian Lasso (Lasso) of park2008bayesian, the Normal-Gamma (NG) prior of griffin2010inference, the Dirichlet-Laplace (DL) prior of bhattacharya2015dirichlet and the Normal- mixture of Inverse Gamma (NMIG) prior of ishwaran2005spike, which is a variant of the stochastic search variable selection (SSVS) prior of george1993variable, george1997approaches. All of these are global-local shrinkage priors and differ from one another only in the choices of $f$ and $g$. In addition, and unless otherwise noted, we use a weakly informative inverted Gamma prior on $\sigma_\varepsilon^2$ with hyperparameters $d_{\sigma}=e_{\sigma}=0.01$.

Using any of these global-local shrinkage priors, MCMC methods can be used to carry out posterior inference and calculate the posterior mean, $\hat{\bm{\beta}}$. This estimate has been shrunk, but not sparsified. { It could be that many elements of $\hat{\bm{\beta}}$ will be close to zero and thus imply a small but negligible effect of $X_{jt}$ on $y_t$. For large $K$, this potentially leads to overfitting issues which is a direct consequence of the fact that shrinkage is limited by a lower bound on the degree of certainty (since we always have prior scaling parameters that might be close to but not exactly equal to zero).} Sparsification solves this by taking $\hat{\bm{\beta}}$ and setting small elements in it to zero.

{Sparsification has been advocated both as a way of improving model interpretability as well as improving forecasts. In regression models with many explanatory variables, various sparsification approaches have been proposed to select the most important variables so as to simplify the task of interpreting the results woody2019. The influential paper of barbieriberger shows that, under certain conditions, the median probability model (i.e. a sparsified model which discards all coefficients with inclusion probabilities below $0.5$) has the best forecast performance. In this paper, we implement sparsification using a recently proposed method where the choice of thresholds is made in an optimal way based on a particular decision theoretic problem.}

We implement sparsification using methods developed in hahncarvalho2015dss and bhattacharya2018signal. We first define the SAVS estimate and then offer some explanation and motivation for it. The SAVS estimate is

with $\bm{X}_j=(X_{j1}, \dots, X_{jT})'$ denoting the $j^{th}$ column of a $T \times K$ matrix $\bm X=(\bm X_1, \dots, \bm X_T)'$, $(x)_+=\text{max}(x, 0)$ and $\text{sign}(x)$ returns the sign of $x$. Note that this is a soft-thresholding approach where all values of $\overline{\gamma}_j $ below a certain value are set to zero and that it only acts on the posterior mean.

The sparsified estimate depends on tuning parameters, $\kappa_j$, which determine the thresholds for each coefficient. Various approaches to selecting these have been proposed in the literature including computationally-intensive approaches such as cross-validation. However, bhattacharya2018signal come up with a surprisingly simple solution. This is to set:

This choice implies a penalty for the $j^{th}$ variable which is ranked in inverse-squared order relative to the magnitude of the $j^{th}$ coefficient. With this choice of thresholds, the SAVS estimate is trivial to calculate.

{To provide some motivation for the SAVS estimate note that ((ref)) can be obtained by first solving an optimization problem closely related to the adaptive Lasso zou2006adaptive:}

Equation ((ref)) tries to find a sparse coefficient vector $\bm \gamma$ that is close to $\hat{\bm \beta}$ while introducing a penalty in case of non-zero elements in $\bm \gamma$.

The typical way to solve this optimization problem is using a coordinate descent algorithm friedman2007pathwise. But, as shown in bhattacharya2018signal, if you initialize this algorithm at $\hat{\bm \beta}$ and then do one iteration you get precisely the simple algorithm described in ((ref)) and ((ref)). It is also noted in bhattacharya2018signal that convergence almost always occurs after one iteration and, hence, stopping after one iteration is a sensible thing to do.

One key shortcoming of computing the SAVS estimate is that uncertainty quantification about $\overline{\bm \gamma}$ is not possible and computing non-linear functions of $\overline{\bm \gamma}$ calls for Monte Carlo integration techniques. bhattacharya2018signal highlight that one potential solution to this issue is to replace $\hat{\bm{\beta}}$ with a draw from the full conditional posterior distribution of $\bm \beta$. This is an insight we build upon in the context of the TVP models which are the focus of this paper.

{A recent paper, woody2019, develops methods for improving estimates of posterior uncertainty in sparsified regression and nonparametric regression models (but not the TVP state space models). This paper provides further theoretical justification for our approach. It sets up an optimization problem where the goal is to find a parsimonious summary of the posterior which minimizes a loss function which combines model fit with a reward for parsimony or a restriction that the posterior summary lies in a more parsimonious class of models. As a simple example consider a regression with large $K$. A loss function could be chosen which would find the optimal regression model with $p<K$ explanatory variables. In the context of Bayesian MCMC estimation of a regression model, the algorithm of woody2019 would perform MCMC on the model with $K$ regressors and project each draw into the sparse posterior for the optimal model with $p$ explanatory variables. This is essentially the same strategy as we will adopt below, although our approach is slightly more general in that our posterior summaries are based on all sparsified MCMC draws. To make this point clear in the context of our simple example, the woody2019 algorithm would choose one specific optimal set of $p$ explanatory variables hahncarvalho2015dss and then project the MCMC draws from the $K$ variable regression into the regression model with the chosen set of $p$ variables. Our algorithm, if used in this simple example, would allow for uncertainty about which specific set of $p$ variables is optimal and, thus, allow for model uncertainty. Apart from this difference, the derivations in woody2019 provide a theoretical justification for our approach and, in particular, the measures of posterior uncertainty it produces.}

Shrinkage and Sparsification in TVP Models

In this section, we develop methods for shrinkage and sparsification in state space models such as the TVP regression and the TVP-VAR. This is achieved using the non-centered parameterization of fs_wagner. We emphasize that the algorithms below do the sparsification at each draw from the MCMC algorithm, allowing for treatment of uncertainty in the shrinkage step. Thus, the algorithms are averaging over different sparsified estimators in a manner similar to Bayesian model averaging.

The TVP Regression Model

The TVP regression model used in this paper takes the form:

where all definitions are the same as in ((ref)) except that $\bm{\beta}_t= ( \beta_{1t}, \dots, \beta_{Kt})'$ are dynamic (time-varying) regression coefficients which follow a random walk with $\bm{w}_t$ being Gaussian innovations with zero mean and variance-covariance matrix $\bm{V}=\text{diag}(v_{1}, \dots, v_{K})$. Each $v_j ~(j=1,\dots, K)$ is a process innovation variance associated with the $j^{th}$ coefficient and thus controls the amount of time-variation in $\beta_{jt}$.

The non-centered parameterization of this model is given by:

with the $j^{th}$ element of $\tilde{\bm{\beta}}_t$ given by $\tilde{\beta}_{jt}=\frac{\beta_{jt}-\beta_{j0}}{ \sqrt{v_{j}}}$, $\sqrt{\bm{V}} = \text{diag}(\sqrt{v_1}, \dots, \sqrt{v_K})$ and $\tilde{\bm \beta}_0 =\bm 0_K$. This equation can be written as:

whereby $\bm{\alpha} = (\bm{\beta}_0', \sqrt{v_1}, \dots, \sqrt{v_K})'$, $\bm{Z}_t = [\bm{X}'_t, ( \tilde{\bm{\beta}}_t \odot \bm{X}_t)']'$ and $\odot$ denotes element-wise multiplication. Conditional on knowing the full history of the states in $\bm \tilde{\bm \beta}_t$, ((ref)) resembles a standard regression model with a (partially) latent covariate vector $\bm Z_t$.

Well-developed MCMC methods exist to carry out Bayesian posterior and predictive inference in state space models such as the TVP regression model under various priors. In this paper, we simulate the full history of the normalized dynamic regression coefficients $\{\tilde{\bm \beta}_t\}_{t=1}^T$ using the forward-filtering backward-sampling algorithm proposed in carterkohn and fs1994. Conditional on $\tilde{\bm \beta}_t$, ((ref)) is a standard regression model, implying that we can simulate $\bm \alpha$ from a Gaussian full conditional posterior distribution and $\sigma^2_\varepsilon$ from an inverted Gamma distribution. The corresponding moments take standard forms and are presented in Appendix (ref).

We propose to do shrinkage on $\bm{\alpha}$ using the global-local mixture priors mentioned in the previous section and described in Appendix (ref). That is, conditional on a draw of the full history of the states, $\{\tilde{\bm{\beta}}_t\}_{t=1}^{T}$, we have the regression model given in ((ref)), and shrinkage can be done exactly as described in the preceding section. For each of the global-local mixture priors we consider, MCMC methods for drawing $\bm{\alpha}$ and $\sigma^2_\varepsilon$, conditional on draws of the states exist. For the Dirichlet-Laplace prior we follow the methods of bhattacharya2015dirichlet. For the NMIG specification, we adopt the algorithm proposed in ishwaran2005spike while for the Horseshoe, the MCMC algorithm developed in makalic2016simple is used. Since the Normal-Gamma prior nests the Bayesian Lasso, we adopt the algorithm put forth in griffin2010inference (see Appendix (ref) for further details).

As highlighted in Section (ref), using shrinkage implies that elements in $\bm \alpha$ are pushed to zero and elements in $\bm Z_t$ might have a small effect on $y_t$. However, in the TVP regression setting, this problem is intensified since the state equation can be written in terms of the sum of the past shocks to the states $\bm w_t$. The corresponding variance of $\bm \beta_t$ thus increases with time and values of $\sqrt{v_j}$ that are close to zero could still induce large aggregate movements in $\beta_{jt}$ over time. In such a situation, sparsification might help since setting $\sqrt{v_j}=0$ directly implies that $\beta_{jt}=\beta_{jt-1}$ for all $t$.

Given a draw from the posterior of $\bm \alpha$, denoted as $\bm \alpha^{(n)}$, from any of the MCMC algorithms is sparsified using SAVS. Applying the SAVS estimator in ((ref)) to each draw from the posterior of $\bm \alpha$ yields:

where $\bm Z_j$ denotes the $j^{th}$ column of $\bm Z = (\bm Z_1, \dots, \bm Z_T)'$, $\kappa_j = |{\alpha}^{(n)}_j|^{-2}$ and $N$ denotes the number of post burn-in MCMC draws. This procedure effectively allows for uncertainty quantification and the computation of potentially non-linear functions of the sparsified parameters such as higher-order forecasts or impulse response functions. Thus, one can think of our proposed procedure as an approximate MCMC algorithm which draws from the sparsified conditional posterior $p(\overline{\bm{\gamma}}| \bm \alpha, \bm Z)$.\footnote{The algorithm is approximate since $\sigma^2_\varepsilon$ does not play a role in the SAVS algorithm. If desired, after each sparsification, one could take a draw of $\sigma_\varepsilon^2$ conditional on the sparsified estimates.} Hence, forecasts produced will average over different sparsified models. That is, one MCMC draw will lead to one particular sparsified model which is used for forecasting, another draw may choose another sparsified model to produce forecasts. Hence, what we are proposing is similar in spirit to Bayesian model averaging. {This feature allows us to calculate posterior inclusion probabilities (PIPs) for each variable. The PIP for a given coefficient is the proportion of MCMC draws for which the coefficient is not set to zero.}

Another possibility would be to use the SAVS algorithm directly on the posterior mean of $\bm{\alpha}$ as is done by hahncarvalho2015dss and bhattacharya2018signal. This procedure yields a point estimate for the time-invariant coefficients and the state innovation variances. {However, one shortcoming of doing this is that $\bm Z_t$ includes latent quantities that need to be integrated out or a plug-in estimate (such as the posterior mean) might be used. However, as puelz2017variable note, this could negatively impact inference since the corresponding uncertainty surrounding $\bm Z_t$ is ignored. Our approach circumvents this by integrating out the latent states contained in $\bm Z_t$. In addition, if the researcher wishes to select a single sparse model, as produced by sparsifying the posterior mean directly, our approach provides an alternate way of choosing the sparsity pattern based on PIPs. }

Another point worth emphasizing about our algorithm is that it is fast. Relative to the computational time required to do MCMC, adding the SAVS step increases the computational burden by a trivial amount. For any empirical specification where MCMC is possible, our proposed algorithm is also possible. Of course, if $K$ is too large, then MCMC methods may be computationally infeasible. In such a case, variational Bayesian methods may be a practical alternative koopkorobilisvb. But with variational Bayes methods, the SAVS algorithm would be applied on the approximate posterior mean and model uncertainty ignored.\footnote{It would be possible to surmount this drawback of variational Bayes by first using variational Bayes to obtain an approximation to the posterior and then applying the SAVS algorithm to draws from this approximation. But this would be computationally demanding, thus undermining the main advantage of variational Bayes.}

The TVP-VAR

The shrink-then-sparsify algorithm we propose for the TVP regression can be extended to handle the TVP-VAR in a straightforward fashion. The idea is to transform the TVP-VAR so that the error covariance matrix in the measurement equation is diagonal. Then the TVP regression algorithm of the preceding sub-section can be applied one equation at a time. Equation-by-equation estimation of VARs is done in several recent papers using transformations similar to the one used here kastnerhuber2017, kpp2019, ccm2016 and the reader is referred to these papers for further details about the computational advantages of this approach. With macroeconomic data it is often important to add stochastic volatility (SV), which leads us to the TVP-VAR-SV specification described in this section.

Let $\bm{y}_{t}$ be an $M \times1$ vector of endogenous variables for $t=1,\ldots,T$. The TVP-VAR-SV can be written as:

where $\bm{X}_t = (\bm{y'}_{t-1}, \dots, \bm{y'}_{t-P}, 1)'$ contains the $P$ lags of $\bm{y}_t$ and an intercept, $\bm{\beta}_t$ is the vector $K=M (MP +1)$ coefficients at time $t$ which is assumed to evolve according to a multivariate random walk. The errors are independent over time with $\bm \varepsilon_t \sim \mathcal{N}(\bm{0}_M, \bm{\Sigma}_t)$. $\bm{\Sigma}_t$ is the time-varying error covariance matrix with

Let $\bm{U}_t$ denote a lower uni-triangular matrix and $\bm{H}_t = \text{diag}(e^{h_{1t}}, \dots, e^{h_{Mt}})$. The $M (M-1)/2$ free elements in $\bm U_t$ follow independent random walks while the $h_{jt}$'s are log-volatilities that evolve according to AR(1) processes,

Here, we let $\mu_j$ denote the unconditional mean, $\rho_j$ the persistence parameter and $\sigma_{\eta, j}^2$ the error variance of the log-volatility process. The initial state $h_0$ is drawn from the stationary distribution of the process. The prior specification on the parameters of the log-volatility equation closely follows kastner2014ancillarity. Specifically, we use a zero mean Gaussian prior with variance $10^2$ on $\mu_j$, a Beta prior on $\frac{\rho_j+1}{2} \sim \mathcal{B}(25, 5)$ and a Gamma prior on $\sigma_{\eta, j}^2 \sim \mathcal{G}(1/2, 1/2)$. This Gamma prior translates into a Gaussian prior on $\pm \sigma_{\eta, j}$ with zero mean and unit variance. In the MCMC algorithm, the full history of $h_{jt}$ as well as the parameters of equation ((ref)) are obtained using the algorithm proposed in kastner2014ancillarity. This algorithm exploits the centered and non-centered parameterization of the non-linear state space model to increase sampling efficiency and samples the full history of the log-volatilities from a $(T-1)$-dimensional multivariate Gaussian distribution.

As noted in ccm2016, kastnerhuber2017, kpp2019, computation is greatly simplified if the model is transformed so that the errors in different equations are independent of one another. This can be achieved by augmenting the $i^{th}$ equation in the system with the contemporaneous values of the first $i-1$ elements in $\bm y_t$. That is, if $y_{it}$ is the $i^{th}$ variable (for $i>1$), we can write the TVP-VAR-SV as a set of $M$ unrelated TVP regressions:

where $\eta_{it}$ and $\eta_{jt}$ are independent for $i\ne j$, $\bm{\beta}_{it}$ denotes the elements of $\bm{\beta}_t$ in the $i^{th}$ equation and $u_{ij,t}$ are the elements of $\bm{U}_t^{-1}$ for $i=2,\dots, M; j=1,\dots,i-1$.

We then write the TVP-VAR-SV using the non-centered parameterization. For equation $i$ we obtain:

Here, we let $\sqrt{\bm{V_{i}^{\beta}}} = \text{diag}\left(\sqrt{v_{i1}^{\beta}}, \dots, \sqrt{v_{iK}^{\beta}}\right)$ and $\sqrt{v_{ij}^{u}}$ denotes the standard deviation of the error in the random walk state equation for the $j^{th}$ VAR coefficient in the $i^{th}$ equation. Similarly, $\sqrt{v_{ij}^{u}}$ is the standard deviation for the random walk state equation for the elements of $\bm{U}_t$. Thus, $\tilde{\bm{\beta}}_{it}$ and $\tilde{u}_{ij,t}$ are the states for equation $i$ and the shocks in the corresponding state equations have unit standard deviation.

Since the errors in the different equations are independent of one another, estimation of one equation at a time using the algorithm of the preceding sub-section, including the SAVS step detailed in ((ref)), can be done. Computation is also sped up since parallelization is feasible. Note also that, since the coefficients in $\bm{U}_t^{-1}$ are appearing as regression coefficients in ((ref)), these can also be shrunk and sparsified. In large TVP-VARs, where there are many such error covariance terms, this is potentially beneficial for forecasting purposes. Notice that we do not only obtain a sparse error covariance matrix but also allow for checking whether the corresponding free elements are time-varying or constant.

Evidence Using Artificial Data

In this section, we present evidence on the performance of the proposed methodology using artificial data generated from different TVP regression models. Across the different data generating processes (DGPs), the covariates are drawn from a Uniform distribution bounded between $-1$ and $1$. The $\bm \beta_t$'s are generated using the non-centered parameterization with $\bm \beta_0 \sim \mathcal{N}(\bm 0_K, 0.1^2 \bm I_K )$ and $\pm \sqrt{v_j} \sim \mathcal{N}(0, 0.1^2), j=1,\dots, K,$ while differing percentages of the elements in $\bm \alpha$ are randomly set to zero. The measurement error variance $\sigma^2_\varepsilon$ is set equal to $0.1^2$.

Before presenting results using repeated samples, the main features of sparsification are illustrated in Figure (ref). The results in the three panels of the figure are based on the Horseshoe prior and use three different single artificial data sets obtained by simulating $T=400$ observations from a large ($K=30$) dynamic regression model. Figure (ref)(a) plots posterior features of $\beta_{jt}$ against time for a case where it is zero (i.e. the DGP is one where $j^{th}$ regressor is not selected) using a non-sparsified and sparsified estimator. Note that the sparsified estimator is precisely correct, it sets $\beta_{jt}=0$ with probability one. Thus, it exactly coincides with the true value and cannot be seen in Figure (ref)(a). The non-sparsified posterior distribution, although the posterior mean is very close to the correct value, has a credible interval that is non-negligible. This reflects estimation uncertainty and could spill over into poor forecast performance using the non-sparsified posterior. The performance of the SAVS algorithm when $\beta_{jt}$ is a non-zero constant (i.e. the DGP is one where $\beta_{jt}=\beta_{jt-1}$ for all $t$) is shown in Figure (ref)(b). In this case, the posterior distributions of the sparsified and non-sparsified models almost coincide. Notice, however, that the credible sets are constant over time for the sparsified model, indicating that the corresponding element in $\sqrt{\bm V}$ is set equal to zero throughout all iterations of the MCMC algorithm. In contrast, Figure (ref)(c) illustrates a case where $\beta_{jt}$ is non-zero and time-varying. Notice that the sparsified and non-sparsified posterior distributions are close to being identical. In this case, it is not desirable to sparsify the corresponding elements in $\bm \alpha$ and the SAVS algorithm is not doing so. Thus, regardless of whether a coefficient is zero, a non-zero constant or time-varying, this figure indicates that our methods estimate it well. They work better than the non-sparsified alternative in cases where there is sparsity and equally well in non-sparse cases.

Table (ref) presents evidence for the importance of sparsification and shrinkage in TVP regression models using different data configurations, priors, numbers of regressors and sample sizes. The DGP described above is modified to reflect varying degrees of sparsity. These different sparsity levels are labeled sparse (with $90$% zeros in $\bm \alpha$), moderate (with $70$% zeros) and dense (with $30$% zeros). To assess how our techniques perform across model dimensions and length of time series involved, we consider variants of each sparsity level with $K=5, 15$ and $30$ explanatory variables and $T=250$ and $400$ observations. The latter are typical values in quarterly and monthly macroeconomic data sets. For each DGP, we generate $100$ artificial data sets and then run each through a sparsified and non-sparsified algorithm using each of the five global-local shrinkage priors listed in Section (ref). We also include a non-informative prior (labelled Flat in the table) which does not do shrinkage. Posterior medians of $\bm \beta_{t}$ are produced and the absolute value of the difference between these and the true value used in the DGP is calculated. The figures in Table (ref) are averages taken over three dimensions: i) the $100$ simulated data sets, ii) time and iii) the $K$ elements of $\bm \beta_t$.

Table (ref) shows the value of sparsification, particularly with sparse DGPs. With the latter, mean absolute errors (MAEs) are lower than their non-sparsified counterpart for every prior and choice for $T$ and $K$. But even in moderately dense specifications, sparsification lowers MAEs in most cases. In the dense specification, sparsification does not improve upon the single best performing non-sparsified model specification. However, in that situation, accuracy differences are found to be negligible.

The benefits of shrinking and sparsifying increase with the number of explanatory variables. Note, for instance, that the unsparsified Flat prior model does not perform that poorly when $K=3$ and $15$, but displays a weak performance when $K=30$. In fact, when $K=3$, the Flat prior works quite well with the sparse specification, provided sparsification is done. This indicates that there are some cases where sparsification is more important than shrinkage.

The choice of $T$ has little impact on the results. In a regression model with constant parameters, we would expect sparsification to be less important as the sample size increases since, with longer time series, the estimation error would decrease. However, with TVP regressions, the number of parameters is also increasing with the sample size which negates this effect. Thus, even with large numbers of observations, the researcher working with TVP models can still benefit from sparsification.

With regards to the different global-local shrinkage priors, no clear pattern emerges where one performs consistently the best across different specifications. When $K=30$ and the DGP is sparse, DL (for $T=250$) and HS (for $T=400$) models that are sparsified are the best performers. When $K=30$ and the DGP is dense, the accuracy of both, the DL and the HS prior deteriorates slightly while the unsparsified NMIG model shows the best performance. Notice that in this situation, accuracy differences across the sparsified and non-sparsified NMIG specification are, however, quite small.

From this discussion it is apparent that identifying a default prior choice is difficult. One key take away from this analysis, however, is that if the DGP is sparse, flexible shrinkage specifications such as the HS, the DL and the NMIG prior in combination with the SAVS algorithm provide accurate parameter estimates. Overall, the table tells a story of the importance of both shrinkage and sparsity, especially in large models, with the precise choice of shrinkage prior being of lesser importance.

In the next step, we assess how well the SAVS algorithm identifies true zeros in $\bm \alpha$. Table (ref) shows average hit rates that measure the percentage of correctly estimated zeros using the SAVS algorithm. From this table, we observe that irrespective of the priors used, our approach works well in identifying the true level of sparsity. For sparse situations, the fraction of correctly identified zeros is often above 95% for most shrinkage priors and model sizes considered. In the case of a Flat prior, we observe values just above 90%, which is remarkable but still well below the percentages observed for the different shrinkage specifications under scrutiny. This slightly weaker performance can be traced back to the fact that without shrinkage, values in $\bm \alpha$ are not pushed to zero and the corresponding penalty $\kappa_j$ is too small. Consistent with the findings in Table (ref), we find no discernible differences in performance across the different shrinkage priors, with all of them displaying a strong performance. In fact, in a sparse setting with $K=30$, the SAVS algorithm identifies almost all zeros correctly, with hit rates being above $99$%.

To sum up, this discussion highlights that sparsification improves estimation accuracy. These improvements tend to increase with model size and the level of sparsity of the DGP. Among the set of competing shrinkage priors, we find no single best performing specification. In terms of correctly predicting zeros in $\bm \alpha$, we found that SAVS works well across all shrinkage priors considered, often correctly identifying above $99$% of the zeros. At this point, and before proceeding to the empirical application, it is worth emphasizing that our analysis only considers whether our shrink-then-sparsify approach improves accuracy of point estimates, ignoring a potential bias-variance tradeoff. One key finding is that applying SAVS never significantly decreases estimation accuracy and correctly predicts a large fraction of true zeros. {In light of Figure (ref), this indicates that, by zeroing out shrunk coefficients, our approach pushes the posterior variance to zero and this could increase predictive accuracy in forecasting applications.}

Forecasting US Macroeconomic Variables

In this section we present results from a forecasting exercise using US quarterly macroeconomic data taken from the FRED-QD database mccracken2016fred that span the period from 1959Q1 to 2017Q4. We focus on forecasting GDP, inflation (based on the GDP deflator) and the Fed Funds rate (henceforth labeled focus variables). Table (ref) provides details on the specific variables included alongside transformations used.

We use the TVP-VAR-SV of Sub-section (ref) combined with the same set of global-local shrinkage priors as in the preceding section. The only specification we do not consider here is the TVP-VAR-SV with a flat prior since this model performs poorly in out-of-sample forecasting and large dimensions.\footnote{The results for the flat prior model are available upon request from the authors.} {More specifically, using weakly informative priors leads to overfitting, which in turn translates into model instability since no penalty is introduced to rule out explosive regions of the parameter space. woody2019 note that in such situations, we fit the noise in the first stage, leading to insufficient posterior variability in the summary.}

For each prior, we use non-sparsified and sparsified versions of the model in order to produce the forecasts. We forecast with small ($M=3$), medium ($M=8$) and large ($M=20$) data sets and set the lag length equal to $2$. Thus, the dimension of the state space in the TVP-VAR-SVs ranges from being moderate to huge. Our forecast evaluation begins in 1997Q1 and runs to the end of the sample. We use root mean squared forecast errors (RMSEs) to evaluate the quality of the point forecasts and average log predictive likelihoods (LPLs) to evaluate the quality of predictive densities. Both are benchmarked relative to a VAR-SV with DL prior, a specification that works well for US macroeconomic data kastnerhuber2017. This is identical to the TVP-VAR-SV with DL prior except that the DL prior now applies directly to the constant VAR coefficients while $\sqrt{\bm V_i^\beta}$ and $\sqrt{v_{ij}^{u}}$ are set equal to zero for all $i, j$. The VAR is transformed to allow for equation-by-equation estimation as described in Sub-section (ref).

Before presenting the results of our forecasting exercise, we present Figure (ref) which sheds light on which variables our algorithm is choosing to predict the focus variables. This figure is produced using the large data set and the HS prior. Previously, we have discussed how doing sparsification for each MCMC draw shares similarities with Bayesian model averaging, allowing us to analyze PIPs. Figure (ref) is a heatmap of these PIPs at the end of the sample. Remember that, in the non-centered parameterization of the TVP-VAR-SV (see equation ((ref))), there are coefficients which appear on the initial states which are constant coefficients. The upper panel of the figure relates to these. The remaining coefficients determine whether there is time-variation relative to the constant coefficients. The lower half of the figure relates to these.

Figure (ref) shows that our methods are inducing a high degree of sparsity in the TVP-VAR-SV in that most of the PIPs are near zero. However, a few of them are not. In terms of the VAR coefficients there is only one coefficient which is always selected (i.e. has a PIP of one). This is the first lag of the 1-year Treasury Bill rate in the equation for the Fed Funds rate. However, an appreciable number of other predictors have PIPs that are substantially above zero but much less than one. In terms of the error covariance matrix, a similar pattern emerges. There is only one error covariance term which is non-zero in every MCMC draw.\footnote{Note that the other green areas refer to the diagonal elements of $\bm U_t$.} This is the covariance between the errors in the equations for two different inflation measures. However, there are several other error covariances with PIPs that are substantially above zero, even if they are below one. We stress that such a finding would not be possible if we were to use the SAVS algorithm directly on the posterior mean as opposed to using it on each MCMC draw. In the former case every PIP would be either zero or one with no values in between.

These patterns are consistent with those found in giannone2017economic who conclude "there seems to be a lot of uncertainty about whether certain predictors should be included in the model, which results into their selection only in a subset of the posterior draws. These findings reflect a substantial degree of collinearity among many predictors that carry similar information, hence complicating the task of structure discovery. In sum, model uncertainty is pervasive and the best prediction is obtained as a weighted average of several models." These features seem to be exactly what our algorithm is uncovering in an automatic fashion.

Finally, it is worth noting that there is evidence of time variation in several of the coefficients and our algorithm is automatically deciding which ones to allow to be time-varying. That is most of the PIPs which are appreciably above zero in the top half of the figure are also above zero in the bottom half. This pattern indicates a non-zero coefficient at time zero which is time varying. But our method also allows for a coefficient to be non-zero but constant. There are some cases which provide evidence of this. For instance, in the GDP growth equation the first lag of S&P500 stock returns has a PIP which is appreciably above zero in the top half of the figure, but is much closer to zero in the bottom half of the figure. This pattern indicates support for a constant coefficient on this predictor.

The evidence in Figure (ref) suggests that shrinking then sparsifying is working in a sensible fashion. But the key test of our methodology is how well it forecasts. Table (ref) presents the results of our forecasting exercise. A comparison of each set of sparsified forecasts to its non-sparsified counterpart shows the benefits of our shrink-then-sparsify strategy, particularly in large models. For $M=8$ and $M=20$, sparsification leads to substantial improvements in both RMSEs and LPLs in almost every case. These improvements are particularly noticeable for GDP forecasting for the one-step-ahead forecasts. In general, the benefits of sparsification are largest when using the DL or Lasso priors. For $M=3$ the benefits of sparsification are less pronounced. In terms of RMSEs, there seems to be no benefits of sparsification, although it does lead to slight improvements in the density forecasts even for this already fairly parsimonious case. This smaller accuracy premium from sparsification can be traced back to the fact that, in small models, {the detrimental influence of irrelevant but non-zero regression coefficients on predictive accuracy is small. In larger models, this effect eventually accumulates, leading to inflated posterior uncertainty and a decreased forecasting accuracy. }

In relation to the benchmark VAR-SV model, it is interesting to note that it is inferior to the TVP-VAR-SV models for the small and medium data sets. Clearly, addition of time-variation in the VAR coefficients helps improve forecasts in these cases. However, in the large data set, the evidence is mixed. In this case, the RMSEs produced by the TVP-VAR-SV are substantially better than those produced by the VAR-SV. However, the density forecasts are not. This could be due to the fact that there is typically a tradeoff between model dimension and parameter change. In small models, there is often a need for a high degree of parameter change to adequately fit patterns in the data and alleviate potential omitted variable biases. But in larger models, the information provided by the additional variables can fit these patterns, leaving less of a role for parameter change. Thus, in high dimensional cases the VAR-SV might be adequate and the extra flexibility provided by a TVP-VAR-SV may not be required. Of course, if the correct specification has a zero coefficient, the non-sparsified approach would try and estimate the time-varying coefficient to be constant over time. But, as illustrated in Figure (ref), estimation uncertainty (although reduced) would still exist which could potentially hurt the forecasting performance of our approach. Sparsification as done in this paper clearly helps, but in the large data set there are still some cases where the VAR-SV is superior. In such cases, a simple extension of our shrink-then-sparsify approach could help. In this paper, we have focused on sparsifying $\bm \alpha$ in equation ((ref)). But any function of the parameters of a model could be sparsified in the same manner and, in particular, sparsifying the change in the states would be possible. This would lead to the constancy of a coefficient over certain periods in time while allowing for movements in other points in time when this kind of sparsification is applied.

{The results in Table (ref) highlight that, when the full hold-out period is considered, sparsification often improves predictive accuracy relative to a non-sparsified model specification. In selected cases, however, the SAVS step also seems to hurt forecasting accuracy. This raises the question whether a practitioner should always use the sparsification step. The findings in the table indicate that relative accuracy gains obtained from using SAVS increase with model size. In large models (with $M=20$), improvements in predictive accuracy for point and density forecasts are often substantial (see, e.g., the improvements for GDP growth and short-term interest rates) while in small models, differences are often negligible and sometimes favor non-sparsified models. This suggests that in large models, using SAVS appears to improve forecasts. However, we would like to stress that applying SAVS yields predictions that are always competitive relative to non-sparsified competitors, even in small-dimensional models. This is corroborated by diebold1995paring tests which suggest that in most cases where the SAVS step lowers predictive accuracy, these decreases are statistically insignificant whereas in the case that sparsification improves forecast accuracy, the differences are often significant. Thus, we can recommend applying the sparsification step in all large-dimensional cases (i.e. for $M\ge 7$) since the computational burden is not increased significantly while the predictive performance is adversely affected in only a few situations in a significant manner. }

The discussion in the previous paragraph provides a simple recommendation that is based on using the full hold-out period. In the next step, we ask whether accuracy differences could also be specific to certain periods in time. To this end, Figure (ref)(a) shows the evolution of the log predictive Bayes factor between the sparsified and non-sparsified large-scale TVP-VAR-SV with the HS prior over the hold-out period.\footnote{Comparable figures for other shrinkage priors reveal similar patterns. Thus. for the sake of brevity, we discuss the results for the HS prior exclusively.} This Bayes factor is obtained by evaluating the one-step-ahead predictive density for the three focus variables jointly after integrating out the remaining variables. To investigate whether the gains in density forecasting performance stem from capturing higher order moments in the predictive distribution or from a more precise point forecast, Figure (ref)(b) shows cumulative squared one-step-ahead forecast errors averaged across the focus variables over time.

Figure (ref)(a) indicates that accuracy premia from sparsification tend to vary significantly over the business cycle. During expansionary stages, sparsification yields modest (in the case of the medium-sized model) to sustained (in the case of the large model) improvements in density forecasting performance relative to the non-sparsified competitor. For the small-scale TVP-VAR-SV, accuracy gains are more muted during expansionary periods. During recessions, in contrast, sparse models tend to be outperformed by their non-sparsified counterparts. Our conjecture is that this stems from the fact that during turbulent times, the sparsified predictive distributions feature a smaller variance, making it harder to capture outlying observations and thus translating in lower log predictive likelihoods.

Our conjecture is confirmed when focusing on point forecasts. In terms of point predictions, we observe that forecast errors are almost identical in the period up to the global financial crisis. During the recession in 2008/2009, forecast errors increase markedly but slightly less so for the sparsified model and for the medium and large dataset. This suggests that the drop in the log predictive Bayes factor is mainly driven by higher order moments, implying that while the accuracy of the point prediction increases, adverse movements in the corresponding predictive variance offset this gain.

Conclusions

Global-local shrinkage priors have enjoyed great popularity in over-parameterized regressions and VARs involving large numbers of variables. And, increasingly, they have been used with TVP versions of these models which are potentially even more over-parameterized. Use of such priors can potentially reduce estimation error and improve forecasts. However, estimation error is not completely eliminated and it is possible that further improvements in forecasting performance can be achieved by adding an additional sparsification step to shrunk estimates to further reduce the lower bound on accuracy associated with shrinkage. In this paper, we have developed methods for doing so. Compared to a recent paper, woody2019, who develop methods for quantifying posterior uncertainty in sparsified models based on using a single sparse model, our approach controls for model uncertainty by sparsifying each MCMC draw. In an artificial data exercise, we have shown that our shrink-then-sparsify approach to TVP regression leads to more accurate estimates for a variety of DGPs. Particularly large gains are found in sparse DGPs. In a macroeconomic forecasting exercise, adding sparsification to shrinkage also leads to substantial improvements in forecast performance if interest centers on using large models.