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.
69,688 characters · 13 sections · 53 citation commands
The Dynamic Triple Gamma as a Shrinkage Process for Time-Varying Parameter Models
\def\spacingset#1{ {#1}} \spacingset{1}
\if11 \fi
\if01 {
} \fi
{\it Keywords:} Bayesian Inference; Hierarchical Priors; Horseshoe Prior; Variable Selection; Markov chain Monte Carlo;
\spacingset{1.9}
It has become evident in recent years that time-series methods must be able to adapt to large, sudden changes in underlying system dynamics. However, allowing too much flexibility can lead to overfitting and poor out-of-sample performance. Therefore, a balance must be struck between local adaptivity and regularization. Time-varying parameter (TVP) models provide a natural framework for exploring this balance, as they are inherently flexible and thus require proper regularization.
In the Bayesian paradigm, regularization is often achieved through shrinkage priors, which are mostly designed with variable selection in mind. However, regularization in TVP models is inherently a problem of variance selection, as it necessitates pulling innovation variances toward zero. Within the TVP framework, fru-wag:sto addressed this limitation by re-casting variance selection as a variable selection problem using a non-centered parameterization of the standard TVP model. In this reparameterized state-space model, shrinkage priors can be adapted to handle variance selection with minimal modification.
This approach to shrinkage is inherently \lq\lq static\rq\rq , as the amount of shrinkage applied to each state never changes over time. Therefore, it may simultaneously suffer from overfitting and underfitting. To see this, consider the well-known Nile data set from kraus1956graphs displayed in Figure (ref). It consists of 100 measurements of the annual flow of the river Nile at Aswan, with an apparent change-point near 1898, most likely due to a change in rainfall patterns. Such a change-point problem necessitates local adaptivity, as it requires a limited number of large jumps (i.e. innovations with high variance) followed and preceded by periods of innovations with low variance. The left-hand panel of Figure (ref) plots the posterior of a time-varying mean model\footnote{Details on the exact model specification and priors used can be found in Appendix (ref).} using an approach with a constant innovation variance. The aforementioned underfitting is apparent around the change-point, as the shift downwards is quite gradual, where it should be fairly sharp. Overfitting, on the other hand, can be seen to the left and right of the change-point, where it seems a lot of noise is tarnishing the estimate of the mean of the time series. Contrast this with the right-hand panel of Figure (ref), which graphs the posterior of a time-varying mean model under the proposed dynamic shrinkage method, which nicely picks up the change point while not suffering from the aforementioned overfitting in the less volatile periods.
This limitation of static shrinkage approaches has been widely recognized, leading to a wealth of literature that aims to develop dynamic shrinkage. An early discussion of dynamic linear models with time-varying variance can be found in pet-etal:dyn. Further pioneering work in this direction was undertaken by koo-kor:for, who developed a time-varying Bayesian model averaging scheme and nak-wes:bay_ana, who used a latent threshold approach to induce time-varying sparsity. Many other approaches are based on the idea of generalizing a shrinkage prior designed for variable selection to handle dynamic shrinkage. Due to its attractive theoretical properties and interpretability, the spike-and-slab approach geo-mcc:var is a popular candidate for extension, with work in this direction done by hub-etal:sho, uri-hed:dyn, and roc-mca:dyn, among others. Another strand of the literature is defined by a continuous shrinkage prior serving as the basis for generalization into a shrinkage process prior, with examples including work by kal-gri:tim, kow-etal:dyn, irie2019bayesian and hub-pfa:dyn. Recent machine learning approaches include the mixture approach by hau-etal:fas or the tree-based contribution of hau-etal:tree.
This paper generalizes the triple gamma shrinkage prior of cad-etal:tri to a stochastic process that also captures dependence in the amount of shrinkage across time. The triple gamma’s spike captures the belief that most innovations are near zero, while its heavy tails still allow for large jumps. Dependence enables volatility information to be shared across time. Since the triple gamma encompasses many popular shrinkage priors as special or limiting cases, this defines a class of stochastic shrinkage processes with diverse shrinkage priors as marginals. The most prominent example is an innovative shrinkage process prior that captures dependence in the amount of shrinkage across time and features the famous horseshoe prior car-etal:hor as marginal distribution; a combination that, to the best of our knowledge, is novel within the existing literature. Despite these features, the full conditional distributions of the posterior remain mostly in closed form, simplifying estimation.
The rest of the paper is structured as follows: Section (ref) lays groundwork by providing a brief overview how global-local (GL) shrinkage priors can be applied to shrinking innovation variances in TVP models toward zero. Section (ref) introduces the key contribution, namely a TVP model with time-varying innovation variances following the innovative shrinkage process defined by a dynamic triple gamma. Section (ref) discusses properties of the dynamic triple gamma prior in more detail, while Section (ref) develops an efficient Markov chain Monte Carlo (MCMC) sampler to estimate TVP models under the dynamic triple gamma. Section (ref) showcases the dynamic triple gamma by fitting a Cholesky stochastic volatility (SV) model to the returns of the components of the EURO STOXX 50 index and demonstrating the out-of-sample forecasting performance via log predictive density scores (LPDS). Finally, Section (ref) concludes.
To set notation for the following discussion, we introduce the standard TVP model:
where the $(1\times d)$-dimensional row vector ${\mathbf x}_t $ contains the regressors and $y_t$ is the regressand, for $t = 1, \dots T$. The innovations' covariance matrix ${\mathbf{{\mathbf{Q}}}} =\mbox{\rm Diag}\left(\theta_1, \ldots, \theta_{d}\right)$ is diagonal, meaning that each coefficient follows a random walk with innovation variance $\theta_j$, i.e. for $j=1,\ldots, d$, ${\beta}_{jt} = {\beta}_{j,t-1} + w_{jt}$, where $w_{jt} \sim \mathcal{N}\left(0 ,\theta_j\right)$. The system is initialized by assuming that the initial value $ {\boldsymbol{\beta}}_{0}$ is randomly drawn from the following normal distribution: $ {\boldsymbol{\beta}}_{0} \sim \mathcal{N} _{d}\left(\boldsymbol{\beta},{\mathbf{{\mathbf{Q}}}} \right)$, where $\boldsymbol{\beta} = (\beta_1, \dots, \beta_d)'$. The error variance can either be assumed to be homoscedastic or follow the SV model of jac-etal:bay_ana where $\sigma^2_t = \mbox{\rm e} ^ {h_t}$ and $h_t$ follows an AR(1) process.
These models are inherently flexible, as the regression coefficients can change at each time point $t \in \{1, \dots, T\}$. One approach to regularizing such models is to shrink the innovation variances $\theta_1, \dots, \theta_d$ toward zero, as smaller variances indicate that innovations are, on average, closer to zero. In the extreme case, where $\theta_j = 0$, we have $\beta_{jt} = \beta_{j,t-1} = \beta_j$ with probability 1. While GL shrinkage priors are primarily designed for variable selection, they can also be used to address variance selection in a TVP model. To this end, fru-wag:sto introduced following useful reparameterization of the standard TVP model:
where $\tilde{\beta}_{j0} \sim \mathcal{N}\left(0 ,1\right)$. The two parameterizations are linked through the simple transformation $\beta_{jt} = \beta_{j} + \sqrt\theta_j \Tilde{\beta}_{jt}$, for all $t = 1, \dots, T$. To see the benefit of this, consider placing a generic GL shrinkage prior on $\sqrt \theta_j$, where $\sqrt\theta_j | \lambda_j, \tau \sim \mathcal{N}\left(0, \tau\lambda_j\right)$, $\tau$ is the global parameter pulling all coefficients toward zero and $\lambda_j \sim p(\lambda_j)$ is the local parameter that allows individual coefficients to be different from zero. When such a prior is placed on $\sqrt\theta_j$, it implies the following GL variance shrinkage prior on $\theta_j$ based on the gamma distribution:
This duality is what makes this reparameterization so useful -- any GL prior originally desigend for variable selection can be placed on $\sqrt \theta_j$ and will automatically imply a shrinkage prior for the innovation variances $\theta_j$ in a TVP model. This is the context in which the triple gamma prior was introduced by cad-etal:tri.
In this paper, we extend the standard TVP model by allowing for dynamic innovation variances. Specifically, each innovation variance \( \theta_j \) is multiplied by a local, time-varying component \( \psi_{jt} \), resulting in \({\mathbf{w}}_{t} \sim \mathcal{N} _{d}\left({\mathbf{0}}, {\mathbf{{\mathbf{Q}}}}_t\right)\) with \({\mathbf{{\mathbf{Q}}}}_t = \mbox{\rm Diag}\left(\theta_1 \psi_{1t}, \ldots, \theta_{d} \psi_{d t}\right)\). It will be shown that the triple gamma prior is also a natural choice as a shrinkage prior for \( \psi_{jt} \) in achieving dynamic shrinkage.
The triple gamma shrinkage prior $TG\left(a, c, \kappa\right)$ is a distribution with three hyperparameters, where $a$ controls the behavior at 0, $c$ controls the tail behavior and $\kappa$ controls global shrinkage. The eponymous representation of the triple gamma is a compound distribution consisting of three gamma random variables, meaning that $ X | a, c, \kappa \sim TG\left(a, c, \kappa\right)$ is equivalent to
It can also be represented in a form closer to the generic GL representation, using the connection between the gamma and the F distribution:
It's variable selection "twin" is equivalent to the normal-gamma-gamma (NGG) prior of gri-bro:hie:
or, in GL form:
This prior is appealing for two reasons. First, the triple gamma prior either is equivalent to or encompasses many popular GL priors as special or limiting cases. As such, it can be considered a fairly generic choice of shrinkage prior, as many popular choices result for specific values of $a$ and $c$, the most prominent example being the horseshoe prior for $a=c=0.5$. This is evident from (ref), since $W \sim \mbox{\rm F}\left(1,1\right)$ implies that $\sqrt{W} \sim t_{1}$ follows a Cauchy distribution. Second, it is mathematically well understood, in particular with respect to it's behavior around the origin and in the tails. This results from the marginal density being available in closed form:
where $\phi = \frac{2c}{\kappa a}$ and $U \left(a,b,z\right)$ is the confluent hypergeometric function of the second kind:
The properties of the confluent hypergeometric function can be used to characterize the pole and tail behavior cad-etal:tri. The triple gamma prior has an infinite pole around $0$ when $a \leq 0.5$ and this pole becomes more extreme the smaller $a$ becomes. Furthermore, the prior has polynomial tails, with smaller values of $c$ equating to heavier tails. This is further reinforced by the fact that the prior moments $\mathbb E(X^k|a, c, \kappa)$ only exist up to $k < 2c$. Interestingly, the horseshoe in particular results as the special case when $a = c = 0.5$. As such, it has the least mass around the origin of all triple gamma specifications with an infinite pole at $0$, while being the most heavy-tailed the triple gamma can be while still still having a finite first moment. If a reader is so inclined, they can therefore think of the triple gamma as a more flexible horseshoe prior that allows for more granular control over the behavior around the origin and in the tails. A visual demonstration of this can be found in Figure (ref), which shows the log density of the variable selection twin of the triple gamma prior, defined on the entire real line. The figure illustrates that specifications with the same $a$ behave very similarly around the origin, while those with the same $c$ exhibit similar behavior in the tails. Finally, $ \tau=2/\kappa$ acts as the global shrinkage component, with larger values of $\kappa$ corresponding to more overall shrinkage.
To introduce dynamic shrinkage through the triple gamma prior, we begin with an exchangeable prior on the innovations and then extend it in Section (ref) to the case that allows for dependence across time. As mentioned above, we extend the TVP model defined in (ref) - (ref), with a triple gamma prior $\theta_j \sim TG ( a^\xi, c^\xi, \kappa_{B}^2)$ on the innovation variance, by multiplying $\theta_j$ with a local component $\psi_{jt}$. If $\psi_{jt}<1$, then global shrinkage induced for coefficient $j$ by $\theta_j$ will be enforced, while $\psi_{jt}>>1$ will apriori allow for large local innovations $w_{jt} $. While other choices for the distribution of $\psi_{jt}$ are feasible, the easiest is to use following exchangeable prior:
This choice leads to $w_{jt}^2| a_j, c_j, \theta_j$ being distributed as \rm i.i.d.\ $TG(a_j, c_j, 2/\theta_j)$, which follows immediately from representation (ref). A similar exchangeable approach to dynamic shrinkage can be seen in hub-etal:ind, where they place exchangeable horseshoe priors on the innovations.
The structure of this setup can be viewed as a doubly hierarchical GL prior on the innovations. The $\theta_j$'s serve as the global component within each state, shrinking all innovations $w_{jt}$ within state $j$ (for $t = 1, \dots, T$), while the $\psi_{jt}$'s act as the local components. Further, the $\theta_j$'s themselves are hierarchical, with following representation involving local components $\xi_j$:
Thus, $\theta_j$ controls the variability within state $j$, while $\kappa_{B}^2$ regulates the range of $\theta_j$ across states.
While the exchangeable triple gamma innovations introduced in Section (ref) are interesting in and of themselves, they do not yet allow information on the amount of shrinkage to be shared across time points. However, the assumption that high volatility in one period is correlated with higher volatility in the next seems fairly reasonable. This requires the local shrinkage factors \( \psi_{jt} \) to follow a stochastic shrinkage process prior, similar to kal-gri:tim and kow-etal:dyn. At the same time, we aim to preserve the triple gamma prior as the marginal distribution of this process, owing to its well-understood mathematical properties, briefly reviewed in Section (ref). This motivates the introduction of the dynamic triple gamma, a novel stochastic shrinkage process prior.
We first begin by stating an alternate, novel representation of the triple gamma prior. We use properties of the $\mbox{\rm F}\left(2 a_j, 2 c_j\right)$-distribution to represent $\psi_{jt}$ in (ref) as\footnote{To see this, note that the F-distribution has a representation as the ratio of two independent gamma distributions, $\mbox{\rm F}\left(2 a_j, 2 c_j\right) \,{\buildrel d \over =}\, \frac{ \mathcal{G}{{\small \left(a_j,1\right)}}/a_j}{ \mathcal{G}{{\small \left(c_j,1\right)}}/c_j} \,{\buildrel d \over =}\, \frac{ \mathcal{G}{{\small \left(a_j,1\right)}}c_j/a_j}{ \mathcal{G}{{\small \left(c_j,1\right)}}}$, which can be written as $\mbox{\rm F}\left(2 a_j, 2 c_j\right) \,{\buildrel d \over =}\, \mathcal{G}^{-1} \left(c_j, \mathcal{G}{{\small \left(a_j,a_j/c_j\right)}}\right)$.}
Interestingly, this alternative representation also leads to a novel sampling scheme for the posterior under triple gamma priors, as can be seen in the MCMC sampler introduced in Section (ref).
The most innovative step in constructing the dynamic triple gamma prior is replacing the exchangeable prior distribution for $ \lambda_{jt} $ in (ref) with a stochastic process that introduces autocorrelation, while preserving the marginal distribution $ \lambda_{jt} \mid a_j, c_j \sim \mathcal{G}{{\small \left(a_j, a_j / c_j\right)}} $ as in the exchangeable case. This ensures the marginal distribution of \( w_{jt} \) remains unchanged. To define such a process, we place the following latent autoregressive gamma (ARG) process on \( \lambda_{jt} \):
which was introduced by gou-jas:aut as an AR-type process for non-negative time series observations. This definition induces the marginal distribution $\lambda_{jt}|a_j, c_j \sim \mathcal{G}{{\small \left(a_j, {a_j}/{c_j}\right)}}$. At the same time, it leads to a process that is autocorrelated due to the dependence of $\lambda_{jt}$ on $\lambda_{j,t-1}$ through the contemporaneous value of $\kappa_{jt}$, with the strength of the autocorrelation governed by the parameter $\rho_j$. More precisely, this construction causes the expectation to be linearly dependent on the previous value, i.e. $\mathbb E[\lambda_{jt}|\lambda_{j,t-1}] = \rho_j \lambda_{j,t-1}$. Note that, as $\lambda_{jt}$ has support on $\mathbb R_{+}$, this only allows for values of $\rho_j$ in $(0, 1)$. To clarify the hierarchical structure of the dynamic triple gamma prior, Figure (ref) provides a visualization, and the complete construction is presented below:
Analogous to the exchangeable case, this construction forms a doubly hierarchical GL prior. $\theta_j$ still acts as the global component within each state $j$, pulling all innovation variances toward zero, while the hierarchical prior placed on $\theta_j$ remains unchanged, inducing a joint prior over all $\theta_j$'s. The key distinction is that the local components, the $\psi_{jt}$'s, now explicitly model dependence across time. Despite this, the marginal distribution of $w_{jt}^2$ remains triple gamma.
From a Bayesian perspective, hyperpriors must be specified for the model parameters $\theta_1, \dots, \theta_d$, ${\beta}_{1}, \dots, {\beta}_{d}$, and, crucially, $\rho_1, \dots, \rho_d$. Since the dynamic triple gamma reduces to the static case when $\psi_{jt} = 1$, we follow cad-etal:tri for the parameters shared with this model, namely $\theta_1, \dots, \theta_d$, ${\beta}_{1}, \dots, {\beta}_{d}$ and the error variance $\sigma^2_t$. A key element of this approach is employing triple gamma priors ${\beta}_{j}^2 \sim TG(a^\tau, c^\tau, \lambda_B^2)$ and $\theta_j \sim TG(a^\xi, c^\xi, \kappa_{B}^2)$ for $j = 1, \dots, d$, allowing the dynamic triple gamma to differentiate between covariates with time-varying, constant, or no effect. Specifically, when $\theta_j = 0$, ${\beta}_{jt} = {\beta}_{j,t-1}=\beta_j$ with probability 1, reducing the effect to a constant value. The triple gamma prior, equivalent to an NGG prior on \( {\beta}_{j} \), also pulls static coefficients toward 0. If, additionally, $\beta_j = 0$, the covariate's effect is constantly zero, effectively excluding it from the model. For additional details on these priors see Appendix (ref).
A prior must also be specified for the $\rho_j$'s, model parameters unique to the dynamic case. Since gou-jas:aut demonstrate that the ARG process degenerates as $\rho_j$ approaches 1, we adopt a generalized beta prior of the first kind (GB1) for $\rho_j$ to mitigate this issue. Its density function is:
valid for $0 < \rho_j < b_\rho$, where $a_\rho \in \mathbb{R}$, $b_\rho > 0$, $\alpha_\rho > 0$, and $\beta_\rho > 0$. This prior generalizes the beta distribution in two ways: it allows the support to range from 0 to $b_\rho$, while the parameter $a_\rho$ controls the power applied to the argument. For bounding the prior away from 1, the first property is especially useful, as it restricts the prior's support. For instance, in the application presented in Section (ref), we set $b_\rho = 0.95$, limiting the range to $(0, 0.95)$.
The dynamic triple gamma shares a particularly close connection to two approaches previously proposed in the literature. First, the normal-gamma autoregressive (NGAR) process of kal-gri:tim uses essentially the same recursive process as defined by (ref)-(ref), albeit with a different parameterization. Through this, they define a strictly stationary process directly on the ${\beta}_{jt}$'s, such that the resulting marginal distribution is a normal-gamma shrinkage prior gri-bro:inf. In our opinion, it seems more natural to induce dynamic shrinkage on the growth process $w_{jt}={\boldsymbol{\beta}}_{jt}- {\boldsymbol{\beta}}_{j,t-1}$, as we do in the present paper, as this allows more flexibility in accommodating small and large jumps in ${\boldsymbol{\beta}}_{jt}$ at the same time.
Second, there exists a (slightly less obvious) connection to the dynamic shrinkage process proposed by kow-etal:dyn, where dependence is induced in the innovations through an AR(1) process on the log of the innovation variance:
where the $\eta_{jt}$ arise from a Z-distribution\footnote{For details on the various non-standard distributions used in this section see Appendix (ref).} and $\phi_j^h$ acts as an autocorrelation parameter. To see how the proposed approach and that of kow-etal:dyn are related, note that the triple gamma prior can also be represented as cad-etal:tri:
where $\mathcal{BP}\left(a_j,c_j\right)$ is the beta-prime distribution, $\tilde{\psi}_{jt} = \frac{a_j} {c_j}\psi_{jt}$, and $\phi_{j} = \theta_j \frac{ c_j}{a_j}$. Based on this, it is possible to verify the following new representation of the triple gamma:
Specifically, from representation ((ref)) we find that $h_{jt}= \log (\phi_{j}\tilde{\psi}_{jt}) = \log \phi_{j} + \log \tilde{\psi}_{jt}$, where $\log \tilde{\psi}_{jt} \sim \mbox{\rm Z}\left(a_j,c_j,0,1\right)$ and therefore $h_{jt} | a_j,c_j,\theta_j \sim \mbox{\rm Z}\left(a_j,c_j,\log \phi_{j},1\right)$.\footnote{Note that $ X \sim \mathcal{BP}\left(a,c\right)$ implies $ Y= \log X \sim \mbox{\rm Z}\left(a,c,0,1\right)$ and therefore $Y + \mu \sim \mbox{\rm Z}\left(a,c,\mu,1\right) $.}
Essentially, the conditional distribution of $w_{jt}^2| a_j,c_j,h_{j,t-1},\mu_j^h, \phi_j^h$ proposed by kow-etal:dyn is equivalent to the triple gamma prior. This result is already interesting on its own, as it allows the properties of the triple gamma to shed light on the characteristics of the dynamic shrinkage process. It further implies that the two approaches coincide when the respective autocorrelation parameters ($\rho_j$ and $\phi_j^h$) are equal to $0$, as they both result in squared innovations that are exchangeable triple gamma, with the same parameters $a_j$ and $c_j$, while $\mu_j^h =\log \phi_j$. The key difference is that the marginal form of the proposed approach is again well-known, i.e. triple gamma. This allows the shrinkage characteristics of the entire process to be better understood, while sacrificing some tractability of the conditional distribution. There are also computational advantages, as the MCMC algorithm proposed in Section (ref) does not need to rely on mixture approximations as in kow-etal:dyn, which tend to be expensive to evaluate.
To gain further understanding of the dynamic triple gamma (DTG), one can rely on alternative representations of the process, which can be achieved by marginalizing out sets of latent variables, such as $\lambda_{j0}, \ldots, \lambda_{jT}$ or $\psi_{j1}, \ldots, \psi_{jT}$. For the first such representation, note that marginalizing ((ref)) over $\{\psi_{jt}\}$ leads to the following representation (see also cad-etal:tri):
Therefore, the innovations conditional on $\lambda_{jt}$ follow a Student-$t$ distribution with $2c_j$ degrees of freedom. This representation also has direct implications for the choice of $c_j$, as the innovations $w_{jt}$ have no finite moments whenever $c_j \leq 0.5$. If $c_j > 0.5$, the innovations have expectation zero. When $c_j > 1$, then the conditional variance of $w_{jt}|c_j, \theta_j$ is finite. From this it can be inferred that a smaller $c_j$ will lead to a process which allows for larger changes in the states a priori, whereas a larger $c_j$ will cause the changes in the states to be more regularized. This interpretation is in line with the effect of $c$ in the static triple gamma discussed in Section (ref).
Further insights into the behavior of the DTG prior can be gleaned by examining the transition density $p(\psi_{jt} | \psi_{j,t-1}, a_j, c_j, \rho_j)$ derived under the assumption that $\lambda_{j,t-1}$ comes from the stationary distribution of the ARG process. This density is available in closed, albeit not well-known, form in Theorem (ref). It allows for impulse response analysis, showcasing how the process responds to a shock while in the steady state. The proof can be found in Appendix (ref).
Theorem (ref) is useful both for building understanding of the properties of the DTG prior, as well as for the MCMC algorithm presented in Section (ref). Figure (ref) plots the transition density for various values of $\rho_j$, $c_j$ and $\psi_{j,t-1}$, highlighting the role both model parameters play in determining the dynamic shrinkage characteristics of the DTG prior. Starting with $\rho_j$, one can see that the degree to which the conditional density deviates from the marginal density for large values of $\psi_{j,t-1}$ depends strongly on the value of $\rho_j$, with values closer to $1$ corresponding to stronger deviation. Whereas all densities in the left-most column, where $\rho_j = 0.1$, hew closely to the marginal distribution, there exists significant variation in the right-most column, where $\rho_j = 0.9$. Despite this, it is interesting to note that fixing $a_j$ at 0.5 still causes all distributions to feature an infinite pole at the origin, irrespective of the value of $\psi_{j,t-1}$, with some distributions displaying multimodal behavior. This effect of $\rho_j$ aligns with the idea of it being a parameter that controls dependence, as, the larger $\rho_j$ is, the more large values of $\psi_{j,t-1}$ increase the probability of observing large values of $\psi_{jt}$.
Shifting focus to the effect of $c_j$ reveals that, as in the static triple gamma, $c_j$ controls the tails of the conditional density. Roughly speaking, smaller values of $c_j$ lead to heavier tails, an idea that is also reflected in the fact that the marginal distribution has no moments for $c_j < 0.5$. This idea is extended to $p(\psi_{jt} | \psi_{j,t-1}, a_j, c_j, \rho_j)$ in Theorem (ref), with the proof again found in Appendix (ref).
The implications of Theorem (ref) align neatly with the intuition gained from Figure (ref). The expectation is a weighted average of the expected value $c_j/(c_j - 1)$ of the marginal $\mbox{\rm F}\left(2a_j, 2c_j\right)$-distribution and an expression depending non-linearly on $\psi_{j,t-1}$. Therefore, it deviates further from the marginal expectation the closer $\rho_j$ is to $1$. The tail controlling behavior of $c_j$ becomes apparent if one lets $c_j$ approach $1$ from above, as $\lim_{c_j \downarrow 1} \mathbb E \left[ \psi_{jt} | \psi_{j,t-1}, a_j, c_j, \rho_j \right] = \infty$. On the other hand, the fact that values of $c_j > 1$ imply heavier regularisation of the tails is also evident, as the expectation, while monotonically increasing in $\psi_{j,t-1}$, is not unbounded. Specifically,
Another benefit of Theorem (ref) is that it allows for the derivation of conditional shrinkage profiles, in the spirit of car-etal:hor. Letting $\theta_j = 1$ for illustrative purposes, one can re-parameterize the conditional density of $w_{jt}$ in the following way, where $\psi_{jt} =1 /{\tau_{jt}} - 1 $:
Using the law of transformation of densities for $\tau_{jt} = \frac{1}{1 + \psi_{jt}}$ yields the following corollary.
$\tau_{jt}$ is confined to the range $(0, 1)$ and characterizes the shrinkage behavior, with $\tau_{jt}$ close to $1$ implying "total" shrinkage, as the variance of $w_{jt}$ then approaches 0 and $\tau_{jt}$ close to $0$ implying the opposite. Figure (ref) graphs the dynamic shrinkage profiles for various values of $\rho_j$, $c_j$ and $\psi_{j,t-1}$. The role of $\rho_j$ as the parameter driving the extent to which the conditional distribution deviates from the marginal distribution is again apparent. Similarly, $c_j$'s importance for the shape of the tails of the distribution is reinforced. In the top row, where $c_j = 0.5$, the result is a dynamic horseshoe, with infinite poles in both corners, irrespective of the magnitude of $\psi_{j,t-1}$, with larger values of $\psi_{j,t-1}$ allocating more prior mass toward $0$ and vice versa. In the bottom two rows, where $c_j \geq 1$, the infinite pole in the left corner disappears, implying some regularization in the tails, even for large values of $\psi_{j,t-1}$. This lines up nicely with the results derived from Theorem (ref), as the expected value of the conditional density only exists for $c_j > 1$.
Performing inference in such a highly parameterized model is challenging both from a theoretical as well as a computational perspective. These challenges stem almost exclusively from the latent non-Gaussian processes $\psi_{jt}$'s and the associated hyperpriors, which becomes evident when noting that the model simplifies to a linear Gaussian state space model when conditioning on the $\psi_{jt}$'s. As this is a well-understood class of models, we can fall back on algorithms that are fairly standard in the literature for the portions of the MCMC algorithm that deal with a conditionally Gaussian state space model. This includes, for example, the simulation smoother of mcc-etal:sim or the ancillarity-sufficiency interweaving strategy (ASIS) of bit-fru:ach to increase the sampling efficiency of $\theta_1, \dots, \theta_d$. Therefore, we focus on the parts of the algorithm that are non-standard here, by conditioning on the innovations $w_{jt} = {\beta}_{jt} - {\beta}_{j,t-1}$. We give a full overview of the MCMC scheme in Appendix (ref). Readers interested in estimating TVP models within the DTG framework can find implementations of the algorithms discussed here in the R package shrinkTVP, available on CRAN.
\paragraph*{Sampling under exchangeable triple gamma innovations}
To give an intuitive understanding of the MCMC sampler, we begin by discussing how to draw samples from a simple model, characterized by a priori exchangeable triple gamma innovations, i.e. $w^2_{jt} |\theta_j, a_j , c_j \mathop{\sim}\limits^{\mathrm{iid}} TG\left(a_j, c_j, 2/\theta_j\right)$. This is a special case of the DTG that results when $\rho_j = 0$, which causes all $d$ sets of latent variables $\kappa_{j1}, \dots, \kappa_{jT}$ in ((ref)) to be equal to 0, while the latent variables $\lambda_{j1}, \ldots, \lambda_{jT}$ in ((ref)) arise independently from $\lambda_{jt}|a_j, c_j \mathop{\sim}\limits^{\mathrm{iid}} \mathcal{G}{{\small \left(a_j, \frac{a_j}{c_j} \right)}}$ as in (ref).
Here, one can employ either the sampler proposed in cad-etal:tri or, based on representation (ref), utilize a generalized and innovative version of the sampler introduced for the horseshoe prior by mak-sch:sim. Applying standard derivations for conditionally conjugate densities involving the inverse gamma and gamma distributions yields:
for $t = 1, \dots, T$ and $j = 1, \dots, d$. This result is valuable in its own right, offering a method for sampling from the posterior of a model under the triple gamma prior without relying on the generalized inverse Gaussian distribution (GIG) as in cad-etal:tri. As the GIG is not as widely implemented in programming languages compared to the gamma and inverse gamma distributions, this outcome facilitates the broader application of the triple gamma prior.
\paragraph*{Sampling for fixed \texorpdfstring{$\rho_j>0$}{rho j >0}}
Moving from the case where the innovations are exchangeable triple gamma, the next simplest case is one where dependence is modeled via a fixed $\rho_j >0$. Naturally, this additionally necessitates the sampling of the latent variables $\kappa_{j1}, \dots, \kappa_{jT}$. Here, the hierarchical structure allows a straightforward Gibbs sampler to be implemented to sample from the posterior of all unknowns. The conditional density (ref) for $\psi_{jt}$ does not require any modification, as conditional on $\lambda_{jt}$, it is independent of $\kappa_{j1}, \dots, \kappa_{jT}$. The full conditional posterior for $\lambda_{jt}$ can, once again, be found through fairly standard derivations based on gamma and inverse gamma densities:\footnote{In several of the sampling steps discussed in this section, adjustments are necessary for the cases where $t \in \{0, 1, T\}$. To maintain brevity, these modifications can be found in Appendix (ref).}
where the case $\rho_j=0$ leads back to ((ref)). The conditional posterior $\kappa_{jt} | \lambda_{jt}, \lambda_{j,t-1}, a_j, c_j, \rho _j $ follows a discrete distribution, independently of all other state variables, with weights given by:
Sampling in this fully conditional setup works well, with good mixing and computations that are able to be performed fairly quickly. However, fixing $\rho_j$ a priori is a fairly strong assumption, particularly because it has a large influence on the properties of the DTG (see Section (ref)).
\paragraph*{Sampling with unknown \texorpdfstring{$\rho_j$}{rho j}}
As opposed to the sampling steps for $\psi_{jt}$, $\lambda_{jt}$ and $\kappa_{jt}$, the full conditional posterior of $\rho_j$ is not available in closed form, necessitating the use of a Metropolis-Hastings-within-Gibbs step. Implementing such a step within a fully conditional Gibbs loop in a straightforward fashion leads to poor mixing. To alleviate this issue, we make use of representations of the hierarchical structure with different sets of random variables marginalized out. For the Metropolis-Hastings (MH) step for $\rho_j$, we construct an approximate likelihood using Theorem (ref):
where $p(\psi_{jt} | \psi_{j,t-1}, a_j, c_j, \rho_j)$ is the density defined in (ref). This, combined with the prior density defined in Section (ref), gives us the posterior up to an unknown constant, allowing the construction of an MH step to generate samples from the conditional posterior $p(\rho_j| a_j, c_j,\psi_{j1}, \dots, \psi_{jT}) $. The important thing to note here is that this approximate likelihood is marginalized both w.r.t. $\lambda_{j0}, \dots, \lambda_{jT}$ and $\kappa_{j1},\dots, \kappa_{jT}$. Thus, to preserve the stationary distribution of the Markov chain, one needs to sample either $\lambda_{j0}, \dots, \lambda_{jT}$ or $\kappa_{j1},\dots, \kappa_{jT}$ from posterior densities that are marginalized w.r.t. to the respective other van-par:par.
It is possible to construct a state space model with the same marginal properties as the one introduced in Section (ref), albeit with $\lambda_{j0}, \dots, \lambda_{jT}$ marginalized out. It is characterized by a generalized beta prime transition density for $\psi_{jt} | \kappa_{jt}, a_j, c_j, \rho _j$, specifically
coupled with the following negative binomial transition density for $\kappa_{jt} | \kappa_{j,t-1}, \psi_{j,t-1} , \rho_j, a_j, c_j$:
with the initial value coming from the stationary distribution of the process governing $\kappa_{jt}$, i.e. $\kappa_{j1} | a_j, \rho_j \sim \mbox{\rm NegBin}\left(a_j , 1-\rho_j \right).$ A representation of the interdependencies of this partially marginalized state space model as a directed graph can be found in Figure (ref). Details on the derivations of transition densities and the stationary distribution can be found in Appendix (ref).
In this alternate state space representation, we can now sample $\kappa_{j1}, \dots, \kappa_{jT}$ conditional on $\psi_{j1}, \dots, \psi_{jT}$ but marginalized w.r.t. $\lambda_{j0}, \dots, \lambda_{jT}$, as required. More specifically
with $z^\psi_{jt} = (1-\pi _{j,t-1})\pi_{jt}a_j \psi_{jt}/(a_j \psi_{jt} + c_j (1-\rho _j) )$ and $\pi_{j,t-1}$ and $\pi_{jt}$ defined as in (ref). Further, it can be shown that the normalizing constant of this distribution is available in closed form. From the definition of the hypergeometric function presented in (ref), it follows immediately that following infinite series can be represented by a hypergeometric function:
Therefore, the normalizing constant of (ref) is equal to
with $a^\kappa_{jt} = a_j + c_j + \kappa_{j,t-1}$ and $b^\kappa_{jt} = a_j + c_j + \kappa_{j,t+1}$.
Based on the normalized weights defined by (ref) and (ref), two main options present themselves to sample $\kappa_{jt}$. One could adapt Algorithm 11.5 in fru:book, to derive a forward-filtering-backward-sampling (FFBS) algorithm. This re-frames $\kappa_{jt}$ as a hidden Markov chain with $K_{\footnotesize \max}$ states, with the rows of the transition matrix given by the normalized version of (ref). The attractiveness of this approach lies in the ability to jointly draw the entire path $\kappa_{j1}, \dots, \kappa_{jT}$, thereby improving mixing. The downside is that the computation of all $K_{\footnotesize \max}^2\times (T-1)$ transition matrices is required before sampling can be performed, which becomes fairly computationally demanding for even moderately sized problems.
An alternate approach consists of sampling $\kappa_{jt}$ directly using inverse transform sampling for the normalized version of (ref). The probabilities have a convenient recursive structure (see Appendix (ref)), allowing for efficient computation, as only the probabilities up to the current realization of $\kappa_{jt}$ need to be evaluated. Thus, particularly when many values of $\kappa_{jt}$ are zero or near zero (as is the case with heavy shrinkage), this method can be executed orders of magnitude faster than FFBS. Further, we found that mixing was fairly good when using this method, despite it being a single move sampler. As a result of its superior speed and generally favorable mixing, we advocate for adopting this method to sample $\kappa_{jt}$.
In summary, the proposed sampler alleviates poor mixing for $\rho_j$ by using the approximate likelihood (ref) marginalized w.r.t. $\lambda_{j0}, \dots, \lambda_{jT}$ and $\kappa_{jt}, \dots, \kappa_{jT}$ to construct a collapsed MH-within-Gibbs step for $\rho_j$. A partially marginalized state space model is utilized to sample $\kappa_{j1}, \dots, \kappa_{jT}$ conditional on $\psi_{j1}, \dots, \psi_{jT}$ but marginalized w.r.t. $\lambda_{j0}, \dots, \lambda_{jT}$. This is done by applying inverse transform sampling to the normalized version of (ref). Next, one can sample $\lambda_{j0}, \dots, \lambda_{jT}$ from (ref) as in the fixed $\rho_j$ case. Finally, conditional on $\lambda_{j0}, \dots, \lambda_{jT}$, realizations of $\psi_{j1}, \dots, \psi_{jT}$ are sampled from (ref) as for exchangeable innovations. It is important to note that the sampling order is important when marginalizing out random variables, as the stationary distribution of the Markov chain is not invariant to permutations of the sampling steps.
The Cholesky SV-model developed by lop-etal:par and further refined by bit-fru:ach provides a method for modeling sparse, time-varying variance-covariance matrices $\bm \Sigma_t$, for $t = 1, \dots, T$, of an $M$-dimensional multivariate time series that follows a conditionally zero mean normal distribution, i.e. $\bm y_t | \bm \Sigma_t \sim \mathcal N_M\left(\bm 0, \bm \Sigma_t\right)$. The key ingredient is the decomposition $\bm \Sigma_t = \bm\Lambda_t\bm D_t \bm \Lambda_t'$, where $\bm \Lambda_t$ is lower unitriangular and $\bm D_t$ is a diagonal matrix. As a consequence, $\bm\Lambda_t^{-1}\bm y_t\sim\mathcal N_M\left(\bm 0, \bm D_t\right)$, which, letting the elements of $\bm\Lambda_t^{-1}$ be denoted as $\gamma_{mj,t}$ for $j < m$, can be expressed as $$ { \left[
\right]\left[
\right] \sim \mathcal{N}_M\left(\mathbf{0}, \mathbf{D}_t\right).} $$ As $\bm D_t$ is a diagonal matrix, $M$ independent TVP models result: {
} To further allow for conditional heteroscedasticity, one can let the $\sigma^2_{mt} = e^{h_{mt}}$ follow a stochastic volatility (SV) law of motion, for $m=1, \ldots, M$, where:
For the priors on the model parameters of the SV process, we follow kas-fru:anc, as these are well-established in the literature, see Appendix (ref) for details.
The equations above can now be estimated as a series of independent TVP regression models with univariate outcomes $y_{mt}$, where the $\gamma_{mj,t}$ (and hence also the elements of $\bm\Lambda_{t}^{-1}$) can be recovered from the regression coefficients as $-\beta_{mj,t}$. The idiosyncratic variance matrix $\bm D_t = \mbox{\rm Diag}\left(e^{h_{1t}}, \dots, e^{h_{Mt}}\right)$ can be recovered from the log volatilities of the individual equations, with the exception of $e^{h_{1t}}$, for which, due to the absence of any regressors, an SV model has to be estimated separately. Together, the draws from $\bm\Lambda_{t}^{-1}$ and $\bm D_t$ define the posterior draws of $\bm \Sigma_t$, binding the draws from the individual equations together. Placing a dynamic triple gamma prior on the innovations of the regression coefficients $\beta_{mj,t}$ allows for covariance modeling that is highly sparse, due to the strong shrinkage imposed by the prior, while still allowing for adaptation to locally different (co)variances, thanks to the dependence induced by the DTG.
We apply the Cholesky SV-model equipped with the dynamic triple gamma prior to a data set of 45 out of 50 of the returns on the companies represented in the EURO STOXX 50 index (with 5 dropped due to data availability issues, see Appendix (ref) for a more detailed description). The order of the data set is alphabetical and it spans 810 time points from the 2nd of January, 2020 to the 21st of February, 2023. Due to the COVID-19 pandemic and the ensuing volatility in financial markets, this represents a particularly challenging data set for econometric modeling. $a_j$ is set to 0.5, as we want to induce heavy shrinkage around the origin, and we test two different choices of $c_j$, namely $c_j = 0.5$ and $c_j = 2.5$, as the former implies a marginal horseshoe prior for the innovations with no moments and very heavy tails, while the latter allows for some regularization of the tails. In both setups, we use the GB1 prior for $\rho_j$ defined in (ref), with hyperparameters $a_\rho = 1$, $b_\rho = 0.95$, $\alpha_\rho = 0.5$, and $\beta_\rho = 0.5$. This results in a prior that has a horseshoe shape ranging from 0 to 0.95, incorporating the prior information that we either expect a state to display very little dependence or a large amount of dependence, but not a middling amount. For the global shrinkage parameters of the innovations as well as the means of the initial values we employ the triple gamma prior distributions ${\beta}_{j}^2 \sim TG(a^\tau, c^\tau, \lambda_B^2)$ and $\theta_j \sim TG(a^\xi, c^\xi, \kappa_{B}^2)$, under the hyperpriors defined in Appendix (ref) with hyperparameters $\alpha_{a^\xi} = \alpha_{a^\tau} = \alpha_{c^\xi} = \alpha_{c^\tau} = 5$, $\beta_{a^\xi} = \beta_{a^\tau} = 10$ and $\beta_{c^\xi} = \beta_{c^\tau} = 2$. For the model parameter of the SV equations, we use the default choices of kna-etal-new:shr_TVP. For comparison, we also estimate a Cholesky SV-model under a static triple gamma prior, using the default hyperparameters laid out in knaus2021shrinkage and under the dynamic shrinkage process of kow-etal:dyn, using the default hyperparameter values suggested by the authors. Each model was run for 100 000 iterations, with a burn-in of 20 000 and a thinning factor of 100.
Generally speaking, most elements of the estimated time varying variance-covariance matrix $\bm\Sigma_t$ are very close to zero for large stretches of time. Notable exceptions are early 2020, which correlates with the onset of the COVID-19 pandemic, early 2022, which coincides with the start of the full-scale Russo-Ukrainian War, and late 2022, where the energy crisis lead to increased volatility in Europe. While all three prior specifications result in posterior estimates with this property, there exist two notable differences. First, the dynamic horseshoe displays large, sudden jumps followed by periods of relative calm, showcasing the flexibility of the dynamic triple gamma prior approach, as well as the heavy tails of the horseshoe prior. An example of this behavior can be seen in the third panel of Figure (ref), in the first four months of 2020. Second, the under-shrinking behavior of the static triple gamma approach becomes apparent in some estimated variances and covariances. To be able to accommodate the drastic changes in the magnitude of the estimates, the associated global shrinkage parameter $\theta_j$ can become too large to effectively shrink noise in calmer periods. In contrast, the dynamic triple gamma can adapt to the locally differing variance requirements of a given state. An example of this can be seen in Figure (ref), which displays the marginal time-varying covariance of two returns, where the estimate from the static approach is much noisier than that under either dynamic specification. High resolution images of all variance-covariance matrices can be found \href{https://imgur.com/a/xoiGkFb}{here.}\footnote{In case the link is not clickable: https://imgur.com/a/xoiGkFb}
To benchmark the out-of-sample forecasting performance of our approach, we compare the one step ahead log-predictive density scores (LPDS) for the last 100 time points in the data set, as in bit-fru:ach. The cumulative mean of this exercise can be found in Figure (ref). Both the slightly regularized dynamic triple gamma as well as the dynamic horseshoe outperform the static triple gamma, as is evidenced by the higher average LPDS throughout virtually the entire sample period. Furthermore, both the static and dynamic triple gamma priors clearly outperform the dynamic shrinkage process of kow-etal:dyn. To gain intuition about how the dynamic triple gamma prior outperforms the static approach, it is instructive to look at Figure (ref). It graphs the point-wise difference in LPDS vis-\`a-vis the static approach, for both the regularized dynamic triple gamma, as well as the dynamic horseshoe. It is interesting to note that the dynamic approaches forecasts better only in around $60 \%$ of time points. However, whenever the forecasting performance of the dynamic approach is worse than the static approach, it is usually not worse by much. In contrast, when the forecasting performance is better, it is often substantially so, leading to forecasting performance that is, on average, much better than under the static approach.
This paper introduced the dynamic triple gamma (DTG) prior as a novel shrinkage process for time-varying parameter (TVP) models. By generalizing the triple gamma prior to a stochastic process, the DTG approach balances local adaptivity and regularization, capturing time-varying innovation variances while remaining computationally tractable. The DTG prior’s theoretical properties were explored, clarifying connections to existing approaches like the NGAR and dynamic shrinkage processes. An application to the EURO STOXX 50 index demonstrated its utility, capturing complex patterns in time-varying (co)-variances and outperforming competing approaches in predictive performance. The DTG prior offers a flexible, interpretable, and computationally efficient approach to modeling time-varying parameters in Bayesian frameworks. Its combination of global shrinkage, local adaptivity, and time-varying dependence makes it a powerful tool for applications in econometrics, financial modeling, and beyond.
Future work could explore sharing information across equations for parameters $\rho_1, \dots, \rho_d$ via a hierarchical prior. For example, clustering these parameters with finite mixture models could encode prior beliefs about groups of covariates with similar dependence structures, enhancing shared learning.
{\bf Disclosure statement:} The authors report there are no competing interests to declare.
\centerline{{\bf SUPPLEMENTARY MATERIAL}} {\bf Appendix:} Appendix with proofs and additional technical details (PDF format).