Copula-GARCH: Time-Varying Volatility Meets Dependence

Two Moods of the Market, Separated

The copula project treated dependence as fixed and the coherent-risk project ended on a diagnosis it could not treat: both static VaR methods failed the independence test because their breaches arrived in clusters. That clustering is volatility clustering, and this project is the model that answers it. Copula-GARCH splits multi-asset returns into two pieces estimated separately — each asset's time-varying volatility from a GARCH marginal, and the dependence of the standardised residuals from a copula — then reassembles them into a risk model whose forecasts breathe with the market.

The marginals, and an honest word about the residuals

The first step fits GARCH(1,1)-$t$ to each of five cross-asset ETFs over 3,772 daily observations. The volatility clustering in the raw returns is overwhelming — Ljung–Box $Q$ statistics in the thousands against a critical value near 31 — and after GARCH it is gone for SPY, TLT and GLD. Not for HYG and EEM, whose residual squares remain autocorrelated ($p=0.07$ and $p=0.00$). The notebooks now say so rather than declaring the residuals i.i.d.: they are close to the invariants the copula step assumes, and every later step inherits that approximation. HYG's persistence also sits pinned at $\alpha+\beta=1$, so its variance process is effectively integrated — harmless for a one-day-ahead VaR, which never needs an unconditional variance.

How much joint-crash risk was really just shared volatility?

Then the result that justifies the whole decomposition. Fit a $t$-copula to the raw returns and to the residuals, and the tail dependence falls: for SPY–HYG, $\lambda_L$ drops from 0.38 to 0.22, with the copula's degrees of freedom rising from 5 to 10. Part of what looked like a tendency to crash together was never dependence at all — it was shared volatility, both assets being turbulent at the same time. Copula-GARCH assigns each effect to its proper source. The Bayesian engine confirms the drop is real rather than sampling noise: fitting bivariate-$t$ models to both series gives $\lambda_L$ posteriors of 0.49 and 0.25, and the entire posterior of the difference sits above zero.

A VaR That Breathes

The payoff is a one-day-ahead VaR built from each day's conditional volatilities scaled against a copula-drawn residual pool. It rises the moment volatility spikes and relaxes as calm returns, where a rolling-window historical VaR only drifts up after its window fills with bad days and then stays high too long. On the 1% backtest the static VaR gets the rate right — 34 violations, 1.04%, Kupiec $p=0.82$ — and fails Christoffersen at $p=0.000$, exactly the clustering the coherent-risk project diagnosed. Copula-GARCH passes both: 39 violations at 1.03%, Kupiec $p=0.84$, Christoffersen $p=0.42$.

Making the comparison fair

That comparison had an asymmetry worth removing, because a backtest is the last place to let one side keep an advantage. As written, copula-GARCH used GARCH parameters estimated on the whole sample while the static baseline saw only a trailing window. So it is also done the hard way — marginals and copula fitted on a 1,000-day burn-in only, variances rolled forward with frozen parameters, both models scored on the remaining days. The verdict survives: copula-GARCH still passes both tests (30 violations, 1.08%, Christoffersen $p=0.34$) and the static VaR still fails independence at $p=0.000$. The advantage is structural, not an artefact of the fitting window.

The One Relationship the Arc Assumed Was a Constant

One number in that residual matrix deserves to stop the reader, and it is the sharpest limitation the whole arc runs into. SPY–TLT — equities against long Treasuries — is the diversification relationship a multi-asset portfolio leans on hardest. Two questions get asked in order. Does the volatility decomposition explain it? No: raw $-0.298$ becomes residual $-0.282$, a change of $+0.016$. Unlike the equity–credit pair, where stripping volatility nearly halved the tail dependence, here the decomposition barely moves the number — the negative equity–bond relationship is not a shared-volatility artefact at all.

Then the second question: is it stable? The rolling 500-day correlation averages $-0.56$ in 2011 and drifts, decade-long, to $+0.16$ in 2024 — crossing zero on 6 January 2023. Split the sample and it is $-0.45$ over 2010–2021 against $+0.09$ over 2022–2024. Eighty-eight per cent of windows are negative, so the textbook story holds for most of the history, which is exactly what makes the reversal easy to miss. This is a problem no earlier project in the arc can solve: shrinkage, Bayesian estimation, Black–Litterman and robust allocation all treat the covariance matrix as a fixed object to be estimated better from noisy data, and no amount of shrinkage repairs an estimate of a parameter that is not constant. An investor who bought the 60/40 hedge on its 2010s correlation held a position whose defining property had quietly reversed.

The same lesson, arrived at from the other direction

The Bayesian residual copula sharpens the same point from the other direction. An LKJ prior over the residual correlation matrix returns every pair with a credible interval, and the intervals are not uniform — SPY–EEM is pinned to within 0.029 while SPY–GLD spans 0.063, a distinction a point-estimate correlation matrix cannot express. For SPY–TLT it returns a posterior mean of $-0.28$ lying entirely below zero, $[-0.31,\,-0.25]$. Emphatic — and about an average over 2010–2024 of a quantity that changed sign inside the window. A tight posterior around a parameter that is not constant is confident about the wrong thing: precise on the average, silent on the drift. Letting volatility and dependence move over time, which is what copula-GARCH does, is the direction the answer lies in.

Notebooks

Downloads

Copula-GARCH Module — Source Code

"""
copulagarch.py -- Copula-GARCH: time-varying volatility + a dependence copula.

Backs the notebooks in  "Copula-GARCH: Time-Varying Volatility Meets Dependence".  It fuses two engines already
in the collection:
  * garch_scratch.py  -- univariate GARCH(1,1)-t for each asset's VOLATILITY,
  * copulas.py        -- a copula for the DEPENDENCE of the standardized residuals.

WHY
---
A static copula fit to raw returns confounds two very different things: the fact
that markets go through calm and turbulent regimes (volatility clustering) and
the fact that assets move together (dependence).  The copula-GARCH decomposition
separates them:
  1. Fit a GARCH marginal to each asset -> the conditional volatility sigma_{j,t}
     and the standardized residual  z_{j,t} = (r_{j,t} - mu_j) / sigma_{j,t}.
     The residuals are close to i.i.d. -- GARCH removes the clustering, so the
     z's are the genuine "invariants" of the horizon project.
  2. Fit a copula to the residuals' dependence.
Re-assembling  r_{j,t} = mu_j + sigma_{j,t} z_{j,t}  with z drawn from the copula
gives a portfolio model whose risk BREATHES with the market: a one-day-ahead VaR
that rises in turbulent regimes and falls in calm ones -- which a static VaR
cannot do, and which is exactly what its clustered backtest violations demanded.
"""

import numpy as np
import garch_scratch as gs
import copulas as cp
from scipy.stats import t as student_t
from scipy.optimize import minimize


# --------------------------------------------------------------------------- #
#  Step 1 -- GARCH marginals                                                    #
# --------------------------------------------------------------------------- #

def fit_marginal(r):
    """Fit a GARCH(1,1) with Student-t innovations to one return series, reusing
    garch_scratch's log-likelihood and variance recursion but a self-contained,
    robust optimiser (avoids the engine's numerical-Hessian step). Returns the
    parameters, the conditional-volatility path, and the standardized residuals z
    (which should be ~ i.i.d.)."""
    r = np.asarray(r, float); mu = r.mean(); rc = r - mu
    r2, s0 = gs.prepare(rc)
    nll = lambda th: -gs.garch_loglik(th, r2, s0, gjr=False)         # 4-vector -> Student-t
    init = np.array([0.05 * np.var(rc), 0.05, 0.90, 7.0])
    best = None
    for start in (init, [0.1 * np.var(rc), 0.10, 0.85, 5.0], [0.02 * np.var(rc), 0.03, 0.95, 10.0]):
        res = minimize(nll, np.asarray(start, float), method="Nelder-Mead",
                       options={"maxiter": 8000, "xatol": 1e-7, "fatol": 1e-7})
        if np.isfinite(res.fun) and (best is None or res.fun < best.fun):
            best = res
    mle = best.x
    s2 = gs.garch_var(mle, r2, s0, gjr=False)
    return dict(theta=mle, mu=mu, sigma=np.sqrt(s2), sigma2=s2,
                z=rc / np.sqrt(s2), nu=float(mle[-1]), se=None,
                persistence=float(mle[1] + mle[2]))


def fit_all(X):
    return [fit_marginal(X[:, j]) for j in range(X.shape[1])]


def std_resid_matrix(marg):
    """Matrix of standardized residuals (T x N)."""
    return np.column_stack([m["z"] for m in marg])


# --------------------------------------------------------------------------- #
#  Step 2 -- copula on the residuals, and re-assembly                           #
# --------------------------------------------------------------------------- #

def residual_pool(marg, R_t, nu_c, n, rng, semiparametric=True):
    """Draw n standardized-residual vectors with the fitted t-copula dependence
    (R_t, nu_c). SEMIPARAMETRIC (default): invert each uniform draw through the
    asset's EMPIRICAL residual distribution -- the standard copula-GARCH choice,
    which calibrates the tails correctly. Set semiparametric=False to use a
    parametric standardized-Student-t marginal instead."""
    N = len(marg)
    U = cp.sim_t_copula(R_t, nu_c, n, rng)
    if semiparametric:
        return np.column_stack([np.quantile(marg[j]["z"], np.clip(U[:, j], 1e-4, 1 - 1e-4))
                                for j in range(N)])
    Z = np.empty((n, N))
    for j in range(N):
        nu = marg[j]["nu"]
        Z[:, j] = student_t.ppf(np.clip(U[:, j], 1e-6, 1 - 1e-6), nu) * np.sqrt((nu - 2) / nu)
    return Z


def dynamic_var(marg, w, Z_pool, alpha=0.01):
    """One-day-ahead copula-GARCH VaR for every day, using each day's conditional
    volatilities to scale a fixed pool of copula residuals. Returns VaR_t (>=0)."""
    w = np.asarray(w, float); N = len(marg)
    sig = np.column_stack([m["sigma"] for m in marg])       # (T, N) conditional vols
    mu = np.array([m["mu"] for m in marg])
    T = sig.shape[0]; var = np.empty(T)
    for t in range(T):
        port = (Z_pool * (w * sig[t])) .sum(1) + float(w @ mu)   # portfolio return sims for day t
        var[t] = -np.quantile(port, alpha)
    return var


def static_var_rolling(r_port, window, alpha=0.01):
    """Rolling historical VaR of the portfolio return series (the static baseline)."""
    var = np.full(len(r_port), np.nan)
    for t in range(window, len(r_port)):
        var[t] = -np.quantile(r_port[t - window:t], alpha)
    return var

References

• Patton, A. J. (2006). Modelling asymmetric exchange rate dependence. International Economic Review 47(2), 527–556. — the copula-GARCH construction used here
• Bollerslev, T. (1986). Generalized autoregressive conditional heteroskedasticity. Journal of Econometrics 31(3), 307–327. — the GARCH marginals
• Jondeau, E. & Rockinger, M. (2006). The copula-GARCH model of conditional dependencies. Journal of International Money and Finance 25(5), 827–853. — two-step estimation and its properties
• Joe, H. (2014). Dependence Modeling with Copulas. Chapman & Hall/CRC. — the semiparametric residual pool and tail-dependence estimation
• Christoffersen, P. F. (1998). Evaluating interval forecasts. International Economic Review 39(4), 841–862. — the independence test the static VaR fails and this model passes
• Kupiec, P. H. (1995). Techniques for verifying the accuracy of risk measurement models. Journal of Derivatives 3(2), 73–84. — the proportion-of-failures test
• Meucci, A. (2005). Risk and Asset Allocation. Springer. — the invariance principle the standardized residuals realise