Multivariate Volatility

Overview

The cross-asset branch of the volatility arc: the earlier examples modelled one series' variance; here the object is the whole conditional covariance matrix $H_t$ — how assets co-move, and how that co-movement changes over time. The central risk-management fact is correlation is not constant — it rises in crises ("diversification meltdown"), exactly when you need diversification most. Four model families are built from scratch and cross-checked, each a different answer to "how do you make $H_t$ both realistic and tractable", then raced in a portfolio application.

Four models, one covariance matrix

Part 1 — CCC (Bollerslev 1990). Splits $H_t = D_t R\, D_t$ into univariate GARCH volatilities $D_t$ and a constant correlation $R$. The decomposition is the keeper — it keeps $H_t$ positive-definite for free and every later model inherits it — but "constant" is the flaw: the data emphatically reject a fixed $R$. CCC is precisely the null hypothesis of Part 2.

Part 2 — DCC (Engle 2002), the workhorse. Lets the correlation become time-varying $R_t$, driven by its own GARCH-like recursion. Its genius is parsimony: no matter how many assets, the entire correlation matrix's dynamics are governed by just two scalars $(a,b)$ — so it scales to hundreds of assets. It beats CCC by 1407 AIC for two extra parameters (a clean nested test). Fit three ways — Engle's two-step QMLE (from scratch), a fully Bayesian NumPyro version that gives an honest posterior over $(a,b)$ and the whole correlation path (the two-step QMLE reports no such uncertainty), and R's rmgarch.

Part 3 — BEKK (Engle & Kroner 1995). Models the covariance directly with a quadratic recursion that guarantees a positive-definite $H_t$. On two assets diagonal BEKK actually wins (AIC 37379 vs DCC's 37439), but it does not scale — the full BEKK's dynamic parameters grow as $2N^2$ (2265 at $N=30$) while DCC stays at two. BEKK for a handful of assets where the covariance structure is paramount; DCC when $N$ is large. From-scratch ML, a Bayesian diagonal BEKK (NumPyro), and R's mgarchBEKK.

Part 4 — Factor Stochastic Volatility (Aguilar & West 2000), the Bayesian capstone. Where MGARCH makes $H_t$ a deterministic function of past returns, factor SV routes co-movement through a low-rank latent factor whose volatilities are their own stochastic processes — delivering a positive-definite $H_t$ (like BEKK) with only $O(N)$ parameters (like DCC), and natively Bayesian. There is no two-step shortcut: MCMC estimates the loadings, the SV parameters, and the entire latent covariance path jointly (clean convergence, $\hat r\approx 1.02$, despite ~32,000 latent states — the heaviest fit in the example at ~44 min). Cross-checked against Kastner's factorstochvol.

Does it pay off?

The portfolio horse race. Every model estimates the same $H_t$ — does a better one build a better portfolio? Feeding each model's covariance into a minimum-variance portfolio out-of-sample, going from no model to a dynamic covariance cuts realized portfolio volatility by ~28% (19.1% → 13.8%), most of it harvested in the 2008 and 2020 crises — with factor SV the lowest-risk. The through-line of the whole volatility arc: co-volatility is persistent, time-varying, and rises together in crises, and modelling that structure — rather than assuming it away — is what turns a covariance matrix into a better decision.

Notebooks

Downloads

Factor SV is fit in NumPyro (no from-scratch module); see the factorsv_* notebooks.

References

• Bollerslev, T. (1990). Modelling the coherence in short-run nominal exchange rates: a multivariate generalized ARCH model. Review of Economics and Statistics 72(3), 498–505. — the CCC model
• Engle, R. F. & Kroner, K. F. (1995). Multivariate simultaneous generalized ARCH. Econometric Theory 11(1), 122–150. — the BEKK model (positive-definite by construction)
• Engle, R. F. (2002). Dynamic conditional correlation: a simple class of multivariate GARCH models. Journal of Business & Economic Statistics 20(3), 339–350. — the DCC model and its two-step estimation
• Aguilar, O. & West, M. (2000). Bayesian dynamic factor models and portfolio allocation. Journal of Business & Economic Statistics 18(3), 338–357. — factor stochastic volatility for the covariance
• Kastner, G., Frühwirth-Schnatter, S. & Lopes, H. F. (2017). Efficient Bayesian inference for multivariate factor stochastic volatility models. Journal of Computational and Graphical Statistics 26(4), 905–917. — the factorstochvol R implementation
• Bauwens, L., Laurent, S. & Rombouts, J. V. K. (2006). Multivariate GARCH models: a survey. Journal of Applied Econometrics 21(1), 79–109. — the multivariate-GARCH landscape

DCC Module — Source Code

"""
Dynamic Conditional Correlation GARCH (Engle 2002) from scratch.

DCC relaxes CCC's constant correlation into a time-varying R_t driven by its own GARCH-like recursion.
Same volatility/correlation split H_t = D_t R_t D_t, but now, on the standardized residuals z_t = D_t^{-1} eps_t,

    Q_t = (1 - a - b) Qbar + a z_{t-1} z_{t-1}' + b Q_{t-1},        Qbar = uncond. corr of z,
    R_t = diag(Q_t)^{-1/2} Q_t diag(Q_t)^{-1/2}      (rescale Q_t to a proper correlation matrix).

Estimated in two steps (Engle 2002): stage 1 = a univariate GARCH per series (reused from ccc.garch11);
stage 2 = maximize the correlation part of the Gaussian log-likelihood over the two scalars (a, b),
    L_c(a,b) = -1/2 sum_t [ log|R_t| + z_t' R_t^{-1} z_t ].
CCC is the special case a = b = 0 (R_t == Qbar), so DCC strictly nests it.
"""
import numpy as np
from scipy.optimize import minimize
from ccc import garch11


def _Rt_path(Z, Qbar, a, b):
    """The full sequence of correlation matrices R_t (T x N x N)."""
    T, N = Z.shape; Q = Qbar.copy(); Rt = np.empty((T, N, N))
    for t in range(T):
        if t > 0:
            Q = (1 - a - b) * Qbar + a * np.outer(Z[t - 1], Z[t - 1]) + b * Q
        d = 1.0 / np.sqrt(np.diag(Q))
        Rt[t] = Q * np.outer(d, d)
    return Rt


def _corr_loglik(Z, Qbar, a, b):
    """Stage-2 correlation log-likelihood (the part that depends on a, b)."""
    T, N = Z.shape; Q = Qbar.copy(); ll = 0.0
    for t in range(T):
        if t > 0:
            Q = (1 - a - b) * Qbar + a * np.outer(Z[t - 1], Z[t - 1]) + b * Q
        d = 1.0 / np.sqrt(np.diag(Q)); Rt = Q * np.outer(d, d)
        _, ld = np.linalg.slogdet(Rt)
        ll += -0.5 * (ld + Z[t] @ np.linalg.solve(Rt, Z[t]))
    return ll


def _full_loglik(Rret, H, Rt, mu):
    """Full Gaussian log-likelihood with H_t = D_t R_t D_t (for AIC vs CCC)."""
    T, N = Rret.shape; e = Rret - np.asarray(mu); D = np.sqrt(H); Z = e / D; ll = 0.0
    for t in range(T):
        _, ld = np.linalg.slogdet(Rt[t])
        logdetH = 2.0 * np.sum(np.log(D[t])) + ld
        ll += -0.5 * (N * np.log(2 * np.pi) + logdetH + Z[t] @ np.linalg.solve(Rt[t], Z[t]))
    return float(ll)


def dcc_fit(Rret):
    """Two-step DCC. Rret is (T x N). Returns dict(fits, H, Z, Qbar, a, b, Rt, loglik)."""
    Rret = np.asarray(Rret, float); T, N = Rret.shape
    fits = [garch11(Rret[:, i]) for i in range(N)]           # stage 1: univariate GARCH per series
    H = np.column_stack([f["h"] for f in fits])
    Z = np.column_stack([f["z"] for f in fits])
    Qbar = (Z.T @ Z) / T                                     # unconditional corr of standardized resids

    def neg(x):
        a, b = x
        if a < 0 or b < 0 or a + b > 0.9999:
            return 1e12
        return -_corr_loglik(Z, Qbar, a, b)

    best = None
    for x0 in ([0.02, 0.97], [0.05, 0.90]):                  # stage 2: optimize (a, b)
        res = minimize(neg, x0, method="Nelder-Mead",
                       options=dict(xatol=1e-7, fatol=1e-7, maxiter=4000))
        if best is None or res.fun < best.fun:
            best = res
    a, b = best.x
    Rt = _Rt_path(Z, Qbar, a, b)
    ll = _full_loglik(Rret, H, Rt, mu=np.array([f["mu"] for f in fits]))
    return dict(fits=fits, H=H, Z=Z, Qbar=Qbar, a=a, b=b, Rt=Rt, loglik=ll)


def pair_corr(Rt, i, j):
    """Extract the conditional-correlation time series between assets i and j."""
    return Rt[:, i, j]