Constant Correlation

In the Constant Conditional Correlation (CCC) model of Bollerslev (1990) the equation for the conditional covariance can be represented as:

where $D_t = diag(\sqrt {h_{11,t}} ,...,\sqrt{h_{nn,t}})$, and R is the positive definite constant correlation matrix.

The conditional variances, h_ii, t, which can be estimated separately, usually admit some type of GARCH(p,q) model, but can be generalized to any admissible model for the variance dynamics.

The conditions for the positivity of the covariance matrix H_t are that R is positive definite and the diagonal elements of D_t, h_ii, t, are positive.

Formally, the model may be represented as follow:

where y_t is a vector representing the series values at time t, μ_t a vector of the conditional mean at time t for each series, ε_t a vector of zero mean residuals which have a multivariate Normal distribution with conditional covariance Σ_t, and LL_t(θ) is the log-likelihood at time t given a parameter vector θ of univariate GARCH (or any other model) parameters. The constant correlation matrix R is usually approximated by its sample counterpart as

where z_t is the vector of GARCH standardized residuals at time t.

Code Snippet

dcc_modelspec(..., dynamics = 'constant')

Dynamic Conditional Correlation

Allowing the correlation to be time varying, whilst retaining the 2-stage separation of dynamics, creates an additional challenge to the positive definiteness constraint at each point in time for the conditional correlation. R. F. Engle (2002) proposed the following proxy process :

where a_j and b_j are non-negative coefficients with the constraint that they sum to less than 1 for stationarity and positive definiteness, and Q̄ is the unconditional matrix of the standardized residuals z_t. To recover R_t, the following re-scaling operation is performed:

Under the assumption that ε_t ∼ N(0, Σ_t), the log-likelihood is:

where LL_V, t(ϕ) is the log-likelihood of the conditional variances given the GARCH parameters ϕ, and LL_R, t(ψ|ϕ) is the log-likelihood of the conditional correlation given the DCC parameters ψ|ϕ (conditioned on the GARCH parameters).

Code Snippet

dcc_modelspec(..., dynamics = 'dcc')

Estimation

Estimation of the CCC and DCC models is done in 2 steps. The first step is to estimate the conditional variances for each series using a univariate GARCH model. The univariate models do not need to conform to the any single type of dynamics, but can be any admissible model. The conditional variances are then used to standardize the residuals which are then used to estimate the conditional correlation matrix. Since the first step is separable, it can be done in parallel for each series. After considerable experimentation, it was chosen to avoid the use of TMB for automatic differentiation and instead use numerical derivatives for the second stage estimation due to speed considerations, and observing that there was little gain in terms of accuracy.

Code Snippet

dcc_modelspec(...) |> estimate(...)

Standard Errors

Since estimation is performed in 2 stages, the log-likelihood is the sum of these 2 stages. More specifically, the second stage log-likelihood is estimated conditioning on the parameters of the first stage. Details of the assumptions behind this and the arguments for the asymptotic distribution of the estimates are put forth in R. F. Engle and Sheppard (2001), who posit that the asymptotic distribution of the 2-stage DCC model estimates with parameters θ = (ϕ, ψ) is:

where:

and

where B₀⁻¹ and M₀ are the Hessian and the covariance matrix of the score vector respectively, also called the bread and meat of the sandwich estimator. In the case of this 2 stage estimation procedure, the Hessian is partitioned into 4 blocks, B₁₁, B₁₂, B₂₂ and M₁₁, M₁₂, M₂₂, where B₁₁ and M₁₁ are the Hessian and covariance matrix of the first stage log-likelihood, and B₂₂ and M₂₂ are the Hessian and covariance matrix of the second stage log-likelihood. B₁₂ and M₁₂ are the cross derivatives of the first and second stage log-likelihoods.

Because the Hessian is partitioned, the calculation of the standard errors is more computationally intensive than the standard OPG (outer product of gradients) method. The package defaults to the OPG method, but the user can choose to calculate the full partitioned Hessian if they so wish, for which parallel functionality has been implemented to provide some speed gains. The package uses numerical derivatives for the second stage estimation, as well as the calculation of the standard errors. When the partitioned Hessian is present, a number of additional estimators become available including the those of the QML and HAC (using the bandwidth of Newey and West (1994)) type.

Code Snippet

dcc_modelspec(...) |> estimate(...) |> vcov(...)

Distributions

The original DCC model was based on the assumption that the residuals are Normally distributed, as shown in equations and .

However, this is not always the case in practice. There are 2 possible approaches to this problem. The first is to use a different multivariate distribution for the DCC model and the second is to use a copula approach.

dcc_modelspec(distribution = 'mvt')

cgarch_modelspec(...)

Generalized Dynamics

A number of extensions to the DCC model have been proposed in the literature. These include the use of asymmetric and generalized dynamics. Cappiello, Engle, and Sheppard (2006) generalize the DCC model with the introduction of the Asymmetric Generalized DCC (AGDCC) where the dynamics of Q_t now become:

where A, B and G are the N × N parameter matrices, z̄_t are the zero-threshold standardized residuals (equal to z_t when less than zero else zero otherwise), Q̄ and N̄ are the unconditional matrices of z_t and z̄_t respectively. The generalization allows for cross dynamics between the standardized residuals, and cross-asymmetric dynamics. However, this comes at the cost of increased computational complexity making it less feasible for large scale systems. Additionally, covariance targeting in such a high dimensional model where the parameters are no longer scalars creates difficulties in imposing positive definiteness constraints. The tsmarch package does not implement this generalized model, but instead a scalar version to capture asymmetric dynamics (adcc):

Code Snippet

cgarch_modelspec(..., dynamics = 'adcc')

Forecasting and Simulation

Because of the non linearity of the DCC evolution process, the multi-step ahead forecast of the correlation cannot be directly solved. R. F. Engle and Sheppard (2001) provide a couple of suggestions which form a reasonable approximation. In the tsmarch package we instead use a simulation approach, similar to the other packages in the tsmodels universe, and in doing so generate a simulated distribution of the h-step ahead correlation and covariance. This also generalizes to the copula model which does not have a closed form solution. Formally, for a draw s at time t, the simulated correlation is given by:

where ϵ^s_t is a vector of standard Normal random variables, Λ is the diagonal matrix of the eigenvalues of R^s_t and E is the matrix of eigenvectors of R^s_t, such that:

For the multivariate Student distribution, the draws for z_t^s are generated as :

where W^s is a random draw from a chi-squared distribution with ν degrees of freedom. For ϵ^s_t, it is also possible to instead use the empirical standardized and de-correlated residuals, an option which is offered in the package.

Code Snippet

dcc_modelspec(...) |> estimate(...) |> predict(...)
dcc_modelspec(...) |> estimate(...) |> simulate(...)

Weighted Distribution

Since all predictions are generated through simulation, the weighted predictive distribution $\bar {y_t}^s$ can be obtained as:

where s is a sample draw from the simulated predictive distribution of the n series, and w_i, t the weight on series i at time t.

Code Snippet

dcc_modelspec(...) |> estimate(...) |> tsaggregate(...)

News Impact Surface

Robert F. Engle and Ng (1993) introduced the news impact curve as a tool for evaluating the effects of news on the conditional variances, and further extended by Kroner and Ng (1998) to the multivariate GARCH case. The news impact surface shows the effect of shocks from a pair of assets on their conditional covariance (or correlation). Formally, the correlation R_i, j given a set of shock (z_i, z_j) is given by:

where z̄, Q̄ and N̄ are defined as in equation .

Code Snippet

dcc_modelspec(...) |> estimate(...) |> newsimpact(...) |> plot(...)

Critique

The DCC model has come under a considerable amount of criticism in the literature, particularly by Caporin and McAleer (2008) who argue that consistency and asymptotic Normality of the estimators has not been fully established, and proposes instead the indirect DCC model (essentially a scalar BEKK) which models directly the conditional covariance matrix and forgoes univariate GARCH dynamics.

Aielli (2013) also points out that the estimation of Q_t as the empirical counterpart of the correlation matrix of z_t is inconsistent since E[z_tz′_t] = E[R_t] ≠ E[Q_t]. They propose a DCC (cDCC) model which includes a corrective step which eliminates this inconsistency, albeit at the cost of targeting which is not allowed.

Estimation

The estimation of the GOGARCH model is performed in 3 steps.

1. Filter the data for conditional mean dynamics (μ_i, t). Similar to the DCC model, this is performed outside the system.
2. Perform Independent Component Analysis (ICA) on the residuals to separate out the additive and statistically independent components.
3. Estimate the conditional distribution of each independent component using an admissible univariate model. In the tsmarch package we use the tsgarch package for estimating the GARCH model for each factor with either a Generalized Hyperbolic (gh) or Normal Inverse Gaussian (nig) distribution

The log-likelihood of the GOGARCH model is composed of the individual log-likelihoods of the univariate independent factor GARCH models and a term for the mixing matrix:

where GH_λi is the Generalized Hyperbolic distribution with shape parameter λ_i for the i-th factor. Since the inverse of the mixing matrix is the un-mixing matrix A⁻¹ = W = K′U, and since the determinant of U is 1 (|U| = 1), the the log-likelihood simplifies to:

In the presence of dimensionality reduction, K is no longer square and we cannot directly compute its determinant. Instead, we use the determinant of KK′ , to reflect the combined effect of whitening and dimensionality reduction:

where m are the number of factors. A key advantage of this approach is the ability to estimate multivariate dynamics using only univariate methods, providing for feasible estimation of large scale systems. The speed advantage of this approach does not come for free, as there is a corresponding need for memory to handle the parallelization offered, particularly when generating higher co-moment matrices.

Code Snippet

gogarch_modelspec(...) |> estimate(...)

Co-Moments

The linear affine representation of the GOGARCH model provides for a closed-form expression for the conditional co-moments. The first 3 co-moments are given by:

where M_f, t², M_f, t³ and M_f, t⁴ are the (N × N) conditional covariance, (N × N²) the conditional third co-moment matrix and the (N × N³) conditional fourth co-moment matrix of the factors, respectively, represented as unfolded tensors

These are calculated from the univariate factor dynamics, and defined as:

where M_k, f, t³, k = 1, …, n and M_{k, l, f, t}⁴, k, l = 1, …, n are the N × N sub-matrices of M_f, t³ and M_f, t⁴, respectively, with elements

Since the factors f_i, t can be decomposed as $z_{i,t}\sqrt{h_{i,t}}$, and given the assumptions on z_i, t then E[f_i, tf_j, tf_k, t|𝔉_t − 1] = 0. It is also true that for i ≠ j ≠ k ≠ l E[f_i, tf_j, tf_k, tf_l, t|𝔉_t − 1] = 0 and when i = j and k = l,

Thus, under the assumption of mutual independence, all elements in the conditional co-moments matrices with at least 3 different indices are zero. Finally, standardizing the conditional co-moments one obtains conditional co-skewness and co-kurtosis of y_t,

where S_{i, j, k, t} represents the series co-skewness between elements i, j, k of y_t, σ_i, t the standard deviation of y_i, t, and in the case of i = j = k represents the skewness of series i at time t, and similarly for the co-kurtosis tensor K_{i, j, k, l, t}. Two natural applications of return co-moments matrices are Taylor type utility expansions in portfolio allocation and higher moment news impact surfaces.

Code Snippet

tscov(...)
tscor(...)
tscoskew(...)
tscoskurt(...)

The weighted moments are calculated as follows:

where w_t is the vector of weights at time t. Using the vectorization operator vec the weighted moments can be re-written compactly as:

One application of this is in the calculation of the Constant Absolute Risk Aversion utility function, using a Taylor series expansion as in Ghalanos, Rossi, and Urga (2015):

where η represents the investor’s constant absolute risk aversion, with higher (lower) values representing higher (lower) aversion, and W̄_t = w′_tE_t − 1[y_t] is the expected portfolio return at time t. Note that this approximation includes two moments more than the mean-variance criterion, with the skewness and kurtosis are directly related to investor preferences (dislike) for odd (even) moments, under certain mild assumptions given in Scott and Horvath (1980).

Code Snippet

gogarch_modelspec(...) |> estimate(...) |> predict(...) |> tsaggregate(...)

Weighted Distribution

The weighted conditional density of the model may be obtained via the inversion of the Generalized Hyperbolic characteristic function through the fast fourier transform (FFT) method as in Chen, Härdle, and Spokoiny (2010) or by simulation as in Section . The tsmarch package generates both simulated predictive distributions as well as the weighted p^*, d^* and q^* functions based on FFT.

It should be noted that while the n-dimensional NIG distribution is closed under convolution, when the distributional higher moment parameters are allowed to be different across assets, as is the case here, this property no longer holds.

Provided that ϵ_t is a n-dimensional vector of innovations, marginally distributed as 1-dimensional standardized GH, the density of weighted asset returns, w_i, ty_i, t is:

where w̄_t = w′_tAH_t^1/2 and hence w̄_i, t is the i-th element of this vector, μ_i, t is the conditional mean of series i at time t, and the parameter vector (μ̄_i, t, δ_i, t, β_i, t, α_i, t, λ_i) represents the generalized hyperbolic distributional parameters which have been converted from their (ρ, ζ) parameterization used in the tsgarch package estimation routine. In order to obtain the weighted density we need to sum the individual log densities of ϵ_i, t:

where, γ = ᾱ_j² − β̄_j², υ = ᾱ_j² − (β̄_j + iu)², and (ᾱ_j, β̄_j, δ̄_j, μ̄_j) are the scaled versions of the parameters (α_i, β_i, δ_i, μ_i) as shown in equation . The density may be accurately approximated by FFT as follows:

Once the density is obtained, we use smooth approximations to obtain functions for the density (d^*), quantile (q^*) and distribution (p^*), with random number generation provided by the use of the inverse transform method using the quantile function.

Code Snippet

gogarch_modelspec(...) |> estimate(...) |> predict(...) |> tsconvolve(...) |> dfft(...)

Factor News Impact Surface

Consider the equation for the construction of the conditional covariance matrix Σ_t:

where H_t is the diagonal conditional covariance matrix of the independent components. The factor news impact covariance surface is a measure of the effect of shocks from the independent factors on the conditional covariance of the returns of assets i and j at time t − 1. In order to construct this we first generate the news impact curves for all the factors using a common grid of shock values, and then calculate the effect of these on the conditional covariance of the returns of assets i and j. Since the factors are independent, we simply loop through all pairwise combinations of these shocks to calculate the covariance surface, keeping the value of the other factors constant:

where :

1. a_k, i is the loading of factor k on asset i,
2. h_k, k|ϵ_k is the conditional variance of factor k from the news impact curve given the shock ϵ_k,
3. a_k, i ⋅ h_k, k|ϵ_k ⋅ a_k, j captures the contribution of factor k to the covariance between assets i and j, based on the shock ϵ_k
   applied to factor k,
4. a_l, i ⋅ h_l, l|ϵ_l ⋅ a_l, j captures the contribution of factor l to the covariance between assets i and j, based on the shock ϵ_l applied to factor l and
5. $\sum_{\substack{m=1 \\ m \neq k, l}}^{n} a_{m,i} \cdot s_{m,m}(0) \cdot a_{m,j}$ captures the contributions of all other factors (except k and l), where their variances are held constant at their unconditional levels (i.e., with no shocks).

This means that the covariance will be affected by both the relative sign and magnitude of the loadings and the relative magnitude of the variance of the factors.

The news impact surface is not restricted to the covariance matrix, but can be extended to the conditional co-skewness and co-kurtosis, but that is only truly meaningful when the higher moment distributional parameters are allowed to vary as in Ghalanos, Rossi, and Urga (2015).

Code Snippet

gogarch_modelspec(...) |> estimate(...) |> newsimpact(...) |> plot(...)

Correlation and Kendall’s τ

Pearson’s product moment correlation R totally characterizes the dependence structure in the multivariate Normal case, where zero correlation also implies independence, but can only characterize the ellipses of equal density when the distribution belongs to the elliptical class. In the latter case for instance, with a distribution such as the multivariate Student, the correlation cannot capture tail dependence determined by the shape parameter. Furthermore, it is not invariant under monotone transformations of original variables making it inadequate in many cases. An alternative measure which does not suffer from this is Kendall’s τ (see Kruskal (1958)), based on rank correlations which makes no assumption about the marginal distributions but depends only on the copula C. It is a pairwise measure of concordance calculated as:

For elliptical distributions, Lindskog, Mcneil, and Schmock (2003), proved that there is a one-to-one relationship between this measure and Pearson’s correlation coefficient ρ given by:

which under certain assumptions (e.g. in the multivariate Normal case) simplifies to $$\frac{2}{\pi} \arcsin \rho_{ij}$$. Kendall’s τ is invariant under monotone transformations making it more suitable when working with non-elliptical distributions. A useful application arises in the case of the multivariate Student distribution where the maximum likelihood approach for the estimation of the correlation matrix R become unfeasible for large dimensions. Instead, estimating the sample counterpart of Kendall’s τ from the transformed margins and then translating that into a correlation matrix as detailed in equation provides for a method of moments type estimator. The shape parameter ν is then estimated keeping the correlation matrix constant, with little loss in efficiency.

Notes on Selected Methods

The table below provides a summary of selected methods, the classes they belong to, the classes they output, whether they are parallelized and the dimensions of the objects involved. For example, the tsaggrate method on predicted and simulated objects will return a list of length 4, with each slot having an object of class tsmodel.distribution of size sim × h representing the predictive simulatated distribution of each of the weighted moments generated by the model. The output for the tsconvolve method which is list(h):list(sim) means that the output is a list of length h where each slot contains a list of length sim.