Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
74,416 characters · 16 sections · 71 citation commands
Dynamic Shrinkage Priors for Large Time-varying Parameter Regressions using Scalable Markov Chain Monte Carlo Methods
\thispagestyle{empty}
\doublespacing
\doublespace The increasing availability of large data sets in economics has led to interest in regressions involving large numbers of explanatory variables. Given the evidence of instability and parameter change in many macroeconomic variables, there is also an interest in time-varying parameter (TVP) regression models and multi-equation extensions such as time-varying parameter Vector Autoregressions (TVP-VARs). This combination of large numbers of explanatory variables with TVPs can lead to regressions with a huge number of parameters. But such regressions are often sparse, in the sense that most of these parameters are zero. In this context, Bayesian methods have proved particularly useful since Bayesian priors can be used to find and impose this sparsity, leading to more accurate inferences and forecasts. A range of priors have been suggested for high-dimensional regression models ishwaran2005spike, park2008bayesian, griffin2010inference, carvalho2010horseshoe, bhattacharya2015dirichlet. There is also a growing literature which extends these methods to the TVP case. Examples include bkk, kg2014, eisenstat2016stochastic, kowal2019dynamic, petrova2019quasi, kalli2019bayesian, bitto2019shrinkage, chan2020reducing, hauzenberger2022fast and fischer2023general.
Most of these papers assume particular forms of parameter change (e.g., it is common to assume parameters evolve according to random walks) and use computationally-demanding Markov Chain Monte Carlo (MCMC) methods. The former aspect can be problematic (e.g., if parameter change is rare and abrupt, then a model which assumes all parameters evolve gradually according to random walks is inappropriate). The latter aspect means these methods are not scalable (i.e., MCMC-based methods cannot handle models with huge numbers of coefficients).
The contributions of the present paper relate to issues of prior elicitation and computation in TVP regressions. With regards to prior elicitation, we develop novel dynamic shrinkage priors for TVP regressions. These modify recent approaches to dynamic shrinkage priors in papers such as kowal2019dynamic. We work with the static representation of the TVP regression model which breaks the coefficients into two groups. One group contains constant coefficients (we call these $\bm \alpha$). The other, which we call $\bm \beta$, are TVPs. In the static representation, the dimension of $\bm \beta$ can be enormous. Our dynamic global-local shrinkage priors are carefully designed to push unimportant elements in $\bm \beta$ to zero in a time-varying fashion. This is done using a global shrinkage parameter that varies over time as well as local shrinkage parameters. The global shrinkage parameter has an interpretation similar to a dynamic factor model with a single factor. This single factor can be used to find periods of time-variation in coefficients and periods when they are constant. Since the assumption of a common volatility factor hampers the use of standard stochastic volatility MCMC algorithms based on a mixture of Gaussians approximation kim1998stochastic, we propose a simple approximation that works particularly well in high dimensional settings.
With regards to computation, we develop a scalable MCMC algorithm. This algorithm is suitable for cases where the posterior for $\bm \beta$, conditional on the other parameters in the model, is Gaussian. This occurs for a wide range of global-local shrinkage priors including the dynamic shrinkage priors used in this paper. In this case, the exact MCMC algorithm of bhattacharya2016fast is the state of the art.\footnote{{kastner2020sparse, hauzenberger2021flexible and korobilis2022new, for example, use this exact algorithm in the context of large VARs to reduce the computational burden of estimating these models.}} However, even it is too computationally slow to handle the huge number of regressors that appear in the static representation of the TVP regression model. Recently, johndrow2017bayes has proposed an approximate algorithm based on this exact algorithm which is computationally much more efficient in sparse models and, thus, is scalable.
In our paper, it is precisely this scalable MCMC algorithm which forms the basis of the algorithm we use. It involves a thresholding step (described below) which we implement in a different manner than johndrow2017bayes. In particular, as opposed to fixing the threshold to a small number, we set it adaptively. Since this would typically imply a number of thresholds that match the dimension of $\bm \beta$, we use a method called Signal Adaptive Variable Selection (SAVS), see bhattacharya2018signal, to determine the thresholds in a novel way. SAVS has the advantage of being computationally fast and easy to implement. Recent papers use SAVS for determining variable relevance hahncarvalho2015dss, portfolio applications puelz2020portfolio or improving macroeconomic forecasts hko2020. We solely use SAVS to identify which variables can be safely set to zero in order to construct an approximate posterior distribution for the TVPs. Thus, the use of SAVS in the context of the algorithm of johndrow2017bayes provides two-fold benefits: computational improvements and more flexibility due to its adaptive nature.
We investigate the use of our methods in artificial and real data. The artificial data exercise demonstrates that our scalable algorithm is a good approximation to exact MCMC and that its computational benefits are substantial. Our application to the eurozone yield curve shows how our methods can effectively pick out small amounts of occasional parameter change in some parameters. Furthermore, allowing for such change in the coefficients improves forecasts.
The remainder of the paper is organized as follows. The second section defines the TVP regression and TVP-VAR models used in this paper. The third section discusses MCMC methods for the regression coefficients and introduces our computationally-efficient approximate method. Section 4 develops different dynamic shrinkage priors and discusses Bayesian estimation. This section also describes a novel method for drawing the volatilities in the context of a multivariate stochastic volatility process with a common factor. Sections 5 and 6 present our artificial data exercise and our empirical application, respectively. Section 7 summarizes and concludes.
The static representation of a TVP regression model involving a $T$-dimensional dependent variable, $\bm y$, and a $T \times K$-dimensional matrix of predictors, $\bm X$ is:
where $\bm \alpha$ is a $K$-dimensional vector of time-invariant coefficients, $\bm \beta_t$ is a $K \times 1$ vector of time-varying coefficients and $\bm L = \text{diag}(\sigma_1, \dots, \sigma_T)$ with $\sigma_t$ denoting time-varying error volatilities. The TVP part of this model arises through the $\bm W \bm \beta$ term. $\bm W$ is a $T \times k (=TK)$ matrix given by:
with $\bm x_t$ denoting a $K$-dimensional sub-vector of $\bm X$. Equation (ref) 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 ((ref)) 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 $\bm W$. For instance, if we are interested in random walk-type behavior in the TVPs, we can set
This specification implies that $\bm \beta$ can be interpreted as the changes in the parameters and multiplication with $\bm W$ yields the cumulative sum over $\bm \beta$. In our empirical exercise, we consider both of these specifications for $\bm 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, clark2011. 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. korobilis2019high uses message passing techniques to estimate large TVP regressions and shows that these large models outperform a range of competing models. Similarly, hkp2020 approximate the TVPs using message passing techniques based on a rotated model representation and sample from the full conditional posterior of $\bm \alpha$ 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, hauzenberger2022fast propose using the singular value decomposition of $\bm W$ in combination with a conjugate shrinkage prior on $\bm \beta$ 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 $\bm \beta$ is extremely sparse. This is the scalable MCMC method, based on posterior perturbations, of johndrow2017bayes.
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 carriero2019large, hko2020. In particular, we can use the following structural representation of the TVP-VAR:
with $\bm y_{t}$ being an $M$-dimensional vector of endogenous variables, $\bm c_t$ denoting an $M$-dimensional vector of intercepts, $\bm A_{pt}$, for $p = 1, \dots, P$, denoting an $M \times M$-dimensional time-varying coefficient matrix that may be stacked in a matrix $\bm A_t = (\bm A_{1t}, \dots, \bm A_{Pt})$. Furthermore, $ \bm \epsilon_t$ is an $M$-dimensional vector of errors and $\bm \Sigma_t = \text{diag}~(\sigma^2_{1t}, \dots, \sigma^2_{Mt})$ refers to its diagonal time-varying covariance matrix. Finally, $\bm A_{0t}$ defines contemporaneous relationships between the elements of $\bm y_t$ and is lower-triangular with zeros on the diagonal.
The $i^{th}~(i=2, \dots, M)$ equation of $\bm y_t$ can be written as a standard TVP regression model:
Here, $\bm x_{it}$ is a $K_i (= MP+i)$-dimensional vector of covariates with $\bm x_{it} = (1, \{y_{jt}\}_{j=1}^{i-1}, \bm y_{t-1}', \dots, \bm y_{t-P}')'$, $\bm \gamma_{it} = (\bm \alpha_i + \bm \beta_{it}) = (c_{it}, \{a_{ij,0t}\}_{j =1}^{i-1}, \bm A_{i\bullet,t})'$ denotes a $K_i$-dimensional vector of time-varying coefficients, with $c_{it}$ referring to the $i^{th}$ element in $\bm c_t$, $a_{ij,0t}$ denoting the $(i,j)^{th}$ element of $\bm A_{0t}$ and $\bm A_{i\bullet,t}$ referring to the $i^{th}$ row of $\bm A_t$. For $i=1$, $\bm x_{1t} = (1, \bm y'_{t-1}, \dots, \bm y'_{t-p})'$ and $\bm \gamma_{1t} = (c_{1t}, \bm A_{1\bullet,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 $\bm W$, this model nests a variety of commonly used specifications in the literature. For instance, if $\bm W$ implies a random walk behavior of the latent states we arrive at a TVP-VAR closely related to the one proposed in primiceri2005time. 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 cogley2010inflation. This model assumes that the variances of the state innovations evolve according to independent stochastic volatility models.
In this section, we explain the MCMC algorithm of johndrow2017bayes and johndrow2020scalable and discuss how we adapt it for our TVP regression model. The parameters in the static representation are $\bm \alpha$ and $\bm \beta$. Since $\bm \alpha$ is typically of moderate size and potentially non-sparse, we use conventional (exact) MCMC methods for it. It is $\bm \beta$ which is high-dimensional and potentially sparse, characteristics the algorithm of johndrow2017bayes is perfectly suited for. Thus, we use this algorithm for $\bm \beta$. Every model used in the empirical application also includes stochastic volatility.
In the following section, we develop an MCMC algorithm to produce draws of $\bm L$. Since there is nothing new in our MCMC algorithm for $\bm \alpha$ and our algorithm for drawing $\bm L$ is discussed later, in this section we will proceed conditionally on them and work with the transformed regression involving dependent variable $\tilde{\bm y} = \bm L^{-1} (\bm y - \bm X \bm \alpha)$ and explanatory variables $\tilde{\bm W} = \bm L^{-1} \bm W$. The appendix provides full details of our MCMC algorithm. In this section, we will also assume that the prior on $\bm \beta$ is (conditional on other parameters) Gaussian with mean zero and a diagonal prior covariance matrix $\bm D_0 = \text{diag}(d_1, \dots, d_{k})$. Many different global-local shrinkage priors have this general form and, in the following section, we will suggest several different choices likely to be well-suited to TVP regressions.
The exact MCMC algorithm of bhattacharya2016fast for drawing $\bm \beta$ proceeds as follows:
bhattacharya2016fast show that this algorithm is fast compared to existing approaches which involve taking the Cholesky factorization of the posterior covariance matrix. However, it can still be slow when $k$ is very large. The computational bottleneck lies in the calculation of $\bm \Gamma = \tilde{\bm W} \bm D_0 \tilde{\bm W}'$ which has computational complexity of order $\mathcal{O}(T^2k)$. In macroeconomic or financial applications involving hundreds of observations, $T^2k = T^3K$ can be enormous.
johndrow2017bayes and johndrow2020scalable propose an approximation to the algorithm of bhattacharya2016fast which, in sparse contexts, will be much faster and, thus, scalable to huge dimensions. The basic idea of the algorithm is to approximate the high-dimensional matrix $\bm \Gamma$ by dropping irrelevant columns of $\tilde{\bm W}$ so as to speed up computation. To be precise, Steps 4 and 5 of the algorithm are replaced with
Here, $\tilde{\bm W}_S$ denotes a $T \times s$-dimensional sub-matrix of $\tilde{\bm W}$ that consists of columns defined by a set $S$ and $\bm D_{0, S}$ is constructed by taking the diagonal elements of $\bm D_0$ also defined by $S$. Let $S = \{j : {\delta_j} =1\}$ denote an index set with $\delta_j$ being the $j^{th}$ element of a $k$-dimensional selection vector $\bm \delta$ with elements $\delta_j = 1$ with probability $p_j$ and $\delta_j = 0$ with probability $(1-p_j)$. johndrow2017bayes approximates $\delta_j$ by setting $\hat{\delta}_j = 0$ if $d_j \in (0, \xi]$ for $\xi$ being a small threshold. Computational complexity is reduced from $\mathcal{O}(T^2k)$ to $\mathcal{O}(T^2 s)$, where $s = \sum_{j=1}^{k} \delta_j$ is the cardinality of the set $S$ or equivalently the number of non-zero parameters in $\bm \beta$. Step 5* yields a draw from the approximate posterior $\hat{p}(\bm \beta | \bullet)$ with the $\bullet$ notation indicating that we condition on the data and the remaining parameters in the model.
The algorithm requires a choice of a threshold for constructing $\bm \delta$. johndrow2017bayes suggest simple thresholding rules that seem to work well in their work with artificial data (e.g., recommendations include setting the threshold to $0.01$ when explanatory variables are largely uncorrelated, but $10^{-4}$ when they are more highly correlated). However, choosing the threshold might be problematic for real data applications and can require a significant amount of tuning in practice. Instead we propose to choose the thresholds in a different way using SAVS.
To explain what SAVS is and how we use it in practice, note first that papers such as hahncarvalho2015dss recommend separating out shrinkage (i.e., use of a Bayesian prior to shrink coefficients towards zero) and sparsification (i.e., setting the coefficents on de-selected variables to be precisely zero so as to remove them from the model) into different steps. First, MCMC output from a standard model (e.g., a regression with global-local shrinkage prior) is produced. Secondly, this MCMC output is then sparsified by choosing a sparse coefficient vector that minimizes the distance between the predictive distribution of the shrunk model and the predictive density of a model based on this sparse coefficient vector plus an additional penalty term for non-zero coefficients. This assumption is critically based on assuming normally distributed shocks. The optimal solution, $\tilde{\bm \beta}$, is then a sparse vector which can be used to construct $\bm \delta$.
The advantages of this shrink-then-sparsify approach are discussed in hahncarvalho2015dss and, in the context of TVP regressions, in hko2020. One important advantage is that estimation error is removed for the sparsified coefficients. When using global shrinkage priors in high dimensional contexts with huge numbers of parameters, small amounts of estimation error can build up and have a deleterious impact on forecasts. By sparsifying, estimation error in the small coefficients is eliminated, thus improving forecasts. This paper differs from the aforementioned papers by using SAVS to approximate the indicators $\bm \delta$ which is then used in our approximate MCMC algorithm.
The SAVS algorithm, developed in bhattacharya2018signal, is a fast method for solving the optimization problem outlined above, making it feasible to sparsify each draw from the posterior of $\bm \beta$. In the present context, our contention is that a strategy which uses SAVS to shrink-then-sparsify our coefficients can be used to provide a sensible estimate of $\bm \delta$ that does not lead to a deterioration in forecast accuracy. Using SAVS, we first produce a sparsified draws $\tilde{\bm \beta}$.\footnote{Precise details for how SAVS works in TVP regressions, along with additional motivation for the approach, are provided in hko2020.} For each draw $\tilde{\bm \beta } = (\tilde{\beta_1}, \dots, \tilde{\beta_{k}})'$, we then set
Each draw of $\hat{\delta}_j$ is used in the construction of $\hat{\bm \Gamma}$ in the MCMC algorithm of johndrow2017bayes described above. We will refer to this algorithm as being approximate to distinguish it from the exact algorithm of bhattacharya2018signal.
For the time-invariant coefficients, $\bm \alpha$, we use a horseshoe shrinkage prior carvalho2010horseshoe. 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 $\bm \beta $. We modify methods outlined in kowal2019dynamic 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 $\bm \beta_t = (\beta_{1t}, \dots, \beta_{Kt})'$ denotes the coefficients at time $t$, $\tau$ denotes a global shrinkage parameter that pushes all elements in $\bm \beta$ towards zero, $\lambda_t$ is a time-specific shrinkage factor that pushes all elements in $\bm \beta_t$ towards zero and $\phi_{jt}$ is a coefficient and time-specific shrinkage term that follows a half-Cauchy distribution.
Thus, the prior covariance matrix of $\bm \beta_t$ is given by:
which implies that $\lambda_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 $\lambda_t$, it is worth summarizing the key innovation of this prior. As opposed to the dynamic horseshoe of kowal2019dynamic, we only introduce persistence in the common shrinkage factor $\lambda_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(\tau \lambda_t)$ and assuming it follows an AR(1) process:
with $\mu = \log \tau$. We consider two possible distributions for $\nu_t $. In the first of these it follows a four parameter $Z$-distribution, $\mathcal{Z}(1/2, 1/2, 0, 0)$, leading to a variant of the dynamic horseshoe prior proposed in kowal2019dynamic (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 cogley2010inflation 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 $\lambda_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 $\kappa_0, \kappa_1$ denote prior variances with the property that $\kappa_1 \gg \kappa_0$. For the Markov switching model, we assume that $d_t$ is driven by a $(2\times2)$-dimensional transition probability matrix $P$ with transition probabilities from state $i$ to $j$ denoted by $p_{ij}$ (with $p_{ii} \sim \mathcal{B}(a_{i,MS}, b_{i,MS})$, for $i = 0,1$, following a Beta distribution a priori). The mixture model assumes that $p(d_t = 1) = \underline{p}$, with $\underline{p} \sim \mathcal{B}(a_{Mix}, b_{Mix})$. In the empirical application we specify $\kappa_1 = 100/K$, $\kappa_0 = 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 $\lambda_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 $\lambda_t$ to evolve according to an AR(1) process), we use a half-Cauchy prior on $\sqrt{\tau} \sim \mathcal{C}^+(0, 1).$
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 $\lambda_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 johndrow2017bayes, can be used. The only modification is that we construct $\bm D_0$ as follows:
with $\lambda_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 $\lambda_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 $\lambda_t$ under the assumption that it evolves according to an AR(1) process, the algorithm proposed in JPR1995 can be used. However, since this algorithm simulates the $\lambda_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 $\bm \beta_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 $\hat{\bm \beta}_t$ be a $K$-dimensional vector of normalized TVPs with typical element $\hat{\beta}_{jt}= {\beta}_{jt}/(\phi_{jt}\tau^{1/2})$. Using ((ref)) and squaring yields:
with $\nu_t = \bm v'_t \bm v_t$ for $\bm v_t \sim \mathcal{N}(\bm 0_K, \bm I_K)$. Notice that $\nu_t$ follows a $\chi^2$ distribution with $K$ degrees of freedom, denoted by $\chi^2_K$. This implies that sampling algorithms that rely on the Gaussian mixture approximation proposed in kim1998stochastic cannot be used. Instead we approximate the $\chi^2_K$ using a well-known limit theorem that implies, as $K \to \infty$,
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 ((ref)) and taking logs yields:
Finally, under the assumption that $(\sqrt{2K}q_t + K)>0$ and by using a Taylor series expansion,\footnote{More precisely, we compute the mean and variance of $\log \hat{\nu_t}$ using a second and first order Taylor series expansion of $\text{E}(\log(K+ \hat{\nu_t}-K))$ and $\text{Var}(\log(K+ \hat{\nu_t}-K))$ around $K$, respectively.} we approximate $\log \hat{v}_t$ with a $\mathcal{N}\left(\log (K) - 1/K, 2/K\right)$ to render ((ref)) 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 \lambda_t$ using the precision sampler outlined, for example, in chan2009efficient and MCCAUSLAND2011199.
The accuracy of this approximation for different values of $K$ is illustrated in Figure (ref). 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 \chi^2_K$ 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 $\lambda_t$, even for small values of $K$.
In this section we illustrate the merits of our approach using synthetic data.
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 ((ref)) 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:
The remaining parameters are set as follows: $\bm \beta_0 = \bm 0$, $\bm L = 0.01 \times \bm I_T$, $\bm \alpha \sim \mathcal{N}(\bm 0, \bm I_K)$ and $\bm X_j \sim \mathcal{N}(\bm 0, \bm I_{T})$ for $j=1,\dots, K$. Based on these, we use the true path of the parameters $\bm \beta_t$ to obtain a realization of $y_t$. In all simulation experiments and for all models considered we simulate $2,500$ draws from the joint posterior of the parameters and latent states and discard the first $500$ draws as burn-in.
We investigate the accuracy of our scalable approximate MCMC methods relative to the exact MCMC algorithm of bhattacharya2016fast (i.e., it is the version of our algorithm which imposes $\delta_j = 1$ for all $j$). Table (ref) shows the ratio of mean absolute errors (MAEs), computed using the posterior mean of $\{\bm \beta_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 likelihoood does not influence the prior and more shrinkage can be achieved.
Thus, Table (ref) 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.
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 $\bm W$. The random numbers refered to below are drawn from the standard Gaussian distribution.
For $K = 1, ...,400$ and $T \in \{100, 200\}$ we randomly draw a $\bm y$ and an $\bm X$. The $\bm W$ is drawn in two ways which correspond to the flexible and random walk specifications of equations ((ref)) and ((ref)), respectively.
In terms of sparsity, we consider four scenarios based on how we choose $\tilde{\bm W}_S$:
Figure (ref) depicts the computational advantages of our approximate MCMC algorithm relative to the exact algorithm of bhattacharya2016fast. It shows the time necessary to obtain a draw of $\bm \beta$. 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.\footnote{The relatively good performance of the exact algorithm in this case is partly due to the fact that we are coding using sparse algorithms. In the flexible specification for $\bm W$, the underlying matrices are block-diagonal and thus exact sampling is already quite fast.}
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 nelson1987parsimonious, diebold2006macroeconomy and assumes that the yield on a security with maturity $\mathfrak{t}$, labeled $r_t(\mathfrak{t})$, features a factor structure:
Here, $L_{t}, S_{t}$ and $C_{t}$ refer to the level, slope and curvature factor, respectively, while $\eta_t(\mathfrak t)$ denotes maturity-specific measurement errors which are independent across maturities and feature variance $\sigma^2_\eta(\mathfrak t)$. $\zeta$ denotes a parameter that controls the shape of the factor loadings. Following diebold2006macroeconomy, we set $\zeta = 0.7308~(12 \times 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 $\bm y_{t} = (L_{t}, S_{t}, C_{t})'$ and estimating the TVP-VAR defined in ((ref)). We use the flexible specification for $\bm W$ in ((ref)) and the approximate algorithm to estimate the model. In addition, we set the lag length to two. After obtaining forecasts for $\bm y_t$, we use ((ref)) to map the factors back to the observed yields. It is worth noting that ((ref)) 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.
To provide some information on the amount of time variation, Figures (ref) and (ref) depicts heatmaps of the posterior inclusion probability (PIPs) for a Nelson-Siegel model with panels a) to d) referring to the four different dynamic priors for $\lambda_t$. These PIPs are the posterior means of the elements of $\bm \delta$.
The main impression provided by Figure (ref) and (ref) 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.
The dataset covers the entire yield curve and includes yields from one-year to thirty-year bonds in one-year steps. We choose $\{1$y, $3$y, $5$y, $7$y, $10$y, $15$y, $30$y$\}$ 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 $\bm W$), we consider the flexible (FLEX) and random walk (RW) specifications defined in ((ref)) and ((ref)). 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 (ref). In addition, we consider as a competitor the conventional TVP-VAR setup of primiceri2005time. We estimate the TVP-VAR only for the Nelson-Siegel model since the original prior overfits in higher dimensions.\footnote{The priors of the conventional primiceri2005time TVP-VAR are informed by OLS estimates using an initial training sample (in our case the initial first $18$ observations). Such an empirical Bayesian calibration strategy is only sensible for models that feature a small number of endogenous variables.}
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 $\bm \beta = \bm 0$ and then using the sampling steps for $\bm \alpha$ 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 gneiting2007strictly 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 (ref) shows the one-month ahead forecast performance of the different models while Table (ref) 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 gneiting2007strictly.
Overall, the evidence in Table (ref) (and Table (ref)) 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 fischer2023general 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 (ref) 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 $\bm 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 $\lambda_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 primiceri2005time 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. {\tiny
}
VARs modelled with many macroeconomic and financial data sets exhibit parameter change and structural breaks. Typically, most parameter change is found in the error covariance matrix. But there can be small amounts of time-variation in VAR coefficients where only some coefficients change and even they only change at points in time. The problem is how to uncover TVPs of this sort. Simply working with a model where all VAR coefficients change can lead to over-fitting and poor forecast performance. In light of this situation, one contribution of this paper lies in our development of several dynamic horseshoe priors which are designed for picking up the kind of parameter change that often occurs in practice. In an application involving eurozone yield data our methods find small amounts of time variation in parameters. In a forecasting exercise we find that appropropriately modeling this time variation leads to forecast improvements.
The second contribution of this paper lies in computation. The approximate MCMC algorithm developed in this paper is scalable in a manner that exact MCMC algorithms are not. Thus, we have developed an algorithm which can be used in the huge dimensional models that are increasingly being used by economists. Finally, we have developed an MCMC algorithm for common stochastic volatility specifications which is particularly well-suited for large $k$ applications such as the one considered in this paper.
{\setstretch{0.85} \addcontentsline{toc}{section}{References}