# Modelling and Forecasting Macroeconomic Risk with Time Varying Skewness
Stochastic Volatility Models ††thanks: This paper is part of my PhD
dissertation at the University of Bologna. I am grateful to Andrea Carriero
for invaluable support and guidance. I thank Luca Fanelli and Umberto
Cherubini for comments and suggestions and participants at the March 2021
Unibo PhD Forum and the July 2022 Unibo PhD Forum.

Andrea Renzetti111Department of Economics, Alma Mater Studiorium Università di
Bologna, Piazza Scaravilli 2, 40126 Bologna, Italy

( First Draft: September 21, 2021
This Draft:
Latest Draft Here)

###### Abstract

In this paper I propose a parametric framework for modelling and forecasting
macroeconomic tail risk based on stochastic volatility models with Skew-Normal
and Skew-t shocks featuring stochastic skewness. The paper develops posterior
simulation samplers for Bayesian estimation of both univariate and VAR models
of this type. In an application, I use the models to predict downside risk to
GDP growth and I show that this approach represents a competitive alternative
to quantile regression. Finally, estimating a medium scale VAR on US data I
show that time varying skewness is a relevant feature of macroeconomic and
financial shocks.

_J.E.L Classification Code: C22, C32, C53_

_Keywords:_ Stochastic Volatility, Stochastic Skewness, Bayesian VARs,
Macroeconomic tail risk

##  1 Introduction

Central banks and policy institutions play a critical role in maintaining
financial stability and fostering economic growth. A key challenge they face
is effectively monitoring the likelihood of severe events that could have
adverse effects on the economy. Failing to adequately assess these risks can
lead to underestimation of potential losses and insufficient policy responses.

To address this challenge, it is essential to develop econometric models that
can accurately predict and assess tail risk in macroeconomic outcomes. In this
paper, I propose an econometric framework specifically designed for modeling
and forecasting macroeconomic tail risk. The framework relies on fully
parametric univariate and multivariate stochastic volatility models with Skew-
Normal and Skew-t shocks. These models aim to capture and predict persistent
time-varying asymmetries in the future distribution of macroeconomic variables
of interest. Capturing these asymmetries is especially relevant given the risk
management nature of the problem of policymaking faced by central banks and
policy institutions [29].

The paper begins by extending the well-known univariate stochastic volatility
model introduced by [25] to explicitly account for time-varying conditional
skewness in the predictive distribution of a given target variable. Then,
building upon this univariate approach, the paper introduces a Bayesian Vector
Autoregressive (VAR) model with stochastic volatility and time-varying
skewness. By allowing to track changes in the shape of the predictive
distribution of multiple time series, this model is suitable for
quantification and forecasting of tail risk to multiple target variables.
Importantly, the model retains all the advantages and familiar toolkit for
policy analysis and scenario analysis associated to the VAR framework.

To demonstrate the effectiveness of the proposed framework, I use the time-
varying skewness stochastic volatility models to assess and predict downside
risk to GDP growth. The findings of this analysis align with the main
conclusions of [2], revealing a nonlinear and asymmetric impact of financial
conditions on the future distribution of GDP growth. Additionally, the models
provide slightly more accurate out-of-sample forecasts of downside risk
compared to quantile regression. Furthermore, by estimating a medium-scale VAR
model, I show that shocks to financial and macroeconomic time series exhibit
significant time-varying volatility and time-varying skewness, suggesting that
taking into consideration both of these features might be of particularly
relevance for accurately assessing upside and downside risk to macroeconomic
indicators.

Related literature A fast-growing body of studies recently used univariate
quantile regression methods for modelling and predicting asymmetries in the
future distribution of macroeconomic variables. For example, [20] used
predictive quantile regression to investigate whether systemic risk indicator
and financial distress indicators predict changes in the lower quantiles of
future macroeconomics shocks. As well, [28] used quantile regression to
examine fluctuations in the risk of a large increase in unemployment. More
recently, [2] used a two step-procedure based on predictive quantile
regression and quantile interpolation to model changes in downside risk to
future GDP growth as a function of current financial and economic
conditions.222The two step approach based on quantile regression of [2] gained
substantial popularity in the literature and has been employed in many other
frameworks to assess and predict tail risk to economic outcomes. Among the
others, [31] used the two step approach of [2] for assessing and predicting
downside and upside risk to inflation while [18] used the same approach for
predicting the probability of large capital out-flows and in-flows to emerging
markets. Despite its popularity, the quantile regression method of [2]
typically fails in the presence of a large information set where fully
parametric models often produce more accurate forecasts of downside risks
[12]. As a matter of fact, when using quantile regression, including multiple
lags of the dependent and independent variables so a to capture the rich
autocorrelation structure of macroeconomic and financial time series becomes
very impractical and often leads to imprecise estimates of the coefficients
and problems such as quantile crossing. Moreover, the entire predictive
distribution of the target variables can only be obtained in two steps by
interpolating the estimated quantiles with a flexible distribution.

In the light of these limitations a new wave of studies have recently brought
some evidences in favour of the use of fully parametric models to assess and
predict tail risk to macroeconomic outcomes. [7] for example, show that
standard GARCH models have superior forecasting performance with respect to
quantile regression methods for forecasting downside risk to GDP growth. As
well, [12] show that a Bayesian VAR with stochastic volatility performs
comparably to quantile regression for estimating and forecasting tail risks.
Here I follow and extend this line of research by considering fully parametric
models featuring both time varying volatility and time varying skewness, as
recently done by [16], [24], [36] and [32]. While the first three
contributions are all univariate 333[16] propose a score driven model with
Skew-t innovations. [23] is the first paper to introduce time varying
conditional skewness in a univariate stochastic volatility model by exploiting
a Noncentral-t distribution for the innovations. [36] exploits the Skew-Normal
distributions but considers a different parametrization for the shocks with
respect to the univariate model that I consider in Section 2.1 relying as well
on a different estimation strategy. in this paper I model time varying
volatility together with time varying skewness both in a univariate and in a
multivariate framework. The main advantages of the multivariate framework is
that it allows to jointly model the dynamic relationship between the target
variables and the risk factors and to explicitly model tail risk to multiple
macroeconomic outcomes of interest. The multivariate model that I propose in
this paper is a VAR model in which Bayesian shrinkage can be conveniently used
to avoid over-fitting when exploiting a potential large information set due
both to the inclusion of larger number of macroeconomic variables and of a
meaningful number of lags needed to properly account for the rich
autocorrelation structure of the macroeconomic and financial time series. The
model features two distinct stochastic processes respectively governing the
time varying volatility and the time varying skewness of the shocks. By
considering distinct stochastic processes for the skewness and the volatility
of the shocks, this model is different from the Bayesian VAR with Skew-Normal
shocks introduced by [32] where the latent stochastic process governing the
shape of the shocks influences not only the conditional skewness, but also the
conditional mean and the conditional variance of the variables in the system.
To my knowledge, this is the first paper that estimates a VAR with two
distinct stochastic processes for the volatility and the skewness of the
shocks. The rest of the article is organized as follows. In Section 2.1 I
present the univariate stochastic volatility models with Skew-Normal and
Skew-t shocks featuring time varying skewness. Then in Section 2.2 I exploit
the same conceptual framework to model time varying skewness together with
time varying volatility in a VAR framework. In both sections I present
posterior simulation samplers used for Bayesian estimation of these models. In
Section 3 I and use the models to predict downside risk to GDP growth and
compare the forecasting performances to the popular two step approach based on
quantile regression by [2]. Finally, in Section 4 I estimate a medium scale
VAR model and show that many macroeconomic and financial variables exhibit
time varying conditional skewness.

##  2 Models

###  2.1 Univariate time varying skewness stochastic volatility model

Stochastic volatility models currently represent the state of the art for
modelling and forecasting macroeconomic and financial time series. The basic
stochastic volatility model of [25] specifies a log-normal auto-regressive
process for the conditional variance with independent innovations in the
conditional mean and conditional variance equation. In a second contribution,
[26] introduce a stochastic volatility model that features correlation between
the volatility and mean innovations (leverage effects) allowing for
conditional skewness, but without modelling it explicitly. [8] present a
stochastic volatility model where the shocks feature a Skew-GED distribution
while [1] introduce a stochastic volatility with Skew-t innovations. Both
contributions explicitly model conditional skewness, but do not allow for time
varying conditional skewness. Here I present a direct extension of the
univariate stochastic volatility model of [25] that instead explicitly allows
for time varying conditional skewness.

In order to model asymmetries in the conditional distribution of the dependent
variable, I assume that the innovations in an otherwise standard stochastic
volatility model follow a potentially asymmetric distribution, being the Skew-
Normal [5] and the Skew-t [6] distribution. The $Skew-
Normal(\zeta,\omega^{2},\lambda)$ is an asymmetric distribution fully
characterized by three parameters: the location parameter $\zeta$, the scale
parameter $\omega^{2}$ and the shape parameter $\lambda$. The shape parameter
$\lambda$ governs the skewness of this distribution. As $\lambda=0$ the Skew-
Normal becomes symmetric and collapses to the Normal. Positive values of
$\lambda$ are associated with a right skewed distribution while negative
values of $\lambda$ are associated with a left skewed distribution. 444See
Appendix A.1 for details on the Skew-Normal and Skew-t. To model time
variation in the shape of the shocks, I treat the shape parameter $\lambda$ as
an additional stochastic process in the model:

$y_{t}=\bm{x_{t}\pi}+\sqrt{h_{t}}\varepsilon_{t}\hskip
56.9055pt\varepsilon_{t}\sim\text{{Skew-
Normal}}(\zeta_{t},\omega_{t}^{2},\lambda_{t})$ (1)
$log(h_{t})=\phi_{h}log(h_{t-1})+\eta_{t}\hskip
71.13188pt\eta_{t}\sim\mathcal{N}(0,\sigma^{2}_{\eta})$ (2)
$\lambda_{t}=\phi_{\lambda}\lambda_{t-1}+\xi_{t}\hskip
113.81102pt\xi_{t}\sim\mathcal{N}(0,\sigma^{2}_{\xi})$ (3)

where $y_{t}$ is the dependent variable observed over the periods
$t=1,\ldots,T$, while $\bm{x_{t}}$ is a row vector of that might contain lags
of the dependent variable and other exogenous regressors and $\bm{\pi}$ is the
column vector of coefficients. I assume that the Skew-Normal shocks have zero
mean and unit variance, that is
$\operatorname*{\mathbb{E}}[\varepsilon_{t}]=0$ and $var(\varepsilon_{t})=1$,
which implies the following constraints on the location and scale parameters:

$\zeta_{t}=-\omega_{t}\delta_{t}\sqrt{\frac{2}{\pi}}\quad\quad\quad\lor t$ (4)
$\omega^{2}_{t}=\left[1-\frac{2}{\pi}\delta_{t}^{2}\right]^{-1}\quad\quad\lor
t$ (5)

where $\delta_{t}=\frac{\lambda_{t}}{\sqrt{1+\lambda_{t}^{2}}}$, with
$-1<\delta_{t}<1$. This parametrization ensures that
$\operatorname*{\mathbb{E}}[y_{t}|\mathcal{I}_{t-1}]=\bm{x_{t}\pi}$. In this
regard, it is important to remark that imposing $\zeta_{t}=0$ instead of (4)
would imply $\operatorname*{\mathbb{E}}[\varepsilon_{t}]\neq 0$, and in
general $\operatorname*{\mathbb{E}}[\varepsilon_{t}|\mathcal{I}_{t-1}]\neq
0$.555 Imposing $\zeta_{t}=0$ instead of (4) leads to a model with a time
varying intercept, shifting the conditional mean of $y_{t}$ proportionally to
$\lambda_{t-1}$. As well, this parametrization ensures that $y_{t}$ features
both time varying conditional volatility and time varying conditional skewness
with the former exclusively governed by the stochastic process in equation (2)
while the latter by the stochastic process in (3).666It is possible to have a
model that features both time varying volatility and time varying skewness by
assuming:  $\begin{array}[]{l}y_{t}=\bm{x_{t}\pi}+\varepsilon_{t}\hskip
56.9055pt\varepsilon_{t}\sim\text{{Skew-
Normal}}(\zeta_{t},\omega_{t}^{2},\lambda_{t})\\\
\lambda_{t}=\phi_{\lambda}\lambda_{t-1}+\xi_{t}\hskip
113.81102pt\xi_{t}\sim\mathcal{N}(0,\sigma^{2}_{\xi})\end{array}$ assuming
$\mathbb{E}(\varepsilon_{t})=0$ (hence (4) still holds) and imposing
$\omega^{2}=1$ which implies $var(\varepsilon_{t})\neq
1=\left(1-\frac{2\delta_{t}^{2}}{\pi}\right)$. However, in this case the
parameter $\lambda_{t}$ would drive both conditional skewness and conditional
volatility. This is not desirable in general, since we might want to model
these two distinct features using different dynamics.

In order to explicitly model heavy-tails, together with time-varying skewness,
I also consider an alternative specification where the innovations are
distributed as a ${Skew-t(\zeta_{t},\omega_{t}^{2},\lambda_{t},\nu)}$ [6]. The
parameter of the degrees of freedom $\nu$ determines the tail thickness of the
Skew-t distribution: as $\nu\rightarrow\infty$ the Skew-t converges to the
Skew-Normal while when $\lambda=0$ the Skew-t collapses to a Student-t with
$\nu$ degrees of freedom. In this case the constraints on the location and
scale parameters that ensure $\operatorname*{\mathbb{E}}[\varepsilon_{t}]=0$
and $var(\varepsilon_{t})=1$ become:

$\zeta_{t}=-\omega_{t}\delta_{t}k_{1}\sqrt{\frac{2}{\pi}}\quad\quad\quad\lor
t$ (6)
$\omega^{2}_{t}=\left(k_{2}-\frac{2}{\pi}k_{1}^{2}\delta_{t}^{2}\right)^{-1}\quad\quad\quad\lor
t$ (7)

where
$k_{1}=\sqrt{\frac{\nu}{2}}\frac{\Gamma(\frac{\nu-1}{2})}{\Gamma(\frac{\nu}{2})}$,
$k_{2}=\frac{\nu}{\nu-2}$ and $\Gamma(.)$ is the Gamma function. This
stochastic volatility model with Skew-t shocks includes as special cases both
the stochastic volatility model with heavy tails without conditional skewness
of [26] and the model with heavy tails and constant conditional skewness of
[1].777The stochastic volatility model with heavy tails of [26] is a
particular version of this model where the shape parameter is constant and
equal to 0, that is $\lambda_{t}=0\hskip 5.69046pt\lor t$. As well, the
stochastic volatility model with skewness and heavy tails of is a particular
version of this model where $\sigma^{2}_{\xi}\rightarrow 0$ and
$\phi_{\lambda}=1$, namely the shape parameter $\lambda_{t}$ is constant. It
is straightforward to modify this specification by assuming a different
dynamics for the log-volatility and the shape parameter in the state equations
(2) and (3). For example if we suspect that some of the variables in
$\bm{x_{t}}$ affect not only the conditional mean, but also the conditional
variance and the conditional skewness of $y_{t}$, we can include them in the
state equations of these two distinct stochastic processes. For example, as it
will be shown in the application to the Growth at Risk framework in Section 3,
motivated by the findings of [2] and subsequent work by [16], [32] and [36] I
consider a specification in which financial condition affect not only the
conditional mean but also the conditional skewness of the future GDP growth
distribution.

####  2.1.1 Priors and estimation of the univariate TVS-SV

For what concerns the specification of the prior distributions for the
parameters of the model, I assume a Normal prior for the regression
coefficients ($\bm{\pi}$) and for the coefficients in the state equations
($\phi_{\lambda}$ and $\phi_{h}$) while I specify an Inverse Gamma Prior for
the variances of the innovations to the log-volatility and to the shape
parameter ($\sigma^{2}_{\eta}$ and $\sigma^{2}_{\xi})$ .

This section develops a posterior sampler which allows for Bayesian estimation
of the univariate models presented above. The estimation strategy leverages on
the fact that $\varepsilon_{t}\sim\text{{Skew
Normal}}(\zeta_{t},\omega^{2}_{t},\lambda_{t})$ has the following stochastic
representation :

$\varepsilon_{t}=\zeta_{t}+\delta_{t}\omega_{t}v_{t}+\sqrt{(1-\delta_{t}^{2})}\omega_{t}z_{t}$
(8)

where $v_{t}\stackrel{{\scriptstyle i.i.d}}{{\sim}}\text{Truncated
Normal}_{[0,\infty)}(0,1)$ and $z_{t}\stackrel{{\scriptstyle
i.i.d}}{{\sim}}\mathcal{N}(0,1)$. Equation (8) implies that conditioning on
the mixing variable $v_{t}$, the random variable $\varepsilon_{t}$ is
distributed as a Normal. This result greatly simplifies the derivation of the
full conditional distributions in the Gibbs Sampler and allows to exploit and
adapt many of the results used for the estimation of the standard stochastic
volatility model with Gaussian innovations [25]. In particular, once I have
obtained a draw from the full conditional posterior distribution of the mixing
variable $v_{t}$, I can exploit and the Gaussian distribution when deriving
the full conditional posterior of the parameters and the unobserved states in
the Gibbs Sampler. As a matter of fact, $\zeta_{t}$, $\omega_{t}$ and
$\delta_{t}$ are neither parameters nor latent states to be estimated.
$\zeta_{t}$ and $\omega_{t}$ satisfy the constraints (4) and (5) and ensure
the correct parameterization of the shocks at each time period $t=1,\ldots,T$,
while $\delta_{t}$ is a one to one map to $\lambda_{t}$, namely
$\delta_{t}=\frac{\lambda_{t}}{\sqrt{1+\lambda_{t}^{2}}}$.

Table 1 presents the details on the Gibbs Sampler while Appendix A.2 reports
the derivations of the full conditional posterior distributions. In Step 1) I
sample the mixing variables $\\{v_{t}\\}_{t=1}^{T}$ from the full conditional
posterior distribution $p(v_{t}|\bm{\Theta},\bm{\lambda},\bm{h})$ which is a
Truncated Normal distribution. Steps 2) 3) 4) 5) 6) are pretty standard: I
draw the regression coefficients $\bm{\pi}$ in the observation equation (1)
and the autoregressive coefficients and the variances in the two state
equations (2) (3) from their respective full conditional posterior
distributions. In Step 7) and Step 9) I draw the initial states for the
volatility $h_{0}$ and the shape parameter $\lambda_{0}$, while in Steps 8)
and 10) I draw the entire history for the volatilities and the shape
parameters. Since it is not feasible to directly sample from the full
conditional distributions of the volatilities
$p(h_{1},\ldots,h_{T}|\bm{\Theta},\bm{v},\bm{\lambda})$ and the shape
parameters $p(\lambda_{i1},\ldots,\lambda_{iT}|\bm{\Theta},\bm{v},\bm{h})$ I
rely on the particle filter to approximate these distributions. I also
considered a version of the algorithm in which I use independence Metropolis
Hastings steps to draw both the volatilities and the shape parameters but I
experienced that the algorithm based on the particle filter has smaller mixing
times.888In the particle steps, in order to alleviate path degeneracy, I
exploit the Ancestor Sampling procedure developed in [30] which enables fast
mixing even when using seemingly few particles. [30] study the properties of
the sampler and provide the formal proof for the convergence of the algorithm.
In the particle approximation, I use the transition equations (2) and (3) as
importance densities and compute the weights accordingly. The details on the
particle steps used to approximate the full conditional posterior distribution
of the volatilities and the shape parameters can be find in Table 4 in the
Appendix A.4. As well, in the Appendix A.4, I report the details on the steps
of the alternative algorithm which relies on the independence Metropolis
Hastings steps to draw the volatilities and the shape parameters.

Table 1: MCMC algorithm for the univariate TVSSV model MCMC for the univariate TVSSV model |
---|---
Initialize $\bm{\Theta}^{(0)},\bm{s}^{(0)}$ |
For $m=0:\text{Total MCMC draws}$ |
1) Draw $\\{v_{t}\\}_{t=1}^{T^{(m+1)}}$ from $p(v_{1}\ldots,{v_{T}}|\bm{\Theta}^{(m)},\bm{\lambda}^{(m)}),\bm{h}^{(m)})$ |
2) Draw $\bm{\pi}^{(m+1)}$ from $p(\bm{\pi}|\bm{\Theta}^{(m)},\bm{v}^{(m)},\bm{\lambda}^{(m)}),\bm{h}^{(m)})$ |
3) Draw $\sigma_{\eta}^{2^{(m+1)}}$ from $p(\sigma^{2}_{\eta}|\bm{\Theta^{(m)},v^{(m)}},\bm{\lambda}^{(m)}),\bm{h}^{(m)})$ |
4) Draw $\sigma_{\xi}^{2^{(m+1)}}$ from $p(\sigma^{2}_{\xi}|\bm{\Theta^{(m)},v^{(m)}},\bm{\lambda}^{(m)}),\bm{h}^{(m)})$ |
5) Draw $\phi_{h}^{(m+1)}$ from $p(\phi_{h}|\bm{\Theta^{(m)},v^{(m)}},\bm{\lambda}^{(m)}),\bm{h}^{(m)})$ |
6) Draw $\phi_{\lambda}^{(m+1)}$ from $p(\phi_{\lambda}|\bm{\Theta^{(m)},v^{(m)}},\bm{\lambda}^{(m)}),\bm{h}^{(m)})$ |
7) Draw $h_{0}^{(m+1)}$ from $p(h_{0}|\bm{\Theta}^{(m)},\bm{v}^{(m)},\bm{\lambda}^{(m)}),\bm{h}^{(m)})$ |
8) Draw $\\{h_{t}\\}_{t=1}^{T^{(m+1)}}$ from $p(h_{1},\ldots,h_{T}|\bm{\Theta^{(m)}},\bm{v^{(m)}},\bm{\lambda}^{(m)}))$ |
Particle Step |
9) Draw $\lambda_{0}^{(m+1)}$ from $\lambda_{0}^{(m+1)}$ from $p(\lambda_{0}|\bm{\Theta}^{(m)},\bm{v}^{(m)},\bm{\lambda}^{(m)}),\bm{h}^{(m)})$ |
10) Draw $\\{\lambda_{t}\\}_{t=1}^{T^{(m+1)}}$ from $p(\lambda_{i1},\ldots,\lambda_{iT}|\bm{\Theta^{(m)}},\bm{v^{(m)}}),\bm{h}^{(m)})$ |
Particle Step |
end |

To estimate the version of the model with Skew-t innovations, I just exploit
the fact that
$\varepsilon_{t}\sim\text{{Skew-t}}(\zeta_{t},\omega^{2}_{t},\lambda_{t},\nu)$
has as well a convenient stochastic representation, namely:

$\varepsilon_{t}=\zeta_{t}+\delta_{t}\omega_{t}o_{t}^{-0.5}v_{t}+\sqrt{(1-\delta_{t}^{2})}\omega_{t}o_{t}^{-0.5}z_{t}$
(9)

$v_{t}\stackrel{{\scriptstyle i.i.d}}{{\sim}}\text{Truncated
Normal}_{[0,\infty)}(0,1)$, $z_{t}\stackrel{{\scriptstyle
i.i.d}}{{\sim}}\mathcal{N}(0,1)$ and $o_{t}\stackrel{{\scriptstyle
i.i.d}}{{\sim}}\mathcal{G}(\frac{\nu}{2},\frac{\nu}{2})$
This is the same same representation of the Skew-Normal except for the
additional mixing variable $o_{t}$. Therefore, conditioning in this case on
the two mixing variables $m_{t}=\\{v_{t},o_{t}\\}$, the random variable
$\varepsilon_{t}$ is distributed as a Normal. Once again, I can exploit and
adapt the derivations of the standard model with Gaussian shocks when deriving
the full conditional posterior distribution in the Gibbs Sampler. In order to
estimate the model it is just needed to consider a further initial step to
draw from $p(o_{1}\ldots,{o_{T}}|\bm{\Theta},\bm{v},\bm{\lambda},\bm{h})$,
namely: 999 In the estimation of the model with heavy tails (Skew-t shocks), I
fix the tail thickness parameters $\nu$ to 5. Given the relative short time
series length of macroeconomic data, it is particularly difficult to make
inference on this parameter. In general, you can draw this parameter adding
another Metropolis Hastings step to draw from
$p(\nu|\bm{\Theta},\bm{v},\bm{o},\bm{\lambda})$

Draw $\\{o_{t}\\}_{t=1}^{T^{(}m+1)}$ from
$p(o_{1}\ldots,{o_{T}}|\bm{\Theta}^{(m)},\bm{v}^{(m)},\bm{\lambda}^{(m)}),\bm{h}^{(m)})$
---

and then adapt Steps 2) to 10) in Table (1) with the new formulas of the full
conditional distributions derived by conditioning on the further mixing
variables $\\{o_{t}\\}_{t=1}^{T}$. In this case, since it is not possible to
directly sample from the full conditional distribution of the mixing variable
$o_{t}$, I use Metropolis Hastings to simulate draws from this distribution.
Appendix A.3 reports the details of this step.

###  2.2 Time varying skewness stochastic volatility VAR model

Given the risk management nature of the problem of policymaking, it is often
the case that the objective of interest is to quantify and predict tail risk
to multiple macroeconomic outcomes [29]. In particular, from a modelling
perspective, we might be interested in a multivariate model that can
characterize asymmetries in the future distribution of multiple macroeconomic
time-series. VAR models [34] emerged as the natural tool to capture the rich
dynamic interrelationship between multiple macroeconomic time series. They
currently represent the workhouse in empirical macroeconomics and are
routinely used for forecasting and policy analyses [35]. In this section I
exploit the conceptual framework presented in the previous section to jointly
model the dynamic behaviour of multiple time series in a Bayesian VAR model
and capture time varying skewness in the conditional distribution of the
variables in the system. The model is given by:

$\bm{y_{t}}=\bm{\Pi_{0}}+\bm{\Pi_{1}y_{t-1}}+\ldots+\bm{\Pi_{p}y_{t-p}}+\bm{A}^{-1}\bm{H_{t}}^{0.5}\bm{\varepsilon_{t}}$
(10)

where $\bm{y}_{t}$ is an $N\times 1$ vector of variables observed over the
periods $t=1,\ldots,T$. $\bm{H_{t}}$ is a diagonal matrix that contains the
volatilities on its main diagonal, namely
${\bm{H_{t}}=diag(h_{1,t}\ldots,h_{N,t})}$ and $\bm{A}^{-1}$ is a lower
triangular matrix with ones on its main diagonal. The log-volatilities evolve
over time according to:

$log(h_{i,t})=\phi_{h,i}log(h_{i,t-1})+\eta_{i,t}\hskip
85.35826pt\eta_{i,t}\sim N(0,\sigma^{2}_{\eta,i})$ (11)

for $i=1,\ldots,N$. In the Gaussian stochastic volatility model of [15] and
[33] it is assumed $\bm{\varepsilon_{t}}\sim N(0,\bm{I})$. In our
specification, $\bm{\varepsilon_{t}}$ is a vector of Skew-Normal shocks,
namely:

$\bm{\varepsilon_{t}}=[\varepsilon_{1t},\ldots,\varepsilon_{Nt}]^{\prime}\hskip
56.9055pt\varepsilon_{it}\sim\text{{Skew-
Normal}}(\zeta_{it},\omega_{it}^{2},\lambda_{it})$ (12)

where the shape parameters $\lambda_{it}$ evolve according to:

$\lambda_{i,t}=\phi_{\lambda,i}\lambda_{i,t-1}+\xi_{i,t}\hskip
113.81102pt\xi_{i,t}\sim N(0,\sigma^{2}_{\xi,i})$ (13)

In order to have $\operatorname*{\mathbb{E}}[\bm{\varepsilon_{t}}]=\bm{0}$ and
$var(\bm{\varepsilon_{t}})=\bm{I}$ the shocks are parameterized imposing the
constraints on the location parameters $\zeta_{it}$ and on the scale
parameters $\omega_{it}$ discussed in the previous section. As in the
univariate framework, I can explicitly model heavy-tails, together with time-
varying skewness, by considering an alternative specification where:

$\bm{\varepsilon_{t}}=[\varepsilon_{1t},\ldots,\varepsilon_{Nt}]^{\prime}\hskip
56.9055pt\varepsilon_{it}\sim\text{{Skew-t}}(\zeta_{it},\omega_{it}^{2},\lambda_{it},\nu)$
(14)

The model nests the constant coefficients version of the popular VAR model
with stochastic volatility introduced by [15] and [33] and considered in
[10].101010As well, the stochastic volatility VAR with fat tails in [14] is
also a special case of this model with $\lambda_{i,t}\hskip 5.69046pt=0\lor
i,t$. [27] stochastic volatility model VAR with Skew-t orthogonal residual is
as well a particular version of this model with $\phi_{\lambda_{i}}=1$ and
$\sigma^{2}_{\xi,i}\rightarrow 0\hskip 5.69046pt\lor i$. In these models, as
long as the short run restrictions implied by the Cholesky ordering are
satisfied, the shocks can be interpreted as structural.111111It is worth to
mention that due to the “Cholesky type” specification of the stochastic
volatility VAR model considered here, the order in which the variables enter
in the VAR matters not only for the identification of the shocks but also for
the estimation of the model. This is fact was stressed first by [33] and more
recently by [4] and [13]. On the lines of the work of [13] I am currently
working on a order invariant version of the model considered in this paper.
This means that, other than for forecasting purposes, the model can be
practically used for policy analysis and structural scenario analyses. Also in
this multivariate framework it is straightforward to modify the specification
of the state equations of the log-volatilities and the shape parameters by
assuming a different dynamics in (11) and (13). For example, as it will be
shown in the empirical application in Section 3, I can capture the nonlinear
relationship between two variables in the VAR by including the lags of one
variable in the state equations of the log-volatility and/or the shape
parameter of the shocks to the other variable.

####  2.2.1 Priors and estimation of the TVS-SV VAR

For what concerns the choice of the prior distributions for the parameters of
the model, I assume a Normal prior for the autoregressive coefficients
$vec(\bm{\Pi})$. As well, following [15], I specify a Normal prior for the
free elements in the matrix $\bm{A}$. Finally, as in the univariate framework,
I specify independent Inverse Gamma priors for the variance of the innovations
to the log-volatilities and to the shape parameters ($\sigma^{2}_{\eta,i}$ and
$\sigma^{2}_{\xi,i})$ and Normal priors for the coefficients in the state
equations ($\phi_{h,i}$ and $\phi_{\lambda,i})$. The estimation strategy for
the VAR model is just a generalization of the one for the univariate model
that again leverages on the stochastic representation of the Skew Normal (8)
and Skew-t (9) shocks. Exploiting this representation, I can write the vector
of Skew-Normal shocks $\bm{\varepsilon_{t}}$ as follows: 121212Note that the
powers on the matrices refer all to diagonal matrices. For example
$(\bm{I}_{N}-\bm{\Delta_{t}}^{2})=diag(\sqrt{1-\delta_{1,t}^{2}},\ldots,\sqrt{1-\delta_{N,t}^{2}})$
or afterwords
$\bm{O_{t}}^{-0.5}=diag\left(\frac{1}{\sqrt{o_{1,t}}},\ldots,\frac{1}{\sqrt{o_{N,t}}}\right)$

$\bm{\varepsilon_{t}}=\bm{\zeta_{t}}+\bm{\Omega_{t}}\bm{\Delta_{t}}\bm{v_{t}}+{\bm{\Omega_{t}}}(\bm{I}_{N}-\bm{\Delta_{t}}^{2})^{0.5}\bm{z_{t}}$
(15)

where:

$\bm{\zeta_{t}}=[\zeta_{1,t},\ldots,\zeta_{N,t}]^{\prime}$

$\bm{\Omega_{t}}=diag(\omega_{1t}\ldots\omega_{Nt})$

$\bm{\Delta_{t}}=diag(\delta_{1t}\ldots\delta_{Nt})$

$\bm{v_{t}}=[v_{1,t},\ldots,v_{N,t}]^{\prime}\hskip 56.9055ptv_{i,t}\sim
TruncatedNormal_{(0,\infty)}(0,1)$

$\bm{z_{t}}=[z_{1,t},\ldots,z_{N,t}]^{\prime}\hskip 56.9055ptz_{it}\sim
N(0,1)$.

As in the univariate framework, I can exploit this result when deriving the
full conditional posterior distributions of the parameters and the unobserved
states in the Gibbs Sampler. As a matter of fact, also in this case,
$\zeta_{it}$ and $\omega_{it}$ respectively stored in the column vector
$\bm{\zeta_{t}}$ and in the diagonal matrix $\bm{\Omega_{t}}$ are neither
parameters nor latent states to be estimated. $\zeta_{it}$ and $\omega_{it}$
are fixed to satisfy the constraints (4) and (5) and ensure the correct
parameterization of the shocks in each equation of the VAR $i=1,\ldots,N$ and
at each time period $t=1,\ldots,T$. As well, the elements in the diagonal
matrix $\bm{\Delta_{t}}$ (that is $\delta_{it}$) are one to one map of the
latent states $\lambda_{it}$.

Table (2) presents the details of the sampler. In Step 1) I draw the mixing
variables $\\{v_{it}\\}_{t=1}^{T}$ for $i=1,\ldots,N$. In Step 2) I draw the
coefficients of the VAR coefficients adapting to my framework the correct
version of the triangular algorithm developed in [10] and corrected in [9].
This approach allows to reduce the computational burden associated to the
system-wide estimation of Bayesian VAR with stochastic volatility and non-
conjugate priors by exploiting a triangularization of the system. In Step 3),
I adapt the approach of [15] to draw the free elements in the matrix $\bm{A}$.
In Step 4) 5) and 6) 7) I draw the variances and the autoregressive
coefficients of the state equations while in Step 8) and 10) I draw the
initial state for the volatilities $h_{i0}$ and the shape parameters
$\lambda_{i0}$ . In Step 9) and 11) I draw the entire path for the
volatilities and the shape parameters, using the Particle Step with Ancestor
Sampling described in Table 4 in the Appendix A.4.

Table 2: MCMC algorithm for the TVSSV VAR model Particle Gibbs Sampler for the TVSSV-VAR model |
---|---
Initialize $\bm{\Theta}^{(0)},\bm{s}^{(0)},\bm{v}^{(0)}$ |
For $m=0:\text{Total MCMC draws}$ |
1) Draw $\\{\bm{v_{it}}\\}_{t=1}^{T^{(m+1)}}$ from $p(\bm{v_{i1}}\ldots,{v_{iT}}|\bm{\Theta}^{(m)},\bm{s}^{(m)})$ | $i=1,\ldots,N$
2) Draw $\bm{\Pi}^{(m+1)}$ from $p(\bm{\Pi}|\bm{\Theta}^{(m)},\bm{v}^{(m)},\bm{s}^{(m)})$ |
3) Draw $\bm{A}^{(m+1)}$ from $p(\bm{A}|\bm{\Theta}^{(m)},\bm{v}^{(m)},\bm{s}^{(m)})$ |
4) Draw $\sigma_{\xi,i}^{2^{(m+1)}}$ from $p(\sigma^{2}_{\xi,i}|\bm{\Theta^{(m)},s^{(m)},v^{(m)}})$ | $i=1,\ldots,N$
5) Draw $\sigma_{\eta,i}^{2^{(m+1)}}$ from $p(\sigma^{2}_{\eta,i}|\bm{\Theta^{(m)},s^{(m)},v^{(m)}})$ | $i=1,\ldots,N$
6) Draw $\phi_{h,i}^{(m+1)}$ from $p(\phi_{h,i}|\bm{\Theta^{(m)},s^{(m)},v^{(m)}})$ | $i=1,\ldots,N$
7) Draw $\phi_{\lambda,i}^{(m+1)}$ from $p(\phi_{\lambda,i}|\bm{\Theta^{(m)},s^{(m)},v^{(m)}})$ | $i=1,\ldots,N$
8) Draw $h_{i,0}^{(m+1)}$ from $p(h_{i,0}|\bm{\Theta}^{(m)},\bm{v}^{(m)},\bm{s}^{(m)})$ | $i=1,\ldots,N$
9) Draw $\\{h_{it}\\}_{t=1}^{T^{(m+1)}}$ from $p(h_{i1},\ldots,h_{iT}|\bm{\Theta^{(m)}},\bm{v^{(m)}},\bm{s}^{(m)})$ | $i=1,\ldots,N$
Particle step: |
10) Draw $\lambda_{i,0}^{(m+1)}$ from $\lambda_{i,0}^{(m+1)}$ from $p(\lambda_{i,0}|\bm{\Theta}^{(m)},\bm{v}^{(m)},\bm{s}^{(m)})$ | $i=1,\ldots,N$
11) Draw $\\{\lambda_{it}\\}_{t=1}^{T^{(m+1)}}$ from $p(\lambda_{i1},\ldots,\lambda_{iT}|\bm{\Theta^{(m)}},\bm{v^{(m)}},\bm{s}^{(m)})$ | $i=1,\ldots,N$
Particle step: |
end |

As in the univariate framework, it is easy to adapt the sampler to a version
of the VAR model with Skew-t shocks. In this case (9) becomes:

$\bm{\varepsilon_{t}}=\bm{\zeta_{t}}+\bm{\Omega_{t}}\bm{\Delta_{t}}\bm{O_{t}}^{-0.5}+{\bm{\Omega_{t}}}(\bm{I}_{n}-\bm{\Delta_{t}}^{2})^{0.5}\bm{O_{t}}^{-0.5}\bm{z_{t}}$
(16)

where $\bm{O_{t}}=diag(o_{1t}\ldots o_{Nt})\hskip 56.9055pto_{it}\sim
Gamma(\frac{\nu}{2},\frac{\nu}{2})$.
It is enough to adapt the Gibbs Sampler by adding another initial step to draw
the mixing variables $\\{o_{it}\\}_{t=1}^{T}$ for $i=1,\ldots,N$.

Draw $\\{o_{it}\\}_{t=1}^{T^{(m+1)}}$ from $p(o_{i1}\ldots,{o_{iT}}|\bm{\Theta}^{(m)},\bm{v}^{(m)},\bm{s_{t}}^{(m)})$ | $i=1,\ldots,N$
---|---

and then to update the formulas of the full conditional posterior
distributions in order to account for the extra terms. Again I use Metropolis
Hastings to simulate draws from
$p(o_{i1}\ldots,{o_{iT}}|\bm{\Theta},\bm{v},\bm{s_{t}})$ for $i=1,\ldots,N$,
since it is not directly possible to sample from these distributions.

##  3 Growth at Risk

The work of [2] (henceforth ABG) pioneered a recently growing body of
research, which examines the main sources of tail risk to GDP growth in
relationship to changes in economic and financial conditions. This section
compares the out of sample Growth-at-Risk (GaR) estimates for the U.S from our
time varying skewness stochastic volatility models to the two step approach
based on quantile regression of ABG . TVSSV models are able to reproduce the
main finding in ABG, which is that deteriorating financial conditions are
associated to shifts of the lower quantiles of the future GDP growth
distribution. As well, in the application, TVSSV models perform comparably if
not even slightly better than quantile regression based methods for
forecasting macroeconomic tail risk. In order to model asymmetric changes in
the conditional distribution of GDP growth as a function of changes in
financial conditions, I consider the following specification of the univariate
TVSSV model:

$\begin{array}[]{l}gdpgrowth_{t}=\pi_{0}+\pi_{1}gdpgrowth_{t-1}+\pi_{2}gdpgrowth_{t-2}+\pi_{3}NFCI_{t-1}+\sqrt{h}_{t}\varepsilon_{t}\\\
\hskip 227.62204pt\varepsilon_{t}\sim Skew-
Normal(\zeta_{t},\omega_{t},\lambda_{t})\\\ \hskip 284.52756pt\text{or}\\\
\hskip 227.62204pt\varepsilon_{t}\sim
Skew-t(\zeta_{t},\omega_{t},\lambda_{t},\nu)\end{array}$ (17)
$log(h_{t})=\phi_{h}log(h_{t-1})+\eta_{t}\hskip
71.13188pt\eta_{t}\sim\mathcal{N}(0,\sigma^{2}_{\eta})$ (18)
$\lambda_{t}=\phi_{\lambda}\lambda_{t-1}+\beta_{1}NFCI_{t-1}+\xi_{t}\hskip
56.9055pt\xi_{t}\sim\mathcal{N}(0,\sigma^{2}_{\xi})$ (19)

In this specification the NFCI directly affects the conditional skewness of
the future GDP growth distribution. More specifically, the coefficient
$\beta_{1}$ captures changes in the skewness of the conditional distribution
of GDP growth as a function of financial conditions. This coefficient is meant
to capture the non-linear relationship between deteriorating financial
conditions and future GDP growth distribution found in ABG. Since our focus is
to model the asymmetric effect of the NFCI on the future GDP growth
distribution, I threat the log-volatilities as exogenous autoregressive
processes, not affected by the NFCI. As a matter of fact, augmenting the state
equation for the log-volatilities with the NFCI index, as it is done in the
state equations of the shape parameters, implies that financial conditions
would affect symmetrically both a upper and the lower quantiles of the future
GDP growth distribution.

Together with the univariate model, I consider as well a bivariate TVSSV-
VAR(2) model where $\bm{y_{t}}=[gdpgrowth,NFCI]^{\prime}$ and:

$\begin{array}[]{l}\bm{y_{t}}=\bm{\Pi_{0}}+\bm{\Pi_{1}y_{t-1}}+\bm{\Pi_{2}y_{t-2}}+\bm{A}^{-1}\bm{H_{t}}^{0.5}\bm{\varepsilon_{t}}\\\
\hskip 142.26378pt\varepsilon_{it}\sim\text{{Skew-
Normal}}(\zeta_{it},\omega_{it}^{2},\lambda_{it})\\\ \hskip
199.16928pt\text{or}\\\ \hskip 142.26378pt\varepsilon_{it}\sim
Skew-t(\zeta_{it},\omega_{it},\lambda_{it},\nu)\end{array}$ (20)
$log(h_{it})=\phi_{h,i}log(h_{it-1})+\eta_{it}\hskip
71.13188pt\eta_{it}\sim\mathcal{N}(0,\sigma^{2}_{i,\eta})\hskip
5.69046pt\text{{i = gdpgrowth, NFCI}}$ (21)
$\lambda_{gdpgrowth,t}=\phi_{\lambda,1}\lambda_{gdpgrowth,t-1}+\beta_{1}NFCI_{t-1}+\xi_{gdpgrowth,t}\hskip
36.98866pt\xi_{i,t}\sim N(0,\sigma^{2}_{\xi,i})$ (22)
$\lambda_{NFCI,t}=\phi_{\lambda,2}\lambda_{NFCI,t-1}+\xi_{NFCI,t}\hskip
170.71652pt\xi_{i,t}\sim N(0,\sigma^{2}_{\xi,i})$ (23)

In this VAR, the dynamic relationship between GDP growth and financial
conditions in modelled jointly. In particular, in this specification, due to
the triangular structure of $\bm{A}^{-1}$ shocks to GDP growth
contemporaneously affect the financial markets, while shocks to NFCI do not
affect GDP growth within the quarter.

I compare the forecasts of downside risk from the TVSSV models to the
forecasts from the quantile regression based method of ABG. Their approach is
based on a two step procedure where in the first step they use predictive
quantile regression to estimate the quantiles of the conditional distribution:

$\hat{Q}_{gdpgrowth_{t+h}|\mathcal{I}_{t}}(\tau)=\bm{\hat{\beta}}^{\tau}X_{t}\hskip
56.9055ptfor\hskip 5.69046pt\tau=0.05,\ldots,0.95\ $ (24)

Then, in the second step, the estimated quantiles are interpolated using a
flexible Skew-t distribution, so as to obtain a complete predictive density
for GDP growth. We specify equation (24) collecting two lags of GDP growth and
one lag of NFCI in the vector $X_{t}$, so as to capture changes in the future
GDP growth distribution as a function of current financial and economic
conditions.

###  3.1 Results

This section presents the results from the estimates of both the univariate
TVSSV models and the VAR TVSSV models with Skew-Normal and Skew-t shocks. The
estimation sample starts in 1971Q1 and the forecasting exercise covers the
period 1995Q1 - 2019Q4. Fig. 1 presents the estimated posterior distribution
for the coefficient $\beta_{1}$ from the univariate time varying skewness
stochastic volatility model. This is the coefficient that in the state
equation of the skewness parameter (22) summarizes how the shape of the
conditional distribution of GDP growth changes as a function of financial
conditions in the previous quarter. As shown in Fig. 1, tighter financial
conditions (increases in the $NFCI$) are on average associated to a decrease
in the skewness of current GDP growth (the posterior mean estimate is
$\hat{\beta_{1}}=-0.26)$. Hence, equation (19), captures the main finding of
ABG, which is that deteriorating financial condition are associated to
movements in the lower quantiles of future GDP growth distribution.

Figure 1: Posterior estimate of $\beta_{1}$

Note: The figure shows the estimated posterior distribution of the coefficient
on the $NFCI_{t-1}$ in the equation of the shape parameter in the TVSSV model
with Skew-t shocks.

Ascertained that the model is able to capture the same asymmetric effect of
financial conditions on the future GDP growth distribution found in ABG, it is
important to understand what is the potential of the model to assess and
predict risk out of sample. Fig. 2 shows the out-of-sample forecasts of Growth
at Risk and Expected Shortfall for the $5^{th}$, $10^{th}$ and $20^{th}$
percentiles while Fig. 3 shows the one quarter ahead estimated recession
probability. We report the results from the stochastic volatility stochastic
skewness model with Skew-t shocks, since the results from the model with Skew-
Normal shocks do not differ qualitatively. As you can notice, during the
Financial Crisis, our parametric models predicts as much downside risk to GDP
growth as the quantile regression method of ABG. As shown in Fig. 3 both the
TVSSV and the TVSSV-VAR models, assign higher probability of recession to the
mild contraction of the U.S. economy following the dotcom bubble in 2000s with
respect to the two step method based on quantile regression.

Figure 2: One quarter ahead Growth at Risk (GaR) and Expected Shortfall
(1995Q1-2019Q4)

Note: The figure shows the estimated $5^{th}$, $10^{th}$, $20^{th}$
percentiles of the one quarter ahead GDP growth predictive distribution (left
panel) and the $5^{th}$, $10^{th}$, $20^{th}$ one quarter ahead expected
shortfall (right panel). In blue estimates from the two step quantile
regression based method by ABG, in red from the TVSSV univariate model with
Skew-t shocks and in yellow the estimates from the TVSSV VAR model.

Figure 3: One quarter ahead recession probability (1995Q1-2019Q4)

Note: The figure shows the estimated one quarter ahead recession
probabilities. In blue estimates from the two step quantile regression based
method by ABG, in red from the TVSSV univariate model with Skew-t shocks and
in yellow the estimates from the TVSSV VAR model.

In terms of forecast accuracy, Table 3 compares the forecasts from our
parametric approaches to the forecasts from the method of ABG. The first two
columns report the results for the average Log Scores and the average
Cumulative Ranked Probability Scores (CRPS), since these two measures are the
most commonly used to evaluate the relative density forecast accuracy of
different models. 131313Defining $y$ the realization of the series to predict,
$f(.)$ the density forecast and $F(.)$ corresponding the cumulative
distribution, Logscores and CRPS are respectively defined as:
$Logscores(f,y)=-log(f(y))$ (25)
$CRPS(f,y)=\int_{-\infty}^{\infty}PS(F(z),\mathbbm{1}\\{y\leq
z\\})dz=\int_{0}^{1}QS_{\alpha}(F^{-1}(\alpha),y)d\alpha$ (26) where
$PS(F(z),\mathbbm{1}\\{y\leq z\\})=(F(z)-\mathbbm{1}\\{y\leq z\\})^{2}$ is the
Brier probability score and $QS_{\alpha}(F^{-1}(\alpha),y)=2(\mathbbm{1}\leq
F^{-1}(\alpha))(F^{-1}(\alpha)-y)$ is the Quantile Score. Looking at average
Logscores, the first column reports the difference between the forecasts from
two step procedure of ABG and the forecasts from the time varying skewness
stochastic volatility models (values greater than zero are associated to more
accurate density forecast w.r.t ABG). According to the average Log-scores, our
parametric models provide more accurate one quarter ahead density forecasts
with respect to ABG. In parenthesis I report the p-values from the Diebold and
Mariano test [17] of equal forecast accuracy and find that for the TVSSV with
Skew-t shocks I am able to reject the null hypothesis of equal forecast
accuracy. For what concerns average CRPS, on the second column, the table
reports the ratio with respect to the model of ABG (values lower than 1 are
associated to more accurate density forecast with respect to ABG). As you can
notice, based on this metrics, the time varying skewness stochastic volatility
models perform as good if not even better than the two step procedure based on
quantile regression. However, in all the cases I am not able to reject the
null of equal forecast accuracy.

Since I aim to assess the ability of the model to correctly characterize
downside risk predictions, on the third column I report the average Quantile
Weighted CRPS introduced by [21] 141414The Quantile Weighted CRPS are computed
as:  $twCRPS=\int_{-\infty}^{\infty}PS(F(z),\mathbbm{1}\\{y\leq
z\\})^{2}w(z)dz=\int_{0}^{1}QS_{\alpha}(F^{-1}(\alpha),y)v(\alpha)d\alpha$
(27) where $v(\alpha)=(1-\alpha)^{2}$ assigns higher weights to the lower
quantiles of the distribution function. and on the fourth, fifth and sixth
column I report the average Quantile Scores for the $5^{th}$, $10^{th}$ and
$20^{th}$ percentiles commonly associated with the tick loss function [19].
Also in this case I report the ratio with respect to the two step approach
based on quantile regression (values lower than 1 are associated to more
accurate density forecast with respect to ABG) and the p-values from the
Diebold-Mariano test in parenthesis. As you can notice, in terms of the
ability of the model to correctly characterize downside risk predictions, I
find that the stochastic volatility models performs comparably if not even
better than ABG. In particular for the TVSSV-VAR with Skew-t shocks I am able
to reject the null of equal forecast accuracy with respect to ABG. The time
series with the CRPS and left Tail Weighted CRPS, can be found in the Appendix
B.3 (Fig. 6). As well, in the Appendix B.3 the histogram with the PITs (Fig.
7) reveals that the forecasts from the TVSSV models, are better-calibrated
with respect to the forecasts from the two-step quantile regression based
method.

Table 3: One quarter ahead out of sample forecasts (1995Q1-2019Q4) | Log scores | CRPS | TwL CRPS | $\bm{QS5^{th}}$ | $\bm{QS10^{th}}$ | $\bm{QS20^{th}}$
---|---|---|---|---|---|---
ABG | 2.4840 | 1.1943 | 0.3623 | 0.2503 | 0.3926 | 0.5842
TVSSV Skew Normal | 0.1946 | 0.9757 | 0.9777 | 1.0334 | 0.9858 | 0.9832
| (0.1788) | (0.7380) | (0.2282) | (0.6027) | (0.3838) | (0.3278)
TVSSV Skew-t | 0.3530 | 0.9609 | 0.9823 | 1.0334 | 0.9911 | 0.9882
| (0.0276) | (0.8659) | (0.3045) | (0.4725) | (0.4472) | (0.3934)
TVSSV VAR Skew Normal | 0.1620 | 0.9805 | 0.9644 | 0.9678 | 0.9718 | 0.9666
| (0.2287) | (0.6781) | (0.1669) | (0.4038) | (0.3209) | (0.2589)
TVSSV VAR Skew-t | 0.0662 | 0.9700 | 0.9610 | 1.0119 | 0.9979 | 0.9633
| (0.1050) | (0.5942) | (0.0366) | (0.5493) | (0.4778) | (0.1052)
Note: For the average Logscores, the first row reports the values from the ABG
method while the
other rows report the difference between the two step procedure and the time
varying skewness
stochastic volatility models. For the other metrics I report the ratio w.r.t
the ABG method.
Inside the parenthesis p-values from the one sided Diebold-Mariano w.r.t the
two step method
of [2]. The bold character indicates rejection of equal forecast accuracy at
5%.

##  4 Time varying skewness in a medium scale VAR

One of the main advantages of the VAR model presented in Section 2.2 is that
it allows to explicitly capture time varying conditional skewness of multiple
time series. In this section I estimate a medium scale VAR model which
includes macroeconomic and financial monthly time series and I investigate the
time varying asymmetric behaviour of the shocks to the variables in the
system. I consider a VAR model with 8 variables being Real personal
consumption expenditures, Industrial Production, Unemployment Rate, average
Weekly Hours Worked, Consumer Price Index, Fed Funds Rate, the spread between
10-Year Treasury and the Fed Funds Rate, the spread between Moody’s Baa
Corporate Bond and the Fed Funds Rate and the Standard and Poors Index. The
variables are in monthly frequency and are taken from the FRED-MD.151515Table
7 in the Appendix reports the variable transformation. I present the results
from the VAR with Skew-t shocks.161616For the VAR with Skew-Normal the
estimated path for the volatilities and shape parameters are almost the same.
I include 13 lags and assume a Minnesota prior structure for the variance
covariance matrix of the regression coefficients. 171717See the Appendix for
the details on the hyper-parameters of the Minnetota Prior.  The estimation
sample is January 1965 - December 2019. Fig. 4 shows the estimated
volatilities while Fig. 5 shows the estimated shape parameters. The dotted
line in blue are the $85^{th}-15^{th}$ credible sets while the red line is the
estimated posterior median.

It is interesting to notice that shocks to the CPI were on average positively
skewed before the 2000s while became left skew for the rest of the sample that
ends on 2019. This switch in the sign of the shape parameter indicates that
conditionally on the past and on the contemporaneous realization of Real
Personal Consumption Expenditures, Industrial Production, Unemployment Rate
and average Weekly Hours Worked, the distribution of CPI was right skewed in
the 1980s, becoming instead left skew from the 2000s. In other words, risk
switched from the upside to the downside. As for the monetary policy shocks,
in the 1980s large positive hikes of the Fed Fund Rate were more frequent,
while from the early 2000s large negative shocks to the Fed Fund Rate become
more likely. Shocks to the average Weekly Hours Worked are skewed to the left
over the entire sample, which means that negative large shocks have been
systematically more frequent than positive large shocks. As well, shocks to
the stock market (SP 500 index) are skewed to the left over the entire sample.
This is in line with the large body of the financial econometrics literature
that studies conditional skewness in asset returns [22]. As for the spread
between 10-Year Treasury and the Fed Funds Rate and the spread between Moody’s
Baa Corporate Bond and the Fed Funds Rate, for most of the sample both the
shocks are skewed to the right meaning that the probability of large positive
shocks has been greater than the probability of large negative shocks. This
finding vanishes starting from 2009 and might be linked to the unconventional
monetary policy following the Great Financial Crisis.

Figure 4: Estimated volatilities

Note: The figure shows the estimated volatilities of the shocks in the TVSSV-
VAR with Skew-t shocks. In red the estimated median, in blue dashed lines the
$85^{th}-15^{th}$ credible sets.

Figure 5: Estimated shape parameters $\lambda_{t}$

Note: The figure shows the estimated shape parameters of the shocks in the
TVSSV-VAR with Skew-t shocks. In red the estimated median, in blue dashed
lines the $85^{th}-15^{th}$ credible sets.

##  5 Conclusion

In this paper I propose a fully parametric framework based on time varying
skewness stochastic volatility models with Skew-Normal and Skew-t shocks for
assessing and forecasting macroeconomic tail risk. First, I consider an
extension of the univariate stochastic volatility model of [25] that
explicitly accounts for time varying skewness in the predictive distribution
of the dependent variable. Then, I introduce a Bayesian VAR model with
stochastic volatility and stochastic skewness to provide an explicit treatment
of conditional skewness when modelling the dynamics of multiple time series. I
compare the time varying skewness stochastic volatility models to the quantile
regression method of [2] to assess and predict tail risk to GDP growth. I find
that the time varying skewness stochastic volatility models considered in this
paper are able to reproduce the main findings of [2], that is the nonlinear
and asymmetric effect of financial conditions on the future GDP growth
distribution. The models predict as much risk as quantile regression during
the Financial crisis while provide slightly more accurate out of sample
forecasts of downside risk over the entire sample. Finally, estimating a
standard medium scale VAR model I find that time varying skewness is a
relevant feature of macroeconomic and financial shocks.

For future research, the VAR model considered in this paper could be used to
study the probability of joint tail events and for constructing structural
scenarios of “at-risk” measures. For example, it could be used to study and
assess stagflation risk, or to analyze scenarios for inflation at-risk and
labour-at risk under different monetary policy paths. As a methodological
extension, particularly interesting would be to consider an order invariant
version of this model.

## References

  * [1] CA Abanto-Valle, VH Lachos and Dipak K Dey  “Bayesian estimation of a skew-student-t stochastic volatility model”  In _Methodology and Computing in Applied Probability_ 17.3 Springer, 2015, pp. 721–738
  * [2] Tobias Adrian, Nina Boyarchenko and Domenico Giannone  “Vulnerable Growth”  In _American Economic Review_ 109.4, 2019, pp. 1263–89  DOI: 10.1257/aer.20161923
  * [3] Christophe Andrieu, Arnaud Doucet and Roman Holenstein  “Particle markov chain monte carlo methods”  In _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 72.3 Wiley Online Library, 2010, pp. 269–342
  * [4] Jonas E. Arias, Juan F. Rubio-Ramirez and Minchul Shin  “Macroeconomic Forecasting and Variable Ordering in Multivariate Stochastic Volatility Models”, 2021  DOI: 10.21799/frbp.wp.2021.21
  * [5] Adelchi Azzalini  “Further results on a class of distributions which includes the normal ones”  In _Statistica_ 46.2, 1986, pp. 199–208
  * [6] Adelchi Azzalini and Antonella Capitanio  “Distributions generated by perturbation of symmetry with emphasis on a multivariate skew t-distribution”  In _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 65.2, 2003, pp. 367–389  DOI: https://doi.org/10.1111/1467-9868.00391
  * [7] Christian Brownlees and André B.M. Souza  “Backtesting global Growth-at-Risk”  In _Journal of Monetary Economics_ 118, 2021, pp. 312–330  DOI: https://doi.org/10.1016/j.jmoneco.2020.11.003
  * [8] Nunzio Cappuccio, Diego Lubian and Davide Raggi  “MCMC Bayesian Estimation of a Skew-GED Stochastic Volatility Model”  In _Studies in Nonlinear Dynamics & Econometrics_ 8.2, 2004  DOI: doi:10.2202/1558-3708.1211
  * [9] Andrea Carriero, Joshua Chan, Todd E. Clark and Massimiliano Marcellino  “Corrigendum to “Large Bayesian vector autoregressions with stochastic volatility and non-conjugate priors” [J. Econometrics 212 (1) (2019) 137–154]”  In _Journal of Econometrics_ 227.2, 2022, pp. 506–512  DOI: https://doi.org/10.1016/j.jeconom.2021.11.010
  * [10] Andrea Carriero, Todd E Clark and Massimiliano Marcellino  “Large Bayesian vector autoregressions with stochastic volatility and non-conjugate priors”  In _Journal of Econometrics_ 212.1 Elsevier, 2019, pp. 137–154
  * [11] Andrea Carriero, Todd E. Clark and Massimiliano Marcellino  “Realtime nowcasting with a Bayesian mixed frequency model with stochastic volatility”  In _Journal of the Royal Statistical Society: Series A (Statistics in Society)_ 178.4, 2015, pp. 837–862  DOI: https://doi.org/10.1111/rssa.12092
  * [12] Andrea Carriero, Todd E. Clark and Massimiliano Marcellino  “Capturing Macroeconomic Tail Risks with Bayesian Vector Autoregressions”, 2020  DOI: 10.26509/frbc-wp-202002r
  * [13] Joshua CC Chan, Gary Koop and Xuewen Yu  “Large order-invariant Bayesian VARs with stochastic volatility”  In _arXiv preprint arXiv:2111.07225_ , 2021
  * [14] Todd E. Clark and Francesco Ravazzolo  “Macroeconomic Forecasting Performance under Alternative Specifications of Time-Varying Volatility”  In _Journal of Applied Econometrics_ 30.4, 2015, pp. 551–575  DOI: https://doi.org/10.1002/jae.2379
  * [15] Timothy Cogley and Thomas J. Sargent  “Drifts and volatilities: monetary policies and outcomes in the post WWII US” Monetary Policy and Learning  In _Review of Economic Dynamics_ 8.2, 2005, pp. 262–302  DOI: https://doi.org/10.1016/j.red.2004.10.009
  * [16] Davide Delle Monache, Andrea De Polis and Ivan Petrella  “Modeling and forecasting macroeconomic downside risk”, 2021  URL: https://ideas.repec.org/p/bdi/wptemi/td_1324_21.html
  * [17] Francis Diebold and Roberto Mariano  “Comparing Predictive Accuracy”  In _Journal of Business & Economic Statistics_ 13.3, 1995, pp. 253–63  URL: https://EconPapers.repec.org/RePEc:bes:jnlbes:v:13:y:1995:i:3:p:253-63
  * [18] Gaston Gelos et al.  “Capital flows at risk: Taming the ebbs and flows”  In _Journal of International Economics_ 134, 2022, pp. 103555  DOI: https://doi.org/10.1016/j.jinteco.2021.103555
  * [19] Raffaella Giacomini and Ivana Komunjer  “Evaluation and Combination of Conditional Quantile Forecasts”  In _Journal of Business & Economic Statistics_ 23.4 [American Statistical Association, Taylor & Francis, Ltd.], 2005, pp. 416–431  URL: http://www.jstor.org/stable/27638838
  * [20] Stefano Giglio, Bryan Kelly and Seth Pruitt  “Systemic risk and the macroeconomy: An empirical evaluation”  In _Journal of Financial Economics_ 119.3, 2016, pp. 457–471  DOI: https://doi.org/10.1016/j.jfineco.2016.01.010
  * [21] Tilmann Gneiting and Roopesh Ranjan  “Comparing Density Forecasts Using Threshold-and Quantile-Weighted Scoring Rules”  In _Journal of Business & Economic Statistics_ 29.3 American Statistical Association, 2011, pp. 411–422  URL: http://www.jstor.org/stable/23243806
  * [22] Campbell R Harvey and Akhtar Siddique  “Conditional skewness in asset pricing tests”  In _The Journal of finance_ 55.3 Wiley Online Library, 2000, pp. 1263–1295
  * [23] Martin Iseringhausen  “The time-varying asymmetry of exchange rate returns: A stochastic volatility – stochastic skewness model”  In _Journal of Empirical Finance_ 58, 2020, pp. 275–292  DOI: https://doi.org/10.1016/j.jempfin.2020.06.008
  * [24] Martin Iseringhausen  “A time-varying skewness model for Growth-at-Risk”  European Stability Mechanism Working Paper, 2021
  * [25] Eric Jacquier, Nicholas G Polson and Peter Rossi  “Bayesian Analysis of Stochastic Volatility Models”  In _Journal of Business & Economic Statistics_ 12.4, 1994, pp. 371–89  URL: https://EconPapers.repec.org/RePEc:bes:jnlbes:v:12:y:1994:i:4:p:371-89
  * [26] Eric Jacquier, Polson and Peter Rossi  “Bayesian analysis of stochastic volatility models with fat-tails and correlated errors”  In _Journal of Econometrics_ 122.1, 2004, pp. 185–212  URL: https://EconPapers.repec.org/RePEc:eee:econom:v:122:y:2004:i:1:p:185-212
  * [27] Sune Karlsson, Stepan Mazur and Hoang Nguyen  “Vector autoregression models with skewness and heavy tails”, 2021  URL: https://ideas.repec.org/p/hhs/oruesi/2021_008.html
  * [28] Michael T Kiley  “Unemployment risk”  In _Journal of Money, Credit and Banking_ Wiley Online Library, 2018
  * [29] Lutz Kilian and Simone Manganelli  “The central bank as a risk manager: quantifying and forecasting inflation risks”, 2003  URL: https://ideas.repec.org/p/ecb/ecbwps/2003226.html
  * [30] Fredrik Lindsten, Michael I Jordan and Thomas B Schon  “Particle Gibbs with ancestor sampling”  In _Journal of Machine Learning Research_ 15 MICROTOME PUBL, 2014, pp. 2145–2184
  * [31] J. López-Salido and Francesca Loria  “Inflation at Risk”, 2020  DOI: 10.17016/FEDS.2020.013
  * [32] Carlos Montes-Galdón and Eva Ortega  “Skewed SVARS: Tracking the structural sources of macroeconomic tail risks”  In _Essays in Honour of Fabio Canova_ 44 Emerald Publishing Limited, 2022, pp. 177–210
  * [33] Giorgio E Primiceri  “Time varying structural vector autoregressions and monetary policy”  In _The Review of Economic Studies_ 72.3 Wiley-Blackwell, 2005, pp. 821–852
  * [34] Christopher A Sims  “Macroeconomics and reality”  In _Econometrica: journal of the Econometric Society_ JSTOR, 1980, pp. 1–48
  * [35] James H. Stock and Mark W. Watson  “Vector Autoregressions”  In _Journal of Economic Perspectives_ 15.4, 2001, pp. 101–115  DOI: 10.1257/jep.15.4.101
  * [36] Elias Wolf  “Estimating growth at risk with skewed stochastic volatility models”  In _Available at SSRN 4030094_ , 2021

##  Appendix A Appendix

###  A.1 Skew Normal and Skew-t: distributions and parameterization

The Skew Normal [5] distribution is:

$\displaystyle
p(\varepsilon_{t}|\zeta,\omega^{2},\lambda)=\frac{2}{\omega}\phi\left(\frac{\varepsilon_{t}-\zeta}{\omega}\right)\Phi\left(\lambda\left(\frac{\varepsilon_{t}-\zeta}{\omega}\right)\right)$

where $\phi(.)$ and $\Phi(.)$ are respectively the pdf and cdf of the standard
Normal. In general, $\varepsilon_{t}\sim Skew-
Normal(\zeta,\omega^{2},\lambda)$ has the following stochastic representation:

$\varepsilon_{t}=\zeta+\delta\omega v_{t}+\sqrt{(1-\delta^{2})}\omega z_{t}$
(28)

where:

$v_{t}\stackrel{{\scriptstyle i.i.d}}{{\sim}}\text{Truncated
Normal}_{[0,\infty)}(0,1)$

$z_{t}\stackrel{{\scriptstyle i.i.d}}{{\sim}}\mathcal{N}(0,1)$

$\delta=\frac{\lambda}{\sqrt{1+\lambda^{2}}}$, with $-1<\delta<1$.

The mean and the variance of $\varepsilon_{t}$ are given by:

$\mathbb{E}[\varepsilon_{t}]=\zeta+\omega\delta\sqrt{\frac{2}{\pi}}$ (29)
$var(\varepsilon_{t})=\omega^{2}\left(1-\frac{2\delta^{2}}{\pi}\right)$ (30)

Assuming $\operatorname*{\mathbb{E}}[\varepsilon_{t}]=0$ and
$var(\varepsilon_{t})=1$ leads to the following constraints on the location
and scale parameters: $\zeta=-\omega\delta\sqrt{\frac{2}{\pi}}$ and
$\omega^{2}=\left(1-\frac{2\delta^{2}}{\pi}\right)^{-1}$. Once we impose these
constraints on the location and scale parameters, with $\lambda=0$ the
distribution collapses to the Standard Normal.

The Skew-t distribution [6] is:

$p(\varepsilon_{t}|\zeta,\omega^{2},\lambda,\nu)=\frac{2}{\omega}t_{\nu}\left(\frac{\varepsilon_{t}-\zeta}{\omega}\right)T_{\nu+1}\left(\lambda\left(\frac{\varepsilon_{t}-\zeta}{\omega}\right)\sqrt{\frac{\nu+1}{\nu\left(\frac{\varepsilon_{t}-\zeta}{\omega}\right)^{2}}}\right)\\\
$ (31)

where $t(.)$ and $T(.)$ are respectively the pdf and cdf of the Student-t with
$\nu$ degrees of freedom.
$\varepsilon_{t}\sim\text{{Skew-t}}(\zeta,\omega^{2},\lambda,\nu)$ has the
following stochastic representation:

$\varepsilon_{t}=\zeta+\delta\omega
o_{t}^{-0.5}v_{t}+\sqrt{(1-\delta^{2})}\omega o_{t}^{-0.5}z_{t}$ (32)

where:

$v_{t}\stackrel{{\scriptstyle i.i.d}}{{\sim}}\text{Truncated
Normal}_{[0,\infty)}(0,1)$

$z_{t}\stackrel{{\scriptstyle i.i.d}}{{\sim}}\mathcal{N}(0,1)$

$o_{t}\stackrel{{\scriptstyle
i.i.d}}{{\sim}}\mathcal{G}(\frac{\nu}{2},\frac{\nu}{2})$

$\delta=\frac{\lambda}{\sqrt{1+\lambda^{2}}}$, with $-1<\delta<1$.

The mean and the variance of $\varepsilon_{t}$ are given by:

$\mathbb{E}[\varepsilon_{t}]=\omega\delta k_{1}\sqrt{\frac{2}{\pi}}$ (33)
$var(\varepsilon_{t})=\omega^{2}\left(k_{2}-\frac{2}{\pi}k_{1}^{2}\delta^{2}\right)$
(34)

with
$k_{1}=\sqrt{\frac{\nu}{2}}\frac{\Gamma(\frac{\nu-1}{2})}{\Gamma(\frac{\nu}{2})}$,
$k_{2}=\frac{\nu}{\nu-2}$.

Assuming $\operatorname*{\mathbb{E}}[\varepsilon_{t}]=0$ and
$var(\varepsilon_{t})=1$ leads to the following constraints on the location
and scale parameters: $\zeta=-\omega\delta k_{1}\sqrt{\frac{2}{\pi}}$ and
$\omega^{2}=\left(k_{2}-\frac{2}{\pi}k_{1}^{2}\delta^{2}\right)^{-1}$. Once we
impose these constraints on the location and scale parameters with $\lambda=0$
the distribution collapses to a Student-t distribution properly re-scaled to
have unit variance (and zero mean).

###  A.2 Full conditional posterior distributions

####  A.2.1 Univariate time varying skewness stochastic volatility model: skew
normal shocks

The full conditional distribution of $\\{v_{t}\\}_{t=1}^{T}$ is given by:

$\scriptsize p(v_{t}|.)\propto
exp\left[-\frac{1}{2}\left(\frac{1}{(1-\delta_{t}^{2})}v_{t}^{2}-\frac{2\delta_{t}h_{t}^{-0.5}}{\omega_{t}(1-\delta_{t}^{2})}(y_{t}-\bm{x_{t}}\bm{\pi}-\zeta_{t}\sqrt{h_{t}})v_{t}\right)\right]\mathbb{I}(0\leq
v_{t}<\infty)$ (35)

this is a Truncated
Normal$\left(\frac{\delta_{t}h_{t}^{-0.5}[y_{t}-\bm{x_{t}}\bm{\pi}]-\delta_{t}\zeta_{t}}{\omega_{t}},1-\delta_{t}^{2}\right)_{[0,\infty)}$

The full conditional distribution of $\bm{\pi}$ is Normal:

$f(\bm{\pi}|.)\sim N(\overline{\mu}_{\pi},\overline{\Sigma}_{\pi})\\\ $ (36)
$\begin{array}[]{l}\overline{\mu}_{\pi}=\overline{\Sigma}_{\pi}\left(\sum_{t=1}^{T}\frac{1}{\sigma^{2}_{t}}\bm{x_{t}}^{\prime}\tilde{y}_{t}+\underline{\Sigma_{\pi}}^{-1}\underline{\mu_{\pi}}\right)\\\
\overline{\Sigma}_{\pi}^{-1}=\underline{\Sigma_{\pi}}^{-1}+\sum_{t=1}^{T}\frac{1}{\sigma^{2}_{t}}\bm{x_{t}}^{\prime}\bm{x_{t}}\\\
\end{array}$

where

$\tilde{y}_{t}\equiv
y_{t}-\sqrt{h_{t}}\zeta_{t}-\sqrt{h_{t}}\omega_{t}\delta_{t}v_{t}$

$\sigma^{2}_{t}\equiv h_{t}\omega^{2}_{t}(1-\delta_{t}^{2})$
while $\underline{\mu_{\pi}}$ and $\underline{\Sigma_{\pi}}$ are the prior
mean and variance covariance matrix.

The full conditional distribution of $\phi_{h}$ is a Normal :

$f(\phi_{h}|.)\sim\mathcal{N}(\bar{\mu}_{\phi_{h}},\bar{\sigma}^{2}_{\phi_{h}})$
(37)
$\bar{\sigma}^{2}_{\phi_{h}}=\left(\sum_{t=1}^{T}\frac{log(h_{t-1})^{2}}{\sigma^{2}_{\eta}}+\frac{1}{\underline{\sigma^{2}_{\phi_{h}}}}\right)^{-1}$
(38)

$\bar{\mu}_{\phi_{h}}=\bar{\sigma}^{2}_{\phi_{h}}\left(\frac{\underline{\mu}_{\phi_{h}}}{\underline{\sigma^{2}_{\phi_{h}}}}+\sum_{t=1}^{T}\frac{log(h_{t-1})log(h_{t})}{\sigma^{2}_{\eta}}\right)$
(39)

where $\underline{\mu}_{\phi_{h}}$ and $\underline{\sigma}^{2}_{\phi_{h}}$ are
prior mean and variance.

The full conditional distribution of $\sigma^{2}_{\eta}$ is an Inverse Gamma :

$\scriptsize
p(\sigma^{2}_{\eta}|.)\propto\left(\frac{1}{\sigma_{\eta}^{2}}\right)^{\frac{T}{2}}exp\left[\sum_{t=1}^{T}-\frac{1}{2\sigma_{\eta,i}^{2}}(ln(h_{t})-\phi_{h}ln(h_{t-1}))^{2}\right]exp\left[-\frac{s_{\sigma_{\eta}^{2}}}{\sigma_{\eta}^{2}}\right]\sigma_{\eta}^{2^{-\nu_{\sigma_{\eta}^{2}}}-1}$
(40)

where $s_{\sigma_{\eta}^{2}}$ and $\nu_{\sigma_{\eta}}^{2}$ are the hyper-
parameters of the Inverse Gamma prior.

The full conditional distribution of $\phi_{\lambda}$ is a Normal :

$f(\phi_{\lambda}|.)\sim\mathcal{N}(\bar{\mu}_{\phi_{\lambda}},\bar{\sigma}^{2}_{\phi_{\lambda}})$
(41)
$\bar{\sigma}^{2}_{\phi_{\lambda}}=\left(\sum_{t=1}^{T}\frac{\lambda_{t-1}^{2}}{\sigma^{2}_{\xi}}+\frac{1}{\underline{\sigma^{2}_{\phi_{\lambda}}}}\right)^{-1}$
(42)

$\bar{\mu}_{\phi_{\lambda}}=\bar{\sigma}^{2}_{\phi_{\lambda}}\left(\frac{\underline{\mu}_{\phi_{\lambda}}}{\underline{\sigma^{2}_{\phi_{\lambda}}}}+\sum_{t=1}^{T}\frac{{\lambda}_{t-1}{\lambda_{t}}}{\sigma^{2}_{\xi}}\right)$
(43)

where $\underline{\mu}_{\phi{\lambda}}$ and
$\underline{\sigma}^{2}_{\phi_{\lambda}}$ are prior mean and variance.

The full conditional distribution of $\sigma^{2}_{\xi}$ is an Inverse Gamma:

$\scriptsize
p(\sigma^{2}_{\xi}|.)\propto\left(\frac{1}{\sigma^{2}_{\xi}}\right)^{\frac{T}{2}}exp\left[\sum_{t=1}^{T}-\frac{1}{2\sigma^{2}_{\xi}}(\lambda_{it}-\phi_{\lambda}\lambda_{it-1})^{2}\right]exp\left[-\frac{s_{\sigma_{\xi}^{2}}}{\sigma_{\xi}^{2}}\right]\sigma_{\xi}^{2^{-}{\nu_{\sigma_{\xi}^{2}}}-1}$
(44)

where $s_{\sigma_{\xi}^{2}}$ and $\nu_{\sigma_{\xi}^{2}}$ are the hyper-
parameters of the Inverse Gamma prior.

The full conditional distribution of $h_{0}$ is
$\mathcal{N}\left(\bar{\mu}_{h0},\bar{\sigma}_{h0}\right)$

$\bar{\mu}_{h0}=\bar{\sigma}_{h0}\left(\frac{\mu_{h0}}{\sigma^{2}_{h0}}+\frac{\frac{log(h_{1})}{\phi_{h}}}{\frac{\sigma^{2}_{\eta}}{\phi_{h}^{2}}}\right)$
(45)

$\bar{\sigma}_{h0}=\frac{\frac{\sigma^{2}_{h_{0}}\sigma^{2}_{\xi}}{\phi^{2}_{h}}}{\sigma^{2}_{h0}+\frac{\sigma^{2}_{\eta}}{\phi^{2}_{h}}}$
(46)

where $\mu_{\lambda_{0}}$ and $\sigma^{2}_{\lambda_{0}}$ are the prior mean
and variance.

The full conditional distribution of $\bm{h}$ is given by:

$\scriptsize p(\bm{h}|.)=\prod_{t=1}^{T}p(h_{t}|h_{t-1},h_{t+1},.)\\\ $ (47)

$\scriptsize p(h_{t}|.)\propto
h_{t}^{-1,5}exp\left[\left(-\frac{1}{2}\left(\frac{y_{t}-\bm{x_{t}}\bm{\pi}-\sqrt{h_{t}}\zeta_{t}-\sqrt{h_{t}}\omega_{t}\delta_{t}v_{t}}{\sqrt{h_{t}}\omega_{t}(1-\delta_{t}^{2})^{0.5}}\right)^{2}-\frac{1}{2}\frac{(lnh_{t}-\mu_{h_{t}})^{2}}{\sigma_{h_{t}}^{2}}\right)\right]$
(48)

$\begin{array}[]{l}\mu_{h_{t}}=\frac{\phi_{h}}{\phi_{h}^{2}+1}\left(lnh_{t+1}+lnh_{t-1}\right)\\\
\sigma_{h}^{2}=\frac{\sigma^{2}_{\eta}}{\phi_{h}^{2}+1}\end{array}$ (49)

The full conditional distribution of $\lambda_{0}$ is
$\mathcal{N}\left(\bar{\mu}_{\lambda 0},\bar{\sigma}_{\lambda 0}\right)$

$\bar{\mu}_{\lambda 0}=\bar{\sigma}_{\lambda
0}\left(\frac{\mu_{\lambda_{0}}}{\sigma^{2}_{\lambda_{0}}}+\frac{\frac{\lambda_{1}}{\phi_{\lambda}}}{\frac{\sigma^{2}_{\xi}}{\phi_{\lambda}^{2}}}\right)$
(50)

$\bar{\sigma}_{\lambda
0}=\frac{\frac{\sigma^{2}_{\lambda_{0}}\sigma^{2}_{\xi}}{\phi^{2}_{\lambda}}}{\sigma^{2}_{\lambda_{0}}+\frac{\sigma^{2}_{\xi}}{\phi^{2}_{\lambda}}}$
(51)

where $\mu_{\lambda_{0}}$ and $\sigma^{2}_{\lambda_{0}}$ are the prior mean
and variance.

The full conditional distribution of $\bm{\lambda}$ is given by:

$\scriptsize
p(\bm{\lambda}|.)=\prod_{t=1}^{T}p(\lambda_{t}|\lambda_{t-1},\lambda_{t+1},.)\\\
$ (52)

$\scriptsize
p(\lambda_{t}|.)\propto\omega_{t}^{-1}(1-\delta_{t}^{2})^{-0,5}exp\left[-\frac{1}{2}\left(\frac{y_{t}-\bm{x_{t}}\pi-\sqrt{h_{t}}\zeta_{t}-\sqrt{h_{t}}o_{t}^{-0.5}\omega_{t}\delta_{t}v_{t}}{\sqrt{h_{t}}o_{t}^{-0.5}\omega_{t}(1-\delta_{t}^{2})^{0.5}}\right)^{2}-\frac{1}{2}\frac{(\lambda_{t}-\mu_{\lambda_{t}})^{2}}{\sigma^{2}_{\lambda_{t}}}\right]$
(53)

$\begin{array}[]{l}\mu_{\lambda_{t}}=\frac{\phi_{\lambda}}{\phi_{\lambda}^{2}+1}\left(\lambda_{t+1}+\lambda_{t-1}\right)\\\
\sigma_{\lambda}^{2}=\frac{\sigma^{2}_{\xi}}{\phi_{\lambda}^{2}+1}\end{array}$
(54)

####  A.2.2 Univariate time varying skewness stochastic volatility model:
skew-t shocks

The full conditional distribution of $\\{v_{t}\\}_{t=1}^{T}$ is given by:

$\scriptsize p(v_{t}|.)\propto
exp\left[-\frac{1}{2}\left(\frac{1}{(1-\delta_{t}^{2})}v_{t}^{2}-\frac{2o_{t}^{0.5}\delta_{t}h_{t}^{-0.5}}{\omega_{t}(1-\delta_{t}^{2})}(y_{t}-\bm{x_{t}\pi}-\zeta_{t}\sqrt{h_{t}})v_{t}\right)\right]\mathbb{I}(0\leq
v_{t}<\infty)$ (55)

this is a truncated normal
$N\left(\frac{\delta_{t}o_{t}^{0,5}h_{t}^{-0.5}[y_{t}-\bm{x_{t}\pi}]-\delta_{t}o_{t}^{0.5}\zeta_{t}}{\omega_{t}},1-\delta_{t}^{2}\right)_{[0,\infty)}$

The full conditional distribution of $\\{o_{t}\\}_{t=1}^{T}$ is given by:

$\scriptsize p(o_{t}|.)\propto
o_{t}^{\frac{\nu+1}{2}-1}exp\left[-\frac{o_{t}}{2}\left(\nu_{t}+\frac{h_{t}^{-1}(y_{t}-\bm{x_{t}}\bm{\pi}-\sqrt{h_{t}}\zeta_{t})^{2}}{\omega_{t}^{2}(1-\delta_{t}^{2})}\right)\right]exp\left[\frac{(y_{t}-\bm{x_{t}}\bm{\pi}-\sqrt{h_{t}}\zeta_{t})(h_{t}^{-0.5}o_{t}^{0.5}\delta_{t}v_{t})}{\omega_{t}(1-\delta_{t}^{2})}\right]$
(56)

The full conditional distribution of $\bm{\pi}$ is Normal:

$f(\bm{\pi}|.)\sim N(\overline{\mu}_{\pi},\overline{\Sigma}_{\pi})\\\ $ (57)
$\begin{array}[]{l}\overline{\mu}_{\pi}=\overline{\Sigma}_{\pi}\left(\sum_{t=1}^{T}\frac{1}{\sigma^{2}_{t}}\bm{x_{t}}^{\prime}\tilde{y}_{t}+\underline{\Sigma_{\pi}}^{-1}\underline{\mu_{\pi}}\right)\\\
\overline{\Sigma}_{\pi}^{-1}=\underline{\Sigma_{\pi}}^{-1}+\sum_{t=1}^{T}\frac{1}{\sigma^{2}_{t}}\bm{x_{t}}^{\prime}\bm{x_{t}}\\\
\end{array}$

where:

$\tilde{y}_{t}\equiv
y_{t}-\sqrt{h_{t}}\zeta_{t}-\sqrt{h_{t}}o_{t}^{-0.5}\omega_{t}\delta_{t}v_{t}$

$\sigma^{2}_{t}\equiv h_{t}\omega^{2}_{t}o_{t}^{-1}(1-\delta_{t}^{2})$

The full conditional distribution of $h_{t}$ is given by:

$\scriptsize p(\bm{h}|.)=\prod_{t=1}^{T}p(h_{t}|h_{t-1},h_{t+1},.)$ (58)

$\scriptsize p(h_{t}|.)\propto
h_{t}^{-1,5}exp\left[\left(-\frac{1}{2}\left(\frac{y_{t}-\bm{x_{t}}\bm{\pi}-\sqrt{h_{t}}\zeta_{t}-\sqrt{h_{t}}o_{t}^{-0.5}\omega_{t}\delta_{t}v_{t}}{\sqrt{h_{t}}o_{t}^{-0.5}\omega_{t}(1-\delta_{t}^{2})^{0.5}}\right)^{2}-\frac{1}{2}\frac{(lnh_{t}-\mu_{h_{t}})^{2}}{\sigma_{h_{t}}^{2}}\right)\right]$
(59)

$\displaystyle{l}$
$\displaystyle\mu_{h_{t}}=\frac{\phi_{h}}{\phi_{h}^{2}+1}\left(lnh_{t+1}+lnh_{t-1}\right)$
(60) $\displaystyle\sigma_{h}^{2}=\frac{\sigma^{2}_{\eta}}{\phi_{h}^{2}+1}$

The full conditional distribution of $\lambda_{t}$ is given by:

$\scriptsize
p(\bm{\lambda}|.)=\prod_{t=1}^{T}p(\lambda_{t}|\lambda_{t-1},\lambda_{t+1},.)$
(61) $\scriptsize
p(\lambda_{t}|.)\propto\omega_{t}^{-1}(1-\delta_{t}^{2})^{-0,5}exp\left[-\frac{1}{2}\left(\frac{y_{t}-\bm{x_{t}}\pi-\sqrt{h_{t}}\zeta_{t}-\sqrt{h_{t}}o_{t}^{-0.5}\omega_{t}\delta_{t}v_{t}}{\sqrt{h_{t}}o_{t}^{-0.5}\omega_{t}(1-\delta_{t}^{2})^{0.5}}\right)^{2}-\frac{1}{2}\frac{(\lambda_{t}-\mu_{\lambda_{t}})^{2}}{\sigma^{2}_{\lambda_{t}}}\right]$
(62) $\displaystyle{l}\scriptsize$
$\displaystyle\mu_{\lambda_{t}}=\frac{\phi_{\lambda}}{\phi_{\lambda}^{2}+1}\left(\lambda_{t+1}+\lambda_{t-1}\right)$
(63)
$\displaystyle\sigma_{\lambda}^{2}=\frac{\sigma^{2}_{\xi}}{\phi_{\lambda}^{2}+1}$

The full conditional distribution of $\bm{v}={v_{1},\ldots,v_{T}}$ is given
by:

$\scriptsize p(\bm{v}|.)=\prod_{t=1}^{T}p(v_{t}|.)$ (64)

$\scriptsize p(v_{t}|.)\propto
exp\left[-\frac{1}{2}\left(\frac{1}{(1-\delta_{t}^{2})}v_{t}^{2}-\frac{2o_{t}^{0.5}\delta_{t}h_{t}^{-0.5}}{\omega_{t}(1-\delta_{t}^{2})}(y_{t}-\bm{x_{t}}\bm{\pi}-\zeta_{t}\sqrt{h_{t}})v_{t}\right)\right]\mathbb{I}(0\leq
v_{t}<\infty)$ (65)

this is is a Truncated Normal
$\left(\frac{\delta_{t}o_{t}^{0,5}h_{t}^{-0.5}[y_{t}-\bm{x_{t}}\bm{\pi}]-\delta_{t}o_{t}^{0,5}\zeta_{t}}{\omega_{t}},1-\delta_{t}^{2}\right)_{[0,\infty)}$

####  A.2.3 VAR with Skew Normal shocks

The full conditional distribution of $vec(\bm{\Pi})$ is
$\mathcal{N}(vec(\bm{\bar{\mu}_{\Pi}}),\bm{\bar{V}_{\Pi}})$, where:

$\scriptsize\bm{\bar{\mu}_{\Pi}}=\bm{\bar{V}_{\Pi}}\left[vec\left(\sum_{t=1}^{T}\bm{X_{t}\tilde{y}_{t}^{\prime}\Sigma_{t}^{-1}}\right)+\bm{\uline{V_{\Pi}}}^{-1}vec(\bm{\uline{\mu_{\Pi}}})\right]$
(66)

with
$\bm{\tilde{y}}_{t}\equiv\bm{y}_{t}-\bm{H_{t}^{0.5}}\bm{A^{-1}}\bm{\zeta}_{t}-\bm{H_{t}^{0.5}}\bm{A^{-1}}\bm{\Omega_{t}}\bm{\Delta_{t}}\bm{v_{t}}$
and:

$\scriptsize\bm{\bar{V}_{\Pi}}=\bm{\uline{V_{\Pi}}}^{-1}+\sum_{t=1}^{T}(\bm{\Sigma_{t}^{-1}}\otimes\bm{X_{t}X_{t}^{\prime}})$
(67)

where
$\bm{\Sigma}_{t}\equiv\bm{A}^{-1}\bm{H}_{t}\bm{\omega_{t}}^{2}(\bm{I}-\bm{\Delta_{t}}^{2})\bm{A^{\prime}}^{-1}$
while $\bm{\uline{\mu_{\Pi}}}$ and $\bm{\uline{V_{\Pi}}}$ are the prior mean
and variance covariance matrix.

The full conditional distribution of the elements in $\bm{A}$ is derived
adapting our framework to the approach of [15]. Considering the system:

$\bm{Au_{t}}=\bm{H_{t}^{0.5}\varepsilon_{t}}$ (68)

since
$\varepsilon_{it}=\zeta_{it}+\omega_{it}\delta_{it}v_{it}+\omega_{it}\sqrt{1-\delta_{it}^{2}}z_{it}$
we have :

$\scriptsize\begin{array}[]{l}u_{1t}=\sqrt{h}_{1t}(\zeta_{1t}+\omega_{1t}\delta_{1t}v_{1t}+\omega_{1t}\sqrt{1-\delta_{1t}^{2}}z_{1t})\\\
u_{2t}=-a_{21}u_{1t}+\sqrt{h}_{2t}(\zeta_{2t}+\omega_{2t}\delta_{2t}v_{2t}+\omega_{2t}\sqrt{1-\delta_{2t}^{2}}z_{2t})\\\
u_{3t}=-a_{31}u_{1t}-a_{32}u_{2t}+\sqrt{h}_{3t}(\zeta_{3t}+\omega_{3t}\delta_{3t}v_{3t}+\omega_{3t}\sqrt{1-\delta_{3t}^{2}}z_{3t})\\\
\hskip 7.11317pt\vdots\hskip 56.9055pt\vdots\hskip 56.9055pt\vdots\\\
u_{Nt}=-a_{N1}u_{1t}-a_{N2}u_{2t}\ldots-
a_{N,N-1}u_{2t}+\sqrt{h}_{Nt}(\zeta_{Nt}+\omega_{Nt}\delta_{Nt}v_{Nt}+\omega_{Nt}\sqrt{1-\delta_{Nt}^{2}}z_{Nt})\end{array}$
(69)

therefore :

$\scriptsize\begin{array}[]{l}u_{1t}-\sqrt{h}_{1t}(\zeta_{1t}+\omega_{1t}\delta_{1t}v_{1t})=\sqrt{h}_{1t}\omega_{1t}\sqrt{1-\delta_{1t}^{2}}z_{1t}\\\
u_{2t}-\sqrt{h}_{2t}(\zeta_{2t}+\omega_{2t}\delta_{2t}v_{2t})=-a_{21}u_{1t}+\sqrt{h}_{2t}\omega_{2t}\sqrt{1-\delta_{2t}^{2}}z_{2t}\\\
u_{3t}-\sqrt{h}_{3t}(\zeta_{3t}+\omega_{3t}\delta_{3t}v_{3t})=-a_{31}u_{1t}-a_{32}u_{2t}+\sqrt{h}_{3t}\omega_{3t}\sqrt{1-\delta_{3t}^{2}}z_{3t}\\\
\hskip 7.11317pt\vdots\hskip 56.9055pt\vdots\hskip 56.9055pt\vdots\\\
u_{Nt}-\sqrt{h}_{Nt}(\zeta_{Nt}+\omega_{Nt}\delta_{Nt}v_{Nt})=-a_{N1}u_{1t}-a_{N2}u_{2t}\ldots-
a_{N,N-1}u_{2t}+\sqrt{h}_{Nt}\omega_{Nt}\sqrt{1-\delta_{Nt}^{2}}z_{Nt}\end{array}$

Since I condition on the parameters, the mixing variables and the latent
states I can define
$\tilde{u}_{it}=u_{it}-\sqrt{h}_{it}(\zeta_{it}+\omega_{it}\delta_{it}v_{it})$
for $i=1,\ldots,N$ and
$\tilde{\sigma}^{2}_{it}=\sqrt{h_{it}}\omega_{it}\sqrt{1-\delta_{it}^{2}}$ and
derive the full conditional posterior for the elements of $\bm{A}$ by
exploiting the system of equations:

$\scriptsize\begin{array}[]{l}\tilde{u}_{1t}=\tilde{\sigma}^{2}_{1t}z_{1t}\\\
\tilde{u}_{2t}=-a_{21}u_{1t}+\tilde{\sigma}^{2}_{2t}z_{2t}\\\
\tilde{u}_{3t}=-a_{31}u_{1t}-a_{32}u_{2t}+\tilde{\sigma}^{2}_{3t}z_{3t}\\\
\hskip 7.11317pt\vdots\hskip 56.9055pt\vdots\hskip 56.9055pt\vdots\\\
\tilde{u}_{Nt}=-a_{N1}u_{1t}-a_{N2}u_{2t}\ldots-
a_{N,N-1}u_{2t}+\tilde{\sigma}^{2}_{Nt}z_{Nt}\end{array}$ (70)

where $z_{it}\sim\mathcal{N}(0,1)$. Assuming a Normal prior for the elements
in $\bm{A}$ and defining $\bm{a}_{i}$ the vector that collects the free
elements in the $i^{th}$ row of the $\bm{A}$ matrix, I can use standard linear
regression results to show that the full conditional posterior of $\bm{a_{i}}$
is given by
$\bm{a_{i}}\sim\mathcal{N}(\bm{\bar{\mu}_{a,i}},\bm{\bar{V}_{a,i}})$ where:

$\begin{array}[]{l}\bm{\bar{\mu}_{a,i}}=\bm{\bar{V}_{a,i}}(\bm{\uline{V_{a,i}}}^{-1}\bm{\uline{\mu_{a}}}+\sum_{t=1}^{T}\tilde{\sigma_{it}}^{2^{-}1}\bm{u_{it}}^{\prime}\tilde{u}_{it})\\\
\bm{\bar{V}_{a,i}}=(\bm{\uline{V_{a,i}}}^{-1}+\sum_{t=1}^{T}\tilde{\sigma_{it}}^{2^{-1}}\bm{u_{it}}^{\prime}\bm{u_{it}})^{-1}\end{array}$
(71)

where $\bm{u_{it}}$ is the vector colleting the right hand variables of the
$i^{th}$ equation in the system above (70) with $i=2,\ldots,N$ and
$\bm{\uline{\mu_{a}}}$ and $\bm{\uline{V_{a,i}}}$ are the prior mean and
variance covariance matrix of the free elements of the $i^{th}$ row of
$\bm{A}$.

####  A.2.4 VAR with Skew-t shocks

The full conditional distribution of $vec(\bm{\Pi})$ is
$\mathcal{N}(vec(\bm{\bar{\mu}_{\Pi}}),\bm{\bar{V}_{\Pi}})$, where:

$\scriptsize\bm{\bar{\mu}_{\Pi}}=\bm{\bar{V}_{\Pi}}\left[vec\left(\sum_{t=1}^{T}\bm{X_{t}\tilde{y}_{t}^{\prime}\Sigma_{t}^{-1}}\right)+\bm{\uline{V_{\Pi}}}^{-1}vec(\bm{\uline{\mu_{\Pi}}})\right]$
(72)

with
$\bm{\tilde{y}}_{t}\equiv\bm{y}_{t}-\bm{H_{t}^{0.5}}\bm{A^{-1}}\bm{\zeta}_{t}-\bm{H_{t}^{0.5}}\bm{A^{-1}}\bm{\Omega_{t}}\bm{\Delta_{t}}\bm{O_{t}}^{-0.5}\bm{v_{t}}$
and:

$\scriptsize\bm{\bar{V}_{\Pi}}=\bm{\uline{V_{\Pi}}}^{-1}+\sum_{t=1}^{T}(\bm{\Sigma_{t}^{-1}}\otimes\bm{X_{t}X_{t}^{\prime}})$
(73)

where
$\bm{\Sigma_{t}}\equiv\bm{A^{-1}}\bm{H}_{t}\bm{\omega_{t}^{2}}(\bm{I}-\bm{\Delta_{t}}^{2})\bm{O_{t}}^{-1}\bm{A^{\prime-1}}$
while $\bm{\uline{\mu_{\Pi}}}$ and $\bm{\uline{V_{\Pi}}}$ are the prior mean
and variance covariance matrix.

The full conditional for $\bm{A}$ is derived following the same steps in the
VAR with Skew-normal shocks just by considering that (69) becomes:

$\scriptsize\begin{array}[]{l}u_{1t}=\sqrt{h}_{1t}(\zeta_{1t}+\omega_{1t}\delta_{1t}o^{-0.5}_{1t}v_{1t}+\omega_{1t}\sqrt{1-\delta_{1t}^{2}}o^{-0.5}_{1t}z_{1t})\\\
u_{2t}=-a_{21}u_{1t}+\sqrt{h}_{2t}(\zeta_{2t}+\omega_{2t}\delta_{2t}o^{-0.5}_{2t}v_{2t}+\omega_{2t}\sqrt{1-\delta_{2t}^{2}o^{-0.5}_{2t}}z_{2t})\\\
u_{3t}=-a_{31}u_{1t}-a_{32}u_{2t}+\sqrt{h}_{3t}(\zeta_{3t}+\omega_{3t}\delta_{3t}o^{-0.5}_{3t}v_{3t}+\omega_{3t}\sqrt{1-\delta_{3t}^{2}}o^{-0.5}_{3t}z_{3t})\\\
\hskip 7.11317pt\vdots\hskip 56.9055pt\vdots\hskip 56.9055pt\vdots\\\
u_{Nt}=-a_{N1}u_{1t}-a_{N2}u_{2t}\ldots-
a_{N,N-1}u_{2t}+\sqrt{h}_{Nt}(\zeta_{Nt}+\omega_{Nt}\delta_{Nt}o^{-0.5}_{Nt}v_{Nt}+\omega_{Nt}\sqrt{1-\delta_{Nt}^{2}}o^{-0.5}_{Nt}z_{Nt})\end{array}$
(74)

###  A.3 Metropolis Hastings Step to draw the mixing variable o

In the time varying skewness stochastic volatility models with Skew-t shocks,
the full conditional distribution of $\\{o_{t}\\}_{t=1}^{T}$ is given by:

$\scriptsize p(o_{t}|.)\propto
o_{t}^{\frac{\nu+1}{2}-1}exp\left[-\frac{o_{t}}{2}\left(\nu_{t}+\frac{h_{t}^{-1}(y_{t}-\bm{x_{t}}\bm{\pi}-\sqrt{h_{t}}\zeta_{t})^{2}}{\omega_{t}^{2}(1-\delta_{t}^{2})}\right)\right]exp\left[\frac{(y_{t}-\bm{x_{t}}\bm{\pi}-\sqrt{h_{t}}\zeta_{t})(h_{t}^{-0.5}o_{t}^{0.5}\delta_{t}v_{t})}{\omega_{t}(1-\delta_{t}^{2})}\right]$
(75)

Since it is not possible to directly sample from this full conditional
distribution, I use Metropolis Hastingss step to draw from this conditional
distribution. I use as proposal

$Gamma\left(\frac{\nu+1}{2},\frac{1}{2}\left[\nu+\frac{h_{t}^{-1}(y_{t}-\bm{x_{t}}\bm{\pi}-\sqrt{h_{t}}\zeta)^{2}}{\omega_{t}^{2}(1-\delta_{t}^{2})}\right]\right)$.
The acceptance probability in the Metropolis Hastings step is:

$p=exp\left[\frac{(y_{t}-\bm{x_{t}\pi}-\sqrt{h_{t}}\zeta_{t})h_{t}^{-0.5}o_{t}^{*^{0.5}}\delta_{t}v_{t}}{\omega_{t}(1-\delta_{t}^{2})}-\frac{(y_{t}-\bm{x_{t}\pi}-\sqrt{h_{t}}\zeta_{t})h_{t}^{-0.5}o^{m_{t}^{0.5}}\delta_{t}v_{t}}{\omega_{t}(1-\delta_{t}^{2})}\right]\\\
$ (76)

where $o^{*}$ is a new draw from the proposal and $o^{m}$ is the previous
draw.

###  A.4 Particle Step in the Gibbs Sampler

Table 4 presents the details on the Particle Step used in the Gibbs Sampler to
draw the volatilities and the shape parameters. $s_{t}$ stands for the generic
unobserved latent state being $log(h_{t})$ in Step 5) and $\lambda_{t}$ in
Step 7) of the Gibbs Sampler in Table 1. As anticipated above a valid particle
approximation to the Gibbs Sampler requires a Conditional Sequential Monte
Carlo update which guarantees that a pre-specified path of the state variables
is ensured to survive all the resampling steps [3]. Hence, if I consider a
generic iteration $m+1$ of the Gibbs Sampler, when using $K$ particles to
approximate $p(h_{1},\ldots,h_{T}|\bm{\Theta},\bm{v},\bm{\lambda})$ and
$p(\lambda_{i1},\ldots,\lambda_{iT}|\bm{\Theta},\bm{v},\bm{h})$ , only $K-1$
particles are generated while the $K^{th}$ particle is set to the pre-
specified path $h_{1:T}^{(m)}$ and $\lambda_{1:T}^{(m)}$. In the particle
approximation, I use the transition equations (2) and (3) as importance
densities $g_{\theta}(s_{t})$ and compute the weights accordingly. I refer to
the original paper, [30] for the details on the Ancestral Sampling step, that
for $t>2$ artificially assign a history to the partial pre-specified path
$s_{t:T}^{(m)}$.

Table 4: Particle Step in the Gibbs Sampler Particle Step with Ancestor Sampling |
---|---
Draw $s_{1}^{k}\sim g_{\theta}(s_{1})$ | for $k=1,\ldots,K-1$
Set $s_{1}^{K}=s_{1}^{(m)}$ |
Compute $w^{k}_{1}=W_{1}(s_{1}^{k})$ and normalize the weights | for $k=1,\ldots,K$
for $t=2:T$ |
Re-sampling step: sample $\\{s_{t-1}^{k}\\}_{k=1}^{K}$ with probabilities given by $\\{w_{t-1}^{k}\\}_{k=1}^{K}$ |
Draw $s_{t}^{k}\sim g_{\theta}(s_{t})$ | for $k=1,\ldots,K-1$
Set $s_{t}^{K}=s_{t}^{(m)}$ |
Ancestral sampling step |
Compute $w^{k}_{t}=W_{t}(s_{t}^{k})$ and normalize the weights | for $k=1,\ldots,K$
end |
Draw $j$ with $Pr(j=k)\propto$ $w_{T}^{k}$ |

In alternative to the particle step, I also considered another version of the
algorithm with independence Metropolis Hastings steps to draw the log-
volatilities and the shape parameters. In particular, I considered a log-
normal proposal density for the volatility (on the lines of [15]) as:

$q(h_{t})\propto
h_{t}^{-1}exp\left[-\frac{(lnh_{t}-\mu_{h_{t}})^{2}}{2\sigma_{h}^{2}}\right]$
(77)

with $\mu_{h_{t}}$ and $\sigma_{h}$ defined in equations (49) for the Skew-
Normal case and (60) for the Skew-t. The acceptance probabilities in the model
with Skew-Normal and Skew-t shocks are respectively given by:

$p=\frac{h_{t}^{*^{-0,5}}exp\left[-\frac{(y_{t}-\bm{x_{t}\pi}-\sqrt{h^{*}_{t}}\zeta_{t}-\sqrt{h^{*}_{t}}\omega_{t}\delta_{t}v_{t})^{2}}{2h_{t}^{*}\omega_{t}^{2}(1-\delta_{t}^{2})}\right]}{h^{m_{t}^{-0,5}}exp\left[-\frac{(y_{t}-\bm{x_{t}\pi}-\sqrt{h_{t}^{m}}\zeta_{t}-\sqrt{h_{t}^{m}}\omega_{t}\delta_{t}v_{t})^{2}}{2h_{t}^{m}\omega_{t}^{2}(1-\delta_{t}^{2})}\right]}$
(78)

$p=\frac{h_{t}^{*^{-0,5}}exp\left[-\frac{(y_{t}-\bm{x_{t}\pi}-\sqrt{h^{*}_{t}}\zeta_{t}-\sqrt{h_{t}^{*}}o_{t}^{-0.5}\omega_{t}\delta_{t}v_{t})^{2}}{2h^{*}_{t}o_{t}^{-1}\omega_{t}^{2}(1-\delta_{t}^{2})}\right]}{h^{m_{t}^{-0,5}}exp\left[-\frac{(y_{t}-\bm{x_{t}\pi}-\sqrt{h_{t}^{m}}\zeta_{t}-\sqrt{h_{t}^{m}}o_{t}^{-0.5}\omega_{t}\delta_{t}v_{t})^{2}}{2h_{t}^{m}o_{t}^{-1}\omega_{t}^{2}(1-\delta_{t}^{2})}\right]}$
(79)

where $h_{t}^{*}$ is the new draw from the proposal distribution, while
$h_{t}^{m}$ is the previous draw. Instead, for the shape parameters, I
considered a Normal proposal:

$q(\lambda_{t})\sim N(\mu_{\lambda_{t}},\sigma^{2}_{\lambda})$ (80)

with $\mu_{h_{t}}$ and $\sigma_{h}$ defined in equations (54) for the Skew-
Normal case and (63) for the Skew-t. The acceptance probabilities in the model
with Skew-Normal and Skew-t shocks are respectively given by:

$p=\frac{\omega_{t}^{*^{-1}}(1-\delta_{t}^{*^{2}})^{-0,5}exp\left[-\frac{(y_{t}-\bm{x_{t}\pi}-\sqrt{h_{t}}\zeta_{t}^{*}-\sqrt{h_{t}}\omega^{*}_{t}\delta^{*}_{t}v_{t})^{2}}{2h_{t}\omega_{t}^{*^{2}}(1-\delta_{t}^{*^{2}})}\right]}{(\omega_{t}^{m})^{-1}(1-\delta_{t}^{m^{2}})^{-0,5}exp\left[-\frac{(y_{t}-\bm{x_{t}\pi}-\sqrt{h_{t}}\zeta^{m}_{t}-\sqrt{h_{t}}\omega^{m}_{t}\delta^{m}_{t}v_{t})^{2}}{2h_{t}\omega_{t}^{m^{2}}(1-\delta_{t}^{m^{2}})}\right]}$
(81)

$p=\frac{\omega_{t}^{*^{-1}}(1-\delta_{t}^{*^{2}})^{-0,5}exp\left[-\frac{(y_{t}-\bm{x_{t}\pi}-\sqrt{h_{t}}\zeta_{t}^{*}-\sqrt{h_{t}}o_{t}^{-0.5}\omega_{t}^{*}\delta^{*}_{t}v_{t})^{2}}{2h_{t}o_{t}^{-1}\omega_{t}^{*^{2}}(1-\delta_{t}^{*^{2}})}\right]}{(\omega_{t}^{m})^{-1}(1-\delta_{t}^{m^{2}})^{-0,5}exp\left[-\frac{(y_{t}-\bm{x_{t}\pi}-\sqrt{h_{t}}\zeta^{m}_{t}-\sqrt{h_{t}}o_{t}^{-0.5}\omega^{m}_{t}\delta^{m}_{t}v_{t})^{2}}{2h_{t}o_{t}^{-1}\omega_{t}^{m^{2}}(1-\delta_{t}^{m^{2}})}\right]}$
(82)

where $\omega_{t}^{*},\zeta_{t}^{*},\delta_{t}^{*}$ are functions of the new
draw from the proposal $\lambda_{t}^{*}$, while
$\omega_{t}^{m},\zeta_{t}^{m},\delta_{t}^{m}$ are $\lambda_{t}^{m}$ are
functions of the previous draw $\lambda_{t}^{m}$.

##  Appendix B Appendix

###  B.1 Priors and hyper-parameters

Table 5 and Table 6 report the specification of the priors and the choice of
the hyper-parameters used for the estimation of the models in the empirical
application.

Table 5: Priors for the parameters of the TVSSV model Parameter | Prior
---|---
$\sigma^{2}_{\xi}$ | InverseGamma$\left(5,0.16\right)$
$\sigma^{2}_{\eta}$ | InverseGamma$\left(5,0.16\right)$
$\phi_{h,\lambda}$ | $\mathcal{N}(1,0.01)$
$\beta_{1}$ | $\mathcal{N}(0,10)$
$\pi_{i}$ | $\mathcal{N}(\underline{\mu_{\pi}},\underline{\sigma_{\pi,i}})$
$log(h_{0})$ | $\mathcal{N}\left(\hat{h}_{0},100\right)$
$\lambda_{0}$ | $\mathcal{N}\left(0,10\right)$

$\hat{h}_{i,0}$ is the estimated variance from an AR(4) model to each series
using an initial sample of 40 observations. In the application in Section 3 I
assume that the elements of $\bm{\pi}$ namely $\pi_{i}$ are centered in zero,
namely $\mu_{\pi,i}=0$ and the variances $\sigma_{\pi,i}$ are set following
[11]. For the VAR I consider the following priors:

Table 6: Priors for the parameters of the VAR TVSSV model Parameter | Prior
---|---
$vec(\bm{\Pi})$ | $\mathcal{N}(vec(\bm{\uline{\mu_{\Pi}}}),\bm{\uline{V_{\Pi}}})$
$a_{ij}$ | $\mathcal{N}(0,100)$

where the elements of $vec(\bm{\uline{\mu_{\Pi}}})$ are equal to zero for the
coefficients on the cross-equation lags and for the intercept. The
coefficients of the own lags are centered in 0 for stationary variables and on
1 for non-stationary variables.

$\uline{V_{\Pi}}$ has the Minnesota type prior:

$v_{ij,l}=\begin{cases*}&$\frac{\theta_{1}}{l^{\theta_{4}}}$ \hskip 71.13188pt
if $\quad i=j$\\\
&$\frac{\sigma^{2}_{i}\theta_{1}\theta_{2}}{\sigma^{2}_{j}l^{\theta_{4}}}$
\hskip 56.9055pt if $\quad i\neq j$\\\ \end{cases*}$ (83)

where I set
$\theta_{1}=0.04\quad\theta_{2}=0.025\quad\theta_{3}=100\quad\theta_{4}=2$. We
estimate $\sigma^{2}_{i}$ from univariate AR(12) regressions.

###  B.2 Variables in the medium scale VAR

Table 7: Variable transformations Variable | Transformation
---|---
Real personal consumption expenditures | $log$
Industrial Production | $log$
Unemployment Rate | $level$
Avg Weekly Hours Worked | $log$
Consumer Price Index | $log$
Fed Funds Rate | $level$
10-Year Treasury Yield - Fed Funds Rate | $level$
Moody’s Baa Corporate Bond Yield - the Fed Funds Rate | $level$
Standard and Poors index | $log$

###  B.3 Other Figures

Figure 6: CRPS and Tail Weighted CRPS (left tail)

Note: The figure above shows the time series of the Cumulative Ranked
Probability Scores (CRPS), while the figure below shows the time series of the
Left Tail Weighted CRPS [21]. In blue estimates from the two step quantile
regression based method by [2], in red from the TVSSV univariate model with
Skew-t shocks and in yellow the estimates from the TVSSV VAR model.

Figure 7: Probability Integral Transforms

Note: Probability Integral Transforms of the forecasts from the quantile
regression based method, the univariate time varying skewness stochastic
volatility models with Skew-Normal and Skew-t shocks, and the VARs with
varying skewness and stochastic volatility with Skew-Normal and Skew-t shocks