Copulas and Tail Dependence

Splitting Dependence From Marginals

Every model in this arc so far — shrinkage, Bayesian estimation, Black–Litterman, robust allocation — has summarised co-movement with a covariance matrix. That is sufficient only if returns are jointly normal, and they are not. A correlation cannot express the thing a risk manager most wants to know: not do these assets move together, but do they crash together. Sklar's theorem separates the two questions — the marginals describe each asset alone, the copula describes only the dependence — so fat-tailed marginals and the dependence structure can be chosen independently.

Same correlation, different disasters

The clearest demonstration in the project takes five copulas calibrated to the same rank correlation $\tau=0.6$ — Kendall's $\tau$, which measures only the tendency of two series to move in the same order rather than by the same amount, so it is unchanged by any monotone transform of either marginal and is therefore the natural correlation to hold fixed when the whole point is to vary the dependence structure and plots their corners. They are indistinguishable by correlation and completely different in the tails: the Gaussian and Frank corners empty out, the Student-$t$ fills both, Clayton fills the lower one (joint crashes), Gumbel the upper (joint rallies). Same number, different disasters. Fitted to five cross-asset ETFs over 521 weeks, the Student-$t$ copula wins on AIC — the Akaike criterion, which scores fit and charges two units per extra parameter, so it can rank non-nested families like these; lower is better — in both engines — Python's grid fit lands on $\nu=10$, R's continuous fit on 8.7, and the data demand tail dependence either way.

A Limit Is Not a Probability

The per-pair results carry a distinction worth stating plainly, because it cuts the other way from the project's argument. The tail-dependence coefficient $\lambda_L$ is a limit as $q\to0$, not a probability at any threshold anyone can measure. Setting an empirical frequency measured at $q=0.10$ against the Gaussian's limit of zero makes the Gaussian look infinitely wrong. Simulate the fitted Gaussian's own conditional-crash frequency at that same $q=0.10$ and it gives 0.47 for SPY–HYG, against 0.52 in the data — nearly right. Both quantities now appear side by side. The Gaussian's failure is real but asymptotic: it opens up only as the threshold shrinks, which is exactly what the notebook's conditional-crash plot always showed.

Where the copula actually bites

The payoff section asks whether any of this changes a number. For the headline 1% VaR of a diversified portfolio, barely — the marginals dominate the aggregate and both models share them. The copula bites in the joint tail, and every figure is now anchored to the 521 observed weeks rather than compared model-to-model in a vacuum. Three assets crashing together happened 2.11% of the time; the Gaussian says 1.49%, the $t$ 1.86%. All five together happened 0.19% of the time; the Gaussian says 0.02%, the $t$ 0.05%. The deeper the event, the more the choice matters — and both copulas still fall short of what actually happened. The $t$ is the better of the two, not a match for the data.

A Tight Posterior Is Not a Known Quantity

The Bayesian notebook produced the sharpest finding, and it is a caution rather than a result. Fitting a Clayton copula to SPY–HYG gives a posterior for $\lambda_L$ with mean 0.63 and a 94% interval of [0.59, 0.66] — tight, and easy to read as a well-known joint-crash probability. It is not. The $t$-copula fitted to the same pair puts $\lambda_L$ at 0.185. The gap between the two families is 0.44; the posterior spans 0.078. Model choice moves the number roughly six times further than parameter uncertainty does, and a credible interval reports only the second.

Asking the data which family has the right shape

Which family is right is a question the data can answer, so it is put to them. The observed conditional-crash frequency decays with depth — 0.64 at $q=0.20$ down to 0.40 at $q=0.02$. Clayton's runs flat at ~0.63, because its single parameter forces the body and the tail to share one number, and fitting 521 mostly-ordinary weeks inflates the tail to compensate. The $t$-copula tracks the decay (mean absolute error 0.057 against Clayton's 0.111) and wins on AIC for this very pair, $-350.2$ against $-321.9$. So the qualitative verdict survives — the data reject "no joint crashes" with near-certainty — but the 0.63 is a Clayton statement, not a data statement. Quoting the interval as if it bounded the crash probability would be precisely the overconfidence copulas exist to cure.

Notebooks

Downloads

Copula Module — Source Code

"""
copulas.py -- Copulas and tail dependence for asset returns.

Backs the notebooks in  "Copulas and Tail Dependence".
Follows Meucci, "Risk and Asset Allocation" (Springer 2005), chapter 2-B.

THE IDEA (Sklar's theorem)
--------------------------
Any joint distribution of returns splits into two independent pieces:
  * the MARGINALS  F_1,...,F_N  (each asset's own distribution), and
  * the COPULA  C  (how the assets move TOGETHER), a distribution on the unit
    cube whose own marginals are uniform.
Formally  F(x_1,...,x_N) = C(F_1(x_1),...,F_N(x_N)).  Modelling the two pieces
separately is powerful: we can fit fat-tailed marginals to each asset and choose
the dependence structure independently.

Why it matters for risk: the popular GAUSSIAN copula has *zero tail dependence* --
it says extreme joint moves (everything crashing at once) are vanishingly rare.
Real markets have strong lower-tail dependence, which the Student-t and Clayton
copulas capture. Mispricing that dependence is, famously, part of what blew up
credit models in 2008.

CONVENTIONS
-----------
X : (T, N) returns.   U : (T, N) 'pseudo-observations' in (0,1) = ranks/(T+1).
Bivariate helpers take two columns; elliptical copulas (Gaussian, t) take a
correlation matrix R.
"""

import numpy as np
from scipy.stats import norm, t as student_t, rankdata, kendalltau


# --------------------------------------------------------------------------- #
#  Pseudo-observations and rank-dependence measures                             #
# --------------------------------------------------------------------------- #

def pseudo_obs(X):
    """Map each column to uniform pseudo-observations via its empirical CDF:
    u = rank(x)/(T+1). This strips away the marginals, leaving only the
    dependence structure -- the copula's data."""
    X = np.asarray(X, float); T = X.shape[0]
    return np.column_stack([rankdata(X[:, j]) / (T + 1) for j in range(X.shape[1])])


def kendall_tau(u, v):
    return kendalltau(u, v)[0]


def spearman_rho(u, v):
    return np.corrcoef(rankdata(u), rankdata(v))[0, 1]


# --------------------------------------------------------------------------- #
#  Gaussian copula                                                              #
# --------------------------------------------------------------------------- #

def fit_gaussian_copula(U):
    """Correlation of the normal scores z = Phi^{-1}(u)."""
    Z = norm.ppf(np.clip(U, 1e-6, 1 - 1e-6))
    return np.corrcoef(Z.T)


def sim_gaussian_copula(R, n, rng):
    Z = rng.multivariate_normal(np.zeros(len(R)), R, size=n)
    return norm.cdf(Z)


def gaussian_copula_loglik(U, R):
    """Sum of log copula density.  c(u) = |R|^{-1/2} exp(-1/2 z'(R^{-1}-I)z)."""
    Z = norm.ppf(np.clip(U, 1e-9, 1 - 1e-9))
    Ri = np.linalg.inv(R); N = len(R)
    sign, logdet = np.linalg.slogdet(R)
    quad = np.einsum("ti,ij,tj->t", Z, Ri - np.eye(N), Z)
    return np.sum(-0.5 * logdet - 0.5 * quad)


# --------------------------------------------------------------------------- #
#  Student-t copula (symmetric tail dependence)                                 #
# --------------------------------------------------------------------------- #

def _t_copula_loglik_given(U, R, nu):
    Z = student_t.ppf(np.clip(U, 1e-9, 1 - 1e-9), nu)
    N = len(R); Ri = np.linalg.inv(R)
    sign, logdet = np.linalg.slogdet(R)
    from scipy.special import gammaln
    q = np.einsum("ti,ij,tj->t", Z, Ri, Z)
    # log multivariate-t density minus sum of log univariate-t densities
    log_num = (gammaln((nu + N) / 2) + (N - 1) * gammaln(nu / 2)
               - N * gammaln((nu + 1) / 2) - 0.5 * logdet)
    log_c = log_num - (nu + N) / 2 * np.log1p(q / nu) \
        + (nu + 1) / 2 * np.sum(np.log1p(Z ** 2 / nu), axis=1)
    return np.sum(log_c)


def fit_t_copula(U, nu_grid=None):
    """Fit correlation R (from Kendall's tau, robust for elliptical copulas) and
    degrees of freedom nu by profile maximum likelihood over a grid."""
    N = U.shape[1]
    R = np.eye(N)
    for i in range(N):
        for j in range(i + 1, N):
            tau = kendalltau(U[:, i], U[:, j])[0]
            R[i, j] = R[j, i] = np.sin(np.pi * tau / 2)      # tau -> linear corr
    if nu_grid is None:
        nu_grid = np.array([3, 4, 5, 6, 8, 10, 15, 20, 30, 50, 100])
    lls = [_t_copula_loglik_given(U, R, nu) for nu in nu_grid]
    nu_hat = nu_grid[int(np.argmax(lls))]
    return R, float(nu_hat)


def sim_t_copula(R, nu, n, rng):
    Z = rng.multivariate_normal(np.zeros(len(R)), R, size=n)
    W = rng.chisquare(nu, size=n) / nu
    Tvar = Z / np.sqrt(W)[:, None]
    return student_t.cdf(Tvar, nu)


# --------------------------------------------------------------------------- #
#  Archimedean copulas (bivariate): Clayton, Gumbel, Frank                       #
# --------------------------------------------------------------------------- #

def clayton_theta(tau):  return 2 * tau / (1 - tau)          # tau = theta/(theta+2)
def gumbel_theta(tau):   return 1 / (1 - tau)                # tau = 1 - 1/theta


def frank_theta(tau):
    """Solve tau = 1 - 4/theta (1 - D1(theta)) for theta (Debye function D1)."""
    from scipy.optimize import brentq
    from scipy.integrate import quad
    def D1(a):
        return quad(lambda t: t / (np.exp(t) - 1), 0, a)[0] / a
    if abs(tau) < 1e-4:
        return 1e-4
    f = lambda th: (1 - 4 / th * (1 - D1(th))) - tau
    lo, hi = (1e-4, 60) if tau > 0 else (-60, -1e-4)
    return brentq(f, lo, hi)


def sim_clayton(theta, n, rng):
    """Conditional-inversion sampling (closed form). Lower-tail dependent."""
    u = rng.uniform(size=n); w = rng.uniform(size=n)
    v = ((w ** (-theta / (1 + theta)) - 1) * u ** (-theta) + 1) ** (-1 / theta)
    return np.column_stack([u, v])


def sim_frank(theta, n, rng):
    u = rng.uniform(size=n); w = rng.uniform(size=n)
    v = -1 / theta * np.log1p(w * (1 - np.exp(-theta)) /
                              (w * (np.exp(-theta * u) - 1) - np.exp(-theta * u)))
    return np.column_stack([u, v])


def _rpstable(alpha, size, rng):
    """Positive alpha-stable subordinator (Laplace transform exp(-t^alpha)),
    alpha in (0,1], via Chambers-Mallows-Stuck."""
    U = rng.uniform(0, np.pi, size=size)
    W = rng.exponential(1.0, size=size)
    return (np.sin(alpha * U) / np.sin(U) ** (1 / alpha)) * \
           (np.sin((1 - alpha) * U) / W) ** ((1 - alpha) / alpha)


def sim_gumbel(theta, n, rng):
    """Marshall-Olkin frailty sampling. Upper-tail dependent."""
    alpha = 1 / theta
    M = _rpstable(alpha, n, rng)
    E = rng.exponential(1.0, size=(n, 2))
    return np.exp(-(E / M[:, None]) ** (1 / theta))


# --------------------------------------------------------------------------- #
#  Tail dependence coefficients                                                 #
# --------------------------------------------------------------------------- #

def tail_dep_t(rho, nu):
    """Lower = upper tail dependence of the bivariate Student-t copula."""
    x = -np.sqrt((nu + 1) * (1 - rho) / (1 + rho))
    return float(2 * student_t.cdf(x, nu + 1))


def tail_dep_clayton(theta):  return dict(lower=2 ** (-1 / theta), upper=0.0)
def tail_dep_gumbel(theta):   return dict(lower=0.0, upper=2 - 2 ** (1 / theta))
def tail_dep_gaussian():      return dict(lower=0.0, upper=0.0)


def empirical_tail_dep(u, v, q=0.05, lower=True):
    """Nonparametric estimate of tail dependence: P(V<q | U<q) (lower) or
    P(V>1-q | U>1-q) (upper), from pseudo-observations."""
    u = np.asarray(u); v = np.asarray(v)
    if lower:
        cond = u < q
        return float(np.mean(v[cond] < q)) if cond.sum() else np.nan
    cond = u > 1 - q
    return float(np.mean(v[cond] > 1 - q)) if cond.sum() else np.nan


# --------------------------------------------------------------------------- #
#  Model comparison                                                             #
# --------------------------------------------------------------------------- #

def aic(loglik, k):
    return 2 * k - 2 * loglik

References

• Sklar, A. (1959). Fonctions de répartition à n dimensions et leurs marges. Publications de l'Institut de Statistique de l'Université de Paris 8, 229–231. — the theorem the whole project rests on
• Joe, H. (2014). Dependence Modeling with Copulas. Chapman & Hall/CRC. — copula families, tail-dependence coefficients and their estimation
• Nelsen, R. B. (2006). An Introduction to Copulas (2nd ed.). Springer. — Archimedean families and the conditional-inversion samplers used here
• McNeil, A. J., Frey, R. & Embrechts, P. (2015). Quantitative Risk Management (2nd ed.). Princeton University Press. — elliptical copulas, tail dependence and the Student-t formula
• Marshall, A. W. & Olkin, I. (1988). Families of multivariate distributions. JASA 83(403), 834–841. — the frailty construction behind the Gumbel sampler
• Hofert, M., Kojadinovic, I., Mächler, M. & Yan, J. (2018). Elements of Copula Modeling with R. Springer. — the copula package used by the independent R engine
• Lewandowski, D., Kurowicka, D. & Joe, H. (2009). Generating random correlation matrices based on vines and extended onion method. Journal of Multivariate Analysis 100(9), 1989–2001. — the LKJ prior on the Bayesian copula's correlation matrix
• Meucci, A. (2005). Risk and Asset Allocation. Springer, ch. 2-B. — copulas in the allocation framework of this arc