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.
99,993 characters · 17 sections · 83 citation commands
Sparse Bayesian State-Space and Time-Varying Parameter Models
Time-varying parameter (TVP) models and, more generally, state space models are widely used in time series analysis to deal with model coefficients that change over time. This ability to capture gradual changes is one of state space models greatest advantages. The flipside of this high degree of flexibility{, however,} is that they run the risk of overfitting with a growing number of coefficients, as many of them might, in reality, be constant over the entire observation period. This will be exemplified in the present chapter with an economic application. We will model US inflation through a TVP Phillips curve, where, out of 18 potentially time-varying coefficients, only a single one actually changes over time. We will show that allowing static coefficients to be time-varying leads to a considerable loss of statistical efficiency, both in uncertainty quantification for the parameters and forecasting future time series observations. We will also show that substantial statistical efficiency can be gained by applying a Bayesian estimation strategy that is able to single out parameters that are indeed constant or even insignificant.
Identifying constant coefficients in a TVP model amounts to a {\em variance selection} problem, involving a decision on whether the variances of the shocks driving the dynamics of a time-varying parameter are equal to zero. Variance selection in latent variable models is known to be a non-regular problem within the framework of classical statistical hypothesis testing har:for. The introduction of shrinkage priors for the variances of a TVP model within a Bayesian framework has proven to be a very useful strategy which is capable of automatically reducing time-varying coefficients to static ones if the model overfits.
In pioneering work, fru-wag:sto reformulated the {\em variance selection} problem for state space models as a {\em variable selection} problem in the so-called non-centered parametrization of the TVP model. This insight established a general strategy for extending shrinkage priors from standard regression analysis to this more general framework. For variance selection in \lq\lq sparse\rq\rq\ state space and TVP models, fru-wag:sto employed discrete spike-and-slab priors, bel-etal:hie relied on the Bayesian Lasso prior, bit-fru:ach applied the normal-gamma prior of gri-bro:inf and cad-etal:tri introduced the triple gamma prior, which is related to the normal-gamma-gamma prior gri-bro:hie and contains the horseshoe prior car-etal:hor as a special case.
The present chapter reviews this literature, starting in Section (ref) with univariate time-varying parameter models. In particular, we will demonstrate that the commonly used inverse gamma prior on the process variances prevents variance selection. Using a ridge prior in the non-centered TVP model instead of an inverse gamma prior provides a simple, yet useful alternative. The ridge prior can be translated into a gamma prior for the variances and leads to more reliable uncertainty quantification in parameter estimation and forecasting for sparse state space models. Starting from the ridge prior, continuous shrinkage priors for variance selection are discussed in Section (ref), whereas Section (ref) discusses discrete spike-and-slab priors. In both sections, we also review strategies for efficient Markov chain Monte Carlo (MCMC) estimation, which is even more challenging for state space models than for standard regression models. Section (ref) discusses extensions to multivariate time series, including TVP Bayesian vector autoregressive models and TVP Cholesky stochastic volatility models, shows how to compare various shrinkage priors through log predictive density scores and addresses the issues of classifying coefficients into dynamic or constant ones. Section (ref) concludes with a brief discussion.
In this section, we consider time-varying parameter (TVP) models for a univariate time series $y_t$. For $t = 1, \ldots, T$, we have that
where $\boldsymbol{\beta}_{t}=(\beta_{1t}, \ldots, \beta_{pt}) ^{\top}$ is a latent state variable and the covariance ${\mathbf{Q}}=\mbox{\rm Diag}\left(\theta_1, \ldots, \theta_p\right)$ of the innovations ${\mathbf{w}}_{t}$ is diagonal. $ {\mathbf x}_t = (x_{1 t}, \ldots, x_{p t})$ is a $p$-dimensional row vector containing the explanatory variables at time $t$. The variables $x_{j t}$ can be exogenous {(i.e. determined outside the model)} control variables and/or be equal to lagged values of $y_{t}$. Usually, one of the variables, say $x_{1 t}$, corresponds to the intercept, but an intercept need not be present. In Section (ref), this approach is extended to multivariate time series ${\mathbf y}_t$.
To fully specify the model, a distribution has to be defined for the initial value $\boldsymbol{\beta}_{0}$ of the state process, with a typical choice being a normal distribution, e.g. $\boldsymbol{\beta}_{0} \sim \mathcal{N} _{p}\left(\boldsymbol{\beta}, {\mathbf{Q}}\right)$, with initial expectation $\boldsymbol{\beta} = (\beta_1, \ldots, \beta_p)^\top$. An alternative choice is to assume a diffuse prior with fixed initial expectation and a very uninformative prior covariance matrix, e.g. $\boldsymbol{\beta}_{0} \sim \mathcal{N} _{p}\left({\mathbf{0}}, 10^5 \cdot{I}_{p} \right)$ where ${I}_{p}$ is the $p$-dimensional identity matrix. However, such a choice is not recommended for TVP models where overfitting presents a concern.
The goal is to recover the unobserved state process $\boldsymbol{\beta}_{0}, \ldots, \boldsymbol{\beta}_{T}$ given the observed time series ${\mathbf y}=(y_1, \ldots,y_T)$. If $\boldsymbol{\beta}$, ${\mathbf{Q}}$ and $\sigma^2$ were known, this is easily achieved by the famous Kalman filter and smoother kal:new. For illustration, a time series $y_t$ is generated from model ((ref)) with $T=200$, $p=3$, $x_{1 t}= 1$, $x_{j t} \sim \mathcal{N}\left(0,1\right)$, $j=2,3$, $\sigma^2=1$, $(\beta_1, \beta_2,\beta_3)=(1,-0.5,0)$ and $(\theta_1, \theta_2,\theta_3)=(0.02,0,0)$. The paths of the hidden process $ \boldsymbol{\beta}_{t}$ are reconstructed using the Kalman filter and smoother based on the true values of $\boldsymbol{\beta}$, $\theta_1$ and $\sigma^2$ and very small values for $\theta_2=\theta_3=10^{-6}$ and compared to the true paths in the left-hand side of Figure (ref). Since the marginal posterior of $\boldsymbol{\beta}_{t}|{\mathbf y}$ is a Gaussian distribution for each $t$, point-wise credible regions for $\boldsymbol{\beta}_{t}$ are easily obtained which are very helpful for uncertainty quantification. Although the TVP model used for estimation overfits, the Kalman smoother is rather accurate in recovering the true paths and clearly indicates that the last coefficients are constant, {\em assuming} that $\theta_2$ and $\theta_3$ are very close to 0.
However, in real-world applications, the variances $\theta_j$ are unknown and estimated from the observed time series, together with the entire path $\mathbf{z}=(\boldsymbol{\beta}_{0}, \ldots, \boldsymbol{\beta}_{T})$. As evident from the Kalman filter, the variances ${\mathbf{Q}}$ of the innovations ${\mathbf{w}}_t $ play an important role in quantifying the loss from propagating the filtering density $ \boldsymbol{\beta}_{t-1}| {\mathbf y}^{t-1} \sim \mathcal{N} _{p}\left({\mathbf m}_{t-1|t-1}, \mathbf{P}_{t-1|t-1}\right),$ given ${\mathbf y}^{t-1}=(y_1, \ldots, y_{t-1})$, into the future to forecast $\boldsymbol{\beta}_{t}$: $$ \boldsymbol{\beta}_{t}|{\mathbf y}^{t-1} \sim \mathcal{N} _{p}\left({\mathbf m}_{t-1|t-1}, \mathbf{P}_{t-1|t-1}+{\mathbf{Q}}\right).$$ A comparably minor change of ${\mathbf{Q}}$ can have a strong effect on uncertainty quantification. For instance, assuming $\theta_2=\theta_3=0.001$ (instead of $10^{-6}$) for the simulated data has a huge effect on the recovered paths, as shown in the right-hand side of Figure (ref). Not only are the credible intervals much broader, we can also no longer be sure if the two coefficients $\beta_{2t}$ and $\beta_{3t}$ are time-varying or constant.
In a maximum likelihood framework, the Kalman filter is used to compute the likelihood function, which is maximized to obtain estimates of $\theta_1, \ldots, \theta_p$, $\sigma^2,$ and $\beta_1, \ldots, \beta_p$ (if the initial means are unknown). Reconstructing $\mathbf{z}=(\boldsymbol{\beta}_{0}, \ldots, \boldsymbol{\beta}_{T})$ then operates conditional on these estimates, see e.g. har:for.
For Bayesian inference, priors are chosen for $\theta_1, \ldots, \theta_p$, $\sigma^2$, and $\beta_1, \ldots, \beta_p$. Given time series observations ${\mathbf y}=(y_1, \ldots, y_T)$, the joint posterior distribution $p(\mathbf{z},\boldsymbol{\beta}, {\mathbf{Q}},\sigma^2|{\mathbf y})$ is the object of interest from which marginal posteriors $p(\boldsymbol{\beta}_{t}|{\mathbf y})$ are derived for each $t$. These can be used for uncertainty quantification as in Figure (ref), while also taking uncertainty in the model parameters into account. Different algorithms have been developed to sample from the joint posterior $p(\mathbf{z},\boldsymbol{\beta}, {\mathbf{Q}},\sigma^2|{\mathbf y})$, in particular two-block Gibbs samplers that alternate between drawing from $p(\mathbf{z}|\boldsymbol{\beta}, {\mathbf{Q}},\sigma^2,{\mathbf y})$ using forward-filtering, backward-sampling (FFBS) car-koh:ong,fru:dat and drawing from $p(\boldsymbol{\beta}, {\mathbf{Q}},\sigma^2| \mathbf{z},{\mathbf y})$.
Both maximum likelihood (ML) and Bayesian inference work well for TVP models where all state variables $\beta_{jt}$ are dynamic. If one of the variances $\theta_j$ is equal to 0, ML estimation leads to a non-regular testing problem, since the true value lies on the boundary of the parameter space har:for. As opposed to this, Bayesian inference is able to deal with such sparse TVP models and, more generally, sparse state space models. The two main challenges from the Bayesian perspective are the choice of an appropriate prior for the variances $\theta_j$ and computational challenges with regards to efficient MCMC estimation.
A popular prior choice for the process variance $\theta_j$ is the inverse gamma distribution,
which is often applied with very small hyperparameters, e.g. $s_0=S_0=0.001$ pet-etal:dyn. Given the latent process $(\beta_{j0}, \ldots, \beta_{jT})$, this prior is conditionally conjugate in the so-called centered parametrization ((ref)), since the density $p(\beta_{j0}, \ldots, \beta_{jT}|\theta_j)$ is the kernel of an inverse gamma distribution. Hence, prior ((ref)) leads to an inverse gamma posterior distribution $p(\theta_j| \beta_{j0}, \ldots, \beta_{jT})$. However, this prior performs poorly when dealing with a sparse TVP model, it is bounded away from zero, making it incapable of inducing strong shrinkage fru-wag:sto.
The effect of choosing a specific prior becomes more apparent when we rewrite model ((ref)) in the non-centered parametrization introduced in fru-wag:sto:
with initial distribution $ \tilde{\boldsymbol{\beta}}_{0} \sim \mathcal{N} _{p}\left({\mathbf{0}}, {I}_{p}\right) $. A linear transformation connects the two parametrizations:
Evidently, both representations are equivalent, and we can specify a prior either on the variances $\theta_j$ in ((ref)) or on the scale parameters $\sqrt \theta_j$ in ((ref)). Since the conjugate prior for $\sqrt \theta _j$ in the non-centered parametrization ((ref)) is the normal distribution, the scale parameter $ \sqrt \theta _j$ is assumed to be Gaussian:
Here, $\sqrt \theta _ j \in \mathbb{R}$ is allowed to take on both positive and negative values. This implies that $\theta_j =(\sqrt \theta_ j)^2$ follows a re-scaled $\chi^2_1$-distribution.introduced such a shrinkage prior (with fixed scale parameter $\tau$) for the process variance in a univariate TVP model (that is $p=1$), and fru-wag:sto extended this idea to state space models with $p>1$. Alternatively, it can be assumed that the prior scale is independent of $\sigma^2$, i.e.
As shown by mor-etal:var, such a prior has certain advantages compared to ((ref)) and allows for the introduction of stochastic volatility in model ((ref)), see kna-etal:shr_tim and Section (ref).
From the viewpoint of variable selection, prior ((ref)) is a ridge prior in a standard regression model, conditional on the hidden path $\mathbf{z}=(\tilde \boldsymbol{\beta}_{0}, \ldots, \tilde \boldsymbol{\beta}_{T})$. Many variable selection priors have been introduced for standard regression models (albeit with known rather than latent regressors), see bha-etal:las for a recent review. Given the non-centered parametrization ((ref)), any of these priors can be, in principle, applied in the context of sparse TVP and state space models for variance selection. And, indeed, the literature has seen an increasing number of papers following this lead bel-etal:hie,bit-fru:ach,cad-etal:tri,fru:com_eff,fru-wag:sto.
Shrinking $\theta_{j}$ toward the boundary value is achieved by shrinking $ \sqrt \theta _j $ toward 0 (which is an interior point of the parameter space in the non-centered parametrization). For a sparse state space model, prior ((ref)) substitutes the inverse gamma prior ((ref)) with a gamma prior. This change in the prior specification is negligible for truly dynamic models, where the posterior distribution $p(\mathbf{z},\boldsymbol{\beta}, {\mathbf{Q}},\sigma^2|{\mathbf y})$ is fairly robust to prior choices $p(\theta_j)$, but has a considerable effect on uncertainty quantification for the unknown path $\mathbf{z}$ for a sparse state space model. This is illustrated in Figure (ref), where the gamma prior $\theta_j \sim \mathcal{G}\left(0.5,10\right)$ is compared to the inverse gamma prior $\theta_j \sim \mathcal{IG} \left(0.001,0.001\right)$ for the inflation data that will be discussed in detail in Section (ref).
A two-block Gibbs sampler is available to sample the latent variables $\mathbf{z}=(\tilde \boldsymbol{\beta}_{0}, \ldots, \tilde \boldsymbol{\beta}_{T})$ and the model parameters $\boldsymbol{\alpha}=(\beta_1, \ldots, \beta_p , \theta_1, \ldots, \theta_p)$ and $\sigma^2$ in the non-centered parametrization, see Algorithm (ref). In step (b), if the prior scale in ((ref)) depends on $\sigma^2$ and, similarly, $\beta_j | \sigma^2 \sim \mathcal{N}\left(0,\sigma^2 \tau\right)$nditional on $\mathbf{z}$, the non-centered parametrization ((ref)) is a standard Bayesian regression model for $\boldsymbol{\alpha}$ with a conjugate prior.
The ridge prior ((ref)) for $ \sqrt \theta _j$ can be rewritten in the following way,
where $\psi^2_j = 1$ is a fixed scale parameter and $\tau$ controls the global level of shrinkage of $\theta_j$, since $\mbox{\rm E}[\theta_j| \tau] = \tau$. In a sparse state space model, we expect that only a fraction of the coefficients are indeed dynamic, while the remaining coefficients are (nearly) constant. This prior perception should be reflected in the choice of the prior distribution of the unknown variances $\theta_1, \ldots, \theta_p$. In this section, we discuss how to incorporate this information through continuous shrinkage priors. In Section (ref), we discuss mixture priors, also called spike-and-slab priors, in the context of variance selection.
Under the ridge prior ((ref)), $\psi^2_j \sim \delta_ 1$ follows a point mass prior on 1, which does not allow for any local adaptation. Continuous shrinkage priors take the form of global-local shrinkage priors in the sense of pol-sco:shr, where $\psi^2_j$ follows a prior $p(\psi^2_j)$ that encourages many small values, representing coefficients that are nearly constant, while at the same time some of the $\psi^2_j$'s are allowed to take on larger values to represent coefficients that are indeed time-varying.
For univariate sparse state space and TVP models, bel-etal:hie introduced the Bayesian Lasso prior par-cas:bay, where $\psi^2 _j$ follows an exponential distribution:
This prior is extended by bit-fru:ach to the normal{-}gamma prior gri-bro:inf, where the exponential prior for $p(\psi^2_j)$ is generalized to a gamma prior:
For both priors, $\tau$ acts as a global shrinkage parameter in a similar manner as for the ridge prior ((ref)), however each innovation variance $\theta_j$ is mixed over its {\em own} (local) scale parameter $\psi^2 _j$, each of which follows an independent exponential ((ref)) or a gamma distribution ((ref)). Hence, the $\psi^2_j$'s play the role of local (component specific) shrinkage parameters. ((ref)) obviously reduces to the Bayesian Lasso prior for $a^\xi=1$, but encourages more prior shrinkage toward small values and, at the same time, more extreme values than the Bayesian Lasso prior for $a^\xi< 1$.
The normal{-}gamma prior ((ref)) for $\sqrt \theta_j $ can be represented in the following way as a \lq\lq double gamma\rq\rq\ on $\theta_j$ bit-fru:ach:
where $\kappa^2_B = 2/\tau$. cad-etal:tri proposed an extension of the double gamma prior ((ref)) to a triple gamma prior, where another layer is added to the hierarchy:
The main difference to the double gamma prior is that the prior scale of the $\xi^2_j$'s is not identical, as each local parameter $\xi^2_j$ depends on yet another local scale parameter $\kappa_j^2$. A similar prior is applied to the initial expectations $\beta_j$:
cad-etal:tri show that the triple gamma prior ((ref)) can be represented as a global-local shrinkage prior in the sense of pol-sco:loc, with the local shrinkage parameter $\psi^2 _j$ arising from an $\mbox{\rm F}\left(2a^\xi, 2 c^\xi\right)$ distribution:
with global shrinkage parameter $\tau=2/\kappa_B^2$. An interesting special case of the triple gamma is the horseshoe prior car-etal:hor which results for $a^\xi=c^\xi =1/2$, since $\psi^2 _j \sim \mbox{\rm F}\left(1,1\right)$ implies that $ \psi _j \sim t _{1}$.
cad-etal:tri show that many other well-known shrinkage priors introduced in a regression context are special cases of the triple gamma, which itself can be regarded as an application of the normal-gamma-gamma prior gri-bro:hie to variance selection in the non-centered parametrization ((ref)).
Among other representations, the triple gamma prior has a representation as a generalized beta mixture prior introduced by arm-etal:gen_bet for variable selection in regression models:
where $\phi^{\xi}=2 c^\xi /(\kappa_B^2 a^\xi)= \tau c^\xi /a^\xi $ and $\mathcal{TPB}\left(a^\xi, c^\xi, \phi^{\xi}\right)$ is the three-parameter beta distribution.
This relationship makes it possible to investigate the shrinkage profile $p(\rho_j )$ of the triple gamma prior. Figure (ref) contrasts a triple gamma prior with $a^\xi=c^\xi=0.1$ with a few of its special or limiting cases, showing the behaviour around the origin, in the tails, as well as the shrinkage profiles.
The graphical representation of the triple gamma prior in Figure (ref) is based on cad-etal:tri who prove the following closed form expression for the marginal prior $p(\sqrt \theta_j |\phi^{\xi}, a^\xi, c^\xi)$:
where $U \left(a,b,z\right)=\int_0^\infty e^{-zt} t^{a-1}(1+t)^{b-a-1} dt$ is the confluent hyper-geometric function of the second kind.
The parameter $a^\xi$ and $c^\xi$ control, respectively, the behaviour of this shrinkage prior at the origin and in the tails. cad-etal:tri prove that the triple gamma prior has an infinite spike at the origin, if $a^\xi \leq 0.5$, where for $a^\xi < 0.5$ and for small values of $\sqrt \theta_j $:
Hence, the infinite spike is more pronounced, the closer $a^\xi $ is to 0. As $\sqrt \theta_j \rightarrow \infty$, the triple gamma prior has polynomial tails, with the shape parameter $c^\xi$ controlling the tail index:
\paragraph*{Choosing the hyperparameters}
A challenging question is how to choose the parameters $a^\xi$, $ c^\xi$ and $\phi^{\xi}$ of the triple gamma prior in the context of variance selection for TVP models.
In high-dimensional settings it is appealing to have a prior that addresses two major issues: first, high concentration around the origin to favor strong shrinkage of small variances toward zero; second, heavy tails to introduce robustness to large variances and to avoid over-shrinkage. For the triple gamma prior, both issues are addressed through the choice of $a^\xi$ and $ c^\xi$.
$a^\xi$ and $c^\xi$ can be fixed{,} as for the Lasso and the horseshoe prior{,} or estimated from the data under a suitable prior. cad-etal:tri, e.g., assume that
restricting the support of $a^\xi $ and $ c^\xi$ to $(0, 0.5)$, ensuring that the triple gamma prior is more aggressive than the horseshoe prior.
Ideally, one should place a hyperprior distribution on the global shrinkage parameter $\phi^{\xi} $. Such a hierarchical triple gamma prior introduces dependence among the local shrinkage parameters $\xi^2_1, \ldots, \xi^2_p$ in ((ref)) and, consequently, among $\theta_1, \ldots, \theta_p$ in the joint (marginal) prior $p(\theta_1, \ldots, \theta_p)$. Introducing such dependence is desirable in that it allows {the prior} to adapt the degree of variance sparsity in a TVP model to the data at hand. For a triple gamma prior with arbitrary $a^\xi$ and (finite) $c^\xi$, cad-etal:tri assume {the} following prior on $\phi^{\xi} $:
Prior ((ref)) reduces to $\phi^{\xi} | a^\xi \sim \mbox{\rm F}\left(2a^\xi ,2 a^\xi\right)$ for $a^\xi = c^\xi $. Hence, for the horseshoe prior, $\phi^{\xi} \sim \mbox{\rm F}\left(1,1\right) $ and the global shrinkage parameter $\tau = \sqrt{\phi^{\xi}}$ follows a Cauchy prior as in bha-etal:las,car-etal:han. As shown by cad-etal:tri, under this hyperprior, the triple gamma prior exhibits behaviour similar to Bayesian Model Averaging (BMA), with a uniform prior on an appropriately defined model size, see Section (ref).
For infinite $c^\xi$, hierarchical versions of the Lasso and the double gamma prior in TVP models are based on a gamma prior for the global shrinkage parameter $ \kappa_B^2$, $ \kappa_B^2 \sim \mathcal{G}\left(d_1, d_2\right)$ bel-etal:hie,bit-fru:ach. This leads to a heavy-tailed extension of both priors, where each marginal density $p(\sqrt \theta_j |d_1, d_2)$ follows a triple gamma prior with the same parameter $a^\xi$ (being equal to one for the Bayesian Lasso) and tail index $c^\xi =d_1$. In this light, very small values of $d_1$ had to be applied in these papers to ensure heavy tails of $p(\sqrt \theta_j |d_1, d_2)$.
The two-block Gibbs sampler outlined in Section (ref) can be extended to perform MCMC inference for continuous shrinkage priors by exploiting the normal scale mixture representation underlying any global-local shrinkage prior.
Assume, for illustration, that we want to apply a normal{-}gamma prior for the initial expectations $\beta_j$ and a double gamma prior for $\theta_j$:
with fixed global shrinkage parameters $a^\tau$, $\lambda_B^2$, $a^\xi$ and $\kappa_B^ 2$. In this case, we can run a three-block Gibbs sampler to draw (a) the latent state process from $p(\mathbf{z}| \boldsymbol{\alpha}, \sigma^2,{\mathbf y})$, (b) the model parameter $\boldsymbol{\alpha}=(\beta_1, \ldots, \beta_p ,\sqrt \theta_1, \ldots, \sqrt \theta_p)$ from $p(\boldsymbol{\alpha}, \sigma^2| \boldsymbol{\lambda}, \boldsymbol{\xi}, \mathbf{z},{\mathbf y})$ conditional on knowing the local scale parameters $\boldsymbol{\lambda}=(\lambda_1, \ldots, \lambda_p) $ and $\boldsymbol{\xi}=(\xi_1^2, \ldots, \xi^2_p)$, and (c) the local scale parameters from $p( \lambda_j| \beta_j,\lambda_B^2)$ and $p(\xi^2_j|\theta_j, \kappa_B^ 2)$ for $j=1, \ldots,p$.
Let us consider step (c) in more {detail}, since sampling the local shrinkage parameter from $\xi^2_j| \theta_j, \kappa_B^ 2$ (and similarly from $ \lambda_j| \lambda_B^2,\beta_j$) is less standard. The double gamma prior $\theta_j|\xi_j^2$ in ((ref)) leads to a density for $\xi^2_j$ given $\theta_j$ which is the kernel of an inverse gamma density. In combination with the gamma prior for $\xi^2_j|a^\xi, \kappa_B^2$ also appearing in ((ref)), this leads to a posterior distribution arising from a generalized inverse Gaussian (GIG) distribution: $ \xi^2_j| \theta_j, a^\xi, \kappa_B^2 \sim \mathcal{GIG}\left(a^\xi-1/2,a^\xi \kappa_B^2,\theta_j\right)$.
A very stable generator from the GIG distribution is implemented in the R-package {\tt GIGrvg} hoe-ley:gig.
\paragraph{To center or to non-center?}
In step (a) and (b) of the three--block sampler described above, we have the option to either work with the centered parametrization ((ref)) or the non-centered parametrization ((ref)). Regardless of the parametrization, sampling the state process is straightforward, using either FFBS car-koh:ong,fru:dat or a one-block sampler such as “all without a loop” (AWOL) bit-fru:ach,kas-fru:anc.
In the centered parametrization, the conditional posterior $\theta_j| \beta_{j0}, \ldots, \beta_{jT},\beta_j$ is again a GIG distribution, since the gamma prior for $\theta_j$ in ((ref)) is combined with the density $p( \beta_{j0}, \ldots, \beta_{jT}| \theta_j, \beta_j)$, which is the kernel of an inverse gamma density. However, like many MCMC schemes which alternate between sampling from the full conditionals of the latent states and the model parameters, the resulting sampler suffers from slow convergence and poor mixing if some of the true process variances are small or even zero.
ally normal prior, $\sqrt \theta_j | \xi_j^2 \sim \mathcal{N}\left(0,\xi_j^2\right)$}, we obtain a joint Gaussian prior for $\boldsymbol{\alpha}=(\beta_1, \ldots, \beta_p ,\sqrt \theta_1, \ldots, \sqrt \theta_p)$, where the local shrinkage parameters $\boldsymbol{\lambda}$ and $\boldsymbol{\xi}$ change the prior scale in a dynamic fashion during MCMC sampling. Hence, in the non-centered parametrization ((ref)), conditional on $\boldsymbol{\lambda}, \boldsymbol{\xi}$ and the latent process $\mathbf{z}$, we are dealing with a Bayesian regression model under a non-conjugate analysis and sampling from $p(\boldsymbol{\alpha}, \sigma^2| \boldsymbol{\lambda}, \boldsymbol{\xi}, \mathbf{z},{\mathbf y})$ can be implemented as in Algorithm (ref).
fru:com_eff discusses the relationship between the various parametrizations for a simple TVP model with $p=1$ and the computational efficiency of the resulting MCMC samplers, see also pap-etal:gen. For TVP models with $p>1$, MCMC estimation in the centered parametrization is preferable for all coefficients that are actually time-varying, whereas the non-centered parametrization is preferable for (nearly) constant coefficients. For practical time series analysis, both types of coefficients are likely to be present and choosing a computationally efficient parametrization in advance is not possible.
section{Application to US inflation modelling}
In our application we model quarterly US inflation (1964:Q1 - 2015:Q4) as a generalized Philips curve with time-varying parameters in the spirit of koo-kor:for. This means that inflation at time $t$ is modeled as
where $y_t$ is inflation at time $t$, ${\mathbf x}_{t-1}$ is a set of $p=18$ predictors including an intercept, exogenous variables from the previous time period and $y_{t-1}$ to $y_{t-r}$, a series of lagged observations of inflation. For the application at hand we assume that $r = 3$. The exogenous predictors included are broad and represent many different potential determinants of inflation. Table (ref) offers an overview of the data, the sources used and the transformations applied to achieve (approximate) stationarity. For this application we assume the error variance $\sigma^2_t$ follows a stochastic volatility specification as in Section (ref).
Three different priors are placed on the expected initial values $\beta_1, \ldots, \beta_p$ and on the variances of the innovations $\theta_1, \dots, \theta_p$, namely the ridge prior, as defined in equation (ref), the Lasso prior, as defined in equation (ref), and the triple gamma prior, as defined in equation (ref). In the case of the Lasso prior, the global shrinkage parameters $\lambda_B^2$ and $\kappa_B^2$ are learned from the data under a gamma prior, specifically $\lambda_B^2 \sim \mathcal{G}\left(0.001, 0.001\right)$ and $\kappa_B^2 \sim \mathcal{G}\left(0.001, 0.001\right)$. In the triple gamma case, the hyperparameters are also learned from the data, under the priors defined in equations (ref) and (ref), with hyperparameter values $\alpha_{a^\xi} = \alpha_{a^\tau} = 5$ and $\beta_{a^\xi} = \beta_{a^\tau} = 10$.
Figure (ref) shows how the three prior setups recovered the same states that were already presented in Figure (ref). While all three are noticeably smaller in scale than the states recovered under the inverse gamma prior, they still differ in this regard as a consequence of the degree of shrinkage imposed, with the triple gamma prior imposing the most, followed by the Lasso prior and the ridge prior, in that order. This can be seen in the parameter for the Dow Jones - the median is virtually zero under the triple gamma prior while displaying much more movement under the other two priors. The parameter of the commodity prices index turns out to be significant, but practically constant under the triple Gamma prior, while the two other priors {also} assign considerable posterior mass to negative values. In the case of the parameter for the treasury bill, the most pronounced movement comes from the state estimated under the triple gamma prior, indicating that truly time-varying states are more likely to be picked up in such a sparse environment if the non time-varying parameters are effectively shrunken towards fixed ones.
Another way to examine the effect that various levels of shrinkage have on the inference that follows is to look at the model implied predictions. Figure (ref) plots the posterior predictive density of the three different models and contrasts these with the true levels of inflation. Two things are noteworthy: first, the stronger the shrinkage imposed by the prior, the less closely the median follows the true observation. This can be seen as shrinkage preventing the model from overfitting. Second, the error variance appears to be larger for the models with more shrinkage, as the spurious time variation in some parameters is dampened, leaving more of the variance to be soaked up by the error term. That this is beneficial for prediction can be seen in Section (ref).
A spike-and-slab prior is a finite mixture distribution with two components, where one component (the {\em spike}) has much stronger global shrinkage than the second component (the {\em slab}). Such mixture shrinkage priors were introduced by geo-mcc:var,geo-mcc:app for variable selection for regression models and aim to identify zero and non-zero regression effects. However, they are useful far beyond this problem and allow, for instance, parsimonious covariance modelling for longitudinal data smi-koh:par, covariance selection in random effects models fru-tue:bay and robust random effects estimation fru-wag:bay.
Discrete spike-and-slab priors were introduced in state space modeling by fru-wag:sto to achieve shrinkage of time-varying state variables toward fixed components. In TVP models, such a prior is introduced for the variance $\theta_j$ and reads
$ \theta_j \sim (1-\pi_{\gamma}) \delta_ 0 + \pi_{\gamma} p_{\mbox{\tiny slab}} ( \theta_j)$, with the spike being a point measure at 0 and $p_{\mbox{\tiny slab}} ( \theta_j)$ being the distribution in the slab. fru-wag:sto introduced the following prior for the scale parameter $\sqrt \theta_j $ in the non-centered parametrization ((ref)), with a ridge prior in the slab:
With $\gamma_j$ being a binary indicator that separates the spike from the slab, $\pi_{\gamma}$ controls the prior occurrence of dynamic coefficients:
Again, this prior can be seen as an extension of the ridge prior, this time with a binary local scale parameter $\psi^2_j=\gamma_j$ taking either the value 0 or 1: $ \sqrt \theta_j |\psi^2_j=\gamma_j \sim \mathcal{N}\left(0, \sigma^2 B_{\gamma} \gamma_j\right)$. A discrete spike-and-slab prior is also applied to the initial expectation $\beta_j$:
with a corresponding binary indicator $\delta_j$ to separate the spike from the slab. The dependence of the prior scale on the error variance $\sigma^2$ in both priors $ p(\beta_j| \sigma^2) $ and $p(\sqrt \theta_j | \sigma^2)$ allows sampling the indicators $\gamma_j$ and $\delta_j$ without conditioning on any model parameters, see Section (ref).
For a TVP model, the initial expectation $\beta_j$ is not identified if the parameter is actually time-varying. Therefore it is not possible to discriminate between $\delta_j=0$ and $\delta_j=1$, if $\gamma_j=1$. For this reason, the following conditional prior for $\delta_j$ given $\gamma_j$ is assumed:
which rules out the possibility that $\delta_j=0$, while $\gamma_j=1$. Combining this conditional prior with ((ref)) leads to a joint prior for each pair $(\delta_j,\gamma_j)$ which has three possible realizations:
As opposed to continuous priors, discrete spike-and-slab priors allow explicit classification of the variables in a TVP model, based on $\delta_j$ and $\gamma_j$:
The probabilities given in ((ref)) are the prior probabilities for classifying coefficients into these three categories. Based on this prior, in a fully Bayesian inference, the joint posterior distribution of $p(\boldsymbol{\delta}, \boldsymbol{\gamma} |{\mathbf y})$ of all indicators $\boldsymbol{\delta}= (\delta_1, \ldots, \delta_p)$ and $\boldsymbol{\gamma}= (\gamma_1, \ldots, \gamma_p)$ is derived and can be used for posterior classification, e.g. by deriving the model most often visited, or the median probability model.
\paragraph*{Choosing hyperparameters for discrete spike-and-slab priors}
First, the prior probabilities $ \pi_{\gamma}$ and $ \pi_{\delta}$ to observe a dynamic or a constant parameter, respectively, have to be chosen. As for standard variable selection, the strategy to fix $ \pi_{\gamma}$ and $\pi_{\delta}$ is very informative on the model sizes. The numbers $p_{d}$, $p_{f}$ and $p_{0}$ of dynamic, constant and zero coefficients, respectively, are given by
Hence, apriori, $p_{d}|\pi_{\gamma} \sim \mathcal{B}in\left(p,\pi_{\gamma}\right)$, $p_{f}|\pi_{\delta}, p_{d} \sim \mathcal{B}in\left(p - p_{d},\pi_{\delta}\right)$, while $p_{0}$ given $p_{d}$ and $p_{f}$ is deterministic, $ p_{0}=p - (p_{d}+ p_{f})$.
Alternatively, a hyperprior can be assumed for both probabilities in order to learn the desired degree of sparsity from the data. Such a hierarchial prior allows more adaptation to the required level of sparsity and assumes that the prior probabilities $\pi_{\delta}$ and $\pi_{\gamma}$ are unknown, each following a beta distribution:
Choosing $a_0^\delta=b_0^\delta=1$ and $a_0^\gamma=b_0^\gamma=1$ implies that the prior on $p_{d}$ is uniform on $\{0, \ldots, p \}$, while $p_{f}| p_{d}$ is uniform on $\{0, \ldots, p - p_{d} \}$.
Second, the prior in the slab has to be specified. For a discrete spike-and-slab prior, all $\theta_j$s with $\gamma_j=0$ and all $\beta_j$s with $\delta_j=0$ are switched off in the non-centered model ((ref)). Hence, a prior has to be chosen for the parameter $\boldsymbol{\beta}_{\delta,\gamma}$ collecting all remaining non-zero $\beta_j$s and $\sqrt \theta_j $s. Under a Gaussian slab distribution, such a prior reads
where $k=p_{f}+2 p_{d}$. However, as for variable selection in regression models, the choice of $ \tau$ is influential in a higher-dimensional setting. A certain robustness is achieved by choosing a hierarchial Student-$t$ slab, where
with hyperpriors $ \lambda ^2 \sim \mathcal{G}\left(a^\lambda, a^\lambda\right)$ and $ \kappa ^2 \sim \mathcal{G}\left(a^\kappa, a^\kappa\right)$ with small degrees of freedom, e.g. $a^{\tau}=a^{\xi}=a^{\lambda}=a^\kappa=0.5$.
Alternatively, fru-wag:sto consider a fractional prior which is commonly used in model selection, as it adapts the prior scale automatically in a way that guarantees model consistency oha:fra. For TVP models, fru-wag:sto defined a fractional prior for $\boldsymbol{\beta}_{\delta,\gamma}$ conditional on the latent process $\mathbf{z}$ as $ p(\boldsymbol{\beta}_{\delta,\gamma}|b, \cdot) \propto \displaystyle p({\mathbf y}|\boldsymbol{\beta}_{\delta,\gamma},\sigma^2, \mathbf{z} )^b $. This prior can be interpreted as the posterior of a non-informative prior combined with a small fraction $b$ of the complete data likelihood $p({\mathbf y}|\boldsymbol{\beta}_{\delta,\gamma},\sigma^2, \mathbf{z} )$.tion{Model space MCMC}
MCMC inference under discrete spike-and-slab priors is challenging, since the sampler is operating in a very high-dimensional model space. Each of the $p$ covariates defines three types of coefficients, hence the sampler needs to navigate through $3^p$ possible models. The various steps of model space MCMC are summarized in Algorithm (ref) for the conjugate slab distribution ((ref)).
Naturally, the most challenging part is Step (a). If $p$ is not too large, then Step (a) can be implemented as a full enumeration Gibbs step by computing the marginal likelihood $p( {\mathbf y}| \boldsymbol{\delta},\boldsymbol{\gamma} ,\mathbf{z})$ for all $3^p$ possible combinations of indicators, as illustrated by fru-wag:sto for unobserved component state space models. Note that{,} conditional on the latent process $\mathbf{z}$, $p( {\mathbf y}| \boldsymbol{\delta},\boldsymbol{\gamma} , \mathbf{z})$ is the marginal likelihood of a constrained version of regression model ((ref)) under the conjugate prior ((ref)) and therefore has a simple closed form. To derive the posterior $p(\boldsymbol{\delta},\boldsymbol{\gamma}| \mathbf{z}, {\mathbf y}) \propto p( {\mathbf y}| \boldsymbol{\delta},\boldsymbol{\gamma} , \mathbf{z}) p(\boldsymbol{\delta},\boldsymbol{\gamma}) $, these marginal likelihoods are combined with the prior $p(\boldsymbol{\delta},\boldsymbol{\gamma})$ {for all models,} which is available in closed form even under the hierarchical prior ((ref)).
In cases where such a full enumeration Gibbs step becomes unfeasible because $p$ is simply too large, Step (a) can be implemented as a single move sampler: loop randomly over all pairs of indicators $(\delta_j,\gamma_j), j=1, \ldots, p,$ and propose to move from the current model $s=(\delta_j,\gamma_j)$ to a new model $s^{\rm \tiny new}=(\delta_j ^{\rm \tiny new},\gamma_j ^{\rm \tiny new})$ with probability $q_{s \rightarrow s^{\rm \tiny new}}$. Accept $(\boldsymbol{\delta},\boldsymbol{\gamma})^{\rm \tiny new}$ with probability $\min(1,\alpha)$ where
The art here is to design sensible moves. One strategy is to move with equal probability to one of the two alternative categories. For instance, if currently $\delta_j=\gamma_j=1$ defines a dynamic coeffcient, then propose, respectively, with probability 0.5 to either move to a fixed coeffcient, where $\delta_j^{\rm \tiny new} = 0$ (while $\gamma_j^{\rm \tiny new} =\gamma_j=1$) or to a zero coeffcient, where $\delta_j^{\rm \tiny new}= \gamma_j^{\rm \tiny new} = 0$. In general, moves involving a change from a fixed to a dynamic coefficient are not easily accepted. Given that $\gamma_j=0$, the current latent path $\mathbf{z}_j =(\tilde \beta_{j0}, \ldots, \tilde\beta_{jT}) $ was sampled from the prior $p(\mathbf{z}_j)$ which can be very different from the smoothed posterior $p(\mathbf{z} _j| \gamma^{\rm \tiny new}_j=1, {\mathbf y}) $, in particular if $T$ is large.
Having updated the vector of indicators $(\boldsymbol{\delta},\boldsymbol{\gamma})$, a modified version of Algorithm (ref) is applied in Step (b) and (c) of Algorithm (ref) to sample the unconstrained model parameters $\boldsymbol{\beta}_{\delta,\gamma} ,\sigma^2 | \mathbf{z}, {\mathbf y}$ and $\mathbf{z} | {\mathbf{\boldsymbol{\vartheta}}} , {\mathbf y}$ in the restricted version of the non-centered parametrization. In particular, the sampling order is interchanged to obtain a valid sampler, since $(\boldsymbol{\delta},\boldsymbol{\gamma})$ are updated without conditioning on the parameter ${\mathbf{\boldsymbol{\vartheta}}}=(\beta_1, \ldots, \beta_p , \sqrt\theta_1, \ldots, \sqrt\theta_p, \sigma^2 )$.
The analysis in Section (ref) is extended, using discrete spike-and-slab priors for $\beta_j$ and $\sqrt \theta_j$ with following slab distributions: (1) Gaussian with $\tau=1$, (2) fractional priors with $b=10^{-4}$ and (3) hierarchial Student-$t$ with $a^{\tau}=a^{\xi}=a^{\lambda}=a^\kappa=0.5$. The hierarchical prior $\sigma^2|C_0 \sim \mathcal{IG} \left(0.5,C_0\right)$, $C_0 \sim \mathcal{G}\left(5,10/3\right)$ is assumed for the (homoscedastic) variance $\sigma^2$. The prior of $\pi_{\delta}$ and $\pi_{\gamma}$ is chosen as in ((ref)) with $a_0^\delta=b_0^\delta=a_0^\gamma=1$ and $b_0^\gamma=2$.
Model space MCMC sampling was run for 100.000 iteration after a burn-in of 10.000. Under the hierarchial Student-$t$ slab, the sampler exhibits acceptance rates of around 20% for all classes of moves. This indicates relatively good performance, given that the latent variables $\mathbf{z}$ are unobserved and imputed under the old indicators. For Gaussian and fractional slabs, the average acceptance rate of moves between fixed and dynamic components was less than 5%. To verify convergence, the sampler was run twice, starting either from a full TVP model with all $\gamma_j$s equal to 1 or from a standard regression model with all $\gamma_j$s equal to 0. Under the Student-$t$ slab, we found high concordance between the models sampled by both chains after burn-in. Under Gaussian and fractional slabs, however, the two chains were sampling totally different models, depending on the starting value.
The time-varying parameters recovered under the hierarchical Student-$t$ slab are shown in Figure (ref). We see a similar discrimination between a dynamic path (treasury bills), a constant path (commodity prices index) and a zero path (Dow Jones) as we saw in Figure (ref) under the triple gamma prior. A more formal discrimination based on the sampled indicators $\delta_j$ and $\gamma_j$ will be performed in Section (ref).tion{Extensions }
Assuming a homoscedastic error variance $\sigma^2$ in the observation equation of the TVP model ((ref)) may create spurious time-variation in the coefficients, as discussed by sim:mac. To be robust against conditional heteroscedasticity, $\sigma^2_t$ is often assumed to be time-varying over $t=1,\dots,T$:
For TVP models, it is common to assume a stochastic volatility (SV) specification jac-etal:bay_ana, where the log volatility $h_t = \log \sigma^2_t $ follows an AR(1) process:
The unknown model parameters $\mu$, $\phi$, and $\sigma_\eta^2$ in (ref) and the entire latent volatility process $\{h_0, h_1,\ldots, h_T\}$ are added to the set of unknown variables. MCMC estimation is easily extended using the very efficient sampler developed by kas-fru:anc and implemented in the R-package {\tt stochvol} hos-kas:mod.
The TVP model ((ref)) introduced in Section (ref) for univariate time series can be easily extended to TVP models for multivariate time series. Consider, as illustration, the following TVP model for a $q$-dimensional time series ${\mathbf y}_{t}$,
where ${\mathbf x}_{t}$ is a column vector of $p$ regressors, and ${\mathbf B}_t $ is a time-varying $(q\times p)$ matrix with coefficient $\beta_{ij,t}$ in row $i$ and column $j$, potentially containing structural zeros or constant values.
\paragraph*{Sparse TVP Cholesky SV models}
One example is the sparse TVP Cholesky SV model bit-fru:ach, which reads for $q = 3$:
where the log volatilities $h_{it}$, $i=1, \ldots, q,$ follow $q$ independent SV processes as defined in ((ref)), with row specific parameters $\mu_i$, $\phi_i$, and $\sigma^2_{\eta,i}$. System ((ref)) consists of three independent univariate TVP models, where no intercept is present. In the first row, no regressors are present either and only the log volatility $h_{1t}$ has to be estimated. In the $i$-th equation, $i-1$ regressors are present and $i-1$ time-varying regression coefficients $\beta_{ij,t}$ as well as the time-varying volatility $h_{it}$ need to be estimated. System ((ref)) can be written as
where ${\mathbf B}_t$ is a $q \times q$ matrix with time-varying coefficients $\beta_{ij,t}$, which are 0 for $j\geq i$. ${\mathbf D}_t=\mbox{\rm Diag}\left(\mbox{\rm e}^{h_{1t}},\ldots,\mbox{\rm e}^{h_{q t}}\right)$ is a diagonal matrix and the $q$-dimensional vector ${\mathbf x}_{t}=(y_{1t},\ldots,y_{q t})^\top$ is equal to ${\mathbf y}_t$.
It is possible to show that this system is equivalent to the assumption of a dynamic covariance matrix, $ {\mathbf y}_t\sim \mathcal{N} _{q}\left({\mathbf{0}},\boldsymbol{\Sigma}_t\right)$, where $\boldsymbol{\Sigma}_t = {\mathbf A}_t {\mathbf D}_t {\mathbf A}_t ^{\top}$ and the dynamic Cholesky factor ${\mathbf A}_t$ is lower triangular with ones on the main diagonal and related to ${\mathbf B}_t$ through ${\mathbf A}_t=({I}_{q}- {\mathbf B}_t)^{-1}$.
Both in ((ref)) as well as in the more general system ((ref)), the unconstrained time-varying coefficients $\beta_{ij,t}$ are assumed to follow independent random walks as in the univariate case:
with initial values $\beta_{ij,0} \sim \mathcal{N}\left(\beta_{ij},\theta_{ij}\right) $. Each of the time-varying coefficients $\beta_{ij,t}$ is potentially constant, with the corresponding process variance $\theta_{ij}$ being 0. A constant coefficient $\beta_{ij,t} = \beta_{ij}$ is potentially insignificant, in which case $\beta_{ij}=0$. Hence, as for the univariate case, discrete spike-and-slab priors as introduced in Section (ref) or continuous shrinkage {priors} as introduced in Section (ref) are imposed on the fixed regression coefficients $\beta_{ij}${,} as well as the process variances $\theta_{ij}${. This defines} a sparse multivariate TVP model for identifying which of these scenarios holds for each coefficient $\beta_{ij,t}$.
It is advantageous to introduce (hierarchical) shrinkage priors which are independent row-wise. For instance, bit-fru:ach, introduce a hierarchical double gamma prior for $\theta_{ij}$ and a hierarchical normal{-}gamma prior for $\beta_{ij}$ for each row $i$ of the TVP Cholesky SV model. Alternatively, independent discrete spike-and-slab priors with row-specific inclusion probabilities can be specified. Any of these choices leads to prior independence across the $q$ rows of the system ((ref)) and both model space MCMC as well as boosted MCMC can be applied row-wise to perform posterior inference.
\paragraph*{Sparse TVP-VAR-SV models}
Another important example are time-varying parameter vector autoregressive models of order $r$ with stochastic volatility (TVP-VAR-SV), where the $q$-dimensional time series ${\mathbf y}_{t}$ is assumed to follow
where $\mathbf{c}_t$ is the $q$-dimensional time-varying intercept, $\boldsymbol{\Phi}_{j,t}$, for $j = 1, \ldots, r $ is a $q \times q$ matrix of time-varying coefficients, and $\boldsymbol{\Sigma}_t$ is the time-varying variance covariance matrix of the error term. Since the influential paper of pri:tim, this model has become a benchmark for analyzing relationships between macroeconomic variables that evolve over time, see cha-eis:bay,eis-etal:sto,fel-etal:sop,koo-kor:lar,nak:tim, among many others.
Since all $q$ equations share the same predictor ${\mathbf x}_t = ( 1, {\mathbf y}_{t -1}^\top, \ldots,{\mathbf y}_{t-r}^\top)^\top$ (a vector of length $p=q r+1$), the TVP-VAR-SV model can be written in a compact notation exactly as in ((ref)) with matrix
All coefficients $\beta_{ij,t}$ in ${\mathbf B}_t$ follow independent random walks as in ((ref)) with initial expectation $\beta_{ij}$ and process variance $\theta_{ij}$. Due to the high dimensional nature of the time-varying matrix ${\mathbf B}_t$, shrinkage priors are instrumental for efficient inference, even for moderately sized systems. For instance, cad-etal:tri introduce independent hierarchical triple gamma priors for $\beta_{ij}$ and $\theta_{ij}$ in each row $i=1, \ldots,q $ of the TVP-VAR-SV model and demonstrate considerable efficiency gain compared to other shrinkage priors, such as the Lasso.
Since $\boldsymbol{\Sigma}_{t}$ is typically a full covariance matrix, the rows of the system (ref) are not independent, as the various components in $\boldsymbol{\varepsilon}_{t}$ are correlated. Following fru-tue:bay, cad-etal:tri use the Cholesky decomposition $\boldsymbol{\Sigma}_t = {\mathbf A}_t {\mathbf D}_t {\mathbf A}_t ^\top $ to represent the TVP-VAR-SV model as a triangular system with independent errors ${\boldsymbol{\eta}}_t \sim \mathcal{N} _{q}\left(\bm 0, {\mathbf D}_t\right) $. ${\mathbf A}_t$ is lower triangular with ones on the main diagonal and the unconstrained elements $a_{ij, t}$ in the $i$-th row and $j$-th column of ${\mathbf A}_t$ again follow random walks, with their own set of shrinkages priors on the corresponding variances and initial expectations.
The TVP-VAR-SV model then has a representation as a system of $q$ univariate TVP models, e.g. for $q=3$:
where $\bm \beta_{t}^i$ is the $i$th row of ${\mathbf B}_t$. For $i>1$, the $i$th equation is a univariate TVP model with the residuals $\eta_{1t}, \ldots, \eta_{i-1,t} $ of the preceding $i - 1$ equations serving as explanatory variables. Nevertheless, the time-varying parameters $\bm \beta_{t}^i$ in each row can be estimated equation by equation cad-etal:tri.
It should be noted that both models might be sensitive to the ordering of the variables of the multivariate outcome ${\mathbf y}_t$, see kil-lue:str for a thorough discussion.
While the discussion of this chapter {is centered around} Gaussian time series, all methods can be extended to non-Gaussian time series, as demonstrated in fru-wag:sto{,} who also considered time series of small counts based on the Poisson distribution. The main idea is to augment auxiliary latent variables $\boldsymbol{\omega}$ such that conditional on $\boldsymbol{\omega}$ a Gaussian TVP model results. Variable and variance selection is then performed conditional on $\boldsymbol{\omega}$, while an additional step in the MCMC scheme imputes $\boldsymbol{\omega}$ given the remaining variables.
Examples include the representation of student-$t$ errors as scale mixtures of Gaussians and binary time series, where the representation $d_t=\mathbb{I} (y_t>0)$ leads to the conditionally Gaussian state space model ((ref)). A similar strategy is pursued in fru-wag:bay,wag-dul:bay for non-Gaussian random effects models and in wag:bay_est for dynamic survival models, see also bha-etal:hor_reg_mac for a recent review on regularisation in complex and deep models.
Log predictive density scores ({\rm LPDS} ) are a widely used scoring rule to compare models; see, e.g., gne-raf:str. As shown by bit-fru:ach, log predictive density scores are also a useful means of evaluating and comparing different shrinkage priors for TVP models. It is common in this framework to use the first $t_0$ time series observations ${\mathbf y}^{{\tiny \mbox{\rm tr}}}=({\mathbf y}_{1},\ldots, {\mathbf y}_{t_0})$ as a \lq\lq training sample\rq\rq , while evaluation is performed for the remaining observations ${\mathbf y}_{t_0+1}, \ldots, {\mathbf y}_T$.
For univariate time series $y_t$, ${\mbox{\rm LPDS}}$ is defined as:
For each point in time, ${\mbox{\rm LPDS}}^{\star}_{t}$ analyzes the performance separately for each $y_t$ and is obtained by evaluating the one-step ahead predictive density $ p(y_{t}| {\mathbf y}^{t-1}) $ given observations $ {\mathbf y}^{t-1}=(y_1,\ldots, y_{t-1})$ up to $t-1$ at the {\em observed} value $y_t$. ${\mbox{\rm LPDS}} $ is an aggregated measure of performance for the entire time series. As shown by fru:bay in the context of selecting time-varying and fixed components for a basic structural state space model, ${\mbox{\rm LPDS}} $ can be interpreted as a log marginal likelihood based on the training sample prior $p({\mathbf{\boldsymbol{\vartheta}}}| {\mathbf y}^{{\tiny \mbox{\rm tr}}})$, since
where ${\mathbf{\boldsymbol{\vartheta}}} = (\beta_1, \ldots, \beta_p$, $\sqrt \theta_1 ,\ldots, \sqrt \theta_p , \sigma^2)$ summarises the unknown model parameters. Hence, log predictive density scores provide a coherent foundation for comparing the predictive power of different types of shrinkage priors.
Determining ${\mbox{\rm LPDS}}^{\star}_{t}$ for each $t =t_0+1, \ldots, T$ can be challenging computationally. In bit-fru:ach, a Gaussian mixture approximation, called the conditionally optimal Kalman mixture approximation, is introduced to determine $p(y_{t}|{\mathbf y}^{t-1})$ independently for each $t$, based on $M$ draws ${\mathbf{\boldsymbol{\vartheta}}}^{(m)}, m=1, \ldots, M$ from the posterior distribution $p({\mathbf{\boldsymbol{\vartheta}}}|{\mathbf y}^{t-1})$.
The whole concept can be extended to multivariate time series by defining
In a triangular system such as the TVP Cholesky SV model and the TVP-VAR-SV model discussed in Section (ref), errors are uncorrelated and we can exploit that
Since we condition on {\em observed} values $y_{1t}, \ldots, y_{i-1,t}$ in equation $i$, ${\mbox{\rm LPDS}}^{\star}_{it}$ can be determined independently for each $t$ {and} for each equation $i$. This allows one to fully exploit the computational power of modern parallel computing facilities.
\paragraph*{Application to inflation modelling}
To demonstrate the benefit that shrinkage provides with regards to out-of-sample prediction, we calculate one-step ahead LPDSs for the last 100 time points of the inflation dataset introduced in Section (ref) and compute the cumulative sum. Six different prior choices are considered here: (1) the triple gamma prior, (2) the horseshoe prior, (3) the double gamma prior, (4) the Lasso prior, (5) the ridge prior and, finally, (6) the inverse gamma prior. Figure (ref) displays the results, with higher numbers equating to better out-of-sample prediction. It is immediately obvious that the inverse gamma prior does not appear to be competitive in this regard. While it displays a high degree of in-sample fit (as evidenced by Figure (ref), Section (ref)), the forecasting performance severely lags behind the other prior choices. Similarly, if not quite as drastically, the ridge prior does not forecast as well as the more strongly regularized approaches. The three priors with the most shrinkage, the triple gamma, the horseshoe and the double gamma, all perform comparably, while the Lasso prior initially lags behind, only to gain ground during the subprime mortgage crisis between 2007 and 2009.
An interesting insight of cad-etal:tri is that the triple gamma prior shows behaviour very similar to a discrete spike-and-slab prior as both $a^\xi$ and $c^\xi$ approach zero. This induces BMA-type behaviour on the joint shrinkage profile $p(\rho_1, \ldots, \rho_p)$, with an infinite spike at all corner solutions, where some $\rho_j$ are very close to one, whereas the remaining ones are very close to zero. For illustration, Figure (ref) compares bivariate shrinkage profiles of various continuous shrinkage priors. The BMA-type behaviour of the triple gamma becomes evident through the large amount of mass placed in the four corners, with the overlayed 500 samples from the prior following suit and clustering in those areas.
Following car-etal:han, a natural way to perform variable selection in the continuous shrinkage prior framework is through thresholding. Specifically, when $ (1-\rho_j) > 0.5$, or $\rho_j< 0.5$, the variable is included, otherwise it is not. Notice that thresholding implies a prior on the model dimension $p_{d}$ defined as
where $\rho_j \sim \mathcal{TPB}\left(a^\xi, c^\xi, \phi^\xi\right)$, see ((ref)). The choice of the global shrinkage parameter $\phi^\xi$ strongly impacts the prior on $p_{d}$. For a symmetric triple gamma prior with $a^\xi=c^\xi$ and $\phi^{\xi}=1$ fixed, for instance, $\pi_{\gamma}=0.5$ and we obtain $p_{d} \sim \mathcal{B}in\left(p, 0.5\right)${,} regardless of $a^\xi${. This leads to }similar problems as with fixing $\pi_{\gamma}=0.5$ for a discrete spike-and-slab prior.a hyperprior on $\phi^\xi$ as discussed in Section (ref) is as vital for variance selection through continuous shrinkage prior as making $\pi_{\gamma}$ random is for a discrete spike-and-slab prior. car-etal:han show that the hyperprior for $\phi^\xi$ defined in ((ref)) leads to a uniform prior distribution on the model dimension $p_{d}$, since $\pi_{\gamma} \sim \mathcal{U}\left[0,1\right]$ is uniformly distributed.
\paragraph*{Application to US inflation modelling}
For illustration, we compare discrete spike-and-slab priors and hierarchical continuous shrinkage priors with regard to classification of the time-varying parameters for the inflation data set introduced in Section (ref). The posterior probabilities of each coefficient to be either zero, fixed or dynamic are estimated from the $M$ posterior draws of $(\delta_j ^{(m)},\gamma^{(m)})$:
and $\mbox{\rm P} (\beta_{jt} \mbox{ \small zero} |{\mathbf y} )= 1- \mbox{\rm P} (\beta_{jt} \mbox{ \small dynamic}|{\mathbf y} )- \mbox{\rm P} (\beta_{jt} \mbox{ \small fixed} |{\mathbf y} )$. The indicators $(\delta_j ^{(m)},\gamma^{(m)})$ are an immediate outcome of the model space MCMC sampler for the discrete spike-and-slab prior and are derived for continuous shrinkage priors using thresholding as explained above.
According to this procedure, none of the coefficients is classified other than zero for the Lasso prior, which is not surprising in light of Figure (ref). Somewhat unexpectedly, the same classification results for the triple gamma prior, for which a clear visual distinction can be made in Figure (ref) between the relatively dynamic coefficient of treasury bills and the other two coefficients which are shrunken toward a fixed coefficient.
As opposed to this, the discrete spike-and-slab prior shows more power to discriminate between the different types of coefficients for this specific data set. The corresponding classification probabilities are reported for each coefficient in Table (ref) and match the behaviour of the recovered time-varying coefficients in Figure (ref). More specifically, treasury bills is clearly classified as dynamic, the commodity prices index is classified as a having positive, but fixed effect on inflation, and the Dow Jones is clearly classified as insignificant.
This chapter illustrates the importance of variance selection for TVP models. If the true model underlying a time series is sparse, with many coefficients being constant or even zero, then a full-fledged TVP model might quickly overfit. To avoid loss of statistical efficiency in parameter estimation and forecasting that goes hand-in-hand with the application of an overfitting model, we generally recommend to substitute the popular inverse gamma prior for the process variances by suitable shrinkage priors. As demonstrated in this chapter, shrinkage priors are indeed able to automatically reduce time-varying coefficients to constant or even insignificant ones.
Within the class of continuous shrinkage priors, flexible priors such as hierarchical versions of the double gamma, the triple gamma or the horseshoe prior typically turn out to be preferable to less flexible priors such as the hierarchical Lasso. These priors often show a comparable behaviour in terms of model comparison through log predictive density scores and they beat the inverse gamma prior by far. This was illustrated with an application to US inflation modelling using a TVP Phillips curve.
Discrete spike-and-slab priors are an attractive alternative to continuous shrinkage priors as they allow explicit classification of the time-varying coefficients into dynamic, constant and zero ones. For continuous shrinkage priors, such a classification can be achieved only indirectly through thresholding and the appropriate choice of the truncation level is still an open issue for TVP models. However, convergence problems with model space MCMC algorithms are common with discrete spike-and-slab priors and the sampler might get stuck in different parts of the huge model space, depending on where the algorithm is intialized. In our illustrative application, discrete spike-and-slab priors were more successful in classifying obviously time-varying coefficients than any continuous shrinkage prior, but only in combination with a Student-$t$ slab distribution. For other slab distributions, in particular Gaussian ones, severe convergence problems with trans-dimensional MCMC estimation were encountered.
A key limitation of any of the approaches reviewed in this chapter is that they can only differentiate between parameters that are constantly time-varying or not time-varying at all. One could think of many scenarios in which a parameter may be required to be time-varying over a stretch of time and be constant elsewhere. The design of suitable {\em dynamic shrinkage priors} that are able to handle such a situation is cutting-edge research in the area of state space and TVP models. Very promising approaches toward dynamic shrinkage priors were put forward by a number of authors, including cas-etal:bay_non,kal-gri:tim,kow-etal:dyn,roc-mca:dyn.
A script to replicate select results from this chapter and instructions on how to download software routines in R is made available as part of the online supplement of this edited volume.
\printindex \cleardoublepage