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.
59,558 characters · 16 sections · 74 citation commands
Vector autoregression models with skewness and heavy tails
Since the seminal work of Sims1980, the vector autoregression (VAR) model has become one of the key macroeconomic models for policy makers and forecasters, see Karlsson2013. The utility of the basic VAR model of Sims1980 has been greatly enhanced by extensions allowing for time varying parameters Primiceri2005,Cogley2005 and stochastic volatility Uhlig1997,Clark2011, Clark2015. These can, however, not fully account for features in the data such as heavy tailed or skewed distributions.
Acemoglu2017 gives a theoretical motivation for the non-Gaussian distribution of macroeconomic variables and the presence of heavy tails and asymmetries is well documented in the literature. For example, Christiano2007 found evidence against Gaussianity by inspecting the skewness and kurtosis properties of residuals from a VAR model and Fagiolo2008 find that the distribution of the output growth rates of OECD countries can be approximated by symmetric exponential-power densities with Laplace tails even after accounting for outliers, autocorrelation and heteroscedasticity. To model the heavy tails Ni2005 propose a VAR model with a multivariate Student's $t$ distribution, while Curdia2014 and Chib2014 impose a similar heavy tailed structural shock in Dynamic Stochastic General Equilibrium (DSGE) models. Karlsson2020, on the other hand, propose a general class of multivariate heavy tailed distributions which includes the normal, $t$ and Laplace distributions as well as their mixture for the error term in the VAR model. Stochastic volatility can also lead to a heavy tailed marginal distribution and Cross2016, Chiu2017 and Liu2019 show that ignoring the stochastic volatility of the shocks will overestimate the fatness of the tails.
As noted by, among others, Curdia2014 the largest shocks occur during recessions, and the skewness of the distribution should be taken into account. Skew-normal and skew-$t$ distributions are common choices for modelling data with skewed distributions. An early application in the VAR literature is Panagiotelis2008 who proposed the use of a multivariate skew-$t$ distribution. Recently Liu2019 estimates different asymmetric and heavy tailed distributions for macroeconomic variables, even though the symmetric Student's $t$ distribution is preferred for monthly data. Delle2020 model the conditional distribution of GDP using a skew-$t$ distribution with time-varying location, scale and shape parameters. Carriero2020 apply a VAR model with a common factor volatility in conditional mean and find evidence of skewness in the unemployment rate and the financial indicator. Carriero2021 account for extreme Covid-19 observations using a VAR model with outlier-augmented stochastic volatility. They show that the model performs on par with a VAR with Student's $t$ distribution.
In this paper, we contribute to the literature by extending the VAR model to account for more realistic assumptions on the multivariate distribution of the variables. We propose a general class of skewed distributions with heavy tails and stochastic volatility for the error terms in the VAR. In doing so we take the generalized hyperbolic skew Student $t$ distribution as our starting point and we refer to this as the GHSkew-$t$-SV class of VAR models. The GHSkew-$t$-SV distribution can be represented as a normal variance-mean mixture and lends itself to straightforward Bayesian inference using a Gibbs sampler with a few Metropolis-Hastings steps. Model comparison and marginal likelihood calculations can be done using the cross entropy method of Chan2018.
In an application to monthly US macro data we compare the in-sample and out of sample forecast performance of 14 VAR models with different assumptions on the tail distribution and stochastic volatility. We find strong support for the VAR models with skewness and heavy tails. Stochastic volatility, heavy tails and skewness all contribute to the in-sample fit. In general, the VAR model with stochastic volatility improves the point and density out-of-sample forecasts. Furthermore, allowing for heavy tailed distributions enhances the out-of-sample forecast which is in agreement with current findings in the literature, see Chiu2017 and Liu2019. An interesting finding is that ignoring the stochastic volatility of the error terms not only overestimates the fatness of the tail distribution but also underestimates the skewness. We recommend that skewness as well as heavy tails should be taken into account for better predictions and in-sample fit.
The rest of the paper is organized as follows. Section 2 introduces the GHSkew-$t$-SV models. Section 3 presents the Bayesian algorithm for inference and the cross entropy methods to calculate the marginal likelihood. Section 4 illustrates some empirical results with monthly macroeconomic data. Section 5 provides some international evidence with the GHSkew-$t$-SV VAR model. Finally, conclusions are reached in Section 6.
We consider different specifications of the skewness and heavy tails in the distribution of the error terms in the VAR model. We also allow for stochastic volatility in addition to the skewness and heavy tails in the distribution of the error terms. Starting with the Gaussian VAR model, we extend the Gaussian structural shock to the generalized hyperbolic skew Student's $t$ distribution using the fact that generalized hyperbolic skew Student's $t$ distribution can be obtained as a variance-mean mixture of the Gaussian distribution.
The Gaussian VAR model with stochastic volatility (Gaussian-SV) is given by
where $\mathbf y_{t}$ is a $k$-dimensional vector of endogenous variables; $\mathbf c$ is a $k$-dimensional vector of constants; $\mathbf B_j$ is a $k \times k$ variate matrix of regression coefficients with $j = 1, \ldots, p$; $\mathbf A$ is a $k \times k$ lower triangular matrix with ones on the diagonal that describes the contemporaneous interaction of the endogenous variables; $\mathbf H_{t}$ is a $k \times k$ diagonal matrix that captures the heteroskedastic volatility; $\boldsymbol {\epsilon}_t$ is a $k$-dimensional vector of error terms that follows a multivariate Gaussian distribution with zero mean vector and identity covariance matrix, i.e. $\boldsymbol {\epsilon}_t \sim \mathcal N_k(\mathbf 0,\mathbf I)$. We assume that the heteroskedastic volatility follows a random walk for $\mathbf H_{t} = diag (h_{1t}, \ldots, h_{kt})$ with
where $\eta_{it} \sim \mathcal N(0,1)$. The Gaussian VAR model without stochastic volatility can be obtained by fixing zero values of $\boldsymbol \sigma^2 = (\sigma_1^2, \ldots, \sigma^2_k)'$ and assuming that $\log h_{it} = \log h_{i0} $ for $i = 1, \ldots, k$ and $t = 1, \ldots, T$.
For notational ease we rewrite the Gaussian VAR model with stochastic volatility in ((ref)) as
where $\mathbf B = (\mathbf c, \mathbf B_1, \ldots, \mathbf B_{p})$ is a $k \times (1+kp)$ variate matrix, $\mathbf x_t = (1, \mathbf y_{t-1}^{'}, \ldots, \mathbf y_{t-p}^{'} )^{'}$ is $(1+kp)$-dimensional vector and $\mathbf u_{t} = \mathbf A^{-1} \mathbf H_t^{1/2} \boldsymbol \epsilon_t $ is a $k$-dimensional vector of heteroskedastic shocks associated with the VAR equations.
The $t$ distribution has been extended in several different ways to allow for skewness and asymmetric behaviour. Among these, Ferreira2007 propose a multivariate skew-$t$ distribution via an affine linear transformation of independent skew-$t$ variables while Sahu2003 use a hidden truncation model to construct a multivariate skew-$t$ distribution where the heavy tail behavior is captured by only one parameter.
We will, however, take the multivariate generalized hyperbolic skew Student's $t$ distribution as our starting point. It is commonly used as it is a general class of distribution which nests the Gaussian distribution and the Student's $t$ distribution as special cases, see Mcneil2015. In addition, Aas2006 showed that the generalized hyperbolic skew Student's $t$ distribution has exponential/polynomial tail behaviors which can handle substantial skewness in comparison to other skew Student's $t$ distributions. As the generalized hyperbolic skew Student's $t$ distribution can be written in term of a Gaussian variance-mean mixture, it is also straightforward to extend the Gibbs sampler for the Gaussian VAR model with stochastic volatility to accommodate heavy tails and skewness.
In the next sections, we extend the distribution of the error term $\mathbf u_{t}$ to more flexible multivariate distributions using two approaches. In the first approach, we rewrite the reduced form VAR into a structural VAR and allow for skewness and heavy tails in the orthogonal shocks. In the second approach, we consider the skewness and heavy tails for each marginal distribution directly by assuming a Gaussian variance-mean mixture for the reduced form errors.
We account for the heavy tailed and asymmetric heteroskedastic shocks in the orthogonal residuals $\mathbf A \mathbf u_{t}$ in the VAR structural form by letting
where $\boldsymbol \gamma = (\gamma_{1}, \ldots, \gamma_{k})^{'}$ is a $k$-dimensional vector of the skewness parameters, the mixing matrix $\mathbf W_{t} = diag(\xi_{1t}, \ldots, \xi_{kt})$ is a $k \times k$ diagonal matrix with $\xi_{it}$ that follows inverse Gamma distribution with shape parameter $\nu_i/2$ and rate parameter $\nu_i/2$, i.e. $\xi_{it} \sim \mathcal {IG} (\frac{\nu_i}{2}, \frac{\nu_i}{2})$, and $\boldsymbol \nu = ( \nu_1, \ldots, \nu_k )^{'}$ is a $k$-dimensional vector that consists of the degrees of freedom; $\mathbf W_t$ and $\boldsymbol \epsilon _t$ are independently distributed. Equation ((ref)) represents the marginal distribution of the orthogonal shock $\mathbf A \mathbf u_{t}$ as a vector of independent univariate generalized hyperbolic skew Student's $t$ distributions (OST). By setting $\boldsymbol\gamma$ to zero we obtain the symmetric and orthogonal $t$ distribution (OT) used by Curdia2014, Clark2015 and Chiu2017. As usual, letting the degree of freedom $\nu_i \rightarrow \infty$ for $i = 1, \ldots, p$, the OT VAR model becomes a Gaussian VAR.
Given the mixing matrix $\mathbf W_{t}$, it holds that $$ \mathbf u_t| \mathbf W_t \sim \mathcal N_k \left( \boldsymbol \mu_t = \mathbf A^{-1} \mathbf W_{t} \boldsymbol\gamma, \boldsymbol \Sigma_t = \mathbf A^{-1} \mathbf W_t^{1/2} \mathbf H_t \mathbf W_t^{1/2}\mathbf A^{-1'} \right). $$ Chiu2017 interprets the mixing matrix $\mathbf W_t$ as capturing the high-frequency shocks in mean and volatility while the stochastic volatility accounts for the low-frequency shocks. The data will determine whether the extreme time variation comes from the volatility shift or from the idiosyncratic heavy tail shocks.
The OST VAR builds the distribution of the error terms from the ground up in terms of the structural form innovations. This makes for a straightforward structural interpretation but also means that the model is sensitive to the identifying assumptions, in this case the triangular structure of $\mathbf A$ and the ordering of the variables. To overcome this we can model the reduced form errors directly as a correlated vector of univariate skew $t$ distributions. We propose a class of multi skew-$t$ (MST) VAR models by assuming that the residuals $\mathbf u_{t}$ are given by
The MST VAR model imposes the skewness and heavy tails directly on the reduced form error in each VAR equation rather than in the idiosyncratic shock. Hence, the marginal distribution of the endogenous variable $y_{it}$ is generalized hyperbolic skew Student's $t$ for $i = 1, \ldots, k$ and $t = 1, \ldots, T$. On the other hand, the OST VAR model considers the reduced form error as a linear combination of orthogonal skew-$t$ shocks. Restricting the mixing variables to be equal for the different equations, $\xi_{1t} = \ldots = \xi_{kt}$, induces a common tail behaviour and the marginal distribution of $\mathbf u_{t}$ is a multivariate generalized hyperbolic skew Student $t$ (Skew-$t$) distribution Mcneil2015. Additionally, setting $\gamma_{1} = \ldots = \gamma_{k} = 0$ results in a multivariate Student $t$ (Student-$t$) distribution.
The defining feature of the multi skew-$t$ in (ref) and what differentiates it from the usual (skew) multivariate $t$ is the equation specific mixing variables $\xi_{it}$. There is thus no commonality in the tail behaviour of the equations.
Also, if $\nu_i \rightarrow \infty$, the MT model becomes a Gaussian VAR with stochastic volatility in spirit of Cogley2005 and Primiceri2005.
In the next section, we illustrate the Bayesian inference and model selection criteria for different specifications of VAR models with/without stochastic volatility and different assumptions on the distribution of the error terms such as Gaussian, Student-$t$, Skew-$t$, orthogonal Student's $t$ (OT), multi Student's $t$ (MT), orthogonal skew Student's $t$ (OST), multi skew Student's $t$ (MST).
In this section, we discuss the prior distributions for the parameters in the MST-SV VAR models and the corresponding general Gibbs sampler scheme. Inference procedures for other model variations are described when needed. We also show how to compute the marginal likelihood based on the cross entropy method proposed by Chan2018.
Denote the set of the MST-SV VAR model parameters and latent variables by \newline $\boldsymbol \theta = \{ vec(\mathbf B)^{'}, \mathbf a^{'}, \boldsymbol \gamma^{'}, \boldsymbol \nu^{'}, \boldsymbol \sigma^{2'}, \mathbf \xi_{1:K,1:T}^{'}, \mathbf h_{1:K,0:T}^{'}\}^{'}$, where $\mathbf a = (a_{2,1}, a_{3,1}, a_{3,2}, \ldots, a_{k,k-1})'$ is the set of elements of the lower triangular matrix $\mathbf A$ and $\mathbf h_{1:K,0}$ is the vector of initial values for the stochastic volatilities. We employ the Minnesota priors for the prior distributions of $\mathbf B$ with the overall shrinkage $l_1 = 0.2$ and the cross-variable shrinkage $l_2 = 0.5$, see Koop2010, and vague prior distributions for other parameters. In details, the Minnesota-type priors assume a Gaussian prior for $vec(\mathbf B)$, i.e. $vec(\mathbf B) \sim \mathcal N ( \mathbf b_{0}, \mathbf V_{\mathbf b_{0}})$, that shrinks the regression coefficients towards univariate random walks with a tighter prior around zero for longer lags. The prior for $\mathbf a$ is also Gaussian, $\mathbf a \sim \mathcal N_{0.5k(k-1)} (0, 10 \mathbf I)$, which implies a weak assumption of no interaction among endogenous variables. The parameters which account for the heavy tails are endowed with Gamma priors, $\nu_i \sim \mathcal G (2,0.1)$ for $i = 1, \ldots, k$ and the skewness parameters are given a normal prior, $\boldsymbol \gamma \sim \mathcal N_k ( \mathbf 0, \mathbf I)$. That is the prior mean of the degrees of freedom of the $t$-distribution is 20 and the skewness has zero prior mean. Finally, the prior for the variance of the shock to the volatility is $\sigma_i^2 \sim \mathcal{G} ( \frac{1}{2}, \frac{1}{2V_{\sigma}})$ which is equivalent to $\pm \sqrt{\sigma_i^2} \sim \mathcal N (0,V_{\sigma})$ see Kastner2014, this prior is less influential in comparison to the conjugated inverse Gamma prior especially when the true value is small. In all cases of VAR model with and without stochastic volatility $\log h_{i0} \sim \mathcal N \left( \log \hat{\Sigma}_{i,OLS}, 4 \right)$ where $\hat{\Sigma}_{i,OLS}$ is the estimated variance of the AR(p) model using the ordinary least square method, see Clark2015.
Given the latent variables $\mathbf \mathbf \xi_{1:K,1:T}$ and the skewness parameters $\boldsymbol \gamma$, the conditional posterior distributions of the remaining parameters in the MST-SV VAR model are similar to those in the Gaussian-SV VAR model. Hence, the MST model can be estimated using a seven-step Metropolis-within-Gibbs Markov chain Monte Carlo (MCMC) algorithm. Let's $\boldsymbol \Psi$ be a set of conditional parameters except the one that we sample from.
The marginal likelihoods of the GHSkew-$t$-SV VAR models require the high-dimensional integration
Following the adaptive importance sampling approach of Chan2018, we divide the model parameters into two groups with the static parameters $\boldsymbol \theta_1 = \{ vec(\mathbf B)^{'}, \mathbf a^{'}, \boldsymbol \gamma^{'}, \boldsymbol \nu^{'}, \boldsymbol \sigma^{2'}, \mathbf h_{0}^{'}\}^{'}$ and the latent states $\boldsymbol \theta_2 = \{ \mathbf \xi_{1:K,1:T}^{'}, \mathbf h_{1:K,1:T}^{'} \}^{'}$. We first use the cross-entropy methods to a the proposal distribution for $\boldsymbol \theta_1$, $f(\boldsymbol \theta_1)$, from the posterior samples. Then, the integrated likelihood $p(y_{1:T}|\boldsymbol \theta_1)$ is calculated using an inner importance sampling loop based on a sparse matrix representation. Algorithm 1 summarizes the marginal likelihood calculation using the adaptive importance sampling approach of Chan2018.
Algorithm 1. (Marginal likelihood estimation via the cross-entropy method)
The number of samples $N$ is chosen such that the variance of the estimated quantity using important sampling is less than one. The parametric families of $f(\boldsymbol \theta_1)$ are the multivariate Gaussian distribution for $(\mathbf B,\mathbf a,\boldsymbol \gamma_{1:k})$, the independent Gamma distribution for $\nu_{1:k}$ and independent inverse Gamma distribution for $\sigma_{1:k}^2$ and $h_{1:k,0}^2$.
In addition to the in-sample model comparison we also assess the forecasting performance of the different specifications of the error distribution in a recursive out of sample forecasting exercise. We compare the forecast accuracy using the mean square forecast error (MSFE) for the point forecast, the log predictive density (LP), and the continuous rank probability score (CRPS) of the posterior predictive distribution for the density forecast.
Let $T_0$ be the last observation in the first estimation sample and $T_1$ the last observation on variable $i$. The MSFE of variable $i$ at $h$ step ahead, for $h = 1, \ldots, H$, is then obtained as,
where $\bar{y}_{i,t+h|t}$ is the mean of the posterior predictive samples using all data up to time $t$ and $y_{i,t+h}^{o}$ is the observed outcome of variable $i$ at $h$ steps ahead. The model with a smaller MSFE is preferred.
The LP of the posterior predictive distribution is computed as,
where $p(y_{i,t+h}^{o} | \mathbf y_{1:t})$ is the $h$-step ahead posterior predictive density function evaluated at the realization of the variable. Following Andersson2008, the LP of the posterior predictive distribution is computed using the Rao-Blackwellization idea which is more stable than the kernel density estimator for extreme observations. In particular, it is evaluated as,
where $\boldsymbol \theta^{(1)}, \ldots, \boldsymbol \theta^{(R)}$ are the posterior samples of the VAR model. The possibly high dimensional integral over intermediate observations implicit in $p(y_{i,t+h}^{o} | \boldsymbol \theta^{(r)}, \mathbf y_{1:t})$ can be approximated by the Monte Carlo approach. For each sample from the posterior we simulate a new path $\mathbf y_{(t+1):(t+h-1)|t}^{(r)}$ using the data generating process for the model and calculate $p(\mathbf y_{i,t+h|t} | \boldsymbol \theta^{(r)}, \mathbf y_{1:t}, \mathbf y_{(t+1):(t+h-1)|t}^{(r)})$. A higher LP value indicates a better density forecasting performance of the model.
The continuous rank probability score (CRPS) is also commonly used to rank the density forecasts. CRPS is obtained as the quadratic difference between the predictive cumulative distribution function and the empirical distribution of the variable Gneiting2007. As Clark2015 noted the CRPS is less sensitive to outliers than the LP and rewards more for values of the predictive density that are close to the outcome.
where $f$ is the predictive density of the variable $y_{i,t+h|t}$, and $(y_{i,t+h|t}, y_{i,t+h|t}^{'}) $ are independent random draws from the predictive density $f$. We apply the Monte Carlo method to simulate 10,000 draws from the predictive density $f$ and compute the expectation.
We estimate a four-variable VAR with industrial production, inflation rate, unemployment rate, Chicago board options exchange's volatility index (VIX) to illustrate the performance and empirical relevance of the different specifications of the error distribution. We use monthly data for the period 01/1970 to 12/2019 from the Federal Reserve Bank of St. Louis, see McCracken2016. Industrial production is included as a growth rate (first difference of the logarithm of the index), the inflation rate is calculated as the first difference of the log of the CPI and the logarithm of the VIX is used. The variables enter with $p=4$ lags.
We compare 14 different specifications for the error distribution: Gaussian; multivariate Student's $t$ and multivariate skew-$t$; orthogonal $t$ (OT) and orthogonal skew-$t$ (OST); multi $t$ (MT) and multi skew-$t$ (MST). All with and without stochastic volatility.
We first estimate the 14 VAR models with and without stochastic volatility using the in-sample dataset. Then we perform an out-of-sample forecasting exercise to measure the forecast accuracy of each VAR model.
The left-hand side of Figure (ref) shows the growth rate of industrial production, inflation rate, unemployment rate and the VIX. Extreme values of the variables are often observed during recession periods based on the NBER indicators. Industrial production growth decreased by more than 4% during the financial crisis in 2008, while the unemployment rate peaked at 10% and the VIX reached as high as 4.2. These unconditional skewed behaviors can be generated by a time-varying variance shock and/or a skewed shock. The right-hand side of Figure (ref) plots the estimated stochastic volatility in the log scale. The volatilities are occasionally higher during recessions which illustrates the relation of the VIX and other macroeconomic variables. We compare the volatilities of the low frequency shocks obtained from the OST-SV model in the solid lines and that of the Gaussian-SV model in the dashed lines. Using the Gaussian-SV model, the volatility of macroeconomic variables might be overestimated during the recessions and crises which is in agreement with the finding of Curdia2014, Chiu2017, among others.
Table (ref) estimates the log marginal likelihood of the VAR models of Section (ref) with and without stochastic volatility. In the class of VAR models without SV, allowing for heavy tails leads to a substantial improvement in the marginal likelihood while the addition of skewness is less useful. Allowing for stochastic volatility leads to a dramatic improvement in the marginal likelihood for all seven specifications. Allowing for heavy tails improves on the Gaussian-SV and the more flexible OST and MST specifications of skewness perform best with a log Bayes factor of 4.8 (MST) and 3.9 (OST) against the third best Student-$t$ specification. It is interesting that skewness plays a more important role in the VAR model with SV than in the VAR model without SV. The flexibility of the tail behaviour in the OST and MST is, however, important as evidenced by the relatively poor performance of the skew-$t$ VAR models where only one mixing variable is used to model the heavy tails.
Next, we take a closer look at the effect of the SV assumption on the skewness and heavy tail parameters in the VAR models. Figure (ref) focuses on the posterior samples of the skewness parameters and the degree of freedom parameters in the MST VAR model with and without stochastic volatility. The left hand side figures show that ignoring the time-varying volatility of the structural shocks overestimates the fatness of the tails. This finding is inline with Chiu2017 and Liu2019, among others. Moreover, we find that the asymmetric distribution is more important in the VAR models with SV. The magnitudes of the skewness parameters for industrial production growth and inflation rate are higher in the case of SV models while the fatness of the tail distribution are smaller. Appendix (ref) also confirms this finding for the OST-SV VAR model. Hence, even considering that the stochastic volatility has captured the time varying uncertainty, there are still asymmetric structural shocks between recession and moderation periods. The heavy tails appear in both industrial production and VIX, however, VIX shows some degree of skewness in the right tail. Although the conditional distributions of inflation rate and unemployment rate are closer to the Gaussian case with SV, there is a possibility that inflation rate is right-skewed. To summarize, when the stochastic volatility assumption reduces the degree of heavy tailed shock, it also increases the level of skewness shock.
To assess the out-of-sample predictive accuracy of the different specifications, we conduct a recursive forecast exercise using the 01/2017 to 12/2019 period as our evaluation sample. We calculate the MSFE for the point forecasts, and the LP and CRPS for the density forecasts. As the VAR models can be nested based on the assumptions on the tail distributions, they are divided into two model groups without and with stochastic volatility for ease of comparison. Using the Gaussian VAR as a benchmark in each group, we test the forecast accuracy using the one-sided Diebold1995 test where the standard errors of the test statistics are computed with the Newey–West estimator, see the discussion in Clark2011. We also compare the Gaussian model with SV to the Gaussian model without SV.
Table (ref) reports the relative improvements in MSFE over the Gaussian VAR models. Each panel reports the MSFE of each variable using the Gaussian VAR model with (and without) stochastic volatility. The relative improvements over the Gaussian models are computed as the ratio of the MSFE of alternative specifications over the Gaussian models during 2007-2019. Entries less than 1 indicate that the given model is better. In the non stochastic volatility VAR group, skewness and heavy tail models improve the point forecast of Industrial production growth up to 12 months ahead but are only statistically significant up to 6 months ahead. While the prediction of unemployment rate is improved in the long run and that of inflation is enhanced in the one step ahead forecast. In the stochastic volatility VAR group, the advantage of skewness and heavy tail models over Gaussian model diminishes. For unemployment there is an improvement overall when allowing for heavy tails and or skewness, significantly so for the one month ahead forecasts.
Table (ref) reports the relative improvements in LP over the Gaussian VAR models. Here entries greater than 0 indicate that the given model is better than the Gaussian non-SV/SV model. The findings for the density forecasts are consistent with the findings for the point forecasts for the non stochastic volatility group. In the VAR with stochastic volatility group, the density forecast of the heavy tail and skewness models is improved at all horizons for the industrial production growth, however it is only statistically significant for the 6 and 12 months ahead forecast. On the other hand, the short run forecast accuracy of VIX increases when including the heavy tail and skewness parameters. The density forecasts in OST and MST VAR models for unemployment rate and inflation rate remain similar to the Gaussian benchmark. In agreement with Clark2011 and Clark2015, stochastic volatility in the Gaussian VAR model is a crucial characteristic to improve out-of-sample density forecast.
Table (ref) reports the relative improvements in CRPS over the Gaussian VAR models where entries greater than 0 indicate that the given model is better. We confirm the previous conclusion by comparing the CRPS among models. However, the effect of heavy tails and skewness is smaller as the CRPS is less sensitive to outliers Clark2015. Skewness and heavy tailed VAR models with stochastic volatility improve significantly on the Gaussian model with SV in the medium term forecast of industrial production and the short term for the unemployment rate. The assumption of stochastic volatility is still essential for the density forecasts.
Next, we concentrate on the effect of skewness parameters in VAR models with stochastic volatility. Figure (ref) shows the probability integral transforms (PITs) of the three month ahead forecast horizon from OT, MT, OST, MST VAR models. Without significant evidence of skewness (see Figure (ref)), the PIT plots for industrial production, inflation and unemployment are very similar. For the VIX the OT-SV model clearly overestimates the left tail of VIX with only a few observations classified as extreme according to the density forecast. The MT VAR model on the other hand has a slight tendency to overestimate the left tail of the VIX. The OST and MST specifications perform better here and we see that allowing for skewness can make a difference. The result is similar for the PIT plots at the other horizons.
Figure (ref) shows the cumulative log Bayes factors of the predictive density for 3 month forecast horizon between the Gaussian-SV and MST-SV models, see the computational details in Geweke2010. Positive values (red) means that MST-SV predicts better than the Gaussian-SV. A common feature across the variables is that the MST-SV performs better than or roughly on par with the Gaussian-SV up to the middle of the recession and performs worse close to the end. For industrial production the MST-SV performs better overall and regains it advantage after the recession. The Gaussian-SV performs better overall for inflation and unemployment but it is noteworthy that, for unemployment, the MST-SV consistently outperforms the Gaussian-SV by a small margin over the expansion. For the VIX the MST-SV also improves its performance during the expansion and does significantly better overall. This pattern suggests that allowing for heavy tails and skewness is not just about accounting for large deviations but also about modelling the more central parts of the distribution well and that the latter can be equally important. Recalling that Chiu2017 interpreted the mixing variables as accounting for high frequency shocks we can also see this as a factor explaining the overall better performance of the MST-SV during the expansion.
Skewness and heavy tails are empirically relevant features in many application areas -- not only the macroeconomic and financial application we consider in this paper. While these features to some extent can be accommodated or masked by time-varying heteroskedasticity modelled as GARCH-type or stochastic volatility processes there is a need for models that explicitly account for skewness and heavy tails in the data. We contribute to this by proposing flexible skew and heavy tailed distributions with the symmetric normal distribution as a special case. Specifically, we introduce a general class of Generalized Hyperbolic Skew Student's $t$ distributions with stochastic volatility for VAR models. The stochastic representation of the GHSkew-$t$ can be written in term of a variance-mean mixture which leads to a straightforward implementation of a Gibbs sampler for posterior inference. We also take advantage of the cross entropy methods by Chan2018 to calculate the model marginal likelihood and compare the in-sample fit among different specifications. In an application to US data we find support for VAR models with skewness and heavy tails. The VAR models with skewness and heavy tails gives better point forecasts and density forecasts compared to Gaussian VAR models for many, but not all, variables we model. We recommend that skewness should be taken into account for improving forecasting performance during recessions and crises.
We thank Pär Österholm for helpful comments. The authors acknowledge financial support from the project “Models for macro and financial economics after the financial crisis" (Dnr: P18-0201, BV18-0018) funded by the Jan Wallander and Tom Hedelius Foundation. Stepan Mazur also acknowledges financial support from the internal research grants of Örebro University. The computations were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC) at HPC2N partially funded by the Swedish Research Council through grant agreement no. 2018-05973.