Efficient Bayesian inference for nonlinear state space models with univariate autoregressive state equation
Alexander Kreuzer Claudia Czado






	
			

			Efficient Bayesian inference for nonlinear state space models with univariate autoregressive state equation
    Alexander Kreuzer [Corresponding author: E-mail: a.kreuzer@tum.de]  Claudia Czado
    
	Zentrum Mathematik, Technische Universität München

	December 30, 2023
===========================================================================================================






Latent autoregressive processes are a popular choice to model time varying parameters. These models can be formulated as nonlinear state space models for which inference is not straightforward due to the high number of parameters. Therefore maximum likelihood methods are often infeasible and researchers rely on alternative techniques, such as Gibbs sampling. But conventional Gibbs samplers are often tailored to specific situations and suffer from high autocorrelation among repeated draws. We present a Gibbs sampler for general nonlinear state space models with an univariate autoregressive state equation. For this we employ an interweaving strategy and elliptical slice sampling to exploit the dependence implied by the autoregressive process. 
Within a simulation study we demonstrate the efficiency of the proposed sampler for bivariate dynamic copula models.
Further we are interested in modeling the volatility return relationship. Therefore we use the proposed sampler to estimate the parameters of stochastic volatility models with skew Student t errors and the parameters of a novel bivariate dynamic mixture copula model.
This model allows for dynamic asymmetric tail dependence.
Comparison to relevant benchmark models, such as the DCC-GARCH or a Student t copula model, with respect to predictive accuracy shows the superior performance of the proposed approach. 






§ INTRODUCTION


There are many situations where statistical models with constant parameters are no longer sufficient to appropriately represent certain aspects of the economy. For example it is well known that volatility of financial assets changes over time (<cit.>). This is why many models that allow for variation in the parameter have been proposed. There are time varying vector autoregressive models (<cit.>, <cit.>), stochastic volatility models (<cit.>), GAM copula models (<cit.>, <cit.>) and many more. Stochastic volatility models and the bivariate dynamic copula model of <cit.> assume that the parameter follows a latent autoregressive process of order 1 (AR(1) process). These two models belong to the class of models that we will study.




In a general time varying parameter framework we consider a d dimensional random variable  at time t, Y_t∈ℝ^d, which is generated from a d dimensional density f(·|s_t). We are interested in models, where the density f(· | s_t) has a univariate dynamic parameter s_t ∈ℝ following an AR(1) process. These models can be formulated as state space models with observation equation  

    Y_t|s_t    ∼ f(y_t|s_t)  independently,

for t=1, …, T. The state equation describes an AR(1) process with mean parameter μ∈ℝ, persistence parameter ϕ∈ (-1,1)  and standard deviation parameter σ∈ (0,∞) and is given by

    s_t = μ + ϕ (s_t-1 - μ) + σϵ_t,

where ϵ_t ∼ N(0,1)  iid for t=1, …, T and s_0|μ,ϕ,σ∼ N(μ, σ^2/1 - ϕ^2). In the state equation we assume Gaussian innovations ϵ_t but in the observation equation we do not put any restrictions on the density f. Thus we allow for nonlinear and non Gaussian state space models. Several established models can be analyzed within this framework. 



By choosing f(·|s_t) as the univariate normal density with mean 0 and variance exp(s_t), denoted by φ(·|0,exp(s_t)), we obtain the stochastic volatility model (<cit.>) given by

    Y_t|s_t ∼φ(y_t|0,exp(s_t))  independently, 
        s_t = μ + ϕ (s_t-1 - μ) + σϵ_t,

for t=1, …, T. By modeling the log variance as a latent AR(1) process this model allows for time varying volatility. Time varying volatility is a stylized fact of financial time series. In this context, the stochastic volatility model  has also shown better performance than the frequently used GARCH models (<cit.>,  <cit.>) for several data sets (<cit.>, <cit.>). 
To allow for heavy tails and skewness, other distributions have been considered in the observation equation. One example is the stochastic volatility model with skew Student t errors (<cit.>), which can also be analyzed within our framework.



Dependence modeling is another research area, where models that allow for time varying parameters have been introduced. Vine copulas (<cit.>, <cit.>, <cit.>) are widely used models to capture complex dependence structures. To name a few, <cit.> and <cit.> employ vine copulas for forecasting the value at risk of a protfolio,  <cit.> gives an overview of applications of vine copulas in finance including asset pricing, credit risk management and portfolio optimization and  <cit.> model the association pattern between gap times with D-vine copulas to study asthma attacks.
A vine copula model is made up of different bivariate copulas with corresponding dependence parameters. Since dependencies may change over time, extensions that allow for variation in the dependence parameter have been proposed. <cit.> introduce a bivariate copula model, where the dependence parameter follows a generalized additive model. 
Another approach is the dynamic bivariate copula model proposed by <cit.> and <cit.>, 
which we will analyze within our state space framework.
For this model, we consider single parameter copula families for which there is a one-to-one correspondence between the copula parameter, denoted by θ, and Kendall's τ. So we can express Kendall's τ as a function of the copula parameter θ and we write τ(θ). The restriction of Kendall's τ to the interval (-1,1) is removed by applying the Fisher' Z transformation F_Z(x) = 1/2log(1+x/1-x). This transformed time varying Kendall's τ is then modeled by an AR(1) process. More precisely, we consider T bivariate random vectors, (U_t1, U_t2)_t=1, …, T∈ [0,1]^T× 2, corresponding to T time points. We assume for t=1, …, T that

    (U_t1, U_t2)|θ_t ∼ c(u_t1,u_t2; θ_t)  independently
       s_t = μ + ϕ (s_t-1 - μ) + σϵ_t,  for   s_t = F_Z(τ(θ_t))),

where c(u_t1,u_t2; θ_t) is a bivariate copula density with parameter θ_t. 

Nonlinear state space models as specified with (<ref>) and (<ref>) are typically difficult to estimate, since there is a large number of parameters and likelihood evaluation requires high dimensional integration. This often makes maximum likelihood approaches infeasible. Gibbs sampling (<cit.>) is a frequently used Bayesian approach to infer parameters of such nonlinear state space models (<cit.>). But the posterior samples, resulting from conventional Gibbs sampling, often suffer from high autocorrelation. Furthermore, the availability of the full conditional distributions or at least an efficient MCMC approach to sample from them is required. This is often tailored to specific situations. We present a Gibbs sampling approach that is designed to handle models with a latent AR(1) process and general likelihood functions as specified by the state space formulation in Equations (<ref>) and (<ref>). 
To sample from the associated posterior distribution we rely on elliptical slice sampling (<cit.>) and on an ancillarity-sufficiency interweaving strategy (<cit.>). Elliptical slice sampling is used to sample the latent states. This allows us to exploit the Gaussian dependence structure, that is implied by the AR(1) process. But even if we provide efficient methods to sample from the full conditionals, the sampler may still suffer from the dependence among the parameters in the posterior distribution. Additionally, its performance may vary for different model parameterizations (<cit.>, <cit.>). This problem is tackled with the ancillarity-sufficiency interweaving strategy, where the parameters of the  latent AR(1) process are sampled from two different parameterizations. The decision between two parameterizations is avoided by using both.   This approach has already shown good results for several models, including univariate and multivariate stochastic volatility models (<cit.>, <cit.>).
The efficiency of our proposed sampler is illustrated with a simulation study. 

The second part of this work has a more applied focus and deals with modeling the volatility return relationship, i.e. the dependence  between an index and the corresponding volatility index. More precisely, we investigate the American index S&P500 and its volatility index the VIX as well as the German index DAX and the VDAX.
 It is important to provide appropriate models for this relationship, since it has influence on hedging and risk management decisions (<cit.>). 
 For our analysis we make use of a two step approach commonly used in copula modeling    (<cit.>) motivated by Sklar's Theorem (<cit.>). We first model the marginal distribution with a univariate skew Student t stochastic volatility model. In the second step we model the dependence for which we propose a dynamic copula model allowing for asymmetric tail dependence. This model is a dynamic mixture of a Gumbel and a Student t copula and can be seen as an alternative to the symmetrized Joe-Clayton copula of <cit.>.
Estimation is carried out through a two step approach, where we first estimate the marginal stochastic volatility models, fix their parameters at point estimates and then estimate the dynamic mixture copula. At both steps, estimation is straightforward with the proposed sampler.
Our model is able to capture several characteristics of the joint distribution of volatility and return. With respect to the marginal distribution we observe positive skewness for volatility indices compared to slight negative skewness for the return indices. In the dependence structure we identify asymmetry and time variation. Finally, we compare the proposed model to several restricted models with constant or symmetric dependence and to a bivariate DCC-GARCH model (<cit.>). Model comparison with respect to predictive accuracy shows the superiority of our approach.


 in which inference is straightforward with the proposed sampler.
 
  estimate the dependence among the marginals with the dynamic mixture copula model allowing for time varying asymmetric tail dependence. 




Further, we propose a dynamic copula model allowing for asymmetric tail dependence in which inference is straightforward with the proposed sampler. The model is a dynamic mixture of a Gumbel and a Student t copula and can be seen as an alternative to the symmetrized Joe-Clayton copula of <cit.>. 





The necessity of such models, that allow for dynamic asymmetric tail dependence, was already observed by <cit.> for exchange rates data. It will here be illustrated with an application to different financial data. 




  



To summarize, the main contribution of this paper is an approach to efficiently sample from the posterior distribution of general nonlinear state space models as specified by Equations (<ref>) and (<ref>). In addition, we propose a dynamic mixture copula for time varying asymmetric tail dependence. We discuss Bayesian inference for this model class and demonstrate how it can be utilized to model the volatility return relationship.

The outline of the paper is as follows: After the introduction, we discuss the proposed MCMC approach in Section <ref>. In Section <ref> we investigate the efficiency of the sampler for bivariate dynamic copula models through an extensive simulation study. Section <ref> deals with modeling the volatility return relationship and Section <ref> concludes.









§ BAYESIAN INFERENCE



In the following, we denote an observation of the d-dimensional random vector Y_t  by y_t and the associated data matrix is given by Y=(y_1, …, y_T)^⊤∈ℝ^T × d.
To subset vectors and matrices we make use of the following notation:
For sets of indices A and B we set
x_A = (x_i)_i ∈ A for a vector x and X_A;B = (x_ij)_i ∈ A, j ∈ B for a matrix X. We use a capital letter to refer to a matrix and small letters to refer to its components. The set {n,…, k} of integers will be abbreviated by n:k.



We consider the state space model as specified by Equations (<ref>) and (<ref>).
To obtain a fully specified Bayesian model we equip the parameters μ, ϕ and σ with prior distributions. We follow <cit.>, who propose the following prior distributions for the latent AR(1) process of the stochastic volatility model:

    μ∼ N(0,σ_μ^2),   ϕ+1/2∼ Beta(a_ϕ, b_ϕ) ,  σ^2 ∼ Gamma(1/2, 1/2B_σ),

where σ_μ, a_ϕ, b_ϕ, B_σ >0. Our standard choice for the prior hyperparameters is σ_μ=100,  a_ϕ = 5,  b_ϕ = 1.5 and B_σ=1 as in <cit.>.
With these prior distributions our Bayesian model is complete.  
For sampling from the posterior distribution of this model we should take into account that sampling efficiency may highly depend on the model parameterization (<cit.>, <cit.>). <cit.> differentiate between two parameterizations: A sufficient augmentation and an ancillary augmentation scheme. In our case a sufficient augmentation is characterized by an observation equation that is free of the parameters μ, ϕ and σ and only depends on the latent states s_1:T. In this case s_1:T is a sufficient statistics for the parameters μ, ϕ and σ. In an ancillary augmentation the state equation is independent of the parameters μ, ϕ and σ, then s_1:T is an ancillary statistics for the parameters μ, ϕ and σ.
 The standard parameterization of our model is already a sufficient augmentation and we refer to this parameterization as given by Equations (<ref>) and (<ref>) as (SA).

    (SA):      Y_t|s_t ∼ f(y_t|s_t)  independently
        s_t = μ + ϕ (s_t-1 - μ) + σϵ_t.

  An ancillary augmentation is obtained by the following parameterization

    s̃_t = s_t - μ -ϕ(s_t-1 - μ)/σ,  with inverse  s_t = μ + ϕ(s_t-1 - μ) + σs̃_t,

for t=1 ,…, T. This reparameterization is obtained by solving Equation (<ref>) for ϵ_t and implies that the state space model is given by

    (AA):      Y_t|s̃_1:T,s_0, μ, ϕ, σ∼ f(y_t|s_t(s̃_1:T,s_0, μ, ϕ, σ)) independently
       s̃_t ∼ N(0,1)  independently,

where s_t(s̃_1:T,s_0, μ, ϕ, σ) is the function that calculates s_t according to (<ref>).  We refer to this model representation as (AA). Instead of deciding between (SA) and (AA), we combine them in an ancillarity-sufficiency interweaving strategy (<cit.>), given by


  * a)  Sample  s_0:T in (SA) from
s_0:T|Y,  μ, ϕ, σ
.

  * b) Sample (μ, ϕ, σ) in (SA) from 
μ, ϕ, σ|Y,s_0:T.

  * c) Move to (AA) via
s̃_t = s_t - μ -ϕ(s_t-1 - μ)/σ,
for t=1, …, T.

  * d) Sample (μ, ϕ, σ) in (AA) from 
μ, ϕ, σ|Y, s_0,s̃_1:T .

  * e) Move back to (SA) via the recursion
s_t = μ + ϕ(s_t-1 - μ) + σs̃_t
for t=1, …, T.


<cit.> employed interweaving for the stochastic volatility model and showed its superior performance with an extensive simulation study. For the stochastic volatility model, <cit.> propose to  move between a sufficient augmentation and a reparameterization of the latent states s_1:T  given by s^K_t = s_t - μ/σ. Within this reparameterization parameters can be sampled conveniently from its full conditional distribution by recognizing a linear regression model. This is possible for the standard stochastic volatility model, but not in our case since our sampler is designed to handle more general likelihood functions. Therefore we have chosen the reparameterization of (SA) such that it is optimal in the sense of <cit.>, i.e. we move between a sufficient and a ancillary augmentation.


Note that reducing the sampler to the first two steps a) and b) results in a standard Gibbs sampler in (SA). This sampler typically suffers from the dependence among the parameters μ, ϕ, σ and the latent states s_1:T in the posterior distribution.








  §.§.§ Step a: Sampling of the latent states in the sufficient augmentation 


To sample the latent states s_0:T from its full conditional in (SA) we make use of elliptical slice sampling as proposed by <cit.>. It was developed  for models, where dependencies are generated through a latent multivariate normal distribution.
In (SA), the AR(1) structure implies, that the vector s_0:T|μ, ϕ, σ has a (T+1) dimensional multivariate normal distribution with mean vector μ^AR and covariance matrix Σ^AR given by

    μ^AR = ( [ μ; μ; ⋮; μ ])
    ∈ℝ^T+1,    Σ^AR = σ^2/1-ϕ^2( [     1     ϕ   ϕ^2     …   ϕ^T;     ϕ     1     ϕ     … ϕ^T-1;     ⋮     ⋮     ⋮           ⋮;   ϕ^T ϕ^T-1 ϕ^T-2     …     1;       ])
       ∈ℝ^(T+1) × (T+1).

(see e.g. <cit.>, Chapter 2.2).
The posterior density is proportional to

    ( ∏_t=1^T f(y_t|s_t) )  φ(s_0:T|μ^AR, Σ^AR) π(μ) π(ϕ) π(σ),

where μ^AR and Σ^AR are given in (<ref>) and π(·) denotes the corresponding prior density as specified in (<ref>).
The initial state can be sampled from its full conditional density given by

    f(s_0|s_1:T, μ, ϕ, σ) = φ(s_0| μ+ϕ(s_1-μ), σ^2).


The full conditional density of the latent states s_1:T is given by

    f(s_1:T|Y, s_0, μ, ϕ, σ) ∝(∏_t=1^T f(y_t|s_t)) φ(s_1:T|μ_1:T|0, Σ_1:T|0),

with a corresponding mean vector μ_1:T|0 and covariance matrix Σ_1:T|0. The mean vector and the covariance matrix of the conditional distribution are derived in a more general way in Appendix A.
By reparameterizing the model with
s_1:T' = s_1:T - μ_1:T|0,
we impose a multivariate normal prior with zero mean. We obtain the situation elliptical slice sampling was designed for. However updating the whole T dimensional vector s_1:T with elliptical slice sampling at once will lead to high autocorrelation in the posterior draws. This is illustrated in Section <ref> and was also observed by <cit.>, where elliptical slice sampling was used for linear regression models. <cit.> circumvent this problem by partitioning the vector s_1:T into smaller blocks. 

We follow this approach and partition the set {1, …, T} into m different blocks B_1, …, B_m ⊂{1, …, T}. Let a_i = min_s ∈ B_i s denote the minimal and b_i = max_s ∈ B_i s denote the maximal index in the i-th block. The blocks are chosen such that  B_i = {t ∈{1, …, T}: a_i ≤ t ≤ b_i }, for i = 1, …, m.
The full conditional density for the i-th block can be expressed as
f(s_B_i|Y, s_0, s_-B_i, μ, ϕ, σ) ∝(∏_t ∈ B_i f(y_t|s_t)) f(s_B_i|s_0, s_-B_i, μ, ϕ, σ),
where -B_i = {1, …, T}∖ B_i. The vector s_B_i|s_0, s_-B_i, μ, ϕ, σ is multivariate normal distributed with mean denoted by μ_B_i| and covariance matrix Σ_B_i| (see Appendix A) and therefore the full conditional density can be written as

    f(s_B_i|Y, s_0, s_-B_i, μ, ϕ, σ) ∝(∏_t ∈ B_i f(y_t|s_t)) φ(s_B_i|μ_B_i|, Σ_B_i|).

To sample the latent states of the i-th block s_B_i from its full conditional we proceed as follows


  * Set
s'_B_i = s_B_i - μ_B_i|

  * Draw s'_B_i from the density 

    f(s'_B_i|Y, s_0, s_-B_i, μ, ϕ, σ) ∝(∏_t ∈ B_i f(y_t|s_t)) φ(s'_B_i|0, Σ_B_i|),

using elliptical slice sampling, where φ(s'_B_i|  0,Σ_B_i|) is interpreted as the prior density for s'_B_i.


  * Set
s_B_i = s'_B_i + μ_B_i| 






  §.§.§ Step b: Sampling of the constant parameters in the sufficient augmentation


In (SA) the observation equation only depends on s_1:T and is independent of the parameters μ, ϕ and σ. The parameters μ, ϕ and σ only depend on s_0:T. This allows to use the same approach  as in <cit.> to sample the parameters μ, ϕ and σ in (SA). We reparameterize the model such that proposals can be found using Bayesian linear regression. We define γ = μ (1-ϕ) and the state equation is given by

    s_t = γ + ϕ s_t-1 + ση_t,

where η_t ∼ N(0,1). For fixed s_0:T, this is a linear regression model with regression parameters γ, ϕ and variance σ^2. Proposals for (μ, ϕ, σ) are found and accepted or rejected as described in <cit.> Section 2.4 (two block sampler).





  §.§.§ Step d: Sampling of the constant parameters in the ancillary augmentation


To sample μ, ϕ and σ in (AA) we deploy a random walk Metropolis-Hastings scheme with Gaussian proposal, where the proposal variance or covariance matrix is adapted during the burn-in period. For the adaptions we use the Robbins Monro process (<cit.>) as suggested by <cit.>. More details are given in Appendix A.



  §.§.§ Implementation

For the implementation of the sampler we use  (<cit.>) which allows to embed C++ code into . In addition we make use of  (<cit.>) to evaluate copula densities and of  (<cit.>). For sampling (μ, ϕ, σ) in (SA) we use corresponding parts of the implementation of the  package  (<cit.>). The  package  (<cit.>) is used to compute effective sample sizes in the following section.





§ ILLUSTRATION OF THE PROPOSED SAMPLER FOR BIVARIATE DYNAMIC COPULA MODELS



We illustrate the MCMC sampler we proposed in the previous section for the bivariate dynamic copula model of <cit.>.
<cit.> have already shown that interweaving improves sampling efficiency a lot for the stochastic volatility model.
We investigate if this is also the case for the bivariate dynamic copula model.
Further we study how the sampling efficiency is affected by the chosen block size and by the data generating process (DGP).
Therefore we perform an extensive simulation study. 
We consider different modifications of the sampler. A sampler is specified by a vector (b,i) which indicates its blocksize (b) and if interweaving is used (i=I) or not (i=NI). We consider ten different sampler specifications (b,i) ∈{1,5,20,100,T}×{I, NI}, where T is the length of the time series. By using blocks of size T, we obtain the sampler which updates the parameters s_1:T jointly with elliptical slice sampling. If we turn off interweaving (i=NI) we obtain a standard Gibbs sampler updating parameters in the sufficient augmentation. These samplers are run for different simulated data sets. A data set is simulated from the bivariate dynamic copula model (see (<ref>)) with parameters: Family, T, μ, ϕ, σ. The parameters are chosen from the following grid  
(Family, T, μ, ϕ, σ) ∈{Gauss, eClayton} ×{500, 1000, 1500}×{0,1}×{0,0.1,0.5,0.9,0.99}×{0.05,0.1,0.2}.
Here, eClayton denotes the extended Clayton copula, which extends the Clayton copula to allow for negative Kendall's τ values. More precisely, the extended Clayton copula 
 has the following density

    c(u_1,u_2;θ) =
       
        c_Clayton(u_1, u_2; θ)  if θ≥ 0 
    
         c_Clayton(1-u_1, u_2; -θ)  if θ < 0,

where c_Clayton(·, ·;θ) is the density of the bivariate Clayton copula with parameter θ (see <cit.>, Chapter 4).
So the extended Clayton copula density is equal to the Clayton copula density for non negative Kendall's τ and equal to a 90 degree rotation of the Clayton copula density for negative Kendall's τ.
Among the different DGPs, most distinct values are considered for  ϕ. We expect its choice to be influential since it controls the dependence among the latent variables. 
With this grid we obtain 180 different DGPs. For each of the different DGPs we generate 100 simulated data sets and for each data set we run the 10 different samplers with the correctly specified copula family for 25000 iterations and discard the first 5000 iterations for burn-in. So, in total we obtain 18000 simulated data sets and each of the 10 samplers is run 18000 times.










Figure <ref> shows trace plots of different parameters (μ, ϕ, σ, s_300) based on one simulated data set for three different sampler specifications. We consider specification (5,I), and the same specification with interweaving turned off, i.e. (5,NI) and the specification (T,I).
We observe that all three samplers produce posterior samples covering the true values. Further, the trace plots suggest that the (5,I) sampler achieves better mixing than the other two considered samplers.  


The runtime of the sampler is mainly affected by the choice of T and the sampler specification. From Table <ref> we see that the runtime is increasing in T and that interweaving adds considerable additional runtime. 










To measure efficiency of the samplers we consider the effective sample size per minute which we call effective sampling rate, similar to <cit.>. We average the effective sampling rate of 100 runs, where the same sampler and the same DGP  was used. Then we obtain 180 average effective sampling rates (AESR) per sampler.
Since the AESR decreases for higher values of T for every sampler, we compare the ten different samplers among DGPs with the same value for T.  For a fixed T we have 60 AESR values per sampler. For each sampler, the minimum of these 60 values (mAESR) is given in Table <ref>.
We consider the minimum since we are interested in samplers that are reliable for all DGPs. 
We see that for the parameters μ, ϕ and σ sampler specifications with interweaving have higher mAESR values, while for the latent states specifications without interweaving perform better.
Although interweaving adds additional runtime, it still increases the mAESR of μ, ϕ and σ considerably.
Further we observe that choosing the blocksize too big results in very low mAESR values for the latent states. In this case many parameters are updated jointly with elliptical slice sampling which results in high autocorrelation among consecutive draws. The performance of samplers with blocksize T is especially poor. For these samplers the length of the time series (T= 500, 1000, 1500) also has strong effects. For example the mAESR for μ for specification (T,NI) decreases by 85 % from 78 to 11 when the length of the time series is increased from 500 to 1500. 
The most inefficient sampler is (T,NI), the sampler without interweaving and with the largest blocksize.
In our opinion the best results are obtained for
sampler specification (5,I). It provides the highest mAESR values for μ, ϕ and σ for all choices of T and also provides rather high mAESR values for the latent states.








In addition to the previous analysis, we investigate how different DGPs affect the samplers. Therefore we consider the best sampler according to Table <ref>, i.e. sampler specification (5,I). In addition we have a look at the same sampler specification without interweaving (5,NI) and the same sampler with joint updates of the latent states, i.e. (T,I). We consider the AESR for σ, which is usually the parameter which causes most problems.
Table <ref> shows for each of these three samplers the DGPs with T=1000 which resulted in the lowest and the highest AESR for σ.
For the (5,I) specification satisfactory AESR values, ranging from 175 to 456, are obtained for all DGPs.
 In the (T,I) specification the latent states are updated jointly. In scenarios with strong dependence among the latent states (ϕ = 0.99) this sampler performs poorly whereas the best performance of the sampler was seen for a DGP with low dependence among the latent states (ϕ=0). The (5,NI) specification is a standard Gibbs sampler in the sufficient augmentation. We see that for this specification the AESR values vary a lot, ranging from 25 to 558. The best performance of this sampler specification was seen for a DGP with high persistence (ϕ=0.99). <cit.>  studied different sampling schemes for the stochastic volatility model and have also seen that a standard Gibbs sampler in the sufficient augmentation performs well for scenarios with high persistence. 








Lastly we compare our sampler to the coarse grid sampler employed by <cit.>. We already covered six DGPs that were also analyzed by <cit.>. Instead of running their sampler we make the comparison with respect to these six cases. There are several points which make the comparison slightly less reliable. First <cit.> did not report exact computation times but they note that 100 000 iterations of their sampler take about 15 minutes. We use this number to calculate the AESR values from the effective sample sizes they report in their paper. Second their calculations were performed on a different computer and third they use different prior distributions for ϕ and σ^2. But we think that this comparison should still give us a rough idea of how the sampling efficiencies compare to each other. From Table <ref> we see that the (5,I) specification considerably outperforms the coarse grid sampler (CG). For every parameter (μ, ϕ, σ) we obtain way higher AESR values. 


























































§ APPLICATION: MODELING THE VOLATILITY RETURN RELATIONSHIP



We investigate the volatility return relationship through the bivariate joint distribution of a stock index and the corresponding volatility index. 
The joint distribution of return and volatility incorporates all the marginal information as well as information about the dependence, which are both relevant for hedging and risk management (<cit.>). Since there has already been evidence for asymmetry in the joint distribution of volatility and return (<cit.>, <cit.>) models which are able to handle these characteristics are necessary.


The two step copula modeling approach motivated by Sklar's theorem (<cit.>) provides a very flexible method for the construction of multivariate distributions. We can combine arbitrary marginal distributions with any copula. Here, we propose a bivariate model that combines the skew Student t stochastic volatility model (<cit.>) for the margins with a novel dynamic mixture copula.  This model allows for asymmetry and heavy tails in the marginal distribution as well as for time varying asymmetric tail dependence in the dependence structure. Both, the marginal as well as the copula model can be estimated with the proposed sampler.




 §.§ Marginal model

The stochastic volatility model with skew Student t errors is obtained by replacing the normal distribution of the stochastic volatility model (<cit.>)  by a skew Student t distribution. This allows for heavy tails and skewness. 
The stochastic volatility model with skew Student t errors as considered by <cit.> is given by

    Y_t = exp(s_t/2) ϵ_t  
       s_t = μ + ϕ (s_t-1 - μ) + σϵ_t,

where ϵ_t|α,df ∼ sst(ϵ_t|α, df) independently for t=1, …, T. We denote by sst(ϵ_t|α, df) the density of the standardized skew Student t distribution with parameters α∈ℝ and df>2 (cf. Appendix B).
Compared to our framework this model has two additional parameters α and df. For these additional parameters we choose the following prior distributions

    α∼ N(0,100),        df ∼ N_>2(5,25),

where N_>2 denotes the normal distribution truncated to (2, ∞). We need to ensure that df > 2 since the standardized skew Student t distribution would not be well defined otherwise. 





    π(δ) = π_nu(exp(δ)+2) * exp(δ)



 Conditional on α and df our sampler can be applied directly to sample (μ, ϕ, σ, s_0:T) from its full conditional. Another approach, which lead to better mixing, is to include the parameters α and df in the interweaving strategy. The sampler is slightly modified in the following way:


  * a)  Sample s_0:T from
s_0:T|Y,  μ, ϕ, σ, α, df
.

  * b) Sample (μ, ϕ, σ, α, df) in (SA) from 
μ, ϕ, σ, α, df|Y,s_0:T .

  * c) Move to (AA) via
s̃_t = s_t - μ -ϕ(s_t-1 - μ)/σ,
for t=1, …, T.

  * d) Sample (μ, ϕ, σ, α, df) in (AA) from 
μ, ϕ, σ, α, df|Y, s_0,s̃_1:T.

  * e) Move back to (SA) via the recursion
s_t = μ + ϕ(s_t-1 - μ) + σs̃_t
for t=1, …, T.



For step a) we proceed as described in Section <ref>. 
For step b) we draw α and df from its univariate full conditional distributions using Metropolis-Hastings, similar to Step d) in Section <ref>.
The parameters (μ, ϕ, σ) are drawn from its full conditional as described in Section <ref>. 

 For step d) we investigated different blocking strategies for the parameters (μ, ϕ, σ, α, df). We compared the different strategies with respect to effective sample sizes and decided to use the following three blocks: (μ,df), (ϕ, σ) and α. Each block is updated using Metropolis-Hastings as in Step d) in Section <ref>.










 §.§ Dependence model

Dependence among financial assets is often modeled with a Student t copula. 
This copula allows for tail dependence symmetric in the upper and lower tail. Evidence against the assumption of symmetric tail dependence has been provided and models to handle this characteristic have become necessary (<cit.>, <cit.>, <cit.>). <cit.> proposes the symmetrized Joe-Clayton  copula. This is a modification of the BB7 copula (<cit.>, Chapter 4) that is symmetric if upper and lower tail dependence coincide, which he describes as a desirable property.
In the application of <cit.> the Student t copula provides the best fit in terms of the likelihood. But they argue that if the focus is on the tails a BB1 or BB7 copula might be more appropriate. The BB1 and BB7 copulas have two parameters which might not be enough, if we want to model three characteristics in a flexible way: upper tail dependence, lower tail dependence and overall dependence as measured with Kendall's τ.  
  We provide another approach to relax the symmetric tail dependence assumption. We propose a mixture of a Student t and a extended Gumbel copula with parameters τ∈ (-1,1), ν >2 and p ∈ [0,1] given by

    C^M(u_1,u_2; τ, ν, p) =  p C^t(u_1,u_2;τ,ν) + (1-p) C^G(u_1,u_2;τ))  ,

where  C^t is the bivariate Student t copula specified by Kendall's τ and the degree of freedom ν and C^G is the bivariate extended Gumbel copula specified by Kendall's τ. The extended Gumbel copula is defined similarly to the extended Clayton copula in Section <ref>, i.e. its density is equal to the Gumbel copula density for positive values of Kendall's τ and equal to a 90 degree rotation of the Gumbel copula density for negative Kendall's τ values. Both copulas C^t and C^G share the dependence parameter τ and we expect the mixture copula to have a similar strength of dependence. 
The corresponding Kendall's τ of the mixture copula is given by

    τ^M
    =   ∫_(0,1)^2 C^M(u_1,u_2;τ, ν, p) c^M(u_1,u_2;τ, ν, p) du_1du_2 =(p^2 + (1-p)^2) τ + 
       +∫_(0,1)^2 p(1-p) ( C^G(u_1,u_2;τ)) (c^t(u_1,u_2;τ,ν) + C^t(u_1,u_2;τ,ν) c^G(u_1,u_2;τ))) du_1du_2.


We evaluated the integral in (<ref>)  numerically for different values of τ, p and ν and observed only negligible difference between τ and τ^M. The upper and lower tail dependence coefficients λ_M^L and λ_M^U of the mixture copula can, for τ>0, be obtained as

    λ_M^L(τ, p, ν) =   lim_u→ 0C^M(u,u)/u = lim_u→ 0p C^t(u_1,u_2;τ,ν) + (1-p) C^G(u_1,u_2;τ)/u
    
    =    p 2 T_ν+1(-√((ν+1)(1-sin(πτ/2))/1+sin(πτ/2)))  + 0
    λ_M^U(τ, p, ν) =       p 2 T_ν+1(-√((ν+1)(1-sin(πτ/2))/1+sin(πτ/2)))  + (1-p)(2 - 2^1-τ),

where we used the well known formulas for the tail dependence coefficients of the Student t and the Gumbel copula (<cit.>, Chapter 4).
Whereas the upper and lower tail dependence coefficients measure dependence in the upper right and lower left corner, we are also interested in the dependence in the upper left and the lower right corner when τ<0.
We consider the following tail dependence coefficients in the upper left corner λ_M^UL and in the lower right corner λ_M^LR if τ<0

    λ_M^LR   = λ_M^L(-τ, p, ν) ,    λ_M^UL = λ_M^U(-τ, p, ν),

analogous to the definition of quarter tail dependence in <cit.>. 



The tail dependence coefficient of the mixture copula is a linear combination of the tail dependence coefficients of its two components, the Student t and the Gumbel copula. The Student t copula has symmetric tail dependence, whereas the Gumbel copula has upper but no lower tail dependence. So we expect upper tail dependence to be higher than lower tail dependence in the mixture copula. The amount of asymmetry in the tails is controlled by p, whereas the copula is symmetric in the tails for p=1 and the level of asymmetry increases as we decrease p.
So this copula allows for great flexibility: The overall dependence can be described by Kendalls's τ, the degrees of freedom parameter controls the upper and lower tail dependence coefficient and p controls the difference between upper and lower tail dependence.  This is visualized in Figure <ref> in Appendix C.
Note that the desirable property according to <cit.> of symmetry in case of coinciding upper and lower tail dependence is here fulfilled.  If we expected higher lower than upper tail dependence we can replace the Gumbel copula by a survival Gumbel copula which has the density c^SG(u_1,u_2) = c^G(1-u_1,1-u_2). 
To allow for time variation we use the mixture copula C^M of (<ref>) within the dynamic bivariate copula model of <cit.>. A nonlinear state space model for T bivariate random vectors (U_t1, U_t2)_t=1, …, T∈ [0,1]^T× 2, corresponding to T time points, is given by

    (U_t1, U_t2) ∼ c^M(u_t1, u_t2; τ_t, ν, p)  independently
       s_t = μ + ϕ (s_t-1 - μ) + σϵ_t,    s_t = F_Z(τ_t)

for t=1, …, T.
We assign a uniform prior on [0,1] for p, a normal prior with mean 5 and standard deviation 20 truncated to the interval (2,∞] for ν and the same priors as in (<ref>) for the remaining parameters. Sampling is done in the following way.


  * Draw log(p/1-p) and log(ν -2 ) from its univariate full conditionals with random walk Metropolis-Hastings with Gaussian proposal (proposal standard deviation: 0.3).

  * Draw μ, ϕ, σ, s_0:T conditioned on p and ν as in Section <ref>.






















































 §.§ Two step estimation

We consider the S&P500 (SPX)  and its volatility index the VIX as well as the DAX and its volatility index the VDAX. The daily log returns from 2006 to 2013 of these indices are obtained from Yahoo finance (https://finance.yahoo.com). With approximately 250 trading days per year this results in 2063 observations, visualized in Figure <ref> in Appendix D. The corresponding data matrix with 2063 rows and 4 columns is denoted by Y.


Combining the marginal and the dependence model we obtain that for T bivariate random vectors (Y_t1,Y_t2)_t=1, …, T∈ℝ^T × 2 the following holds



    (Y_t1,Y_t2) ∼
       C^M(ssT( y_t1/exp(s_t1^st/2) | α_1^st, df_1^st), ssT( y_t2/exp(s_t2^st/2) | α_2^st, df_2^st);F_Z^-1( s_t^cop),  ν^cop,  p^cop)





    c^M(ssT(x_t1| α_1^st, df_1^st), ssT(x_t2| α_2^st, df_2^st);F_z^-1( h_t^cop),  ν^cop,  p^cop) ·
       · sst(x_1| α_1^st, df_1^st) sst(x_2| α_2^st, df_2^st),


independently, where

    s_tj^st   = μ_j^st + ϕ_j^st (s_t-1;j^st - μ_j^st) + σ_j^stϵ_tj^st
    
    s_t^cop   = μ^cop + ϕ^cop (s_t-1^cop - μ^cop) + σ^copϵ_t^cop

and ϵ_tj^st, ϵ_t^cop∼ N(0,1) iid, α_j^st, df_j^st as in (<ref>), ν^cop, p^cop as in (<ref>) and μ_j^st, μ^cop, ϕ_j^st, ϕ^cop, σ_j^st, σ^cop, s_0j^st, s_0^cop as in (<ref>) and (<ref>) for j=1,2 and t=1, … T. Here, ssT denotes the distribution function of the standardized skew Student t distribution (cf. Appendix B). We refer to the probability integral transforms ssT( y_t1exp(-s_t1^st/2)| α_1^st, df_1^st) and  ssT( y_t2exp(-s_t2^st/2)| α_2^st, df_2^st) for t=1, …, T as copula data.







For inference we rely on a two step approach. We first estimate marginal distributions and based on the resulting estimated copula data we estimate the copula parameters. This approach is also called inference for margins (<cit.>) and is commonly used in (Bayesian) copula modeling (<cit.>, <cit.>, <cit.>, <cit.>, <cit.>).  


We first estimate marginal distributions and based on these estimates we estimate the copula parameters. This approach is also called inference for margins (<cit.>) and is commonly used in copula modeling. Due to its advantages, also many Bayesian researchers rely on a two step approach for copula models (<cit.>, <cit.>, <cit.>, <cit.>, <cit.>). 
Here, the likelihood consists of a product of the marginal densities and the copula density c^M(u_1, u_2, θ). 
In our case all three arguments u_1, u_2, θ of the copula are modeled dynamically. We obtained high values for the copula density by choosing high dependence and scale u_1, u_2 appropriately with their dynamical marginal variances. This scaling of the variances may increase the copula contribution to the likelihood, while decreasing the marginal contribution. Thus different parameter choices may lead to similar values of the likelihood. Such multimodal posterior distributions are difficult to estimate in general. Furthermore this lead to results that do not seem reasonable. For example we observed copula data deviating far from being marginally uniformly distributed and very high oscillations in the latent log variances. Although oscillations are supposed to be controlled by the AR(1) prior, we have seen that the effect in the likelihood outweighs the effect of the prior. This shows that a joint Bayesian estimation approach for this setup might run into problems. In a two step approach, s_1:T;j^st is estimated only based on data y_1:T;j (j=1,2) and the above-mentioned problems can be circumvented. Furthermore the two step approach allows for faster inference by dividing the problem into smaller parts, which converge faster. In addition to that, the calculation of copula data several times in every step, as required by a fully Bayesian approach, is here avoided. This is particularly important here since we require expensive evaluations of the standardized skew  Student t distribution function. Given these reasons we restrict ourselves to a two step approach. 












First we fit a skew Student t stochastic volatility model for each of the indices. For each index we run the sampler (5,I) for 31000 iterations and discard the first 1000 draws as burn-in.
As it is typical for financial data all indices show a high persistance parameter ϕ (Posterior mode estimates for ϕ: SPX: 0.99, VIX: 0.90, DAX: 0.99, VDAX: 0.96). A notable difference is that for stock indices we observe negative skewness, whereas for the volatility indices positive skewness is observed (Posterior mode estimates for α: SPX: -0.51, VIX: 1.33, DAX: -0.48, VDAX: 0.96). Evidence for negative skewness has also been observed for the log returns of other stock indices, as e.g., for the NASDAQ by <cit.>.  Posterior mode estimates, posterior quantiles and effective samples sizes for several parameters of the four marginal models are summarized in Table <ref> in Appendix D. The estimated daily log variances are shown in Figure <ref>. In the end of 2009, the estimated variances are high for all indices due to the financial crisis. 








In the next step we obtain data on the [0,1] scale by applying the probability integral transform using the posterior mode estimates of the marginal parameters. We refer to this data as pseudo copula data and it is obtained as

    û_tj = ssT(y_tjexp(-ŝ_tj^st/2);α̂_j^st, d̂f̂_j^st),
 
where ŝ_tj^st, α̂_j^st, d̂f̂_j^st are the posterior mode estimates of the corresponding marginal skew Student t stochastic volatility model for t=1, …, T, j=1,…, 4. 
In the copula data marginal characteristics are removed and what is left is information about the dependence structure.
Based on the pseudo copula data two dynamic mixture copula models are fitted, one corresponding to the pair (SPX,VIX) and one corresponding to the pair (DAX,VDAX). For each pair we obtain 31000 iterations with the sampler specification (5,I) and discard the first 1000 draws as burn-in.  The posterior mode estimates for p are 0.29 for the model for (SPX,VIX) and 0.66 for the model for (DAX,VDAX), respectively. (Further, posterior statistics for the model parameters μ, ϕ, σ, p and ν are shown in Table <ref> in Appendix D). So both fitted models allow for asymmetric tail dependence, whereas the asymmetry is stronger for the (SPX,VIX) model. For these models tail dependence in the upper left corner λ_M^UL is stronger than the one in the lower right corner λ_M^LR. This means that joint extreme comovements, where the stock index decreases and the volatility index increases are more likely to occur than vice versa, which agrees with the statement that the market reacts more extreme in bad market situations (<cit.>). 





 

















The time varying estimates for Kendall's τ and the tail dependence coefficients are shown in Figure <ref>. The figure visualizes the asymmetry in tail dependence. We also observe changes in tail dependence as time evolves: for (SPX,VIX), λ_M^UL ranges from 0.41 to 0.71 and λ_M^LR from 0.004 to 0.10 . For the pair (DAX,VDAX), λ_M^UL ranges from 0.12 to 0.67  and λ_M^LR from 0.003 to 0.35. This variation over time in tail dependence goes hand in hand with variation in Kendall's τ. For (SPX,VIX), Kendall's τ ranges from -0.76 to -0.48 and for (DAX,VDAX) from  -0.80 to -0.30. 










 §.§ Out of sample predictions

We aim to further support the findings we obtained through the dynamic copula model, i.e. that the dependence structure is asymmetric and varies over time. Therefore we consider several restrictions with respect to the dependence structure. Giving up time variation leads a constant mixture copula, giving up asymmetry leads to a dynamic Student t copula and giving up time variation and asymmetry leads to a constant Student t copula. 
In addition, we compare our model to the frequently used dynamic conditional correlation (DCC) GARCH model of <cit.>. The DCC-GARCH allows for time varying symmetric dependence.
So we take five different models into consideration. These models are summarized in Table <ref>. The model ℳ_dyn^mix is the model given in (<ref>), which was also used for the previous analysis.


We compare the models with respect to cumulative pseudo log predictive scores (<cit.>), which are obtained by evaluating the corresponding density at point estimates, instead of averaging over all posterior draws.

 In comparison to other multivariate scoring rules, such as the energy score or the variogram score (<cit.>), pseudo log predictive scores have here the advantage that they can be computed fast since they require only one evaluation of the density per observation. 
We consider T+K observations of dimension two, stored in the data matrix Y_1:(T+K);1:2, where the first T observations are used to train the model and the last K are used for testing.
 


  * ℳ_dyn^mix: A mixture copula  model with dynamic dependence combined with skew Student t stochastic volatility margins as given in (<ref>) which was also used in the previous analysis.

  * ℳ_const^mix: A  mixture copula model with constant dependence combined with skew Student t stochastic volatility margins, i.e we replace F_Z^-1(h_t,cop) in (<ref>) by a parameter τ which does not depend on t. In comparison to ℳ_dyn^mix the parameters h_0:T,cop, μ_cop, ϕ_cop and σ_cop are replaced by a single parameter τ.

  * ℳ_dyn^t: A bivariate dynamic copula model as in (<ref>) based on a Student t copula, i.e. based on a copula with symmetric tail dependence, combined with skew Student t stochastic volatility margins. This model is obtained by replacing the copula density c^M in (<ref>) by the density of bivariate Student t copula, whereas the parameter p is no longer needed.

  * ℳ_DCC: The DCC(1,1)-GARCH(1,1) model of <cit.>.








  
The DCC-GARCH model is estimated based on T data points in the training period. Based on this model, we obtain rolling one-day ahead estimates of the covariance matrix for each day in the test set. The pseudo log predictive scores are obtained  by  evaluating the corresponding multivariate normal log densities at the observations. To estimate the DCC-GARCH models we used the R package  of <cit.>.
 

Similarly, the model ℳ_dyn^mix is estimated with the training data. Instead of computing daily updates for all model parameters we fix the constant parameters at their posterior mode estimates to save computation time. The dynamic parameters are updated daily and the one-day ahead forecasts are obtained by evolving the AR(1) process. To obtain the pseudo log predictive scores we evaluate the density implied by (<ref>) at the corresponding observations. Appendix D contains a detailed description of this procedure. The pseudo log predictive scores for ℳ_const^mix, ℳ_dyn^t and ℳ_const^t are obtained similarly. 


This procedure for calculating the pseudo log predictive scores is applied for both data sets corresponding to the pairs (SPX,VIX) and (DAX,VDAX). Using the last two years (2012 - 2013) of our data set as test data yields K=517. As training period we use T=1000 which corresponds to a training period of approximately four years. 

Table <ref> summarizes the cumulative pseudo log predictive scores.
In both cases, for the (SPX,VIX) as well as for the (DAX,VDAX) data, the best model is provided by the dynamic mixture copula model ℳ_dyn^mix. 
Furthermore, we see that in both cases the constant mixture copula model ℳ_const^mix is preferred over the constant and dynamic Student t copula models ℳ_const^t and ℳ_dyn^t. For the (DAX,VDAX) data, the second best model is provided by the constant mixture copula model ℳ_const^mix. For this data the rolling window estimates of Kendall's τ in Figure <ref> vary less than for the (SPX,VIX) data. For the (SPX,VIX) data, the DCC-GARCH ℳ_DCC yields the second best cumulative pseudo log predictive scores.
 









§ CONCLUSION



We propose a sampler, applicable to general nonlinear state space models with univariate autoregressive state equation. Sampling efficiency is demonstrated for bivariate dynamic copula models within a simulation study. Furthermore we use the sampler to estimate the parameters of a dynamic bivariate mixture copula model. This mixture copula model turns out to be a good candidate to model the volatility return relationship, since in our application it produces more accurate forecasts than a bivariate DCC-GARCH model or a Student t copula model.


In this work, there are two objectives that might be extended: The sampler and the bivariate mixture copula model. The sampler could be extended  to allow for a broader class of models. For example we might consider autoregressive processes of higher order in the state equation. In this case we can still rely on elliptical slice sampling and on an interweaving strategy.  Another extension could relax the assumption of a Gaussian dependence structure in the state equation by replacing the  autoregressive process by  a D-vine copula model. In this case elliptical slice sampling can no longer be applied to sample the latent states and an alternative sampling method is required.

The bivariate dynamic mixture copula could serve as a building block for regular vine copula models. Thus we could extend the bivariate model to arbitrary dimensions. This is interesting if we study not only the bivariate  volatility return relationship, but for example the dependence structure among several exchange rates. 







spbasic    






§ APPENDIX A. DETAILS ON THE SAMPLING PROCEDURE




 §.§ Sampling of the latent states in the sufficient augmentation


Here we derive μ_B_i| and Σ_B_i|. By the conditional independence assumptions of the AR(1) process and the way we defined the blocks B_1, …, B_m we obtain 

    f(s_B_i|s_0, s_-B_i, μ, ϕ, σ)    = f(s_B_i|  s_a_i-1,s_b_i+1, μ, ϕ, σ),  for  i = 1, …, m-1,  and 
    
    f(s_B_m|s_0, s_-B_m, μ, ϕ, σ)    = f(s_B_m|  s_a_m-1, μ, ϕ, σ).

Conditional on μ, ϕ and σ, the vector (s_a_i-1, s_B_i, s_b_i+1) is multivariate normal distributed with mean vector μ^AR_(a_i-1, B_i, b_i+1)∈ℝ^c_i+2 and covariance matrix Σ^AR_(a_i-1, B_i, b_i+1);(a_i-1, B_i, b_i+1)∈ℝ^(c_i+2)×(c_i+2), where c_i is the cardinality of 
 B_i.
Thus the vector
s_B_i|s_a_i-1, s_b_i+1, μ, ϕ, σ follows a multivariate normal distribution with mean vector μ_B_i| and covariance matrix Σ_B_i| given by

    μ_B_i|   = μ^AR_B_i + Σ^AR_B_i;(a_i-1,b_i+1)1-ϕ^2/(1-ϕ^2(c_i+1))σ^2( [        1 -ϕ^c_i+1; -ϕ^c_i+1        1;          ])
    ( [ s_a_i-1 - μ; s_b_i+1 - μ;             ]), 
    Σ_B_i|   =Σ^AR_B_i;B_i -Σ^AR_B_i;(a_i-1,b_i+1)1-ϕ^2/(1-ϕ^2(c_i+1))σ^2( [        1 -ϕ^c_i+1; -ϕ^c_i+1        1;          ])
    Σ^AR_(a_i-1,b_i+1);B_i.


The vector s_B_m|s_a_m-1,  μ, ϕ, σ corresponding to the last block is multivariate normal distributed with mean vector μ_B_m| and covariance matrix Σ_B_m| obtained as

    μ_B_m|   = μ^AR_B_m + Σ^AR_B_m;a_m-11-ϕ^2/σ^2 (s_a_m-1 - μ), 
    Σ_B_m|   =Σ^AR_B_m;B_m -Σ^AR_B_m;a_m-11-ϕ^2/σ^2Σ^AR_a_m-1;B_m.


We need to sample from N( 0, Σ_B_i|) several times during elliptical slice sampling. Instead of working with the c_i × c_i covariance matrix Σ_B_i| we can more efficiently sample from the c_i dimensional normal distribution  by using the conditional independence assumptions of the AR(1) process. It holds that


    f(s'_B_i|  h_a_i-1,h_b_i+1, μ, ϕ, σ)    = 
    ∏_t=0^s_i-1
    f(h'_a_i+t|  h_a_i-1:a_i+t-1  ,h_b_i+1, μ, ϕ, σ)
       = 
    ∏_t=0^s_i-1
    f(h'_a_i+t|  h_a_i+t-1,h_b_i+1, μ, ϕ, σ),

where  f(h'_a_i+t|  h_a_i+t-1,h_b_i+1, μ, ϕ, σ) is the univariate normal density with mean

    μ'_a_i+t|a_i+t-1,b_i+1 = 1/1-ϕ^2(s_i+1-t))( (ϕ-ϕ^2s_i+1-2t) (h_a_i+t-1-μ) + (ϕ^s_i-t - ϕ^s_i+2-t (h_b_i+1-μ)) ),

and variance

    σ^2_a_i+t|a_i+t-1,b_i+1 = σ^2/1-ϕ^2(1 - 1/1-ϕ^2(s_i+1-t)(ϕ^2 -2ϕ^2(s_i-t+1)+ϕ^2(s_i-t))).


So we can sample h'_B_i = (h'_a_i+t)_t=0, …, s_i-1 conditioned on h_a_i-1,h_b_i+1, μ, ϕ, σ recursively by

    ϵ_t ∼ N(0,1) 
       h'_a_i+t = μ'_a_i+t|a_i+t-1,b_i+1 + √(σ^2_a_i+t|a_i+t-1,b_i+1)ϵ_t,

for t = 0, …, s_i-1.





    f(s_B_i|  s_a_i-1,s_b_i+1, μ, ϕ, σ)    = 
    ∏_t=0^c_i-1
    f(s_a_i+t|  s_a_i-1:a_i+t-1  ,s_b_i+1, μ, ϕ, σ)
       = 
    ∏_t=0^c_i-1
    f(s_a_i+t|  s_a_i+t-1,s_b_i+1, μ, ϕ, σ),

where  f(s_a_i+t|  s_a_i+t-1,s_b_i+1, μ, ϕ, σ) is the univariate normal density with mean

    μ_a_i+t|a_i+t-1,b_i+1 = 1/1-ϕ^2(c_i+1-t))( (ϕ-ϕ^2c_i+1-2t) (s_a_i+t-1-μ) + (ϕ^c_i-t - ϕ^c_i+2-t (s_b_i+1-μ)) ),

and variance

    σ^2_a_i+t|a_i+t-1,b_i+1 = σ^2/1-ϕ^2(1 - 1/1-ϕ^2(c_i+1-t)(ϕ^2 -2ϕ^2(c_i-t+1)+ϕ^2(c_i-t))).


So we can sample s_B_i = (s_a_i+t)_t=0, …, s_i-1 conditioned on s_a_i-1,s_b_i+1, μ, ϕ, σ recursively by

    s_a_i+t∼ N(μ_a_i+t|a_i+t-1,b_i+1, σ^2_a_i+t|a_i+t-1,b_i+1),

for t = 0, …, c_i-1 and then s_B_i - μ_B_i| is a sample from N( 0,Σ_B_i|). 








 §.§ Sampling of the constant parameters in the ancillary augmentation


To sample μ, ϕ and σ in (AA) we deploy an adaptive random walk Metropolis-Hastings scheme as suggested by <cit.>, where tuning parameters are selected automatically using the Robbins Monro process (<cit.>). For sampling, it is convenient to move to unconstrained parameter spaces which is achieved by the following transformations

    ψ = ln(σ),   ξ = F_Z(ϕ).

Here F_Z(x) = 1/2log(1+x/1-x) is Fisher's Z transformation. The from (<ref>) implied log prior densities for ξ and ψ are given by

    ln(π(ξ))    = (a_ϕ-1)ln(F_Z^-1(ξ)+1)+(b_ϕ-1) ln(1-F_Z^-1(ξ))+ln(1-(F_Z^-1(ξ))^2) + c_1
    ln(π(ψ))    = -ψ - 1/2B_σexp(2ψ) + 2ψ + c_2.
    ,


where c_1 ∈ℝ and c_2 ∈ℝ are constants. The log posterior density in (AA) is obtained as

    lp_(AA)(μ,    ξ, ψ, s_0, s̃_1:T|Y) = ∑_t=1^T ln(f(y_t|s_t(s̃_1:T, μ, F_Z^-1(ξ), exp(ψ))))- 1/2∑_t=1^T s̃_t^2 
       + ln(φ(s_0|μ, exp(ψ)^2/1-F_Z^-1(ξ)^2))+ln(π(μ)) + ln(π(ξ)) + ln(π(ψ)) + c_3,

where c_3 ∈ℝ is a constant. We sample (μ, ϕ, σ) in two blocks, one block for μ and one block for (ϕ, σ).



  §.§.§ Update for μ

To sample the mean parameter μ from its full conditional we propose a new state μ_prop in the r-th iteration of the MCMC procedure by

    μ_prop∼ N(μ_cur, σ_MH,μ^r-1),

where μ_cur is the current value for μ. The proposal μ_prop is accepted with probability

    R = exp(lp_(AA)(μ_prop, ξ, ψ, s_0, s̃_1:T|Y)-lp_(AA)(μ_cur, ξ, ψ, s_0, s̃_1:T|Y))

and the scaling parameter σ_MH,μ^r is updated according to <cit.> by

    ln(σ_MH,μ^r) = ln(σ_MH,μ^r-1) + 4.058(R - 0.44)/r-1.

The scaling parameter is reduced, if the acceptance probability is larger than 0.44 and increased if the acceptance probability is smaller than 0.44. We target an average acceptance probability of 0.44, as recommended by <cit.> for univariate random walk Metropolis-Hastings. The constant 4.058 controls the step size and is chosen as suggested by <cit.>. 



  §.§.§ Joint update for ϕ and σ


In the r-th iteration, a two dimensional proposal (ξ_prop, ψ_prop) for (ξ_, ψ_)   is obtained by

    (ξ_prop, ψ_prop)^⊤∼ N_2( (ξ_cur, ψ_cur)^⊤, Σ_MH,ξ,ψ^r-1 ),

where (ξ_cur, ψ_cur) are the current values. The proposal is accepted with probability

    R = exp(lp_(AA)(μ, ξ_prop, ψ_prop, s_0, s̃_1:T|Y)-lp_(AA)(μ, ξ_cur, ψ_cur, s_0, s̃_1:T|Y)).


For adapting the covariance matrix we follow a suggestion of <cit.>. Let I_n denote the n-dimensional identity matrix. We set Σ_MH,ξ,s^r = I_2 if r<100 and 

    Σ_MH,ξ,s^r = (σ_MH,ξ,s^r)^2(Σ̂^r + (σ_MH,ξ,s^r)^2/r I_2)          if  r ≥ 100.

Here Σ̂^r is the empirical covariance matrix of (ξ^i,ψ^i)_i=1, … r, the first r samples for (ξ,ψ), and


    ln(σ_MH,ξ,s^r) = ln(σ_MH,ξ,s^r-1) + 6.534(R - 0.234)/r-1.
 

The matrix Σ̂^r + (σ_MH,ξ,ψ^r)^2/r I_2 is a positive definite estimate of the covariance matrix. This covariance estimate is scaled by (σ_MH,ξ,ψ^r)^2 to obtain the covariance matrix for the proposal in the next iteration. The scaling (σ_MH,ξ,ψ^r)^2 is tuned to achieve an average acceptance probability of 0.234 as suggested by <cit.> for multivariate random walk Metropolis-Hastings.
To reduce computational cost the empirical covariance matrix Σ̂^r can be updated in every step by the following recursion (see e.g. <cit.>)

    Σ̂^r = r-2/r-1Σ̂^r-1 + 1/r ((ξ^r,ψ^r)^⊤ - μ̂^r-1)((ξ^r,ψ^r)^⊤ - μ̂^r-1)^⊤,

where μ̂^r-1 is the sample mean of (ξ^i, ψ^i)_i=1, …, r-1. We also update the sample mean recursively by

    μ̂^r=1/r((r-1)μ̂^r-1 + (ξ^r,ψ^r)^⊤).



We have seen that the adaptions for the μ and the (ϕ,σ) updates tend to be very small after burn-in and therefore we only adapt during the burn-in period. This also ensures a correct sampling procedure without the need to verify the validity of an adaptive MCMC scheme.



§ APPENDIX B. THE STANDARDIZED SKEW STUDENT T DISTRIBUTION



According to <cit.>, the density of the univariate skew Student t distribution with parameters ξ∈ℝ, ω∈ (0,∞), α∈ℝ and df ∈ (0,∞)   is given by

    st(x|ξ, ω, α, df )= 2/ωt(x |df)T(αx-ξ/ω√(df +1/(x-ξ/ω)^2+df)|df+1),

where t(·|df) is the density function of the univariate Student t distribution with df degrees of freedom and T(·|df) the corresponding distribution function.
The expectation and variance of a random variable X following a skew Student t distribution with parameters ξ, ω, α as above and df > 2 are given by 

    E(X) = ξ +  ω b_dfδ,  and  Var(X) = ω^2(df/df -2  - b_df^2 δ^2  ),

where δ = α^2/√(1+ α^2) and b_df  = √(df/π)Γ(df-1/2)/Γ(df/2).
If we set

    ω   = √(1/(df/df -2  - b_df^2 δ^2  ))  and  ξ = - ω b_dfδ = - √(1/(df/df -2  - b_df^2 δ^2  )) b_dfδ,

only the parameters α and df remain unknown and the random variable has zero mean and a variance of one. We refer to the corresponding distribution as the standardized skew Student t distribution. Its density is denoted by sst and is obtained as

    sst(x|α, df) = st(x |  - √(1/(df/df -2  - b_df^2 δ^2  )) b_dfδ, √(1/(df/df -2  - b_df^2 δ^2  )), α, df).







§ APPENDIX C. ADDITIONAL MATERIAL FOR THE BIVARIATE DYNAMIC MIXTURE COPULA (SECTION <REF>)
































§ APPENDIX D. ADDITIONAL MATERIAL FOR THE APPLICATION (SECTION <REF>)






 §.§ Daily log returns






 §.§ Posterior statistics









































 §.§ Calculating the log predictive score

We describe in detail how we proceed for model ℳ_dyn^mix. We consider T+K observations of dimension two, stored in the data matrix Y_1:(T+K);1:2, where the first T observations are used to train the model and the last K are used for evaluation. 



Step 1: (Model fitting based on the training period)


  * We fit two marginal skew Student t stochastic volatility models to y_1:T;1 and y_1:T;2. This yields R_train draws of the parameters denoted by s_1:T;j^st,r, μ_j^st,r, ϕ_j^st,r, σ_j^st,r, α_j^st,r and df_j^st,r, r=1, …, R_train and corresponding posterior mode estimates  ŝ_1:T;j^st, μ̂_j^st, ϕ̂_j^st, σ̂_j^st, α̂_j^st and d̂f̂_j^st for j=1,2.
 

  * We estimate the copula data 

    û_tj = ssT(y_tjexp(-ŝ_tj^st/2) | α̂_j^st, d̂f_j^st)

 for t=1, …, T, j=1,2.



  * We fit the dynamic bivariate mixture copula model introduced in (<ref>) based on the pseudo copula data Û_1:T;1:2 and obtain posterior draws s_1:T^cop,r, μ^cop,r, ϕ^cop,r, σ^cop,r, ν^cop,r, p^cop,r for r=1, …, R_train and corresponding posterior mode estimates ŝ_1:T^cop, μ̂^cop, ϕ̂^cop, σ̂^cop, ν̂^cop, p̂^cop.





Step 2: (The one-day ahead predictive density at time point T+1)


  * For j=1,2,we obtain an estimate for the log variance at time point T+1, ĥ_T+1;j^st, as the posterior mode estimate of h^st,r_T+1j = μ_j^st,r + ϕ_j^st,r (h^st,r_Tj - μ_j^st,r) , r=1, …, R_train.  


  * We obtain an estimate for Fisher's Z transform of Kendall's τ at time point T+1, ĥ_T+1^cop, as the posterior mode estimate of h^cop,r_T+1 = μ^cop,r + ϕ^cop,r ( h^cop,r_T - μ^cop,r) , r=1, …, R_train.


  * The one-day ahead predictive density at time point T+1 evaluated at (z_1,z_2) is given by 

    f_T+1^p(z_1,z_2)=c_T+1^p(z_1,z_2) g_T+1^p(z_1,z_2),

where

    c_T+1^p(z_1,z_2)=c^M(ssT(x_1 |α̂_1^st, d̂f̂_1^st), ssT(x_2|α̂_2^st, d̂f̂_2^st);F_z^-1(ĥ_T+1^cop), ν̂^cop, p̂^cop)

and

    g_T+1^p(z_1,z_2)=sst(x_1 |α̂_1^st, d̂f̂_1^st) sst(x_2 |α̂_2^st, d̂f̂_2^st) exp(-ĥ_T+1;1^st/2) exp(-ĥ_T+1;2^st/2)

with x_j = z_jexp(-ĥ_T+1;j^st/2).



Step 2: (The one-day ahead predictive density)



Estimating the one-day ahead predictive density at time T+k, 1≤ k≤ K would usually require to fit daily models with observations up to time T+k-1 for k=1, …, K. In order to save computational resources we use another approach where we only update the dynamic parameters, i.e. the log variances and Kendall's τ. For the constant parameters we use the estimates from the training period 1, …, T. In this case we found that it is enough to only consider a time horizon of 100 time points, i.e. to estimate a dynamic parameter at time T+k we consider data in the period T+k-100, …, T+k-1.  We proceed as follows  to obtain the one-day ahead predictive density at time point T+k with 1≤ k≤ K.


  * We consider a skew Student t stochastic volatility model as in (<ref>), where we keep the parameters μ, ϕ, σ, α and df fixed and only update the latent log variances. Therefore we draw the latent log variances  s_(T+k-100):(T+k-1);j^st conditional on y_(T+k-100):(T+k-1);j, μ̂_j^st, ϕ̂_j^st, σ̂_j^st, d̂f̂_j^st and α̂_j^st for j=1,2. We denote the draws by s^st,r_(T+k-100):(T+k-1);j, r=1, …, R_test for j=1,2. Corresponding posterior mode estimates are denoted by ŝ_(T+k-100):(T+k-1);j^st, j=1,2.

  * We estimate the copula data via the probability integral transform, i.e.
for j=1,2 and t = T+k-100, … T+k-1 we calculate

    û_tj=ssT(y_tjexp(-ŝ_tj^st/2)  | α̂_j^st, d̂f̂_j^st).


  * We fit the dynamic mixture copula model to the pseudo copula data Û_(T+k-100):(T+k-1);(1:2) where we keep the constant parameters fixed. We only update s_(T+k-100):(T+k-1)^cop conditional on 
Û_(T+k-100):(T+k-1);(1:2), μ̂^cop, ϕ̂^cop, σ̂^cop, ν̂^cop, p̂^cop. The corresponding draws are denoted by s^cop,r_(T+k-100):(T+k-1), r=1, …, R_test and the posterior mode estimates by ŝ_(T+k-100):(T+k-1)^cop.



  * For j=1,2, we obtain an estimate for the log variance at time point T+k as ŝ_T+k;j^st = μ̂_j^st + ϕ̂_j^st (ŝ^st_T+k-1;j - μ̂_j^st).  


  * We obtain an estimate for Fisher's Z transform of Kendall's τ at time point T+k,  as  ŝ_T+k^cop = μ̂^cop + ϕ̂^cop (ŝ_T+k-1^cop - μ̂^cop).  


  * The predictive density evaluated at (y_1,y_2) is given by 

    f_T+k^p(y_1,y_2)=c_T+k^p(y_1,y_2) g_T+k^p(y_1,y_2),

with

    c_T+k^p(y_1,y_2)=c^M(ssT(x_1 |α̂_1^st, d̂f̂_1^st), ssT(x_2|α̂_2^st, d̂f̂_2^st);F_Z^-1(ŝ_T+k^cop), ν̂^cop, p̂^cop),

where c^M is the density of the mixture copula defined in (<ref>) and

    g_T+k^p(y_1,y_2)=sst(x_1  |α̂_1^st, d̂f̂_1^st) sst(x_2 |α̂_2^st, d̂f̂_2^st) exp(-ŝ_T+k;1^st/2) exp(-ŝ_T+k;2^st/2),

with x_j = y_jexp(-ŝ_T+kj^st/2) for j=1,2. 
 

Step 3: (The cumulative pseudo log predictive score)



The cumulative pseudo log predictive score is obtained as

    LP    = ∑_k=1^Klog(f^p_T+k(y_T+k;1,y_T+k;2)).



During the training period we run R_train = 31000 iterations with a burn-in of 1000, while for updating only the dynamic parameters 11000 iterations with a burn-in of 1000 is enough, i.e. we use R_test=11000.






Supplementary material






§ ELLIPTICAL SLICE SAMPLING


Elliptical slice sampling was proposed by <cit.>. It was developed  for models where dependencies are generated through a latent multivariate normal distribution.
 
We assume that the posterior density for a parameter vector θ given data D is proportional to

    f(θ|D) ∝ℓ(θ|D) φ(θ |  0 , Σ),
 
where ℓ(θ|D) is the likelihood function and φ(θ |  0 , Σ) is the multivariate normal density with zero mean and covariance matrix Σ. 
<cit.> consider the Metropolis-Hastings sampler of <cit.> where a proposal θ' is obtained from the following stochastic representation

    θ' = √(1 - α^2)θ + α v,  v ∼ N( 0,Σ).

Here α∈ [-1,1] is a fixed step size parameter. The proposal is accepted with probability 

    min(1,ℓ(θ')/ℓ(θ)).

Elliptical slice sampling adapts the step size parameter α during sampling. This eliminates the need to select the parameter before sampling and it may be a better approach for situations where good choices of the step size parameter depend on the region of the state space. 
<cit.> first suggest alternatively to propose a new state by

    θ' = cos(ω)θ  + sin(ω)  v,  v ∼ N( 0,Σ).

Here the angle ω corresponds to the step size. As we move ω towards zero the proposal gets closer to the initial value θ.  
<cit.> argue that (<ref>) provides a more flexible choice for the proposals compared to (<ref>), if the parameter ω is also updated, which is here the case. In elliptical slice sampling we first draw an angle ω from the uniform distribution on [0,2π] and obtain a proposal as outlined in (<ref>). This proposal is accepted according to (<ref>). If the proposal is not accepted a new angle is selected with a slice sampling approach (<cit.>) such that the angle approaches zero as more samples are rejected. This ensures that at some point the proposal will be accepted. The approach is outlined in  Algorithm <ref>. 
<cit.> show that this samples from a Markov chain, where (<ref>) is the corresponding stationary distribution.








§ THE DCC-GARCH MODEL

The DCC-GARCH model was introduced by <cit.>. Mathematical propoerties were developed by <cit.>. In the DCC-GARCH model  it is assumed that a d dimensional random vector ϵ_t∈ℝ^d at time t is multivariate normal distributed with zero mean and dynamic covariance matrix H_t ∈ℝ^d × d, i.e. 

    ϵ_t|H_t ∼ N( 0,H_t),

for t=1, …, T. The covariance matrix can be written as 

    H_t = D_tC_tD_t,

where  C_t is a d × d correlation matrix and D_t is a d × d diagonal matrix. The diagonal matrix D_t contains the marginal standard deviations. We assume GARCH(P_j,Q_j) innovations for the j-th diagonal entry of D_t, i.e.

    d_t,j^2 = ω_j  + ∑_p=1^P_jα_pjϵ_t-p,j^2 +  ∑_q=1^Q_jβ_qj d_t-q,j^2,

where


  * ω_j > 0,

  * d_t,0>0,

  * ∑_p=1^P_iα_jp + ∑_q=1^Q_iβ_qj<1 and the roots of 1 - ∑_p=1^P_iα_jpZ^p + ∑_q=1^Q_iβ_qjZ^q lie outside the unit circle,

  * α_jp for all p ∈{1, … P_j} and β_jp for all q ∈{1, … Q_j}  are such that d_tj^2 is positive. 
 
  The correlation matrix is decomposed as follows

    C_t = diag(L_t)^-1/2L_tdiag(L_t)^-1/2,

where L_t is a positive definite matrix. Further we denote by
 
    r_t = D_t^-1ϵ_t

 the standardized return and obtain L̅ as

    L̅ = 1/T∑_t=1^T r_tr_t^⊤.
 
 For L_t we assume the following dynamic structure

    L_t = (1 - ∑_m=1^M a_m - ∑_n=1^Nb_n) L̅ + ∑_m=1^M a_m r_t-mr_t-m^⊤ + ∑_n=1^N b_n L_t-n,

where 


  * α_m ≥ 0 for all m ∈{1, … M} and b_n ≥ 0 for all n ∈{1, …, N},

  * ∑_m=1^M a_m + ∑_n=1^N b_n < 1,

  * L_0 is positive definite.

<cit.> show that under the above conditions H_t is a proper covariance matrix. The DCC(M,N)-GARCH(P,Q) model is obtained by setting P_j = P and Q_j=Q for j = 1, …, d.












§ FURTHER RESULTS FOR THE APPLICATION



 §.§ Results for the marginal models











 §.§ Results for the dependence models













