Skew selection for factor stochastic volatility models
Jouchi Nakajima









Skew selection for factor stochastic volatility models
    Jouchi Nakajima
[
Bank for International Settlements. E-mail: jouchi.nakajima@bis.org. The views expressed herein are those of the author alone and do not necessarily reflect those of Bank for International Settlements.]


January 2019

    
================================================================================================================================================================================================================================================




This paper proposes factor stochastic volatility models with skew error distributions. The generalized hyperbolic skew t-distribution is employed for common-factor processes and idiosyncratic shocks. Using a Bayesian sparsity modeling strategy for the skewness parameter provides a parsimonious skew structure for possibly high-dimensional stochastic volatility models. Analyses of daily stock returns are provided. Empirical results show that the skewness is important for common-factor processes but less for idiosyncratic shocks. The sparse skew structure improves prediction and portfolio performance.



Key words: Factor stochastic volatility; Generalized hyperbolic skew t-distribution; Portfolio allocation; Skew selection; Stock returns. 













§ INTRODUCTION



Multivariate analyses of stochastic volatility (SV) have been widely discussed in both the theoretical and empirical contexts of financial econometrics literature <cit.>. Among them, the class of factor stochastic volatility (FSV) models is a theoretically and practically effective approach to understanding dynamic variances and covariances in multivariate series of financial variables <cit.>. Summarizing common dynamics among the series by latent factors yields an efficient modeling strategy for multivariate volatility, and for high-dimensional applications in particular. Previous studies have proposed generalizations on the basis of a standard FSV model, for example, with time-varying factor loadings <cit.>, and heavy-tailed error distributions <cit.>.


The skewness for a return distribution of the financial variables has been recognized as an important aspect of model fit and predictive performance <cit.>. The skew error distributions have been incorporated into univariate SV models <cit.> and a Cholesky-type multivariate SV <cit.>. Little is known, however, about the skew error distribution of the FSV models.


To fill this gap, this paper introduces an FSV model with skew error distribution. The generalized hyperbolic (GH) skew t-distribution <cit.> is incorporated into common-factor processes and the idiosyncratic shocks in the FSV model. It is simple and straightforward to assume that each of the factor processes and idiosyncratic shocks is associated with a different skewness parameter. Given this full model, a key question is whether all of the skewness parameters are necessary to describe the skewness in the multivariate series of financial variables. Because the factor processes summarize the common dynamics in multivariate SV series, some of the skewness parameters in the idiosyncratic shocks may not be relevant, and furthermore, not all of the latent processes may need the skew distribution. If so, then, assuming all possible skewness parameters can lead to more uncertainty in model-fit estimation and prediction, especially for high-dimensional applications.


To address this point, the paper proposes a skew selection strategy for the proposed FSV model that includes the skew error distribution. It employs Bayesian sparsity modeling, which has been widely exploited in the literature <cit.>. A so-called spike-and-slab prior selects a zero or non-zero skewness parameter for each factor process and each idiosyncratic shock. This approach allows the proposed FSV model to select-out redundant skewness from the full model and to provide a parsimonious multivariate skew structure. The proposed approach can easily scale up, which has practical relevance for high-dimensional problems. The paper is related to the literature on Bayesian sparsity modeling in the context of multivariate SV for reducing the number of parameters. <cit.> propose a Bayesian stochastic search approach; <cit.> develops a global-local shrinkage prior; and <cit.> and  <cit.> incorporate time-varying sparsity with a latent thresholding process.


This paper provides two empirical analyses using multivariate daily stock returns to illustrate the proposed modeling framework. The first analysis illustrates the key characteristics of the proposed model with a concise dataset of five US-sector indices. A forecasting exercise and a Value-at-Risk (VaR) analysis for portfolio are examined. In the second analysis, a larger-scale dataset of 20 US individual stock returns is analyzed, which highlights the advantages of the skew FSV model and the sparsity skewness modeling in a high-dimensional, practical setting. Portfolio performance in an out-of-sample forecasting exercise is provided. Empirical results from these analyses show that the skewness is important for common-factor processes but less for idiosyncratic shocks. The sparse skew structure contributes to improving stock returns and VaR forecasts as well as portfolio performance.


<cit.> develop the same type of the skew FSV model, incorporating the sparsity prior for the skewness parameter, as well as for the leverage-effect parameter. Their empirical analysis using daily stock returns shows that the FSV model with skewness performs better than a no-skew model for VaR forecasts. The focus of the current paper is different; it concentrates on how a sparse skew structure plays a role in improving the forecasting ability of the skew FSV model. The empirical analyses use several diagnostics of forecasting ability in addition to the VaR measure and reveal that the skew FSV model with the skew selection mechanism performs better than the one without it.


Section 2 formulates the FSV models that includes the skew distribution. Section 3 explains a Bayesian estimation strategy for model fitting. Section 4 provides the first empirical analysis using five US-sector stock indices. Section 5 shows the second analysis using 20 individual stock returns. Section 6 concludes.








§ THE MODEL



This section describes the GH skew t-distribution in a univariate case (Section <ref>), and then defines the FSV models with skew distribution (Section <ref>). The identification issue for the factor loadings is addressed (Section <ref>), and the sparsity structure for the skewness is introduced (Section <ref>).








 §.§ GH skew t-distribution
 


The random variable that follows the GH skew t-distribution, defined by w_it (t=1,2,…) for a series i, can be written in the form of normal variance-mean mixture:
	
    w_it =  m_i + β_i z_it + √(z_it)ε_it,
    	  ε_it∼ N(0,1),

where z_it follows the generalized inverse gaussian (GIG) distribution <cit.>. This original form of the GH skew t-distribution is quite flexible for determining a shape of distribution. Previous studies <cit.> discuss that the parameters in the GH skew-t distribution are typically difficult to jointly estimate. Following <cit.>, to derive a parsimonious structure of the skewed distribution for the SV models, the current paper assumes that (i) z_it∼ IG(ν_i/2,ν_i/2), where IG denotes the inverse gamma distribution; (ii) m_i = -β_ic_i, where c_i ≡E(z_it) = ν_i/(ν_i-2), for E(w_it) = 0; and (iii) ν_i>4 for the finite variance of w_it. The resulting density function describes the shape of distribution with two parameters β_i and ν_i. Although the skewness and kurtosis of the distribution are functions of both these two parameters <cit.>, the β_i and ν_i essentially represent the degree of skewness and kurtosis, respectively.


The literature discusses several alternatives of skew t-distributions <cit.>. Specifications differ among the distributions, which yield different shapes of the density, tail behaviors in particular. This paper uses the GH skew t-distribution defined above to exploit the form of the normal variance-mean mixture, as it makes Bayesian computation for model-fitting easy and efficient as described later in the paper. While it is of interest to compare different specifications of the skew t-distributions in terms of model fit and forecasting, this paper sticks to the GH skew t-distribution because the main purpose of the analysis is to investigate the skew structure of financial time series in the context of FSV models. Running a horse race among the different skew t-distributions is left for future work.








 §.§ Factor stochastic volatility with skew distribution
 


Define a k× 1 vector of response, _t=(y_1t,…,y_kt)', for t=1,…,T. A standard FSV model is formulated as
	
    _t    =   _t + _t,   _t ∼ N(, _t), 
    _t    ∼    N(, _t),

where _t=(f_1t,…,f_pt)' is a p× 1 vector of the latent, common factors;  is a k× p matrix of factor loadings; _t = (exp(h_1t),…,exp(h_kt)), and _t = (exp(h_k+1,t),…, exp(h_qt)), with q=k+p. This form of the model arises from the idea of dimension reduction in the unconditional variance of _t, denoted by _t, as _t∼ N(, _t). With a large k, the matrices in eqns. (<ref>) and (<ref>) yield an effective reduction of dimensionality in _t, with _t=_t'+_t. Because the matrices _t and _t are diagonal and a typical choice of p is small, the number of parameters can be reduced effectively.


We build the proposed skew t-distribution on the simple form of the FSV model defined in eqns. (<ref>) and (<ref>). A variety of generalizations can be easily introduced, however. A local trend can be added to the observation equation (<ref>); the factor loading  can be time-varying; and the process of _t can be dynamic, for example, as it follows a vector autoregressive (VAR) model <cit.>. The following discussion on incorporating the skewness and its estimation method can be straightforwardly applied to the more elaborated FSV models.

The new class of the FSV models with the GH skew t-distribution is defined by
	
    _t    =   _t + _t^1/2_t,   _t = (w_1t,… w_kt)', 
    
    	w_it   =   β_i (z_it - c_i) + √(z_it)ε_it,    i=1,…,k, 
    _t    =   _t^1/2_̌t,   _̌t = (v_1t,… v_pt)', 
    
    	v_jt   =   β_k+j (z_k+j,t - c_k+j) + √(z_k+j,t)ε_k+j,t,    j=1,…,p,

where each element of _t and _̌t follows the GH skew t-distribution specified by eqn. (<ref>). Define _t=(h_1t,…,h_qt)', and _t=(ε_1t,…,ε_qt)'. The latent SV process is defined as the standard autoregressive form:
	
    _t+1 =  + (_t-) + _t,

where = (μ_1,…,μ_q)', = (ϕ_1,…,ϕ_q), _t = (η_1t,…,η_qt), and 
	
    ([ _t; _t ])  ∼  N( , [ [  ;   ]] ),

where = (σ_1^2,…,σ_q^2), and = (ρ_1σ_1,…,ρ_qσ_q). This formulation makes each of the processes {h_it}_t=1^T follow the univariate autoregressive process. The parameter ρ measures the so-called leverage effect, which is widely addressed in financial econometrics <cit.>. In a stock market, for example, a large price drop tends to lead to a large increase in the volatility the next day. This behavior can be described using a negative correlation (ρ<0) between the two shocks (ε_t, η_t).



A key aspect of the proposed model structure is that, conditional on (, _t), the model reduces to q-fold state-space representation of the univariate SV model. To see this, define
	
    ỹ_it = {[ y_it - _(i)_t      (i ≤ k),;          f_jt (i>k; j=i-k), ]. i=1,… , q,

where _(i) denotes a 1× p vector of the i-th row of . Because of the diagonal structure of the matrices (_t, _t, , ), we obtain the univariate SV model with leverage effect, that is,
	
    ỹ_it   =   {β_i (z_it-c_i) + √(z_it)ε_it}exp(h_it/2), 
    
    	h_i,t+1   =   μ_i + ϕ_i (h_it - μ_i) + η_it, 
    ( [ ε_it; η_it;      ])
    	   ∼    N(, _i),
    	  and  _i=
    	( [      1 ρ_iσ_i; ρ_iσ_i  σ_i^2 ]),

for i=1,…,q. In a Bayesian estimation framework, this structure enables us to easily scale up the dimension of the responses for the computation of parameter estimates. With the conditionally independent priors described below, we can obtain an efficient and fast parallel computation.









 §.§ Identification
 

A general modeling issue on the latent factors is restriction for identification. The factor models require a proper parameter restriction on the loading matrix or another valid approach. The most common approach, the one taken here, is the upper triangular form of the factor-loading matrix to order factors according to the selected first p time-series variables <cit.>. Specifically, we assume that b_ii=1, for i=1,…,p; and b_ij=0, for i<j, and i≤ p. Although this approach constrains the model because it depends on the selection and ordering of the first p variables, it is practical to track it with a robustness check that examines different sets of the first p series. Alternative strategies can be found in the existing studies <cit.>.








 §.§ Sparsity prior for the skewness
 

The Bayesian sparsity modeling strategy is incorporated into the skewness parameters to develop the skew selection in the FSV model. 
The sparsity prior for β_i is given by
	
    β_i  ∼ κ N(β_i|0, τ_0^2) + (1-κ)δ_0(β_i),      i=1,…,q,

where δ_0 denotes the Dirac delta function at zero, and 0<κ<1. This prior assigns the normal distribution for β_i with probability κ and the point mass at β_i=0 with probability 1-κ. This is the so-called spike-and-slab prior that is often used in the literature <cit.>. This approach yields an effective skew selection for the FSV model, selecting out redundant skewness parameters from the full set of skewness and providing parsimonious multivariate skew structure.


It is not easy to theoretically derive the skewness of y_it in the resulting FSV model. The following brief simulation study addresses how the skewness of y_it depends on the skewness parameters for the idiosyncratic shocks and the common-factor processes. A sample of size T=1,000, k=3, p=2 is simulated from the skew FSV model. The parameter values are set as ϕ_i=0.995, σ_i=0.05, ρ_i=-0.5, and ν_i=8 for i=1,…,q. These values are selected following previous studies <cit.>. For the level of volatility, we set μ_i=-11, for i=1,2,3, while μ_i=-10 for i=4,5; to assume that the factor volatilities are higher than volatilities of idiosyncratic shocks, which are typically observed in practice. The free parameters in the factor loading matrix  are randomly generated from the uniform distribution U[0.5,1.5]. 








Figure <ref> exhibits summaries of the skewness of y_it from 1,000 sets of simulated data for different values of the skewness parameter: ≡(β_1,…,β_q)= (i) (-1,-1,-1,0,0), (ii) (0,0,0,-1,0), (iii) (0,0,0,-1,-1), and (iv) -_1× 5.  Case (i) assumes that only the idiosyncratic shocks follow the skew distribution; the resulting skewness of _t is not so large as Case (iv) which assumes both the idiosyncratic shocks and the factor processes follow the skew distribution. This reflects the structure of the factor model: y_it is the sum of the factor and the idiosyncratic shock. The level of the factor volatility is assumed to be higher than that of idiosyncratic shock in this simulation exercise, the symmetry of the factor processes weakens the degree of skewness caused from the idiosyncratic shock.


Case (ii) assumes that only the first factor process follows the skew distribution. The degree of skewness of the first series (Y1) is larger than that in Case (i), because the skewness in the first factor process dominates the idiosyncratic shock. In contrast, the degree of skewness in the second (Y2) and the third (Y3) series is slightly smaller than the first series, because the second factor, which follows the symmetric distribution but does not appear in the first series, weakens the degree of skewness. The level of volatility for the second factor, however, is typically smaller than that for the first factor in practice, which leads to a smaller difference between the skewness for the first series and the others. The key finding here is that the resulting skewness of _it also depends on . 


Case (iii) assumes that both two factor processes follow the skew distribution, which exhibits a notably similar degree of skewness as Case (iv). From this result, the assumption of β_i≠ 0 for all i=1,…,q, as in Case (iv), can be redundant to describe the skewness of _t. The skew selection structure with the sparsity prior for  is expected to effectively address this point for the skew FSV model, exploring the best parsimonious structure to select out the redundant parameters and improve the prediction performance.










§ BAYESIAN ESTIMATION



The Bayesian analysis and computation are developed for the proposed FSV model based on the Markov chain Monte Carlo (MCMC) method. The sampling estimation scheme can be built on popular efficient conditional samplers for univariate SV models with a leverage effect <cit.>, for the state space dynamic models <cit.>, and for the multivariate FSV models <cit.>. <cit.> utilize particle Gibbs with an ancestor sampling algorithm to sample the volatility process in the skew FSV model. In contrast, the current paper shows that the full posterior analysis can be implemented by using the simpler algorithm of multi-move sampler for the volatility without any particle sampling.


Based on observations = {_1,…,_T}, the conditional samplers are developed for each of the following latent variables and model parameters:
	
	
  * The q-fold univariate SV processes {h_i1,…,h_iT}, their parameters _i ≡{μ_i,ϕ_i,σ_i,ρ_i, ν_i}, and mixing latent processes {z_i1,…,z_iT} for i=1,…,q.
	
  * The factor process _t and the loading matrix .
	
  * The skewness parameter  and the sparsity parameter κ.
	
A key computational strategy here is decoupling the univariate SV processes in the conditional samplers. Generally, a sampler for the SV process requires a relatively heavier computational burden and makes the mixing of the MCMC algorithm slower. Its simple multivariate extension and, further, with the skewness distribution can be more challenging. However, as mentioned before, the representation of the normal variance-mean mixture leads to the efficient computation with those conditional samplers for the SV processes and their parameters. In addition, the conditional independence of the univariate SV leads to a conditional sampler for β_i under the sparsity prior.

Each of the conditional sampler is outlined as follows.




 Generation of SV and mixing latent process 


Based on eqns. (<ref>)–(<ref>), each univariate SV process is generated using the conditional sampler for the SV with leverage effect. <cit.> develop the efficient multi-move sampler for this type of non-linear state space model and <cit.> utilize it on the SV model with the GH skew t-distribution, which is taken in the current analysis. For the generation of the mixing latent process z_it, the conditional posterior distribution is the product of the kernel of the inverse gamma distribution and another part. A standard Metropolis-Hasting algorithm is employed with a candidate generated from the inverse gamma distribution. 

To develop the conditional sampler for the parameters defining each of the SV process, we simply assume traditional priors that are independent across i: normal priors for μ_i, shifted beta priors for ϕ_i and ρ_i, inverse gamma priors for σ_i^2, and truncated gamma priors for ν_i with ν_i>4. Then, the conditional samplers are implemented with standard posterior distributions for μ_i and the MH algorithm for (σ_i,ρ_i) and ν_i <cit.>.


Note that because these latent processes {h_i1,…,h_iT, z_i1,…,z_iT} and the parameter _i are conditionally independent across i in the conditional posterior inference, we can implement parallel computing for the generation of q-fold SV processes.






 Generation of factor process and loadings 


Define η̂_it = (h_i,t+1 - μ_i) - ϕ_i(h_it - μ_i), and σ̂_it^2 = z_it(1-ρ^2_i)e^h_it, for i=1,…,q. Conditional on the mixing latent process, SV process and its parameters, we obtain
	
    y_it   =   _(i)_t + α_it + û_it,   i=1,…,k, 
    
    	f_jt   =   α_k+j,t + û_k+j,t,    j=1,…,p,

where
	
    α_ℓ t   =   {β_i(z_ℓ t-c_ℓ ) + √(z_ℓ t)ρ_ℓη̂_ℓ t/σ_i}e^h_ℓ t/2, 
    û_ℓ t   ∼    N(0,σ̂_ℓ t^2),

for ℓ=1,…,q.

Define _̱i=(B_i1,…,B_ir)' as the r× 1 vector of free parameters in the i-th row of , with r=min(i-1,p), for i=2,…,k. We generate the sample of _̱i, for i=2,…,k, based on the following regression model:
	
    ŷ_it =  d_it + _it'_̱i + û_it,   û_it∼ N(0,σ̂_it^2),

where ŷ_it = y_it - α_it, _it = (f_1t, …, f_rt)', and
	
    d_it = {[             f_it   (for i=2,…,p),;                0 (for i=p+1,…,k). ].

We set a prior as _̱i ∼ N(_̱i0, _i0). Then, the conditional posterior distribution of _̱i is N(_i, _i), where
	
    _i  = (_i0^-1 + ∑_t=1^T _it_it'/σ̂_it^2)^-1,  _i  = _i (_i0^-1_̱i0 + ∑_t=1^T (ŷ_it - d_it)_it/σ̂_it^2).



For the generation of _t, define _t = (ŷ_1t,…,ŷ_kt)', and _t = (û_1t,…,û_kt)'. We obtain the following form:
	
    _t    =   _t + _t, 
    _t    ∼    N(_t, _t),

where _t = (α_k+1,t,…,α_qt)', and _t = (σ̂_k+1,t^2,…,σ̂_qt^2). We generate _t from the conditional posterior distribution, N(_t, _t), where _t = (_t^-1 + '^-1)^-1, and _t = _t (_t^-1_t + '_t^-1_t), with _t = (σ̂_1t^2,…,σ̂_kt^2), for t=1,…,T.





 Generation of skewness parameters 


Before we derive the posterior distribution of the β_i under the sparsity prior, consider a simple normal prior N(0,τ_0^2) for β_i. Then, the conditional posterior distribution of β_i appears to be the normal distribution; we denote it as N(β̂_i,τ̂_i^2). Turning into the sparsity prior, the conditional posterior distribution is derived as conjugate as
	
    β_i | · ∼ κ̂_i N(β_i | β̂_i,τ̂_i^2) + (1-κ̂_i)δ_0(β_i),

where κ̂_i=κγ_i / (κγ_i+1-κ), and γ_i=exp(β̂_i^2/2τ̂_i^2)τ̂_i/τ_0. In this conditional posterior, β_i follows the normal distribution with probability κ̂, and shrunk at zero otherwise.


For the generation of κ, we set a beta prior and derive the conjugate posterior distribution conditional on the number of β_i's such that β_i ≠ 0.









§ A STUDY OF US-SECTOR INDICES



The first real-data analysis applies the proposed FSV model to a series of k=5 daily stock returns of US-sector indices. This analysis demonstrates how the skew error distribution and sparsity modeling for the skewness parameter play the role of improving the model fit and forecasting ability of the FSV model.







 §.§ Data and setup


The data consist of five selected S&P500 Sector indices, listed in Table <ref>. The sample period spans from January 4, 2006 to December 29, 2017, resulting in T=3,019 business days. The returns are computed as the log difference of the daily closing price, and demeaned prior to estimation.








The number of factors is determined as p=1 from the conventional factor analysis on the k=5 series. The first conventional factor explains 83% of variations, and the second one only an additional 8%. A screen-plot shows a significant drop from the first to the second factor. Further, an analysis of predictive performance, which is explained below, confirms that the model with p=1 factors performs better than the one with p=2.


In theory, there is no ideal method to determine a selection of the first p variables in _t, which technically induce the factor processes under the triangular identification structure embedded in the factor-loading matrix . Ideally, those variables are strongly correlated to other series to efficiently extract the common factors. A simple practical approach taken here is to compute a correlation matrix of _t and to sum up  each row of the matrix, which represents a degree of correlation with other series. In the dataset, the Industrials sector is the series with the highest degree of correlation, and is selected as the first variable in _t, as shown in Table <ref>.


For the computation, the following priors are assumed: μ_i ∼ N(-11,1), (ϕ_i+1)/2 ∼ B(20,1.5), σ_i^-2∼ G(20, 0.01), (ρ_i+1)/2 ∼ B(1, 1), ν_i ∼ G(24, 0.8)I[ν_i>4], β_i ∼κ N(β_i|0, 10) + (1-κ)δ_0(β_i), and κ∼ B(2, 2), where B and G denotes the beta and gamma distributions, respectively. These priors mainly reflect those assumed in the existing studies <cit.>. The size of MCMC iteration is 50,000 samples after the burn-in period of 5,000 samples. To check a convergence of the MCMC draws, the convergence diagnostic of <cit.> is computed. In the following analyses, the null that the Markov chain converges is not rejected at 5 percent significance level.


To measure the efficiency of the MCMC algorithm, the relative numerical efficiency measure, or so-called effective sample size is computed <cit.>. The effective sample size for selected parameters appears to be mostly less than 100, which is comparable to the standard univariate SV models in the existing studies. This figure shows that the developed MCMC algorithm mixes well and yields a practically efficient estimation scheme.

The analysis focuses on the following five FSV specifications:
	
	
  * Model S0: Symmetric t-distribution for _t and _t (β_i≡0).
	
  * Model SY: Skew t-distribution for _t, symmetric t-distribution for _t (β_i=0, for i>p), no sparsity on β_i for all i (i.e., κ≡ 1).
	
  * Model SF: Skew t-distribution for _t, symmetric t-distribution for _t (β_i=0, for i≤ p), no sparsity on β_i.
	
  * Model SYF: Skew t-distribution for both _t and _t, no sparsity on β_i.
	
  * Model SSYF: Skew t-distribution for both _t and _t, with the sparsity prior on β_i, for all i.
	
A key focus here is on which skewness structure– for only the idiosyncratic shocks, latent processes, or both of them– is desirable and how the sparsity skew selection works in forecasting the multivariate returns.






 §.§ Posterior estimates










Prior to the forecasting exercise, posterior estimates based on a fit of the SSYF model are reported. Figure <ref> plots posterior means of the SV process for the factor process and idiosyncratic shocks. The SV of the factor process (F1) exhibits the higher level of volatility and the largest variation among the SV processes, which indicates the variation explained by the latent factor dominates the one by each idiosyncratic shock. The volatilities of five stock indices are highly correlated and their common variation is summarized well by the latent factor. The factor SV increases largely at the time of market turmoil in the global financial crisis around 2009 and the European debt crisis around 2011. The SV of idiosyncratic shock to the Energy sector (Y3) exhibits a different trajectory from  other sectors. In particular, it increases relatively higher than the other series in late 2014 and afterwards. This kind of decoupling in the SV process could reflect an industry-specific business environment. The oil price and other commodity price dropped significantly from late 2014 to 2016, which may have driven up the volatility of stock prices in the energy sector.
















Figure <ref> shows a result of posterior estimates for selected parameters obtained from the SSYF model, with posterior median exhibited in dots and solid line, the 50% and 90% credible intervals in dark and light filled area, respectively.  The posterior median of ϕ_i is quite high, around 0.985 to 0.998, which suggests a strong persistence in the SV process. The parameter estimates for the factor is notably different from those for the idiosyncratic shocks. The factor process is characterized by a lower value of ϕ_i and a higher value of σ_i and μ_i than the idiosyncratic shocks, which implies a larger variation and a higher level of the factor SV process. In addition, the posterior median of ρ_i for the factor is about -0.4 with its 90% credible intervals excluding zero. In contrast, those of all the idiosyncratic shocks include zero. This result suggests that the leverage effect is observed only in the factor process, but not in the idiosyncratic shocks. The degree ν_i of freedom in the GH skew-t distribution is lower for the Consumer Staples (Y4) and Information Technology (Y5) sectors than the others, which indicates heavier-tails of idiosyncratic shocks in those sectors.


Figure <ref> plots the posterior estimates of the skewness parameter β_i from Models SYF and SSYF. For both models, the posterior median of β_i for the factor (F1) is about -1.0 and the 90% credible intervals do not include zero. In contrast, β_i's for idiosyncratic shocks are mostly less relevant. In Model SYF, the 90% credible intervals for Y1 to Y3 include zero. The posterior distributions for those sectors are remarkably shrunk towards zero in Model SSYF. While the 90% credible intervals for Y5 excludes zero in Model SYF, the posterior distribution is partly shrunk in Model SSYF with its posterior median at zero. The posterior probability of β_i=0 is about 80–95% for Y1 to Y3 and 70% for Y5. The Consumer Staples sector (Y4) appears to have relevant skewness in its idiosyncratic shocks with the posterior probability of β_i=0 only 1%.







 §.§ Forecasting performance



A recursive forecasting exercise is examined to study the model comparison among the competing specifications. The analysis is based on a sequence of updating data and forecasting returns one to five days ahead, in order to illustrate a practical use of the model in the financial market. Specifically, first, the model parameters are estimated for a period, over t=1…,T_1, using the posterior computation with the MCMC algorithm. Through the MCMC algorithm, samples of returns (y_T_1+1,…,y_T_1+5) are generated by simulating the latent process from one to five days ahead to obtain the posterior predictive density π(y_T_1+h|y_1,…,y_T_1) for h=1,…,5. We then update the returns for five business days and estimate the model for the returns (y_1,…,y_T_1+5) using the full MCMC analysis with the computation of posterior predictive density for the next five days, (y_T_2+1,…,y_T_2+5), where T_2=T_1+5. We repeat this prediction exercise recursively for T_j=T_j-1+5, j=2,…,J, until we obtain forecasts for 500 business days, that is, J=100.


We set T_1=2,519, which corresponds to the first pre-sample period, to January 6, 2016; the forecasting exercise period goes from January 7, 2016 to December 29, 2017. This experiment reflects practical forecasting activity in the financial market, given the historical data over about 10 years at the starting point, recursively forecasting returns in the following week with the updated data. The exercise period is quite long and includes periods of both bull and bear market sentiments.


The forecasting performance of each model is compared based on two diagnostics. The first one is out-of-sample predictive density. Specifically, the log predictive density ratio (LPDR) for Model i's h-day ahead forecast is computed as
	
    LPDR_h(M_i)  = ∑_j=1^J log{ p_M_i(_T_j+h|_1:T_j)/p_M_0(_T_j+h|_1:T_j)},

where p_M denotes the predictive density under model M, and T_1 denotes the set of time points for the recursive forecasting. This quantity measures how high the likelihood is, when evaluated at the realized return (_T_j+h), that is, how accurate the predictive density is, relative to the base model M_0.

Table <ref> reports the result of the LPDR for five models relative to Model S0, showing the LPDR at each horizon h=1,…,5 and the sum of all the horizons. Because all the LPDRs are positive, the forecasting performance of the newly introduced skew-t FSV models is superior to the standard FSV model with the symmetric t-distribution. Model SF performs better than Model SY, which suggests that the skewness in the factor process is more important for prediction than the ones in the idiosyncratic shocks. Interestingly, Model SYF performs worse than Model SF. This result suggests that the statistically irrelevant skewness parameters for the idiosyncratic shocks, as observed in the left panel of Figure <ref>, worsens the predictive ability of the skew FSV model, unless the shrinkage method is applied. The LPDRs of Model SSYF, which is the highest among the competing models for all the horizons, confirms that the shrinkage method for the skewness parameters effectively selects statistically relevant ones and consequently improves the forecasting performance.
 








The second forecasting exercise is based on Value-at-Risk (VaR). The VaR is a standard tool to measure a risk of possible shortfalls in portfolio. We compute the VaR for portfolios that consist of the five stock indices using the forecasts from each competing model, and evaluate the accuracy of the VaR forecast by comparing it with the realized loss of the portfolio return. This exercise focuses on the prediction ability of the skew-t FSV models for a tail risk of the stock returns, which has practical importance in portfolio management.


We consider a standard portfolio allocation strategy based on a mean-variance optimization following <cit.>. The analysis uses standard Bayesian mean-variance optimization. Let _t and _t denote the forecast mean vector and variance matrix of _t, respectively. Note that these variables are computed from the posterior predictive distribution obtained in the recursive forecasting, as illustrated in the previous exercise. Define _t as a vector of portfolio weights, and m̅ as a target level of return. We solve the optimization problem with respect to the portfolio weights _t as it minimizes the forecast variance of the portfolio return among the restricted portfolios whose expectation is equal to m̅. We also assume that the amount of money invested into the portfolio is fixed during the sample period, and that the resource can be freely reallocated to arbitrary long or short positions without any transaction cost.


The mathematical form of the optimization problem is minimizing _t'_t_t, subject to _t'_t = m̅, and _t' = 1. The solution is _t^(m̅) = _t (' _t _t _t - _t' _t _t ), where _t = (m̅ - _t) / d_t, and d_t = (' _t )(_t' _t _t) - (' _t _t)^2, with _t = _t^-1. Note that the portfolio is reallocated every business day. This experiment assumes a practical situation in which investors allocate their resource every day based forecasts updated weekly. The analysis uses a range of daily target returns of m̅=0.005%, 0.01%, and 0.02%, implying a yearly return of approximately 1.25%, 2.5%, and 5.0%, respectively. In addition, we consider a target-free portfolio that is implemented by minimizing the optimization problem only with _t' = 1.

For a diagnostic of model comparison in this VaR forecast exercise, we employ the test statistic of <cit.>. Let n denote the number of VaR violations, and N the total number of forecasting days. Note that 
the expected number of violations for α quantile is α N. Under the null hypothesis that the expected ratio of violations is equal to α, the likelihood ratio statistic is given by
	
    2log{(n/N)^n (1-n/N)^N-n} - 2log{α^n (1 - α)^N-n}.

As this statistic is asymptotically distributed as χ^2(1), we use it to test whether the VaR forecast is accurate <cit.>. The analysis examines the likelihood ratio test for α=0.5%, 1%, and 5% levels.







Table <ref> reports the number of VaR violations and results of the likelihood ratio test. The value with the asterisk (“å”) shows that the null hypothesis is rejected at a 10% significance level. For Model S0, the target-free portfolio does not pass the test. For Models SY, SF, and SYF, some of the target-return portfolios exhibit a rejection of the test. In particular, the number of violations in Model SY indicates its optimistic VaR forecast, possibly due to the lack of skewness parameter in the factor process. Finally, Model SSYF provides no portfolio that is rejected in the selected range of portfolio allocation rules and VaR levels. This confirms that the proposed skew selection strategy works remarkably well for forecasting VaR in portfolio risk management.






§ A HIGHER-DIMENSIONAL STUDY: US INDIVIDUAL STOCKS










The second empirical analysis uses daily returns of 20 individual stocks to illustrate a high-dimensional application of the skew FSV models. The stocks are randomly selected from the compositions of S&P500 index, and listed in Table <ref>. The selected set of securities covers major industrial sectors. The sample period is T=751 business days from January 3, 2014 to December 30, 2016. As in the previous analysis, the returns are computed as the log difference of the daily closing price.

We examine the number of factors p=1 and 2 from the conventional factor analysis. The first and second conventional factors explain more than 50% of variations in total. In addition, an additional factor does not appear to improve the return performance in the following forecasting exercise. The Morgan Stanley and Boeing are selected for the first two series of _t for the p=2 factors. These two countries from the Financials and Industrials sectors are considered to parsimoniously summarize common movements among the individual stock series. The model with p=1 loads Morgan Stanley for the first factor. Note that a robustness check with randomly selected stocks for the first two series revealed that this selection of the first two series do not significantly affect the main result and its implication in this analysis.


A model comparison is examined based on out-of-sample forecasting performance in the recursive forecasting exercise in the same manner as the previous analysis. The posterior predictive distributions are computed for one- to five-day horizons. The first set of data spans to January 5, 2016 and the forecasting period is from January 5, 2016 to December 30, 2016, to obtain the forecasts over 250 business days. Cumulative returns from the optimized portfolios allocations with the target returns (m̅) are compared among the competing models. Model SY is excluded from the model comparison because it performs relatively less well in the previous analysis. For the posterior computation, all the prior settings are the same as the previous analysis.








Table <ref> reports the cumulative returns over the forecasting period for portfolios with three target returns. Model SSYF provides the uniformly highest cumulative returns both for p=1 and p=2 among the competing models for each target return. The return from the p=1 SSYF model performs slightly better than the p=2 model. This result confirms that the skew FSV model with parsimonious skew selection structure improves the forecasting performance. Interestingly, Model SF is outperformed by Model S0 for p=2, while Model SYF outperforms Model S0, which indicates that the skewness for idiosyncratic shocks is partly relevant to describe dynamics of the high-dimensional multivariate returns.









Figure <ref> shows posterior estimates of selected parameters from Model SSYF with p=1. The posterior estimates of ϕ_i indicates quite high persistence for the SV processes of both the factors and idiosyncratic shocks, while some of the idiosyncratic shocks exhibit relatively wider credible intervals. For the parameter ρ_i, the posterior median of the factor is about -0.7, suggesting that the leverage effect exists as a common behavior among the individual stocks. Additionally, several idiosyncratic shocks (Y1, Y7, Y17) exhibit relevant leverage effect. The posterior estimates of the degree ν_i of the freedom are quite low for several stocks, with some of the posterior medians being below 10. This estimate indicates that the price movement of individual stocks is relatively more volatile and the return distribution is more heavily-tailed than the sector indices examined in the previous analysis.


The estimates of β_i show considerable shrinkage for most of the idiosyncratic shocks. While some of them exhibit that the posterior medians or credible intervals depart from zero, 16 of 20 idiosyncratic shocks exhibit a posterior median at zero. For Y18 to Y20, the skewness is relatively more relevant than the other series with the 50% credible intervals excluding zero. The result in the forecasting performance partly reflects this relevant skewness for idiosyncratic shocks. The β_i for the factor is around -1.0, which indicates the importance of the skewness in the factor process. 








§ CONCLUDING REMARKS



This paper develops the FSV model with GH skew-t distribution. The Bayesian sparsity prior is embedded to the skewness parameters, which provides a parsimonious skew structure for the FSV model. Two empirical analyses using daily stock returns of US-sector indices and individual stocks show that the advantage of the skew selection structure that effectively shrinks irrelevant skewness parameters and consequently improves forecasting performance.


Possible extensions can be explored in several directions. The factor loadings can be time-varying and the factor process can be more complicated dynamics, such as VAR and VARMA processes. From another perspective, it is of interest to model a time-varying skewness that take non-zero values when relevant and zero otherwise using time-varying sparsity techniques <cit.>.




1.4

chicago




