Triple the Gamma—A Unifying Shrinkage Prior for Variance and Variable Selection in Sparse State Space and TVP Models

Abstract

Time-varying parameter (TVP) models are very flexible in capturing gradual changes in the effect of explanatory variables on the outcome variable. However, in particular when the number of explanatory variables is large, there is a known risk of overfitting and poor predictive performance, since the effect of some explanatory variables is constant over time. We propose a new prior for variance shrinkage in TVP models, called triple gamma. The triple gamma prior encompasses a number of priors that have been suggested previously, such as the Bayesian Lasso, the double gamma prior and the Horseshoe prior. We present the desirable properties of such a prior and its relationship to Bayesian Model Averaging for variance selection. The features of the triple gamma prior are then illustrated in the context of time varying parameter vector autoregressive models, both for simulated dataset and for a series of macroeconomics variables in the Euro Area.

1. Introduction

Model selection in a high-dimensional setting is a common challenge in statistical and econometric inference. The introduction of Bayesian model averaging (BMA) techniques in the statistical literature (Brown et al. 2002; Cottet et al. 2008; Raftery et al. 1997) has led to many interesting applications, see, among others, (Frühwirth-Schnatter and Tüchler 2008; Kleijn and van Dijk 2006; Koop and Potter 2004; Sala-i-Martin et al. 2004) for early references in econometrics.

Selecting explanatory variables for possibly very high-dimensional regression problems though shrinkage priors is an attractive alternative to BMA which relies on discrete mixture priors, see Bhadra et al. (2019) for an excellent review. There is a vast and growing literature on shrinkage priors for regression problems that focuses on the following aspects. First, how to choose sensible priors for high-dimensional model selection problems in a Bayesian framework, second, how to design efficient algorithms to cope with the associated computational challenges and third, to investigate, both from a theoretical and a practical viewpoint, how such priors perform in high-dimensional problems.

A striking duality exists in this very active area between Bayesian and traditional approaches. For many shrinkage priors, the mode of the posterior distribution obtained in a Bayesian analysis can be regarded as a point estimate from a regularization approach, see Fahrmeir et al. (2010) and Polson and Scott (2012a). One such example is the popular Lasso (Tibshirani 1996) which is equivalent to a double-exponential shrinkage prior in a Bayesian context (Park and Casella 2008). However, the two approaches differ when it comes to selecting penalty parameters that impact the sparsity of the solution. One advantage of the Bayesian framework in this context is that the penalty parameters are considered to be unknown hyperparameters which can be learned from the data. Such “global-local” shrinkage priors (Polson and Scott 2011) adjust to the overall degree of sparsity that is required in a specific application through a global shrinkage parameter and separate signal from noise through local, individual shrinkage parameters.

While the inclusion of potentially many explanatory variables though shrinkage priors in regression models is addressed in a vast literature, the use of shrinkage priors for more general econometric models in time series analysis, such as state space models and time-varying parameter (TVP) models is, in comparison, less well-studied. Sparsity in the context of such models refers to the presence of a few large variances among many (nearly) zero variances in the latent state processes that drive the observed time series data. A common goal in this setting is to recover a few dynamic states, driven by such a state space model, among many (nearly) constant coefficients. As shown by Frühwirth-Schnatter and Wagner (2010), this variance selection problem can be cast into a variable selection problem in the non-centered parametrization of a state space model. Once this link has been established, shrinkage priors that are known to perform well in high-dimensional regression problems can be applied to variance selection in state space models, as demonstrated for the Lasso (Belmonte et al. 2014) and the normal-gamma (Bitto and Frühwirth-Schnatter 2019; Griffin and Brown 2017).

Despite this already existing variety, we introduce a new shrinkage prior for variance selection in sparse state space and TVP models in the present paper called triple gamma prior, as it has a representation involving three gamma distributions. This prior can be related to various shrinkage priors that were found to be useful for high-dimensional regression problems, such as the generalized beta mixture prior (Armagan et al. 2011), and contains the popular Horseshoe prior (Carvalho et al. 2009, 2010) as a special case. Furthermore, the half-t and the half Cauchy (Gelman 2006; Polson and Scott 2012b), suggested as robust alternatives to the inverse gamma distribution for variance parameters in hierarchical models, as well as the Lasso and the double gamma, are special cases of the triple gamma. In this context, the triple gamma can also be regarded as an extension of the scaled beta2 distribution (Pérez et al. 2017).

Among Bayesian shrinkage priors, usually a clear distinction is made between two-group mixture or spike-and-slab priors and continuous shrinkage priors, of which the triple gamma is a special case. An important contribution of the present paper is to show that the triple gamma provides a bridge between these two approaches and has the following property which is favourable both in sparse and dense situations. One of the hyperparameters allows high concentration over the region in the shrinkage profile that is relevant for shrinking noise, while the other hyperparameter allows high concentration over the region that prevents overshrinking of signals. This allows the triple gamma prior to exhibit behavior that very much resembles Bayesian model averaging based on discrete spike-and-slab priors, with a strong prior concentration at the corner solutions where some of the variances are nearly zero. While this is reminiscent of the Horseshoe prior, the shrinkage profile induced by the triple gamma is more flexible than that of a Horseshoe. Thanks to the estimation of the hyperparemters, it is not constrained to be symmetric around one half, enabling adaption to varying degrees of sparsity in the data.

The triple gamma prior also scores well from a computational perspective. While exploring the full posterior distribution for spike-and-slab priors leads to computational challenges due to the combinatorial complexity of the model space, Bayesian inference based on Markov chain Monte Carlo (MCMC) methods is straightforward for continuous shrinkage priors, exploiting their Gaussian-scale mixture representation (Bitto and Frühwirth-Schnatter 2019; Makalic and Schmidt 2016). An extension of these schemes to the triple gamma prior is fairly straightforward.

We will study the empirical performance of the triple gamma for a challenging setting in econometric time series analysis, namely for time-varying parameter vector autoregressive models with stochastic volatility (TVP-VAR-SV models). Since the influential paper of Primiceri (2005) (see Del Negro and Primiceri (2015) for a corrigendum), this model has become a benchmark for analyzing relationships between macroeconomic variables that evolve over time, see Nakajima (2011), Koop and Korobilis (2013), Eisenstat et al. (2014), Chan and Eisenstat (2016), Feldkircher et al. (2017) and Carriero et al. (2019), among many others. Due to the high dimensionality of the time-varying parameters, even for moderately sized systems, shrinkage priors such as the triple gamma prior are instrumental for efficient inference.

The rest of the paper is organized as follows—in Section 2, we define the triple gamma prior and discuss some of its properties. The close relationship between the triple gamma and spike-and-slab priors applied in a BMA context is investigated in Section 3.2. Section 4 introduces an efficient MCMC scheme and Section 5 provides applications to TVP-VAR-SV models. Section 6 concludes the paper.

2. The Triple Gamma as a Prior for Variance Parameters

2.1. Motivation and Definition

To motivate the triple gamma prior, consider the state space form of a TVP model for a univariate time series $y_t$. For $t=1,\dots ,T$, we have that

$$\begin{array}{cc}\hfill & {\mathit{\beta }}_t={\mathit{\beta }}_{t-1}+{\mathbf{w}}_t,{\mathbf{w}}_t\sim {\mathcal{N}}_d\left(\mathbf{0},\mathbf{Q}\right),\hfill \\ \hfill & y_t={\mathbf{x}}_t{\mathit{\beta }}_t+{\epsilon }_t,{\epsilon }_t\sim \mathcal{N}\left(0,{\sigma }_t^2\right),\hfill \end{array}$$

(1)

where $\mathbf{Q}=\mathrm{Diag}\left({\theta }_1,\dots ,{\theta }_d\right)$ and the initial value of the state process follows a normal distribution, ${\mathit{\beta }}_0\sim {\mathcal{N}}_d\left(\mathit{\beta },\mathbf{Q}\right)$, with initial mean $\mathit{\beta }={({\beta }_1,\dots ,{\beta }_d)}^{\top }$. ${\mathbf{x}}_t=(x_{t1},\dots ,x_{td})$ is a d-dimensional row vector containing the explanatory variables at time t. The variables $x_{tj}$ can be exogenous control variables and/or be equal to lagged values of $y_t$. Usually, one of the variables, say $x_{t1}$, corresponds to the intercept, but an intercept need not be present. This approach can be straightforwardly adapted to the multivariate case as for the TVP-VAR-SV model that will be considered in Section 5.

The error variance ${\sigma }_t^2$ in the observation equation is either homoscedastic (${\sigma }_t^2\equiv {\sigma }^2$ for all $t=1,\dots ,T$) or follows a stochastic volatility (SV) specification (Jacquier et al. 1994), where the log volatility $h_t=log{\sigma }_t^2$ follows an AR(1) process. Specifically,

$$\begin{array}{c}\hfill h_t|h_{t-1},\mu ,\varphi ,{\sigma }_{\eta }^2\sim \mathcal{N}\left(\mu +\varphi (h_{t-1}-\mu ),{\sigma }_{\eta }^2\right).\end{array}$$

(2)

For Bayesian inference, priors have to be chosen for the unknown variances ${\theta }_1,\dots ,{\theta }_d$ and the unknown initial means ${\beta }_1,\dots ,{\beta }_d$. In order to shrink dynamic coefficients to static ones and, in this way, avoid overfitting, a shrinkage prior is placed on ${\theta }_j$ that puts a lot of prior mass close to zero. One such prior is the double gamma prior, employed recently by Bitto and Frühwirth-Schnatter (2019). The double gamma prior can be expressed as a scale-mixture of gamma distributions, with the following hierarchical representation:

$$\begin{array}{c}\hfill {\theta }_j|{\xi }_j^2\sim \mathcal{G}\left(\frac{1}{2},\frac{1}{2{\xi }_j^2}\right),{\xi }_j^2|a^{\xi },{\kappa }_B^2\sim \mathcal{G}\left(a^{\xi },\frac{a^{\xi }{\kappa }_B^2}{2}\right).\end{array}$$

(3)

In the double gamma prior, each innovation variance ${\theta }_j$ is mixed over its own scale parameter ${\xi }_j^2$, each of which has an independent gamma distribution, with a common hyperparameter ${\kappa }_B^2$. Moreover, the ${\xi }_j^2$’s play the role of local (component specific) shrinkage parameters, while the parameter ${\kappa }_B^2$ is a (common) global shrinkage parameter.

We propose an extension of the double gamma prior to a triple gamma prior, where another layer is added to the hierarchy:

$$\begin{array}{c}\hfill {\theta }_j|{\xi }_j^2\sim \mathcal{G}\left(\frac{1}{2},\frac{1}{2{\xi }_j^2}\right),{\xi }_j^2|a^{\xi },{\kappa }_j^2\sim \mathcal{G}\left(a^{\xi },\frac{a^{\xi }{\kappa }_j^2}{2}\right),{\kappa }_j^2|c^{\xi },{\kappa }_B^2\sim \mathcal{G}\left(c^{\xi },\frac{c^{\xi }}{{\kappa }_B^2}\right).\end{array}$$

(4)

The main difference with the double gamma prior is that the prior scale of the ${\xi }_j^2$’s is not identical, but each ${\xi }_j^2$ depends on its component specific scale ${\kappa }_j^2$. We will show in Section 2.2 that the triple gamma prior can be represented as a global-local shrinkage prior in the sense of Polson and Scott (2012a) where the local shrinkage parameters ${\xi }_j^2$ arise from an $\mathrm{F}\left(2a^{\xi },2c^{\xi }\right)$ distribution. Hence, the triple gamma prior contains the Horseshoe prior and many other well-known shrinkage priors as special cases, as will be discussed in Section 2.3.

The shrinkage behaviour of the triple gamma prior becomes even more apparent when we rewrite model (1) in the non-centered parametrization introduced in Frühwirth-Schnatter and Wagner (2010):

$$\begin{array}{cc}\hfill & {\tilde{\mathit{\beta }}}_t={\tilde{\mathit{\beta }}}_{t-1}+{\tilde{\mathbf{w}}}_t,{\tilde{\mathbf{w}}}_t\sim {\mathcal{N}}_d\left(\mathbf{0},{\mathbf{I}}_d\right),\hfill \\ \hfill & y_t={\mathbf{x}}_t\mathit{\beta }+{\mathbf{x}}_t\mathrm{Diag}\left(\sqrt{{\theta }_1},\dots ,\sqrt{{\theta }_d}\right){\tilde{\mathit{\beta }}}_t+{\epsilon }_t,{\epsilon }_t\sim \mathcal{N}\left(0,{\sigma }_t^2\right),\hfill \end{array}$$

(5)

with ${\tilde{\mathit{\beta }}}_0\sim {\mathcal{N}}_d\left(\mathbf{0},{\mathbf{I}}_d\right)$, where ${\mathbf{I}}_d$ is the d-dimensional identity matrix. Both representations are equivalent and we can specify a prior either on the variances ${\theta }_j$ in (1) or the scale parameters ${\sqrt{\theta }}_j$ in (5). Using the fact that ${\theta }_j/{\xi }_j^2\sim {\chi }_1^2$ and the ${\chi }_1^2$-distribution can be represented as ${\chi }_1^2=Z_j^2$, where $Z_j\sim \mathcal{N}\left(0,1\right)$ follows a standard normal distribution, we can match prior (4) to the non-centered parametrization (5). This yields

$$\begin{array}{c}\hfill {\sqrt{\theta }}_j|{\xi }_j^2\sim \mathcal{N}\left(0,{\xi }_j^2\right),{\xi }_j^2|a^{\xi },{\kappa }_j^2\sim \mathcal{G}\left(a^{\xi },\frac{a^{\xi }{\kappa }_j^2}{2}\right),{\kappa }_j^2|c^{\xi },{\kappa }_B^2\sim \mathcal{G}\left(c^{\xi },\frac{c^{\xi }}{{\kappa }_B^2}\right).\end{array}$$

(6)

In (6), we could force ${\sqrt{\theta }}_j$ to take on only positive values, however, we do not impose such a constraint and allow ${\sqrt{\theta }}_j$ to take on negative values. Since the half-normal ${\sqrt{\theta }}_j\sim \mathcal{N}\left(0,{\xi }_j^2\right)I\{{\sqrt{\theta }}_j>0\}$ also implies that ${\theta }_j\sim {\xi }_j^2{\chi }_1^2$, the question arises whether the negative half is of importance. Whenever inference is performed under the non-centered parametrization (5), as is done in Section 4, restricting the prior to the positive half will lead to automatic truncation of the full conditional posterior $p({\sqrt{\theta }}_j|{\tilde{\mathit{\beta }}}_0,\dots ,{\tilde{\mathit{\beta }}}_T,\mathit{y},·)$ to the positive part during MCMC sampling. If the positive and the negative mode of the marginal posterior $p({\sqrt{\theta }}_j|\mathit{y})$ are well-separated, then this will not matter. However, if the true value of ${\theta }_j$ is close or equal to zero and $p({\sqrt{\theta }}_j|\mathit{y})$ is concentrated at zero, this truncation will introduce a bias, because the negative half is not accounted for.

Interestingly, prior (6) is related to the so-called normal-gamma-gamma prior consider by Griffin and Brown (2017) in the context of defining hierarchical shrinkage priors for regression models. This relation is helpful in choosing a prior on the fixed coefficients ${\beta }_1,\dots ,{\beta }_d$. To allow shrinkage of these coefficients toward insignificant ones in a TVP model, we extend Bitto and Frühwirth-Schnatter (2019) further by assuming such a normal-gamma-gamma prior on ${\beta }_1,\dots ,{\beta }_d$:

$$\begin{array}{c}\hfill {\beta }_j|{\tau }_j^2\sim \mathcal{N}\left(0,{\tau }_j^2\right),{\tau }_j^2|a^{\tau },{\lambda }_j^2\sim \mathcal{G}\left(a^{\tau },\frac{a^{\tau }{\lambda }_j^2}{2}\right),{\lambda }_j^2|c^{\tau },{\lambda }_B^2\sim \mathcal{G}\left(c^{\tau },\frac{c^{\tau }}{{\lambda }_B^2}\right).\end{array}$$

(7)

In Section 2.4, we will discuss hierarchical versions of both priors, by putting a hyperprior on the parameters ${\kappa }_B^2$, ${\lambda }_B^2$, $a^{\xi }$, $a^{\tau }$, $c^{\xi }$, and $c^{\tau }$.

2.2. Properties of the Triple Gamma Prior

In this section, we study the mathematical properties of the triple gamma prior. It is shown in Theorem 1 that the triple gamma prior is a global-local shrinkage prior where the local shrinkage parameters arise from the $\mathrm{F}\left(2a^{\xi },2c^{\xi }\right)$ distribution. Furthermore, a closed form of the marginal shrinkage prior $p({\sqrt{\theta }}_j|{\varphi }^{\xi },a^{\xi },c^{\xi })$ is given in Theorem 1, which is proven in Appendix A.

Theorem 1.

For the triple gamma prior defined in (4), with $a^{\xi }>0$ and $c^{\xi }>0$, the following holds:

(a) 

It has following representation as a local-global shrinkage prior:

$$\begin{array}{cc}\hfill & {\sqrt{\theta }}_j|{\psi }_j^2,{\kappa }_B^2\sim \mathcal{N}\left(0,\frac{2}{{\kappa }_B^2}{\psi }_j^2\right),{\psi }_j^2|a^{\xi },c^{\xi }\sim \mathrm{F}\left(2a^{\xi },2c^{\xi }\right).\hfill \end{array}$$

(8)

(b) 

The marginal prior $p({\sqrt{\theta }}_j|{\varphi }^{\xi },a^{\xi },c^{\xi })$ takes the following form with ${\varphi }^{\xi }=\frac{2c^{\xi }}{{\kappa }_B^2{a}^{\xi }}$,

$$\begin{array}{c}\hfill p({\sqrt{\theta }}_j|{\varphi }^{\xi },a^{\xi },c^{\xi })=\frac{\Gamma (c^{\xi }+\frac{1}{2})}{\sqrt{2\pi {\varphi }^{\xi }}B(a^{\xi },c^{\xi })}U\left(c^{\xi }+\frac{1}{2},\frac{3}{2}-a^{\xi },\frac{{\theta }_j}{2{\varphi }^{\xi }}\right),\end{array}$$

(9)

where $U\left(a,b,z\right)$ is the confluent hyper-geometric function of the second kind:

$$\begin{array}{c}\hfill U\left(a,b,z\right)=\frac{1}{\Gamma \left(a\right)}{\int }_0^{\infty }{e}^{-zt}{t}^{a-1}{(1+t)}^{b-a-1}dt.\end{array}$$

In Figure 1, we can see the marginal prior distribution of ${\sqrt{\theta }}_j$ under the triple gamma prior $a^{\xi }=c^{\xi }=0.1$ and under other well-known shrinkage priors which are special cases of the triple gamma, see

Table 1. Priors on ${\sqrt{\theta }}_j$ which are equivalent to (top) or special cases of (bottom) the triple gamma prior.

Prior for ${\sqrt{\mathit{\theta }}}_{\mathit{j}}$		${\mathit{a}}^{\mathit{\xi }}$	${\mathit{c}}^{\mathit{\xi }}$	${\mathit{\kappa }}_{\mathit{B}}^2$	${\mathit{\varphi }}^{\mathit{\xi }}$
$\mathcal{N}\left(0,{\psi }_j^2\right),{\psi }_j^2\sim GG\left(a^{\xi },c^{\xi },{\varphi }^{\xi }\right)$	normal-gamma-gamma	$a^{\xi }$	$c^{\xi }$	$\frac{2c^{\xi }}{{\varphi }^{\xi }{a}^{\xi }}$	${\varphi }^{\xi }$
$\mathcal{N}\left(0,\frac{1}{{\kappa }_j}-1\right),{\kappa }_j\sim \mathcal{TPB}\left(a^{\xi },c^{\xi },{\varphi }^{\xi }\right)$	generalized beta mixture	$a^{\xi }$	$c^{\xi }$	$\frac{2c^{\xi }}{{\varphi }^{\xi }{a}^{\xi }}$	${\varphi }^{\xi }$
$\mathcal{N}\left(0,{\psi }_j^2\right),{\psi }_j^2\sim \mathrm{SBeta}2\left(a^{\xi },c^{\xi },{\varphi }^{\xi }\right)$	hierarchical scaled beta2	$a^{\xi }$	$c^{\xi }$	$\frac{2c^{\xi }}{{\varphi }^{\xi }{a}^{\xi }}$	${\varphi }^{\xi }$
$\mathcal{DE}\left(0,\sqrt{2}{\psi }_j\right),{\psi }_j^2\sim \mathcal{G}\left(c^{\xi },\frac{1}{{\lambda }^2}\right)$	normal-exponential-gamma	1	$c^{\xi }$	$2{\lambda }^2{c}^{\xi }$	$\frac{1}{{\lambda }^2}$
$\mathcal{N}\left(0,{\tau }^2{\psi }_j^2\right),{\psi }_j\sim t_1$	Horseshoe	$\frac{1}{2}$	$\frac{1}{2}$	$\frac{2}{{\tau }^2}$	${\tau }^2$
$\mathcal{N}\left(0,\frac{1}{{\kappa }_j}-1\right),{\kappa }_j\sim \mathcal{B}\left(1/2,1\right)$	Strawderman-Berger	$\frac{1}{2}$	1	4	1
$\mathcal{N}\left(0,{\tau }^2{\tilde{\xi }}_j\right),{\tilde{\xi }}_j\sim \mathcal{G}\left(a^{\xi },a^{\xi }\right)$	double gamma	$a^{\xi }$	∞	$\frac{2}{{\tau }^2}$	-
$\mathcal{N}\left(0,{\tau }^2{\tilde{\xi }}_j\right),{\tilde{\xi }}_j\sim \mathcal{E}\left(1\right)$	Lasso	1	∞	$\frac{2}{{\tau }^2}$	-
$t_{\nu }\left(0,{\tau }^2\right)$	half-t	∞	$\frac{\nu }{2}$	$\frac{2}{{\tau }^2}$	-
$t_1\left(0,{\tau }^2\right)$	half-Cauchy	∞	$\frac{1}{2}$	$\frac{2}{{\tau }^2}$	-
$\mathcal{N}\left(0,B_0\right)$	normal	∞	∞	$\frac{2}{B_0}$	-

. Theorem 1 also allows us to give a closed form for the prior $p({\theta }_j|{\varphi }^{\xi },a^{\xi },c^{\xi })=p({\sqrt{\theta }}_j|{\varphi }^{\xi },a^{\xi },c^{\xi })/{\sqrt{\theta }}_j$.1

Global-local shrinkage priors are typically compared in terms of the concentration around the origin and the tail behaviour. For the triple gamma prior $p({\sqrt{\theta }}_j|{\varphi }^{\xi },a^{\xi },c^{\xi })$, the two shape parameters $a^{\xi }$ and $c^{\xi }$ play a crucial role in this respect, see Theorem 2 which is proven in Appendix A.

Theorem 2.

The triple gamma prior (9) satisfies the following:

(a) 

For $0<a^{\xi }<0.5$ and small values of ${\sqrt{\theta }}_j$,

$$\begin{array}{c}\hfill p({\sqrt{\theta }}_j|{\varphi }^{\xi },a^{\xi },c^{\xi })=\frac{\Gamma (\frac{1}{2}-a^{\xi })}{\sqrt{\pi }{\left(2{\varphi }^{\xi }\right)}^{a^{\xi }}B(a^{\xi },c^{\xi })}{\left(\frac{1}{{\sqrt{\theta }}_j}\right)}^{1-2a^{\xi }}+O\left(1\right).\end{array}$$

(b) 

For $a^{\xi }=0.5$ and small values of ${\sqrt{\theta }}_j$,

$$\begin{array}{c}\hfill p({\sqrt{\theta }}_j|{\varphi }^{\xi },a^{\xi },c^{\xi })=\frac{1}{\sqrt{2\pi {\varphi }^{\xi }}B(0.5,c^{\xi })}\left(-log{\theta }_j+log\left(2{\varphi }^{\xi }\right)-\psi (c^{\xi }+0.5)\right)+O(|{\theta }_{j}log{\theta }_j|),\end{array}$$

where $\psi (·)$ is the digamma function.

(c) 

For $a^{\xi }>0.5$,

$$\begin{array}{c}\hfill {\displaystyle \underset{{\sqrt{\theta }}_j\to 0}{lim}}p({\sqrt{\theta }}_j|{\varphi }^{\xi },a^{\xi },c^{\xi })=\frac{\Gamma (c^{\xi }+\frac{1}{2})\Gamma (a^{\xi }-\frac{1}{2})}{\sqrt{2\pi {\varphi }^{\xi }}\Gamma \left(c^{\xi }\right)\Gamma \left(a^{\xi }\right)}.\end{array}$$

(d) 

As ${\sqrt{\theta }}_j\to \infty$,

$$\begin{array}{c}\hfill p({\sqrt{\theta }}_j|{\varphi }^{\xi },a^{\xi },c^{\xi })=\frac{\Gamma (c^{\xi }+\frac{1}{2}){\left(2{\varphi }^{\xi }\right)}^{c^{\xi }}}{\sqrt{\pi }B(a^{\xi },c^{\xi })}{\left(\frac{1}{{\sqrt{\theta }}_j}\right)}^{2c^{\xi }+1}\left[1+O\left(\frac{1}{{\theta }_j}\right)\right].\end{array}$$

From Theorem 2, Part (a) and (b), we find that the triple gamma prior $p({\sqrt{\theta }}_j|{\varphi }^{\xi },a^{\xi },c^{\xi })$ has a pole at the origin, if $a^{\xi }\le 0.5$. According to Part (a), the pole is more pronounced, the closer $a^{\xi }$ gets to 0. For $a^{\xi }>0.5$, we find from Part (c) that $p({\sqrt{\theta }}_j|{\varphi }^{\xi },a^{\xi },c^{\xi })$ is bounded at zero by a positive upper bound which is finite, as long as $0<c^{\xi }<\infty$. Part (d) shows that the triple gamma prior $p({\sqrt{\theta }}_j|{\varphi }^{\xi },a^{\xi },c^{\xi })$ has polynomial tails, with the shape parameter $c^{\xi }$ controlling the tail index. Prior moments $\mathrm{E}({\left({\sqrt{\theta }}_j\right)}^k|{\varphi }^{\xi },a^{\xi },c^{\xi })$ exist up to $k<2c^{\xi }$. Hence, the triple gamma prior has no finite moments for $c^{\xi }<1/2$.

Finally, additional useful representations of the triple gamma prior as a global-local shrinkage prior are summarized in Lemma 1 which is proven in Appendix A. Representation (a) shows that the triple gamma is an extension of the double gamma prior where the Gaussian prior ${\sqrt{\theta }}_j|{\xi }_j^2\sim \mathcal{N}(0,{\xi }_j^2)$ is substituted by a heavier-tailed Student-t prior, making the prior more robust to large values of ${\sqrt{\theta }}_j$. Representation (b) and (c) will be useful for MCMC inference in Section 4. Representations (c) and (d) show that for a triple gamma prior with finite $a^{\xi }$ and $c^{\xi }$, ${\varphi }^{\xi }$ acts as a global shrinkage parameter, in addition to $2/{\kappa }_B^2$.

Lemma 1.

For $a^{\xi }>0$ and $c^{\xi }>0$, the triple gamma prior (4) has the following alternative representations:

$$\begin{array}{cc}\hfill & \left(a\right){\sqrt{\theta }}_j|{\tilde{\xi }}_j^2,c^{\xi },{\kappa }_B^2\sim t_{2c^{\xi }}\left(0,\frac{2}{{\kappa }_B^2}{\tilde{\xi }}_j^2\right),{\tilde{\xi }}_j^2|a^{\xi }\sim \mathcal{G}\left(a^{\xi },a^{\xi }\right),\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill & \left(b\right){\sqrt{\theta }}_j|{\check{\xi }}_j^2,c^{\xi },{\kappa }_B^2\sim t_{2c^{\xi }}\left(0,\frac{2}{a^{\xi }{\kappa }_B^2}{\check{\xi }}_j^2\right),{\check{\xi }}_j^2|a^{\xi }\sim \mathcal{G}\left(a^{\xi },1\right).\hfill \end{array}$$

(11)

Additional representations for $0<a^{\xi }<\infty$ and $0<c^{\xi }<\infty$ based on ${\varphi }^{\xi }=\frac{2c^{\xi }}{{\kappa }_B^2{a}^{\xi }}$ are

$$\begin{array}{cc}\hfill & \left(c\right){\sqrt{\theta }}_j|{\check{\xi }}_j^2,{\check{\kappa }}_j^2,{\varphi }^{\xi }\sim \mathcal{N}\left(0,{\varphi }^{\xi }{\check{\xi }}_j^2/{\check{\kappa }}_j^2\right),{\check{\xi }}_j^2|a^{\xi }\sim \mathcal{G}\left(a^{\xi },1\right),{\check{\kappa }}_j^2|c^{\xi }\sim \mathcal{G}\left(c^{\xi },1\right),\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill & \left(d\right){\sqrt{\theta }}_j|{\tilde{\psi }}_j^2,{\varphi }^{\xi }\sim \mathcal{N}\left(0,{\varphi }^{\xi }{\tilde{\psi }}_j^2\right),{\tilde{\psi }}_j^2|a^{\xi },c^{\xi }\sim \mathcal{BP}\left(a^{\xi },c^{\xi }\right),\hfill \end{array}$$

(13)

where $\mathcal{BP}\left(a^{\xi },c^{\xi }\right)$ is the beta-prime distribution.2

2.3. Relation of the Triple Gamma to Other Shrinkage Priors

The triple gamma prior can be related to the very active research on shrinkage priors in a Bayesian framework in various ways. On the one hand, popular priors for variance parameters introduced as robust alternatives to the inverse gamma prior are special cases of the triple gamma, see Table 1. For instance, in (8), ${\psi }_j^2$ converges a.s. to 1, as $a^{\xi }\to \infty$ and $c^{\xi }\to \infty$, and the triple gamma reduces to a normal distribution for ${\sqrt{\theta }}_j$, applied for univariate TVP models (Frühwirth-Schnatter 2004) and unobserved component state space model (Frühwirth-Schnatter and Wagner 2010). For $c^{\xi }\to \infty$, $\mathrm{F}\left(2a^{\xi },2c^{\xi }\right)$ converges to the $\mathcal{G}\left(a^{\xi },a^{\xi }\right)$ distribution and the triple gamma reduces to the Bayesian Lasso for $a^{\xi }=1$ (Belmonte et al. 2014) and otherwise to the double gamma (Bitto and Frühwirth-Schnatter 2019) applied in sparse TVP models.

Gelman (2006) introduced the half-t and the half-Cauchy prior for variance parameters in hierarchical models, by assuming that ${\sqrt{\theta }}_j$ follows a “folded” t-distribution, that is, a t-distribution truncated to $[0,\infty )$, see also Polson and Scott (2012b). In (10), ${\tilde{\xi }}_j^2$ converges a.s. to 1 as $a^{\xi }\to \infty$ and the triple gamma reduces to a $t_{2c^{\xi }}$- distribution and to the Cauchy distribution for $c^{\xi }=1/2$, however without being “folded”, since we allow ${\sqrt{\theta }}_j$ to take on negative values, as explained in Section 2.1.

On the other hand, the triple gamma prior is related to popular shrinkage priors in regression models. It extends the generalized beta mixture prior introduced by Armagan et al. (2011) for variable selection in regression models,

$$\begin{array}{c}\hfill {\beta }_j|{\xi }_j^2\sim \mathcal{N}\left(0,{\xi }_j^2\right),{\xi }_j^2\sim \mathcal{G}\left(a^{\xi },{\lambda }_j\right),{\lambda }_j\sim \mathcal{G}\left(c^{\xi },{\varphi }^{\xi }\right),\end{array}$$

to variance selection in state space and TVP models. This is evident from rewriting (4) as ${\xi }_j^2\sim \mathcal{G}\left(a^{\xi },{\lambda }_j\right),{\lambda }_j\sim \mathcal{G}\left(c^{\xi },{\varphi }^{\xi }\right)$. We exploit this relationship in Section 3.1 to investigate the shrinkage profile of a triple gamma prior. Using Armagan et al. (2011, Definition 2), the triple gamma prior can be written as

$$\begin{array}{c}\hfill {\sqrt{\theta }}_j|{\rho }_j\sim \mathcal{N}\left(0,1/{\rho }_j-1\right),{\rho }_j|a^{\xi },c^{\xi },{\varphi }^{\xi }\sim \mathcal{TPB}\left(a^{\xi },c^{\xi },{\varphi }^{\xi }\right),\end{array}$$

(14)

where $\mathcal{TPB}\left(a^{\xi },c^{\xi },{\varphi }^{\xi }\right)$ is the three-parameter beta distribution with density:

$$\begin{array}{c}\hfill p\left({\rho }_j\right)=\frac{1}{B(a^{\xi },c^{\xi })}{\left({\varphi }^{\xi }\right)}^{c^{\xi }}{\rho }_j^{c^{\xi }-1}{(1-{\rho }_j)}^{a^{\xi }-1}{\left(1+({\varphi }^{\xi }-1){\rho }_j\right)}^{-(a^{\xi }+c^{\xi })}.\end{array}$$

(15)

From (14) and (15), it becomes evident that the Strawderman-Berger prior ${\sqrt{\theta }}_j|{\rho }_j\sim \mathcal{N}\left(0,1/{\rho }_j-1\right)$, ${\rho }_j\sim \mathcal{B}\left(1/2,1\right)$ (Berger 1980; Strawderman 1971) is that special case of the triple gamma prior where ${\varphi }^{\xi }=1$, $a^{\xi }=1/2$, and $c^{\xi }=1$.

The special case of a triple gamma, where $a^{\xi }=c^{\xi }=1/2$, corresponds to a Horseshoe prior (Carvalho et al. 2009, 2010) on ${\sqrt{\theta }}_j$ with global shrinkage parameter ${\tau }^2=2/{\kappa }_B^2$, since ${\psi }_j^2\sim \mathrm{F}\left(1,1\right)$ implies that ${\psi }_j\sim t_1$. The Horseshoe prior has been introduced for variable selection in regression models and has been shown to have excellent theoretical properties in this context for the “nearly black” case (van der Pas et al. 2014). The triple gamma is a generalization of the Horseshoe prior, with a similar shrinkage profile, however with much more mass close to the corner solutions. Most importantly, as will be discussed in Section 3.1, this leads to a BMA-type behaviour of the triple gamma prior for small values of $a^{\xi }$ and $c^{\xi }$.

The vast literature on shrinkage priors contains many more related priors. Rescaling ${\xi }_j^2=2/\left({\kappa }_B^2\right){\psi }_j^2$ in (8), for instance, yields a representation involving a scaled beta2 distribution,3

$$\begin{array}{c}\hfill {\sqrt{\theta }}_j|{\xi }_j^2\sim \mathcal{N}\left(0,{\xi }_j^2\right),{\xi }_j^2|a^{\xi },c^{\xi },{\varphi }^{\xi }\sim \mathrm{SBeta}2\left(a^{\xi },c^{\xi },{\varphi }^{\xi }\right),\end{array}$$

(16)

as is easily derived from (A2). The scaled beta2 was introduced by Pérez et al. (2017) in hierarchical models as a robust prior for scale parameters, ${\sqrt{\theta }}_j$, and variance parameters, ${\theta }_j$, alike. Based on (16), the triple gamma can be seen as a hierarchical extension of this prior which puts a scaled beta2 distribution on the scaling parameter ${\xi }_j^2$ of a Gaussian prior for ${\sqrt{\theta }}_j$, see Table 1. Griffin and Brown (2017) termed prior (16) gamma-gamma distribution, denoted by $GG\left(a^{\xi },c^{\xi },\varphi \right)$.

For $a^{\xi }=1$, the triple gamma reduces to the normal-exponential-gamma which has a representation as a scale-mixture of double exponential $\mathcal{DE}\left(0,\sqrt{2}{\psi }_j\right)$-distributions, see Table 1.

It has been considered for variable selection in regression models (Griffin and Brown 2011) and locally adaptive B-spline models (Scheipl and Kneib 2009). The R2-D2 prior suggested by Zhang et al. (2017) for high-dimensional regression models is another special case of the triple gamma. It reads

$$\begin{array}{c}\hfill {\beta }_j\sim \mathcal{N}\left(0,{\sigma }^2{\varphi }_j\omega \right),({\varphi }_1,\dots ,{\varphi }_d)\sim \mathcal{D}\left(a^{\tau },\dots ,a^{\tau }\right),\omega \sim \mathcal{G}\left(a,\tau \right),\tau \sim \mathcal{G}\left(b,1\right),\end{array}$$

where $a=d{a}^{\tau }$ and ${\sigma }^2$ is the residual error variance of the regression model. As shown by Zhang et al. (2017), this implies the following prior for the coefficient of determination: $R^2\sim \mathcal{B}\left(a,b\right)$ which motivates holding a fixed, while $a^{\tau }$ decreases as d increases. Using that ${\varphi }_j\omega \sim \mathcal{G}\left(a^{\tau },\tau \right)$, we can show that the R2-D2 prior is equivalent to the following hierarchical normal gamma prior applied in Bitto and Frühwirth-Schnatter (2019) for TVP models:

$$\begin{array}{c}\hfill {\beta }_j|{\tau }_j^2\sim \mathcal{N}\left(0,{\tau }_j^2\right),{\tau }_j^2\sim \mathcal{G}\left(a^{\tau },a^{\tau }{\lambda }_B^2/2\right),{\lambda }_B^2\sim \mathcal{G}\left(b,2{\sigma }^2/a^{\tau }\right).\end{array}$$

The popular Dirichlet-Laplace prior, ${\sqrt{\theta }}_j|{\psi }_j\sim \mathcal{DE}\left(0,{\psi }_j\right)$, however, is not related to the triple gamma as the prior scale ${\psi }_j$ rather than the prior variance ${\psi }_j^2$ follows a gamma distribution, see again Table 1.

2.4. Using the Triple Gamma for Variance Selection in TVP Models

A challenging question is how to choose the parameters $a^{\xi }$, $c^{\xi }$ and ${\kappa }_B^2$ or ${\varphi }^{\xi }$ of the triple gamma prior in the context of variance selection for TVP models. In addition, in a TVP context, the shrinkage parameters $a^{\tau }$, $c^{\tau }$ and ${\lambda }_B^2$ or ${\varphi }^{\tau }=2c^{\tau }/\left(a^{\tau }{\lambda }_B^2\right)$ for the prior (7) of the initial values ${\beta }_j$ have to be selected.

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 }$, see Theorem 2. First of all, we need values $0<a^{\xi }\le 0.5$ to induce a pole at 0. Second, values of $0<c^{\xi }<0.5$ will lead to very heavy tails. For very small values of $a^{\xi }$ and $c^{\xi }$, the triple Gamma is a proper prior that behaves nearly as the improper normal-Jeffrey’s prior (Figueiredo 2003), where $p\left({\sqrt{\theta }}_j\right)\propto 1/{\sqrt{\theta }}_j$ and $p\left({\rho }_j\right)\propto {\rho }_j^{-1}{(1-{\rho }_j)}^{-1}$.

Ideally, we would place a hyper prior distribution on all shrinkage parameters which would allow us to learn the global and the local degree of sparsity, both for the variances and the initial values. Such a hierarchical triple gamma prior introduces dependence among the local shrinkage parameters ${\xi }_1^2,\dots ,{\xi }_d^2$ in (4) and, consequently, among ${\theta }_1,\dots ,{\theta }_d$ in the joint (marginal) prior $p({\theta }_1,\dots ,{\theta }_d)$. Introducing such dependence is desirable in that it allows to learn the degree of variance sparsity in TVP models, meaning that how much a variance is shrunken toward zero depends on how close the other variances are to zero. However, first naïve approaches with rather uninformative, independent priors on ${\kappa }_B^2$, $a^{\xi }$, $c^{\xi }$ and ${\lambda }_B^2$, $a^{\tau }$, $c^{\tau }$ were not met with much success and we found it necessary to carefully design appropriate hyper priors.

Hierarchical versions of the Bayesian Lasso (Belmonte et al. 2014) and the double gamma prior (Bitto and Frühwirth-Schnatter 2019) in TVP models are based on the gamma prior ${\kappa }_B^2\sim \mathcal{G}\left(d_1,d_2\right)$. Interestingly, this choice can be seen as 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 light of this relationship, it is not surprising that very small values of $d_1$ were applied in these papers to ensure heavy tails of $p({\sqrt{\theta }}_j|d_1,d_2)$. Since a triple gamma prior already has heavy tails, we choose a different hyperprior in the present paper.

For the case $a^{\xi }=c^{\xi }=1/2$, the global shrinkage parameter $\tau$ of the Horseshoe prior typically follows a Cauchy prior, $\tau \sim t_1$ (Bhadra et al. 2017b; Carvalho et al. 2009), see also Bhadra et al. (2019, Section 5). The relationship ${\varphi }^{\xi }=2/{\kappa }_B^2={\tau }^2$ between the various global shrinkage parameters (see Table 1) implies in this case ${\varphi }^{\xi }\sim \mathrm{F}\left(1,1\right)$ or, equivalently, ${\kappa }_B^2/2\sim \mathrm{F}\left(1,1\right)$.

For a triple gamma prior with arbitrary $a^{\xi }$ and $c^{\xi }$, this is a special case of the following prior:

$$\begin{array}{c}\hfill \left.\frac{{\kappa }_B^2}{2}\right|a^{\xi },c^{\xi }\sim \mathrm{F}\left(2a^{\xi },2c^{\xi }\right),\end{array}$$

(17)

which will be motivated in Section 3.2. Under this prior, the triple gamma prior exhibits BMA-like behavior with a uniform prior on an appropriately defined model size (see Theorem 3). Prior (17) is equivalent to following representations:

$$\begin{array}{cc}\hfill & {\kappa }_B^2|a^{\xi }\sim \mathcal{G}\left(a^{\xi },d_2\right),d_2|a^{\xi },c^{\xi }\sim \mathcal{G}\left(c^{\xi },{\displaystyle \frac{2c^{\xi }}{a^{\xi }}}\right),\hfill \\ \hfill & {\varphi }^{\xi }|a^{\xi },c^{\xi }\sim \mathcal{BP}\left(c^{\xi },a^{\xi }\right).\hfill \end{array}$$

(18)

Concerning $a^{\xi }$ and $c^{\xi }$, we choose the following priors:

$$\begin{array}{c}\hfill 2a^{\xi }\sim \mathcal{B}({\alpha }_{a^{\xi }},{\beta }_{a^{\xi }}),2c^{\xi }\sim \mathcal{B}({\alpha }_{c^{\xi }},{\beta }_{c^{\xi }}).\end{array}$$

(19)

Hence, we are restricting the support of $a^{\xi }$ and $c^{\xi }$ to $(0,0.5)$, following the insights brought to us by Theorem 2.

We follow a similar strategy for the parameters $a^{\tau }$, $c^{\tau }$ and ${\lambda }_B^2$ (${\varphi }^{\tau }$) of the prior (7) of the initial values ${\beta }_j$:

$$\begin{array}{c}\hfill \left.\frac{{\lambda }_B^2}{2}\right|a^{\tau },c^{\tau }\sim \mathrm{F}\left(2a^{\tau },2c^{\tau }\right),2a^{\tau }\sim \mathcal{B}({\alpha }_{a^{\tau }},{\beta }_{a^{\tau }}),2c^{\tau }\sim \mathcal{B}({\alpha }_{c^{\tau }},{\beta }_{c^{\tau }}),\end{array}$$

(20)

which is equivalent to ${\lambda }_B^2|a^{\tau }\sim \mathcal{G}\left(a^{\tau },e_2\right)$, $e_2|a^{\tau },c^{\tau }\sim \mathcal{G}\left(c^{\tau },2c^{\tau }/a^{\tau }\right)$, and ${\varphi }^{\tau }|a^{\tau },c^{\tau }\sim \mathcal{BP}\left(c^{\tau },a^{\tau }\right)$.

An interesting special case is the “symmetric” triple gamma, where $a^{\xi }=c^{\xi }$. Despite this constraint, the favourable shrinkage behaviour is preserved and decreasing $a^{\xi }=c^{\xi }$ toward zero simultaneously leads to a high concentration around the origin and a heavy-tailed behaviour. For a symmetric triple gamma prior, the global shrinkage parameter ${\varphi }^{\xi }$ is independent of $a^{\xi }$ and $c^{\xi }$ and is related to the global shrinkage parameters ${\kappa }_B^2$ through ${\varphi }^{\xi }=2/{\kappa }_B^2$. This induces shrinkage profiles that are symmetric around 1/2, see Section 3.1. Interestingly, a symmetric triple gamma resolves the question whether to choose a gamma or an inverse gamma prior for a variance parameter ${\psi }_j^2$. It implies the same symmetric beta-prime distribution on the variance, ${\psi }_j^2\sim \mathrm{F}\left(2a^{\xi },2a^{\xi }\right)=\mathcal{BP}\left(a^{\xi },a^{\xi }\right)$, and the information, ${\left({\psi }_j^2\right)}^{-1}\sim \mathcal{BP}\left(a^{\xi },a^{\xi }\right)$, and can be represented as a gamma prior with the scale arising from an inverse gamma prior or, equivalently, as an inverse gamma prior with the scale arising from a gamma prior:

$$\begin{array}{c}{\psi }_j^2={\check{\xi }}_j^2\times \frac{1}{{\check{\kappa }}_j^2},{\left({\psi }_j^2\right)}^{-1}={\check{\kappa }}_j^2\times \frac{1}{{\check{\xi }}_j^2},{\check{\xi }}_j^2\sim \mathcal{G}\left(a^{\xi },1\right),{\check{\kappa }}_j^2\sim \mathcal{G}\left(a^{\xi },1\right).\hfill \end{array}$$

3. Shrinkage Profiles and BMA-Like Behavior

3.1. Shrinkage Profiles

In the sparse normal-means problem where $\mathit{y}|\mathit{\beta }\sim {\mathcal{N}}_d\left(\mathit{\beta },{\sigma }^2{\mathbf{I}}_d\right)$ and ${\sigma }^2=1$, the parameter ${\rho }_j=1/(1+{\psi }_j^2)$ appearing in (14) is known as shrinkage factor and plays a fundamental role for comparing different shrinkage priors, as ${\rho }_j$ determines shrinkage toward 0.

Also in a variance selection context, it is evident from (14) that values of ${\rho }_j\approx 0$ will introduce no shrinkage on ${\theta }_j$, whereas values of ${\rho }_j\approx 1$ will introduce strong shrinkage of ${\theta }_j$ toward 0. Hence, the prior $p\left({\rho }_j\right)$, also called shrinkage profile, will play an instrumental role in the behaviour of different shrinkage priors. Following Carvalho et al. (2010), shrinkage priors are often compared in terms of the prior they imply on ${\rho }_j$, that is, how they handle shrinkage for small “observations” (in our case innovations) and how robust they are to large “observations”. Note that we ideally want a shrinkage profile that has a pole in zero (heavy tails to avoid over-shrinking signals) and a pole in one (spikiness to shrink noise). The Horseshoe prior, for example, implies ${\rho }_j\sim \mathcal{B}\left(1/2,1/2\right)$ which is a shrinkage profile that takes this much desired form of a “horseshoe”, see Figure 2.

For the triple gamma prior, the shrinkage profile is given by the three-parameter beta prior $p\left({\rho }_j\right)$ provided in (15). For ${\varphi }^{\xi }=1$, ${\rho }_j\sim \mathcal{B}\left(c^{\xi },a^{\xi }\right)$ and ${\kappa }_B^2=2c^{\xi }/a^{\xi }$. Choosing small values $a^{\xi }<<1$ will put prior mass close to 1, choosing small values $c^{\xi }<<1$ will put prior mass close to 0, whereas values for both $a^{\xi }$ and $c^{\xi }$ smaller than one will induce the form of a horseshoe prior for ${\rho }_j$. Evidently, for ${\varphi }^{\xi }=1$, a symmetric triple gamma prior with $a^{\xi }=c^{\xi }$ implies a Horseshoe prior for ${\rho }_j$ that is symmetric around 0.5. This is illustrated in Figure 2 for a symmetric triple gamma with $a^{\xi }=c^{\xi }=0.1$.

In Figure 2 we can also see the shrinkage profile for the Bayesian Lasso and the double gamma, which correspond to a triple gamma where $c^{\xi }\to \infty$.4 For the Bayesian Lasso with $a^{\xi }=1$, it is clear that the shrinkage profile $p\left({\rho }_j\right)$ converges to a constant for ${\rho }_j\to 1$, while there is no mass around ${\rho }_j=0$. This means that this prior tends to over-shrink signals, while not shrinking the noise completely to zero. A double gamma prior with $a^{\xi }<1$ has the potential to shrink the noise completely to zero, as $p\left({\rho }_j\right)$ has a pole at ${\rho }_j=1$, but $p\left({\rho }_j\right)$ has also zero mass around ${\rho }_j=0$, meaning the prior encourages over-shrinking of signals.

When we make ${\kappa }_B^2$ random, we obtain a “prior density” of shrinkage profiles, see Figure 3. We can see that such hierarchical versions of the Lasso and the double gamma have shrinkage profiles that resemble the ones of the Horseshoe and the triple gamma. We have used ${\kappa }_B^2\sim \mathcal{G}\left(0.01,0.01\right)$ for the Lasso and the double gamma, $2/{\kappa }_B^2\sim \mathrm{F}\left(1,1\right)$ for the Horseshoe and $2/{\kappa }_B^2\sim \mathrm{F}\left(0.2,0.2\right)$ for the triple gamma, see Section 2.4.

3.2. BMA-Type Behaviour

Bayesian model averaging (BMA) provides to statisticians and practitioners an essential and coherent tool to account for model uncertainty. In a multiple regression setting, the uncertainty is inherent in the choice of variables to be included. In a TVP framework, there is additional uncertainty about the time-variation of the state parameters, that is, which explanatory variables have a static and which ones a dynamic effect on the response variable. In this section, we show that the triple gamma prior mimics the typical BMA behavior, thus allowing us to incorporate model uncertainty with respect to time variation.

From the perspective of Bayesian model averaging, an ideal approach for handling sparsity in TVP models would be the use of discrete mixture priors as suggested in Frühwirth-Schnatter and Wagner (2010),

$$\begin{array}{c}\hfill p\left({\sqrt{\theta }}_j\right)=(1-\pi ){\delta }_0+\pi ·p_{slab}\left({\sqrt{\theta }}_j\right),\end{array}$$

(21)

with ${\delta }_0$ being a Dirac measure at 0, while $p_{slab}\left({\sqrt{\theta }}_j\right)$ is the prior for non-zero variances. In terms of shrinkage profiles, the discrete mixture prior (21) has a spike at ${\rho }_j=1$, with probability $1-\pi$, and a lot of prior mass at ${\rho }_j=0$, provided that the tails of $p_{\mathrm{slab}}\left({\sqrt{\theta }}_j\right)$ are heavy enough. The mixture prior (21) is considered the “gold standard” in BMA, both theoretically and empirically, see for example, Johnstone and Silverman (2004). However, MCMC inference under this prior is extremely challenging. As opposed to this, MCMC inference for the triple gamma prior is straightforward, see Section 4.

In this section, we relate the triple gamma prior to BMA based on the discrete mixture prior (21). An interesting insight is that the triple gamma prior shows a behaviour very similar to a discrete mixture prior, if both $a^{\xi }$ and $c^{\xi }$ approach zero. This induces BMA-type behaviour on the joint shrinkage profile $p({\rho }_1,\dots ,{\rho }_d)$, with a spike at all corner solutions, where some ${\rho }_j$ are very close to one, whereas the remaining ones are very close to zero.

The bivariate shrinkage profiles shown in Figure 4 give us some intuition about the convergence of a symmetric triple gamma prior with $a^{\xi }=c^{\xi }\to 0$ toward a discrete spike and slab mixture. As opposed to the Lasso and the double gamma prior, the Horseshoe and the triple gamma prior put nearly all prior mass on the “corner solutions”, which correspond to the four possibilities (a) ${\rho }_1={\rho }_2=0$, that is, no shrinkage on ${\theta }_1$ and ${\theta }_2$, (b) ${\rho }_1=1,{\rho }_2=0$, that is, shrinkage of ${\theta }_1$ toward 0 and no shrinkage on ${\theta }_2$, (c) ${\rho }_1=0,{\rho }_2=1$, that is, shrinkage of ${\theta }_1$ toward 0 and no shrinkage on ${\theta }_2$, and (d) ${\rho }_1={\rho }_2=1$, that is, shrinkage of both ${\theta }_1$ and ${\theta }_2$ toward 0.

A very important aspect of BMA is that of choosing a prior for the model dimension, K, see for example, Fernández et al. (2001) and Ley and Steel (2009). In the discrete mixture prior (21), the distribution of K depends on the choice of $\pi$. Fixing $\pi$ corresponds to a very informative prior on the model dimension, for example $\pi =0.5$ assigns more prior probability to models of dimension $d/2$ and lower prior probability to empty or full models. In fact, let ${\delta }_j$ be the indicator that tells us if the j-th coefficient is included in the model, then we have that $K={\sum }_{j=1}^d{\delta }_j\sim \mathrm{Binom}\left(d,\pi \right)$. Placing a uniform prior on $\pi$ has been shown to be a good choice, since it corresponds to placing a prior on K which is uniform on $\{0,\dots ,d\}$. Note that $\pi$ will be learned using information from all the variables. In this sense, $\pi$ is a global shrinkage parameter which will adapt to the degree of sparsity.

Following ideas in Carvalho et al. (2009), we believe that a natural way to perform variable selection in the continuous shrinkage prior framework is though thresholding. Specifically, we say that when $(1-{\rho }_j)>0.5$, or ${\rho }_j<0.5$, the variable is included, otherwise it is not. Notice that this classification via thresholding makes perfect sense in the case of a triple gamma of which the Horseshoe is a special case, but less so for a Lasso or double gamma prior, even if the shrinkage profile shows a Horseshoe-like behaviour for hierarchical versions of these priors (see again Figure 3). Notice that thresholding implies a prior on the model dimension K. Specifically,

$$\begin{array}{c}\hfill K={\displaystyle \sum _{j=1}^d}I\{{\rho }_j<0.5\}\sim \mathrm{Binom}\left(d,{\pi }^{\xi }\right),{\pi }^{\xi }=\mathrm{Pr}({\rho }_j<0.5),\end{array}$$

(22)

where ${\rho }_j|a^{\xi },{\varphi }^{\xi }\sim \mathcal{TPB}\left(a^{\xi },b^{\xi },{\varphi }^{\xi }\right)$, see (15). The choice of ${\varphi }^{\xi }$ (or ${\kappa }_B^2$) will strongly impact the prior on K. For a symmetric triple gamma with $a^{\xi }=c^{\xi }$, for instance, and fixed ${\varphi }^{\xi }=1$, that is ${\kappa }_B^2=2$, we obtain $K\sim \mathrm{Binom}\left(d,0.5\right)$, since ${\pi }^{\xi }=0.5$ regardless of $a^{\xi }$. Hence, we have to face similar problems as with fixing $\pi =0.5$ for the discrete mixture prior (21).

Placing a hyper prior on ${\varphi }^{\tau }$ and ${\varphi }^{\xi }$ (or equivalent ones on ${\lambda }_B^2$ and ${\kappa }_B^2$), as we did in Section 2.4, is as vital for BMA-type variable and variance selection through the triple gamma prior, as making $\pi$ random is for the discrete mixture prior (21). Ideally, we would like to have a uniform distribution on the model size K. We show in Theorem 3 that the hyperprior for ${\kappa }_B^2$ defined in (17) achieves exactly this goal, since ${\pi }^{\xi }$ is uniformly distributed, see Appendix A for a proof.

Theorem 3.

For a hierarchical triple gamma prior with fixed $a^{\xi }>0$ and $c^{\xi }>0$ the probability ${\pi }^{\xi }$ defined in (22) follows a uniform distribution, ${\pi }^{\xi }\sim \mathcal{U}\left[0,1\right]$, under the hyper prior

$$\begin{array}{c}\hfill \left.\frac{{\kappa }_B^2}{2}\right|a^{\xi },c^{\xi }\sim \mathrm{F}\left(2a^{\xi },2c^{\xi }\right),\end{array}$$

(23)

or, equivalently, under the hyper prior

$$\begin{array}{cc}\hfill & {\varphi }^{\xi }|a^{\xi },c^{\xi }\sim \mathcal{BP}\left(c^{\xi },a^{\xi }\right).\hfill \end{array}$$

(24)

Finally, it is important to point out that the thresholding approach allows us to estimate posterior inclusion probabilities, that is the probability that the corresponding variable is included in the model or, in the case of variance selection, that the corresponding parameter is time varying. In our simulations (Section 5.3) and in our application (Section 5.4), we will estimate the posterior inclusion probabilities obtained under different shrinkage priors.

4. MCMC Algorithm

Let $\mathit{y}=(y_1,\dots ,y_T)$ be the vector of time series observations and let $\mathit{z}$ be the set of all latent variables and unknown model parameters in a TVP model. Moreover, let ${\mathit{z}}_{-x}$ denote the set of all unknowns but x. Bayesian inference based on MCMC sampling from the posterior $p(\mathit{z}|\mathit{y})$ is summarized in Algorithm 1. The hierarchical priors introduced in Section 2.4 are employed, where ${\kappa }_B^2$ follows (17), $(a^{\xi },c^{\xi })$ follow (19), and $(a^{\tau },c^{\tau },{\lambda }_B^2)$ follow (20). For certain sampling steps, the hierarchical representation (18) is used for ${\kappa }_B^2$, and similarly for ${\lambda }_B^2$.

Algorithm 1 extends several existing algorithms such as the MCMC schemes introduced for the Horseshoe prior by Makalic and Schmidt (2016) and for the double gamma prior by Bitto and Frühwirth-Schnatter (2019). We exploit various representations of the triple gamma prior given in Lemma 1 and choose representation (12) as the baseline representation of our MCMC algorithm:

$$\begin{array}{cc}\hfill & {\beta }_j|{\check{\tau }}_j^2,{\check{\lambda }}_j^2,{\varphi }^{\tau }\sim \mathcal{N}\left(0,{\varphi }^{\tau }{\check{\tau }}_j^2/{\check{\lambda }}_j^2\right),{\check{\tau }}_j^2|a^{\tau }\sim \mathcal{G}\left(a^{\tau },1\right),{\check{\lambda }}_j^2|c^{\tau }\sim \mathcal{G}\left(c^{\tau },1\right),\hfill \\ \hfill & {\sqrt{\theta }}_j|{\check{\xi }}_j^2,{\check{\kappa }}_j^2,{\varphi }^{\xi }\sim \mathcal{N}\left(0,{\varphi }^{\xi }{\check{\xi }}_j^2/{\check{\kappa }}_j^2\right),{\check{\xi }}_j^2|a^{\xi }\sim \mathcal{G}\left(a^{\xi },1\right),{\check{\kappa }}_j^2|c^{\xi }\sim \mathcal{G}\left(c^{\xi },1\right),\hfill \end{array}$$

where ${\varphi }^{\tau }=2c^{\tau }/\left({\lambda }_B^2{a}^{\tau }\right)$ and ${\varphi }^{\xi }=2c^{\xi }/\left({\kappa }_B^2{a}^{\xi }\right)$. All conditional distributions in our MCMC scheme are available in closed form, except for the ones for $a^{\xi }$, $c^{\xi }$, $a^{\tau }$ and $c^{\tau }$, for which we will resort to a Metropolis-Hastings (MH) step within Gibbs. Several conditional distributions are the same as for the double gamma prior and we apply Algorithm 1 of Bitto and Frühwirth-Schnatter (2019). We provide more details on the derivation of the various densities in Appendix B.

Algorithm 1. MCMC inference for TVP models under the triple gamma prior.

Choose starting values for all global shrinkage parameters $(a^{\tau },c^{\tau },{\lambda }_B^2,a^{\xi },c^{\xi },{\kappa }_B^2)$ and local shrinkage parameters ${\{{\check{\tau }}_j^2,{\check{\lambda }}_j^2,{\check{\xi }}_j^2,{\check{\kappa }}_j^2\}}_{j=1}^d$, and repeat the following steps: (a)  Define for $j=1,\dots ,d$, ${\tau }_j^2={\varphi }^{\tau }{\check{\tau }}_j^2/{\check{\lambda }}_j^2$ and ${\xi }_j^2={\varphi }^{\xi }{\check{\xi }}_j^2/{\check{\kappa }}_j^2$ and sample from the posterior $p({\tilde{\mathit{\beta }}}_0,\dots ,{\tilde{\mathit{\beta }}}_T,{\beta }_1,\dots ,{\beta }_d,\sqrt{{\theta }_1},\dots ,\sqrt{{\theta }_d}|{\{{\xi }_j^2,{\tau }_j^2\}}_{j=1}^d,\mathit{y})$ using Algorithm 1, Steps (a), (b), and (c) in Bitto and Frühwirth-Schnatter (2019). In the homoscedastic case, use Step (f) of this algorithm to sample from ${\sigma }^2|{\mathit{z}}_{-{\sigma }^2},\mathit{y}$. For the SV model (2), sample the parameters μ, ϕ, and ${\sigma }_{\eta }^2$ as in Kastner and Frühwirth-Schnatter (2014), for example, using the R-packagestochvol(Kastner 2016). (b)  Use the prior $p({\sqrt{\theta }}_j|{\check{\kappa }}_j^2,a^{\xi },c^{\xi })$, marginalized w.r.t. ${\check{\xi }}_j^2$, to sample $a^{\xi }$ from $p(a^{\xi }|{\mathit{z}}_{-a^{\xi }},\mathit{y})$ via a random walk MH step on $z=log(a^{\xi }/(0.5-a^{\xi }))$. Propose $a^{\xi ,(*)}=0.5e^{z^{*}}/(1+e^{z^{*}})$, where $z^{*}\sim \mathcal{N}\left(z^{(m-1)},v^2\right)$ and $z^{(m-1)}=log(a^{\xi ,(m-1)}/(0.5-a^{\xi ,(m-1)}))$ depends on the previous value $a^{\xi ,(m-1)}$ of $a^{\xi }$, accept $a^{\xi ,(*)}$ with probability $$\begin{array}{c}\hfill min\left\{1,\frac{q_a\left(a^{\xi ,(*)}\right)}{q_a\left(a^{\xi ,(m-1)}\right)}\right\},q_a\left(a^{\xi }\right)=p(a^{\xi }|{\mathit{z}}_{-a^{\xi }},\mathit{y})a^{\xi }(0.5-a^{\xi }),\end{array}$$ and update ${\varphi }^{\xi }=2c^{\xi }/\left({\kappa }_B^2{a}^{\xi }\right)$. Explicit forms for $p(a^{\xi }|{\mathit{z}}_{-a^{\xi }},\mathit{y})$ and $log{q}_a\left(a^{\xi }\right)$ are provided in (A3) and (A4). Similarly, use the prior $p({\beta }_j|{\check{\lambda }}_j^2,a^{\tau },c^{\tau })$, marginalized w.r.t. to ${\check{\tau }}_j^2$, to sample $a^{\tau }$ via a random walk MH step and update ${\varphi }^{\tau }=2c^{\tau }/\left(a^{\tau }{\lambda }_B^2\right)$. (c)  Sample ${\check{\xi }}_j^2$, $j=1,\dots ,d$, from a generalized inverse Gaussian distribution, see (A5): $$\begin{array}{cc}\hfill & {\check{\xi }}_j^2|{\check{\kappa }}_j^2,{\theta }_j,a^{\xi },{\varphi }^{\xi }\sim \mathcal{GIG}\left(a^{\xi }-\frac{1}{2},2,{\displaystyle \frac{{\check{\kappa }}_j^2{\theta }_j}{{\varphi }^{\xi }}}\right).\hfill \end{array}$$ (25) Similarly, update ${\check{\tau }}_j^2$, $j=1,\dots ,d$, conditional on $a^{\tau }$: $$\begin{array}{c}\hfill {\check{\tau }}_j^2|{\beta }_j,{\check{\lambda }}_j^2,a^{\tau },{\varphi }^{\tau }\sim \mathcal{GIG}\left(a^{\tau }-\frac{1}{2},2,{\displaystyle \frac{{\check{\lambda }}_j^2{\beta }_j^2}{{\varphi }^{\tau }}}\right).\end{array}$$ (d)  Use the marginal Student-t distribution $p({\sqrt{\theta }}_j|{\check{\xi }}_j^2,c^{\xi },{\kappa }_B^2)$ given in (11) to sample $c^{\xi }$ from $p(c^{\xi }|{\mathit{z}}_{-c^{\xi }},\mathit{y})$ via a random walk MH step on $z=log(c^{\xi }/(0.5-c^{\xi }))$. Propose $c^{\xi ,(*)}=0.5{\mathrm{e}}^{z^{*}}/(1+{\mathrm{e}}^{z^{*}})$, where $z^{*}\sim \mathcal{N}\left(z^{(m-1)},v^2\right)$ and $z^{(m-1)}=log(c^{\xi ,(m-1)}/(0.5-c^{\xi ,(m-1)}))$ depends on the previous value $c^{\xi ,(m-1)}$ of $c^{\xi }$, accept $c^{\xi ,(*)}$ with probability $$\begin{array}{c}\hfill min\left\{1,\frac{q_c\left(c^{\xi ,(*)}\right)}{q_c\left(c^{\xi ,(m-1)}\right)}\right\},q_c\left(c^{\xi }\right)=p(c^{\xi }|{\mathit{z}}_{-c^{\xi }},\mathit{y})c^{\xi }(0.5-c^{\xi }),\end{array}$$ and update ${\varphi }^{\xi }=2c^{\xi }/\left({\kappa }_B^2{a}^{\xi }\right)$. Explicit forms for $p(c^{\xi }|{\mathit{z}}_{-c^{\xi }},\mathit{y})$ and $log{q}_c\left(c^{\xi }\right)$ are provided in (A6) and (A7). Similarly, to sample $c^{\tau }$ via a random walk MH step use the marginal distribution of ${\beta }_j|{\check{\tau }}_j^2,a^{\tau },c^{\tau }$ with respect to ${\check{\lambda }}_j^2$ and update ${\varphi }^{\tau }=2c^{\tau }/\left(a^{\tau }{\lambda }_B^2\right)$. (e)  Sample ${\check{\kappa }}_j^2$, for $j=1,\dots ,d$, from following gamma distribution, see (A8): $$\begin{array}{c}\hfill {\check{\kappa }}_j^2|{\theta }_j,{\check{\xi }}_j^2,c^{\xi },{\varphi }^{\xi }\sim \mathcal{G}\left(\frac{1}{2}+c^{\xi },{\displaystyle \frac{{\theta }_j}{2{\varphi }^{\xi }{\check{\xi }}_j^2}}+1\right).\end{array}$$ (26) Similarly, update ${\check{\lambda }}_j^2$, $j=1,\dots ,d$, conditional on $c^{\tau }$: $$\begin{array}{c}\hfill {\check{\lambda }}_j^2|{\beta }_j,{\check{\tau }}_j^2,c^{\tau },{\varphi }^{\tau }\sim \mathcal{G}\left(\frac{1}{2}+c^{\tau },{\displaystyle \frac{{\beta }_j^2}{2{\varphi }^{\tau }{\check{\tau }}_j^2}}+1\right).\end{array}$$ (f)  Sample $d_2$ from $d_2|a^{\xi },c^{\xi },{\kappa }_B^2\sim \mathcal{G}\left(a^{\xi }+c^{\xi },{\kappa }_B^2+\frac{2c^{\xi }}{a^{\xi }}\right)$, see (A9); sample from ${\kappa }_B^2$ from following gamma distribution, $$\begin{array}{c}\hfill {\kappa }_B^2|{\{{\theta }_j,{\check{\kappa }}_j^2,{\check{\xi }}_j^2\}}_{j=1}^d,a^{\xi },c^{\xi },d_2\sim \mathcal{G}\left(\frac{d}{2}+a^{\xi },{\displaystyle \frac{a^{\xi }}{4c^{\xi }}}{\displaystyle \sum _{j=1}^d}\frac{{\check{\kappa }}_j^2}{{\check{\xi }}_j^2}{\theta }_j+d_2\right),\end{array}$$ (27) see (A10), and update ${\varphi }^{\xi }=2c^{\xi }/\left({\kappa }_B^2{a}^{\xi }\right)$. Similarly, sample $e_2$ from $e_2|a^{\tau },c^{\tau },{\lambda }_B^2\sim \mathcal{G}\left(a^{\tau }+c^{\tau },{\lambda }_B^2+\frac{2c^{\tau }}{a^{\tau }}\right)$, sample ${\lambda }_B^2$ from $$\begin{array}{cc}\hfill & {\lambda }_B^2|{\{{\beta }_j,{\check{\lambda }}_j^2,{\check{\tau }}_j^2\}}_{j=1}^d,a^{\tau },c^{\tau },e_2\sim \mathcal{G}\left(\frac{d}{2}+a^{\tau },{\displaystyle \frac{a^{\tau }}{4c^{\tau }}}{\displaystyle \sum _{j=1}^d}\frac{{\check{\lambda }}_j^2}{{\check{\tau }}_j^2}{\beta }_j^2+e_2\right),\hfill \end{array}$$ and update ${\varphi }^{\tau }=2c^{\tau }/\left(a^{\tau }{\lambda }_B^2\right)$.

The MCMC scheme in Algorithm 1 is not a full conditional scheme, as several steps are based on partially marginalized distributions. That means that the sampling order matters. For instance, in Step (b), we marginalize w.r.t. ${\check{\xi }}_1^2,\dots ,{\check{\xi }}_d^2$, hence we need to update ${\check{\xi }}_1^2,\dots ,{\check{\xi }}_d^2$ after sampling $a^{\xi }$, before we update $c^{\xi }$ in Step (d) conditional on ${\check{\xi }}_1^2,\dots ,{\check{\xi }}_d^2$. Similarly, due to marginalization in Step (d), we need to update ${\check{\kappa }}_1^2,\dots ,{\check{\kappa }}_d^2$, before we update $d_2$ in Step (f). Furthermore, both Step (b) and Step (d) are based on the marginal prior of ${\kappa }_B^2$, given in (17). Hence, in Step (f), $d_2$ has to be updated from $d_2|a^{\xi },c^{\xi },{\kappa }_B^2$, before ${\kappa }_B^2$ is updated conditional on $d_2$.

For a symmetric triple gamma prior, where $a^{\xi }=c^{\xi }$, the MCMC scheme in Algorithm 1 has to be modified only slightly. Either $q_a\left(a^{\xi }\right)$ in Step (b) is adjusted and Step (d) is skipped, setting $c^{\xi }=a^{\xi }$, or $q_c\left(c^{\xi }\right)$ in Step (d) is adjusted and Step (b) is skipped, setting $a^{\xi }=c^{\xi }$. In Appendix B, we provide details in (A11) for the first case and in (A12) for the second case. Similar modifications are needed, if $a^{\tau }=c^{\tau }$. All other steps in Algorithm 1 remain the same for $a^{\xi }=c^{\xi }$ and/or $a^{\tau }=c^{\tau }$.

5. Applications to TVP-VAR-SV Models

5.1. Model

In this section, we consider a generalization of the TVP model (1), where ${\mathit{y}}_t$ is a m-dimensional time series, observed for $t=1,\dots ,T$. The time series ${\mathit{y}}_t$ is assumed to follow a time-varying parameter vector autoregressive model with stochastic volatility (TVP-VAR-SV) of order p:

$$\begin{array}{c}\hfill {\mathit{y}}_t={\mathit{c}}_t+{\mathsf{\Phi }}_{1,t}{\mathit{y}}_{t-1}+{\mathsf{\Phi }}_{2,t}{\mathit{y}}_{t-2}+\dots {\mathsf{\Phi }}_{p,t}{\mathit{y}}_{t-p}+{\mathit{\epsilon }}_t,{\mathit{\epsilon }}_t\sim {\mathcal{N}}_m\left(\mathbf{0},{\mathsf{\Sigma }}_t\right),\end{array}$$

(28)

where ${\mathit{c}}_t$ is the m-dimensional time-varying intercept, ${\mathsf{\Phi }}_{j,t}$, for $j=1,\dots ,p$ is an $m\times m$ matrix of time-varying coefficients, and ${\mathsf{\Sigma }}_t$ is the time-varying variance covariance matrix of the error term. The TVP-VAR-SV model can be written in a more compact notation as the following TVP model:

$$\begin{array}{c}\hfill {\mathit{y}}_t=({\mathbb{I}}_m\otimes {\mathbf{x}}_t){\mathit{\beta }}_t+{\mathit{\epsilon }}_t,{\mathit{\epsilon }}_t\sim {\mathcal{N}}_m\left(\mathbf{0},{\mathsf{\Sigma }}_t\right),\end{array}$$

(29)

where ${\mathbf{x}}_t=({\mathit{y}}_{t-1}^{\prime },\dots ,{\mathit{y}}_{t-p}^{\prime },1)$ is a row vector of length $mp+1$ and the time-varying parameter ${\mathit{\beta }}_t$ is defined as ${\mathit{\beta }}_t={({{\mathit{\beta }}_t^1}^{\prime },\dots ,{{\mathit{\beta }}_t^m}^{\prime })}^{\prime }$, where ${\mathit{\beta }}_t^i={({{\mathsf{\Phi }}_{1,t}}_{i•},\dots ,{{\mathsf{\Phi }}_{p,t}}_{i•},c_{t,i})}^{\prime }$. Here, ${{\mathsf{\Phi }}_{j,t}}_{i•}$ denotes the i-th row of the matrix ${\mathsf{\Phi }}_{j,t}$ and $c_{t,i}$ denotes the i-th element of ${\mathit{c}}_t$. Since the influential paper of Primiceri (2005) (see Del Negro and Primiceri (2015) for a corrigendum), this model has become a benchmark for analyzing relationships between macroeconomic variables that evolve over time, see Nakajima (2011), Koop and Korobilis (2013), Eisenstat et al. (2014), Chan and Eisenstat (2016), Feldkircher et al. (2017) and Carriero et al. (2019), among many others.

Following Frühwirth-Schnatter and Tüchler (2008), we use a Cholesky decomposition of the time-varying covariance matrix ${\mathsf{\Sigma }}_t$, that is ${\mathsf{\Sigma }}_t={\mathit{A}}_t{\mathit{D}}_t{\mathit{A}}_t^{\prime }$, where ${\mathit{D}}_t$ is a diagonal matrix and ${\mathit{A}}_t$ is lower unitriangular matrix, see Carriero et al. (2019) and Bitto and Frühwirth-Schnatter (2019) for related models. We denote with $a_{ij,t}$ the element at the i-th row and j-th column of ${\mathit{A}}_t$, and with ${\sigma }_{i,t}^2$ the i-th diagonal element of ${\mathit{D}}_t=\mathrm{Diag}\left({\sigma }_{1,t}^2\dots {\sigma }_{m,t}^2\right)$. In total, we have $m(m-1)/2+m(mp+1)$ (potentially) time-varying parameters. Using the Cholesky decomposition, we can rewrite the system as:

$$\begin{array}{c}\hfill {\mathit{y}}_t=({\mathbb{I}}_m\otimes {\mathbf{x}}_t){\mathit{\beta }}_t+{\mathit{A}}_t{\mathit{\eta }}_t,{\mathit{\eta }}_t\sim {\mathcal{N}}_m\left(\mathbf{0},{\mathit{D}}_t\right),\end{array}$$

(30)

where ${\mathit{\eta }}_t={({\eta }_{1,t},\dots ,{\eta }_{m,t})}^{\top }$. The idiosyncratic shocks ${\eta }_{i,t}\sim \mathcal{N}\left(0,{\sigma }_{i,t}^2\right)$ follow independent SV processes as in (2), with row specific parameters. Specifically, with $h_{i,t}=log{\sigma }_{i,t}^2$, we have that the logarithm of the elements of the diagonal matrix ${\mathit{D}}_t$ follow independent AR(1) processes:

$$\begin{array}{c}\hfill h_{i,t}={\mu }_i+{\varphi }_i(h_{i,t-1}-{\mu }_i)+{\nu }_{i,t},{\nu }_{i,t}\sim \mathcal{N}\left(0,{\sigma }_{\eta ,i}^2\right),\end{array}$$

for $i=1,\dots ,m$. Here, ${\mu }_i$ is the mean, ${\varphi }_i$ is the persistence parameter, and ${\sigma }_{\eta ,i}^2$ is the variance of the ith log-volatility $h_{i,t}$.

It is possible to write the TVP-VAR-SV model (30) as a system of m univariate TVP models as in (1):

$$\begin{array}{cc}\hfill y_{1,t}=& {\mathbf{x}}_t{\mathit{\beta }}_t^1+{\eta }_{1,t},{\eta }_{1,t}\sim \mathcal{N}\left(0,{\sigma }_{1,t}^2\right),\hfill \\ \hfill y_{2,t}=& {\mathbf{x}}_t{\mathit{\beta }}_t^2+a_{21,t}{\eta }_{1,t}+{\eta }_{2,t},{\eta }_{2,t}\sim \mathcal{N}\left(0,{\sigma }_{2,t}^2\right),\hfill \\ \hfill y_{3,t}=& {\mathbf{x}}_t{\mathit{\beta }}_t^3+a_{31,t}{\eta }_{1,t}+a_{32,t}{\eta }_{2,t}+{\eta }_{3,t},{\eta }_{3,t}\sim \mathcal{N}\left(0,{\sigma }_{3,t}^2\right),\hfill \\ \hfill & \cdots \hfill \\ \hfill y_{m,t}=& {\mathbf{x}}_t{\mathit{\beta }}_t^m+a_{m1,t}{\eta }_{1,t}+\dots +a_{m,m-1,t}{\eta }_{m-1,t}+{\eta }_{m,t},{\eta }_{m,t}\sim \mathcal{N}\left(0,{\sigma }_{m,t}^2\right).\hfill \end{array}$$

Note that for $i>1$, the i-th equation of this system is a TVP model where the residuals of the preceding $i-1$ equations are added as explanatory variables:

$$\begin{array}{c}\hfill y_{i,t}={\mathbf{x}}_t{\mathit{\beta }}_t^i+{\displaystyle \sum _{j=1}^{i-1}}{a}_{ij,t}{\eta }_{j,t}+{\eta }_{i,t},{\eta }_{i,t}\sim \mathcal{N}\left(0,{\sigma }_{i,t}^2\right),\end{array}$$

and all time-varying parameters follow a random walk as in the TVP model (1):

$$\begin{array}{cc}\hfill & {\beta }_{j,t}^i={\beta }_{j,t-1}^i+v_{ij,t},v_{ij,t}\sim \mathcal{N}\left(0,{\theta }_{ij}^{\beta }\right),\mathrm{for}i=1,\dots ,m,\mathrm{and}j=1,\dots ,mp+1,\hfill \\ \hfill & a_{ij,t}=a_{ij,t-1}+w_{ij,t},w_{ij,t}\sim \mathcal{N}\left(0,{\theta }_{ij}^a\right),\mathrm{for}i=1,\dots ,m,\mathrm{and}j=1,\dots ,i-1,\hfill \end{array}$$

with initial values ${\beta }_{j,0}^i\sim \mathcal{N}\left({\beta }_{ij}^{\beta },{\theta }_{ij}^{\beta }\right)$ and $a_{ij,0}\sim \mathcal{N}\left({\beta }_{ij}^a,{\theta }_{ij}^a\right)$. Here, ${\beta }_{j,t}^i$ denotes the jth element of the vector ${\mathit{\beta }}_t^i$.

To achieve shrinkage for each VAR coefficient ${\beta }_{j,t}^i$ as well as for each Cholesky factor $a_{ij,t}$, we proceed as in Section 2 and introduce shrinkage priors for the initial expectations ${\beta }_{ij}^{\beta }$ and ${\beta }_{ij}^a$ as well as the variances ${\theta }_{ij}^{\beta }$ and ${\theta }_{ij}^a$. We do this independently for each equation of the system. Within each equation, the ${\beta }_{ij}^{\beta }$s and ${\beta }_{ij}^a$s are assumed to follow independent shrinkage priors to allow for flexibility in the prior structure, and similarly for ${\theta }_{ij}^{\beta }$ and ${\theta }_{ij}^a$:

$$\begin{array}{cc}\hfill & {\beta }_{ij}^x\sim \mathcal{N}\left(0,{\varphi }_i^{\tau ,x}{\check{\tau }}_{ij}^{x,2}/{\check{\lambda }}_{ij}^{x,2}\right),{\check{\tau }}_{ij}^{x,2}\sim \mathcal{G}\left(a_i^{\tau ,x},1\right),{\check{\lambda }}_{ij}^{x,2}\sim \mathcal{G}\left(c_i^{\tau ,x},1\right),{\varphi }_i^{\tau ,x}=2c_i^{\tau ,x}/\left({\lambda }_{B,i}^{x,2}{a}_i^{\tau ,x}\right),\hfill \\ \hfill & {\sqrt{\theta }}_{ij}^x\sim \mathcal{N}\left(0,{\varphi }_i^{\xi ,x}{\check{\xi }}_{ij}^{x,2}/{\check{\kappa }}_{ij}^{x,2}\right),{\check{\xi }}_{ij}^{x,2}\sim \mathcal{G}\left(a_i^{\xi ,x},1\right),{\check{\kappa }}_{ij}^{x,2}\sim \mathcal{G}\left(c_i^{\xi ,x},1\right),{\varphi }_i^{\xi ,x}=2c_i^{\xi ,x}/\left({\kappa }_{B,i}^{x,2}{a}_i^{\xi ,x}\right),\hfill \end{array}$$

(31)

where $x=\beta$ for the VAR-coefficients and $x=a$ for the elements of ${\mathit{A}}_t$. Following Section 2.4, the priors for the global shrinkage parameters in the ith equation read

$$\begin{array}{cc}\hfill & {\lambda }_{B,i}^{x,2}|a_i^{\tau ,x},c_i^{\tau ,x}\sim \mathrm{F}\left(2a_i^{\tau ,x},2c_i^{\tau ,x}\right),2a_i^{\tau ,x}\sim \mathcal{B}\left({\alpha }_{a^{\tau }},{\beta }_{a^{\tau }}\right),2c_i^{\tau ,x}\sim \mathcal{B}\left({\alpha }_{c^{\tau }},{\beta }_{c^{\tau }}\right),\hfill \\ \hfill & {\kappa }_{B,i}^{x,2}|a_i^{\xi ,x},2c_i^{\xi ,x}\sim \mathrm{F}\left(2a_i^{\xi ,x},2c_i^{\xi ,x}\right),2a_i^{\xi ,x}\sim \mathcal{B}\left({\alpha }_{a^{\xi }},{\beta }_{a^{\xi }}\right),2c_i^{\xi ,x}\sim \mathcal{B}\left({\alpha }_{c^{\xi }},{\beta }_{c^{\xi }}\right).\hfill & \hfill \end{array}$$

(32)

5.2. A Brief Sketch of the TVP-VAR-SV MCMC Algorithm

Our algorithm exploits the aforementioned unitriangular decomposition to estimate the model parameters equation-by-equation. Due to the prior structure introduced in (31), the estimation of ${\mathit{\beta }}_t^i$ and the $a_{ij,t}$’s is separated into two blocks, with the algorithm cycling through the m equations, alternating between sampling ${\mathit{\beta }}_t^i$ conditional on ${\mathsf{\Sigma }}_t$ and sampling the $a_{ij,t}$s and $d_{i,t}$s conditional on the VAR coefficients ${\mathit{\beta }}_t^i$. Given a set of initial values, the algorithm repeats the following steps:

Algorithm 2. MCMC inference for TVP-VAR-SV models under the triple gamma prior.

Choose starting values for all global and local shrinkage parameters in prior (31) for each equation and repeat the following steps:

For $i=1,\dots ,m$, update all the unknowns in the ith equation: (a)  Conditional on ${\mathit{A}}_t$ and ${\mathit{D}}_t$, create ${\check{y}}_{i,t}=y_{i,t}-{\sum }_{j=1}^{i-1}{a}_{ij,t}{\eta }_{j,t}$ and define the following TVP model: $${\check{y}}_{i,t}={\mathbf{x}}_t{\mathit{\beta }}_t^i+{\eta }_{i,t},{\eta }_{i,t}\sim \mathcal{N}\left(0,{\sigma }_{i,t}^2\right).$$ Apply Algorithm 1 (sans the step for the variance of the observation equation) to this univariate TVP model, to draw from the conditional posterior distribution of the time-varying VAR-coeffcients ${\mathit{\beta }}_t^i,$ for $t=0,\dots ,T,$ their initial expectations ${\beta }_{ij}^{\beta }$, the process variances ${\theta }_{ij}^{\beta }$, the local shrinkage parameters ${\check{\tau }}_{ij}^{\beta ,2},{\check{\lambda }}_{ij}^{\beta ,2},{\check{\xi }}_{ij}^{\beta ,2},{\check{\kappa }}_{ij}^{\beta ,2}$, as well as the global shrinkage parameters ${\lambda }_{B,i}^{\beta ,2},{\kappa }_{B,i}^{\beta ,2},a_i^{\tau ,\beta },c_i^{\tau ,\beta },a_i^{\xi ,\beta }$, and $c_i^{\xi ,\beta }$. (b)  For $i>1$, create $y_{i,t}^{🟉}=y_{i,t}-{\mathbf{x}}_t{\mathit{\beta }}_t^i$, conditional on ${\mathit{\beta }}_t^i$, and define the following TVP model: $$y_{i,t}^{🟉}={\displaystyle \sum _{j=1}^{i-1}}{a}_{ij,t}{\eta }_{j,t}+{\eta }_{i,t},{\eta }_{i,t}\sim \mathcal{N}\left(0,{\sigma }_{i,t}^2\right),$$ where the residuals from the previous $i-1$ equations, $({\eta }_{1,t},\dots ,{\eta }_{i-1,t})$, are used as explanatory variables and no intercept is present. Apply Algorithm 1 to this univariate TVP model, to sample the volatilities ${\sigma }_{i,t}^2$ and the time-varying coefficients $a_{ij,t}$ in the ith row of ${\mathit{A}}_t$ for $t=0,\dots ,T$ from the respective conditional posteriors, as well as the initial expectations ${\beta }_{ij}^a$, the process variances ${\theta }_{ij}^a$, the local shrinkage parameters ${\check{\tau }}_{ij}^{a,2},{\check{\lambda }}_{ij}^{a,2},{\check{\xi }}_{ij}^{a,2},{\check{\kappa }}_{ij}^{a,2}$ and the global shrinkage parameters ${\lambda }_{B,i}^{a,2},{\kappa }_{B,i}^{a,2},a_i^{\tau ,a},c_i^{\tau ,a},a_i^{\xi ,a}$, and $c_i^{\xi ,a}$.

In the following applications, we run our algorithm for $M=200,000$ iterations, discarding the first $100,000$ iterations as burn-in, and then keeping the output of one every 100 iterations.

5.3. Illustrative Example with Simulated Data

To illustrate the merit of our methodology in the context of TVP-VAR-SV models, we simulate data from two TVP-VAR-SV models with $T=200$ points in time, $p=1$ lags and $m=7$ equations, with varying degrees of sparsity. In the dense regime, approximately 30% of the values of $\mathit{\beta }$ and $\mathit{\theta }$ (here referring to the means of the initial states and the variances of the innovations as defined in Section 2, respectively) are truly zero, while in the sparse regime approximately 90% are truly zero. We show results for the triple gamma prior, the Horseshoe prior, the double gamma and the Lasso.

Regarding the priors on the hyperparameters, we use prior (32) with ${\alpha }_{a^{\tau }}={\alpha }_{c^{\tau }}={\alpha }_{a^{\xi }}={\alpha }_{c^{\xi }}=1$ and ${\beta }_{a^{\tau }}={\beta }_{c^{\tau }}={\beta }_{a^{\xi }}={\beta }_{c^{\xi }}=6$ for the triple gamma. The probability density function of the corresponding beta prior is monotonically increasing, with a maximum at $0.5$. This prior places positive mass in a neighborhood of the Horseshoe, but allows for more flexibility. In practice, placing a prior on the spike and slab parameters of the triple gamma, instead of fixing them to 0.5 as in the Horseshoe, allows us to learn the shrinkage profile from the data, including asymmetric profiles.

We assume that the global shrinkage parameters ${\lambda }_{B,i}^{\beta ,2},{\kappa }_{B,i}^{\beta ,2},{\lambda }_{B,i}^{a,2}$, and ${\kappa }_{B,i}^{a,2}$ follow a $\mathrm{F}\left(1,1\right)$ distribution for the Horseshoe prior which corresponds to the prior in Carvalho et al. (2009) and a $\mathcal{G}\left(0.001,0.001\right)$ distribution for the Lasso and the double gamma prior, as suggested in Belmonte et al. (2014) and Bitto and Frühwirth-Schnatter (2019). Concerning the spike parameters $a_i^{\tau ,a},a_i^{\xi ,a},a_i^{\tau ,\beta }$, and $a_i^{\xi ,\beta }$ of the double gamma, we employ a rescaled beta prior to force them to be smaller than $0.5$. Specifically, we use a $B(4,6)$ prior which places most of its mass between $0.05$ and $0.4$, a range that Bitto and Frühwirth-Schnatter (2019) have found to induce desirable shrinkage characteristics.

Figure 5 shows the posterior path of a permanently non-significant state, that is a state where the true ${\beta }_{j,t}^i=0$ for $t=1,\dots ,T$, in the sparse regime. The entire set of states for the triple gamma prior can be found in Appendix C. Note that, while the zero line is contained in the 95% posterior credible interval for all priors, said interval is thinner under the triple gamma prior and the double gamma prior than under the Lasso and the Horseshoe prior.

We calculate the posterior inclusion probabilities based on the thresholding approach introduced in Section 3.2, comparing the triple gamma prior to widely used special cases. In a variance selection context, the posterior inclusion probabilities reflect the uncertainty on whether a state should be time varying or constant over time. Figure 6 shows the posterior inclusion probabilities for the variance of the innovations (${\theta }_{ij}^{\beta }$’s) under four different shrinkage priors, for the sparse and the dense scenario, respectively. The cells are shaded in gray when the corresponding true state parameter is time-varying (${\theta }_{ij}^{\beta }\ne 0$), while the background is white when the corresponding true state parameter is not time-varying (${\theta }_{ij}^{\beta }=0$). In this simulated example, the posterior inclusion probabilities under the triple gamma prior are consistently higher for the variances that are actually different from 0, even when they are very small. This outcome is in line with the analytical results derived in Section 2.2, which show that the tails of the triple gamma prior are heavier than those of the other priors.

5.4. Modeling Area Macroeconomic and Financial Variables in the Euro Area

Our application investigates a subset of the area wide model of the European Union of Fagan et al. (2005), which comprises quarterly macroeconomic data spanning from 1970 to 2017. We include seven of the variables present in the dataset, namely real output (YER), prices (YED), short-term interest rate (STN), investment (ITR), consumption (PCR), exchange rate (EEN) and unemployment (URX). A more detailed description of the data and the transformations performed to make the time series stationary can be found in

Table A1. Data overview.

Variable	Abbreviation	Description	Tcode
Real output	YER	Gross domestic product (GDP) at market prices in millions of Euros, chain linked volume, calendar and seasonally adjusted data, reference year 1995.	1
Prices	YED	GDP deflator, index base year 1995. Defined as the ratio of nominal and real GDP.	1
Short-term interest rate	STN	Nominal short-term interest rate, Euribor 3-month, percent per annum	2
Investment	ITR	Gross fixed capital formation in millions of Euros, chain linked volume, calendar and seasonally adjusted data, reference year 1995.	1
Consumption	PCR	Individual consumption expenditure in millions of Euros, chain linked volume, calendar and seasonally adjusted data, reference year 1995.	1
Exchange rate	EEN	Nominal effective exchange rate, Euro area-19 countries vis-à-vis the NEER-38 group of main trading partners, index base Q1 1999.	1
Unemployment	URX	Unemployment rate, percentage of civilian work force, total across age and sex, seasonally adjusted, but not working day adjusted.	2

Note: Data was retrieved from https://eabcn.org/page/area-wide-model. Tcode $=1$ indicates that differences of logs were taken, while Tcode $=2$ implies that the raw data was used.

in Appendix D. To stay in line with the literature, for example, Feldkircher et al. (2017), we estimate a TVP-VAR-SV model with $p=2$ lags on all endogenous variables. The hyperparameter choices are the same as in Section 5.3. As in the example with simulated data, we run the algorithm for $M=\mathrm{200,000}$ iterations, discarding the first 100,000 iterations as burn-in, and then keeping the output of one every 100 iterations.

Figure 7 and Figure 8 display the posterior inclusion probabilities for the means of the initial states and the innovation variances of the VAR coefficients, respectively. A few things about Figure 7 are noteworthy. First, the posterior inclusion probabilities on the diagonal, meaning those belonging to the parameter of each equation’s own autoregressive term, appear to be those that are the highest, while off diagonal elements are more likely to be excluded. Second, the equation for the short-term interest rate is characterized by a large amount of parameters with a high inclusion probability, across all priors. Third, the first lag tends to have higher posterior inclusion probabilities than the second lag, which is in line with the literature. In most cases, the triple gamma prior can be seen to have either the largest or the smallest posterior inclusion probability compared to the other priors. This can be seen as a reflection of the fact that the triple gamma prior places more mass on the edges of the shrinkage profile, as illustrated in Section 3.

Now, we shift our focus to the posterior inclusion probabilities for the ${\theta }_{ij}^{\beta }$’s plotted in Figure 8. Compared to the means of the inital states, almost all inclusion probabilities are essentially zero. This lack of variability is unsurprising, as it is well known (see, e.g., Feldkircher et al. (2017)) that stochastic volatility in a TVP-VAR model for macroeconomic variables can explain a large part of the variability in the data. However, the triple gamma prior appears to allow posterior distributions that place slightly more mass on models with some time variation, in particular with respect to the financial variables.

Figure 9 and Figure 10 display the posterior median of ${\beta }_{ij}^{\beta }$ and $\left|{\sqrt{\theta }}_{ij}^{\beta }\right|$, respectively. Here the triple gamma can be seen to be quite conservative, both in terms of which parameters to include, as well as their magnitude. In particular the medians of the $\left|{\sqrt{\theta }}_{ij}^{\beta }\right|$ are interesting, as they are closest to zero under the triple gamma prior, despite having the highest posterior inclusion probabilities among all considered priors.

In Figure A3 and Figure A4 in Appendix D, all the posterior paths of ${\mathsf{\Phi }}_{1,t}$ and ${\mathsf{\Phi }}_{2,t}$ under the triple gamma prior are shown.

6. Conclusions

In the present paper, shrinkage for time-varying parameter (TVP) models was investigated within a Bayesian framework with the goal to automatically reduce time-varying parameters to static ones, if the model is overfitting. This goal was achieved by suggesting the triple gamma prior as a new shrinkage priors for the process variances of varying coefficients, extending previous work using spike-and-slab priors, the Bayesian Lasso, or the double gamma prior. The triple gamma prior is related to the normal-gamma-gamma prior applied for variable selection in highly structured regression models (Griffin and Brown 2017). It contains the well-known Horseshoe prior as a special case, however it is more flexible, with two shape parameters that control concentration at zero and the tail behaviour. This leads to a BMA-type behaviour which allows not only variance shrinkage, but also variance selection.

In our application, we considered time-varying parameter VAR models with stochastic volatility. Overall, our findings suggest that the family of triple gamma priors introduced in this paper for sparse TVP models is successful in avoiding overfitting, if coefficients are, indeed, static or even insignificant. The framework developed in this paper is very general and holds the promise to be useful for introducing sparsity in other TVP and state space models in many different settings.

A number of extensions seem to be worth pursuing. First of all, the triple gamma prior is relevant not only for TVP models, but for any model containing variance parameters such as random-effect models or Bayesian p-splines models (Scheipl and Kneib 2009). Second, in particular, in ultra-sparse settings, modifications of the triple gamma prior seem sensible. Currently, the hyperprior for the global shrinkage parameter of the triple gamma prior is selected in a way that it implies a uniform prior on “model size”. A generalization of Theorem 3 would allow the choice of hyper priors that induce higher sparsity. Furthermore, in the variable selection literature, special priors such as the Horseshoe+ (Bhadra et al. 2017a) were suggested for very sparse, ultra-sparse high dimensional settings. Exploiting once more the non-centered parametrization of a state space model, it is straightforward to extend this prior to variance selection using following hierarchical representation:

$$\begin{array}{c}\hfill {\sqrt{\theta }}_j|{\kappa }_j^2,{\xi }_j^2\sim \mathcal{N}\left(0,\frac{2}{{\kappa }_B^2}{\kappa }_j^2{\xi }_j^2\right),{\kappa }_j\sim t_1,{\xi }_j\sim t_1.\end{array}$$

We leave these extensions for future research.

Finally, an important limitation of our approach is that shrinking a variance toward zero implies that a coefficient is fixed over the entire observation period of the time series. In future research we will investigate dynamic shrinkage priors (Kalli and Griffin 2014; Kowal et al. 2019; Ročková and McAlinn 2020) where coefficients can be both fixed and dynamic.