Copulas and Tail Dependence¶
The dependence-modelling engine of the Risk and Asset Allocation arc — from scratch¶
"Correlation is a single number that pretends to summarise how two assets move together. It cannot tell you the one thing risk managers care about most: do they crash together?"
Self-contained; no prior reading of Meucci's Risk and Asset Allocation (Springer 2005, ch. 2) required.
The problem with correlation¶
Every model so far — shrinkage, Bayesian estimation, Black–Litterman, robust allocation — summarised co-movement with a covariance matrix. That is enough only if returns are jointly normal. They are not. Two facts that a correlation cannot express:
- Marginals are fat-tailed and often skewed — each asset has its own non-normal distribution.
- Dependence is stronger in the tails — assets that are mildly correlated in calm markets plunge together in crashes. This is tail dependence, and the Gaussian model assigns it a value of zero.
Copulas solve both by separating the two questions. Sklar's theorem says any joint distribution factorises as $$F(x_1,\dots,x_N)=C\big(F_1(x_1),\dots,F_N(x_N)\big),$$ where the marginals $F_i$ describe each asset alone and the copula $C$ — a distribution on the unit cube with uniform margins — describes only the dependence. We can therefore fit fat-tailed marginals to each asset and choose the dependence structure independently. Choosing a copula with the right tail dependence is, famously, part of what the 2008 credit models got catastrophically wrong.
Roadmap¶
- Sklar's theorem in practice — pseudo-observations strip away the marginals.
- The copula zoo — Gaussian, Student-$t$, Clayton, Gumbel, Frank, and their tails.
- Tail dependence — the coefficient, and why Gaussian $=0$ is dangerous.
- Fitting to real cross-asset data — the $t$-copula wins.
- The risk payoff — portfolio tail risk under a Gaussian vs a $t$ copula.
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import copulas as cp
plt.rcParams.update({"figure.figsize": (9,4.5), "axes.grid": True, "grid.alpha": .25,
"axes.spines.top": False, "axes.spines.right": False, "font.size": 11})
BLUE, ORANGE, GREEN, RED, GREY, PURP = "#2b6cb0","#dd6b20","#2f855a","#c53030","#718096","#6b46c1"
np.set_printoptions(precision=3, suppress=True)
R = pd.read_csv("crossasset_weekly.csv", index_col=0, parse_dates=True)
ASSETS = list(R.columns); X = R.values; T, N = X.shape
print("Cross-asset panel:", T, "weeks,", ASSETS)
print(R.corr().round(2))
Cross-asset panel: 521 weeks, ['SPY', 'TLT', 'GLD', 'HYG', 'EEM']
SPY TLT GLD HYG EEM
SPY 1.00 0.03 0.17 0.75 0.75
TLT 0.03 1.00 0.39 0.21 0.06
GLD 0.17 0.39 1.00 0.27 0.32
HYG 0.75 0.21 0.27 1.00 0.65
EEM 0.75 0.06 0.32 0.65 1.00
The dataset¶
Weekly log-returns (%) of five liquid ETFs spanning the major asset classes, 2015–2024 (521 weeks):
| Ticker | Exposure |
|---|---|
| SPY | US equities (S&P 500) |
| TLT | Long US Treasuries (bonds) |
| GLD | Gold |
| HYG | US high-yield corporate credit |
| EEM | Emerging-market equities |
We chose a cross-asset set on purpose: it contains very different dependence structures — equities and high-yield credit crash together (strong tail dependence), while bonds and gold behave almost independently of equities, sometimes rallying when stocks fall. Correlation alone (above) hints at this but cannot capture the tail behaviour that copulas reveal.
1. Sklar's theorem in practice: pseudo-observations¶
To isolate the copula from the data we replace each asset's returns by their ranks, scaled to $(0,1)$: $$u_{t,i} = \frac{\text{rank}(x_{t,i})}{T+1}\ \approx\ F_i(x_{t,i}).$$ These pseudo-observations have uniform margins by construction, so any structure left in their joint scatter is pure dependence — the copula's data. The transform is monotone, so it preserves the ranking (and hence rank-based dependence) while erasing the marginal shape.
U = cp.pseudo_obs(X)
i, j = ASSETS.index("SPY"), ASSETS.index("HYG")
fig, ax = plt.subplots(1, 2, figsize=(12, 4.6))
ax[0].scatter(X[:,i], X[:,j], s=8, alpha=.4, color=BLUE)
ax[0].set_xlabel("SPY weekly return (%)"); ax[0].set_ylabel("HYG weekly return (%)")
ax[0].set_title("Raw returns (marginals + dependence mixed)")
ax[1].scatter(U[:,i], U[:,j], s=8, alpha=.4, color=PURP)
ax[1].set_xlabel("SPY rank $u$"); ax[1].set_ylabel("HYG rank $v$")
ax[1].set_title("Pseudo-observations (uniform margins = pure copula)")
plt.tight_layout(); plt.show()
print("Notice the crowding in the bottom-LEFT corner of the copula scatter: when SPY has a very bad")
print("week (u near 0), HYG almost always does too (v near 0). That clustering IS lower-tail dependence,")
print("and it is invisible to a single correlation number (SPY-HYG corr = %.2f)." % np.corrcoef(X[:,i],X[:,j])[0,1])
Notice the crowding in the bottom-LEFT corner of the copula scatter: when SPY has a very bad week (u near 0), HYG almost always does too (v near 0). That clustering IS lower-tail dependence, and it is invisible to a single correlation number (SPY-HYG corr = 0.75).
2. The copula zoo¶
Different copulas encode different shapes of dependence — especially in the corners. We simulate five classic copulas, each calibrated to the same rank correlation (Kendall's $\tau=0.6$), and look at their scatter:
- Gaussian — the dependence implied by a multivariate normal. Elliptical, and crucially no tail dependence: the corners thin out.
- Student-$t$ — elliptical like the Gaussian but with symmetric tail dependence: joint booms and joint crashes cluster.
- Clayton — lower-tail dependence only (joint crashes), the natural model for downside risk.
- Gumbel — upper-tail dependence only (joint rallies).
- Frank — symmetric but with no tail dependence (like Gaussian in the tails, different in the middle).
rng = np.random.default_rng(1); n = 3000; tau = 0.6
rho = np.sin(np.pi*tau/2)
sims = {
"Gaussian": cp.sim_gaussian_copula(np.array([[1,rho],[rho,1]]), n, rng),
"Student-t (nu=3)": cp.sim_t_copula(np.array([[1,rho],[rho,1]]), 3, n, rng),
"Clayton (lower)": cp.sim_clayton(cp.clayton_theta(tau), n, rng),
"Gumbel (upper)": cp.sim_gumbel(cp.gumbel_theta(tau), n, rng),
"Frank (no tails)": cp.sim_frank(cp.frank_theta(tau), n, rng),
}
fig, ax = plt.subplots(1, 5, figsize=(16, 3.4))
for a,(name,S) in zip(ax, sims.items()):
a.scatter(S[:,0], S[:,1], s=4, alpha=.35, color=BLUE)
a.set_title(name, fontsize=10); a.set_xticks([0,1]); a.set_yticks([0,1])
plt.suptitle("Five copulas at the same Kendall's tau = 0.6 — watch the corners", y=1.05)
plt.tight_layout(); plt.show()
print("All five share the same rank correlation, yet their tails differ completely. The Gaussian and")
print("Frank corners are sparse (independence in the extremes); the t fills BOTH corners; Clayton fills")
print("the bottom-left (crashes); Gumbel the top-right (rallies). Same 'correlation', different disasters.")
All five share the same rank correlation, yet their tails differ completely. The Gaussian and Frank corners are sparse (independence in the extremes); the t fills BOTH corners; Clayton fills the bottom-left (crashes); Gumbel the top-right (rallies). Same 'correlation', different disasters.
3. Tail dependence — and why Gaussian $=0$ is dangerous¶
The lower tail-dependence coefficient makes this precise: $$\lambda_L=\lim_{q\to0^+}\Pr\big(U_2<q \mid U_1<q\big)$$ — the probability one asset is in a crash given the other is. Closed forms:
| copula | $\lambda_L$ | $\lambda_U$ |
|---|---|---|
| Gaussian | $0$ | $0$ |
| Student-$t(\rho,\nu)$ | $2\,t_{\nu+1}\!\big(-\sqrt{\tfrac{(\nu+1)(1-\rho)}{1+\rho}}\big)$ | same |
| Clayton$(\theta)$ | $2^{-1/\theta}$ | $0$ |
| Gumbel$(\theta)$ | $0$ | $2-2^{1/\theta}$ |
The Gaussian copula's $\lambda_L=0$ says that no matter how correlated two assets are, simultaneous extreme crashes are asymptotically impossible. That is empirically false and was a load-bearing assumption in pre-2008 credit models. We show it by estimating $\Pr(V<q\mid U<q)$ across shrinking $q$ for the real SPY–HYG pair, against samples from the fitted Gaussian and $t$ copulas.
One caution that governs how the next section's table must be read: $\lambda_L$ is a limit as $q\to0$, not a probability at any threshold you can measure. At a finite $q$ every copula — the Gaussian included — predicts far more joint crashes than its own limit, so an empirical frequency at $q=0.10$ and a $\lambda_L$ of $0$ are simply not comparable quantities. We will therefore report both, side by side.
# empirical vs Gaussian vs t: conditional joint-crash probability as q shrinks
R_t, nu = cp.fit_t_copula(U); R_g = cp.fit_gaussian_copula(U)
rho_ij = R_t[i,j]
big = 60000
Ug = cp.sim_gaussian_copula(np.array([[1,rho_ij],[rho_ij,1]]), big, rng)
Ut = cp.sim_t_copula(np.array([[1,rho_ij],[rho_ij,1]]), nu, big, rng)
qs = np.linspace(0.02, 0.30, 15)
emp = [cp.empirical_tail_dep(U[:,i], U[:,j], q, True) for q in qs]
gaus = [cp.empirical_tail_dep(Ug[:,0], Ug[:,1], q, True) for q in qs]
tcop = [cp.empirical_tail_dep(Ut[:,0], Ut[:,1], q, True) for q in qs]
plt.plot(qs, emp, "o-", color=RED, lw=2, label="SPY–HYG (data)")
plt.plot(qs, tcop, "s--", color=BLUE, lw=2, label="fitted Student-t copula")
plt.plot(qs, gaus, "^--", color=GREY, lw=2, label="fitted Gaussian copula")
plt.xlabel("tail threshold q"); plt.ylabel("P( HYG crash | SPY crash )")
plt.title("Joint-crash probability: the Gaussian copula vanishes, the data does not")
plt.legend(); plt.gca().invert_xaxis(); plt.show()
print("As q -> 0 the Gaussian curve heads to ZERO (its lower tail dependence), while the data stays high")
print("and the t-copula tracks it far better. t-copula tail dependence lambda_L = %.2f vs Gaussian 0.00."
% cp.tail_dep_t(rho_ij, nu))
As q -> 0 the Gaussian curve heads to ZERO (its lower tail dependence), while the data stays high and the t-copula tracks it far better. t-copula tail dependence lambda_L = 0.18 vs Gaussian 0.00.
4. Fitting copulas to the real cross-asset data¶
We fit the two elliptical copulas that scale to many assets — Gaussian and Student-$t$ — to all five ETFs, and compare them by AIC (lower is better). The $t$-copula has one extra parameter (the degrees of freedom $\nu$ controlling tail fatness); if it wins, the data demand tail dependence.
kg = N*(N-1)//2
ll_g = cp.gaussian_copula_loglik(U, R_g)
ll_t = cp._t_copula_loglik_given(U, R_t, nu)
print("Fitted Student-t copula degrees of freedom nu = %.0f (small nu = fat joint tails)\n" % nu)
print(" Gaussian copula : loglik %7.1f AIC %8.1f" % (ll_g, cp.aic(ll_g, kg)))
print(" Student-t copula: loglik %7.1f AIC %8.1f <- %s" %
(ll_t, cp.aic(ll_t, kg+1), "WINS" if cp.aic(ll_t,kg+1) < cp.aic(ll_g,kg) else "loses"))
# Per-pair tail dependence -- comparing like with like.
# lambda_L is a q->0 LIMIT. An empirical frequency measured at a finite q=0.10 cannot be
# set against it: at that depth EVERY copula, the Gaussian included, predicts far more joint
# crashes than its own limit. So we also simulate each fitted copula's OWN P(V<q|U<q) at
# q=0.10 -- the only column directly comparable with the data.
big = 200000
print("\nLower-tail dependence by pair:")
print(" P(V<q | U<q) at q = 0.10 | lambda_L (q -> 0)")
print(" pair data t-cop Gaussian | t-cop Gaussian")
for a in range(N):
for b in range(a+1, N):
rt, rg = R_t[a,b], R_g[a,b]
Ut_ = cp.sim_t_copula(np.array([[1,rt],[rt,1]]), nu, big, rng)
Ug_ = cp.sim_gaussian_copula(np.array([[1,rg],[rg,1]]), big, rng)
print(" %-4s-%-4s %.2f %.2f %.2f | %.2f 0.00" %
(ASSETS[a], ASSETS[b],
cp.empirical_tail_dep(U[:,a], U[:,b], 0.10, True),
cp.empirical_tail_dep(Ut_[:,0], Ut_[:,1], 0.10, True),
cp.empirical_tail_dep(Ug_[:,0], Ug_[:,1], 0.10, True),
cp.tail_dep_t(rt, nu)))
print("\nRead the two blocks against each other. At q=0.10 the Gaussian and the t are nearly")
print("indistinguishable and both sit close to the data: the Gaussian's failure is ASYMPTOTIC,")
print("not something visible at the 10% level. That is exactly what the previous plot showed --")
print("the curves separate only as q shrinks. Setting the empirical q=0.10 figure against the")
print("Gaussian's q->0 limit of zero would make the gap look far larger than it is at any")
print("threshold one can actually measure with %d weeks." % T)
Fitted Student-t copula degrees of freedom nu = 10 (small nu = fat joint tails)
Gaussian copula : loglik 460.9 AIC -901.8
Student-t copula: loglik 477.3 AIC -932.5 <- WINS
Lower-tail dependence by pair:
P(V<q | U<q) at q = 0.10 | lambda_L (q -> 0)
pair data t-cop Gaussian | t-cop Gaussian
SPY -TLT 0.27 0.11 0.10 | 0.00 0.00 SPY -GLD 0.25 0.16 0.15 | 0.01 0.00
SPY -HYG 0.52 0.48 0.47 | 0.18 0.00
SPY -EEM 0.60 0.51 0.49 | 0.20 0.00 TLT -GLD 0.31 0.27 0.25 | 0.04 0.00 TLT -HYG 0.29 0.19 0.16 | 0.02 0.00
TLT -EEM 0.17 0.10 0.09 | 0.00 0.00
GLD -HYG 0.23 0.20 0.18 | 0.02 0.00 GLD -EEM 0.25 0.25 0.21 | 0.03 0.00
HYG -EEM 0.44 0.39 0.39 | 0.11 0.00 Read the two blocks against each other. At q=0.10 the Gaussian and the t are nearly indistinguishable and both sit close to the data: the Gaussian's failure is ASYMPTOTIC, not something visible at the 10% level. That is exactly what the previous plot showed -- the curves separate only as q shrinks. Setting the empirical q=0.10 figure against the Gaussian's q->0 limit of zero would make the gap look far larger than it is at any threshold one can actually measure with 521 weeks.
# fitted t-copula correlation matrix
fig, ax = plt.subplots(figsize=(5.2,4.4))
im = ax.imshow(R_t, vmin=-1, vmax=1, cmap="RdBu_r")
ax.set_xticks(range(N)); ax.set_xticklabels(ASSETS); ax.set_yticks(range(N)); ax.set_yticklabels(ASSETS)
for a in range(N):
for b in range(N): ax.text(b, a, "%.2f"%R_t[a,b], ha="center", va="center", fontsize=9)
ax.set_title("Student-t copula correlation"); fig.colorbar(im, fraction=.046); plt.show()
print("Equities and credit (SPY, HYG, EEM) form a tightly-dependent block: SPY-HYG %.2f, SPY-EEM %.2f." % (R_t[0,3], R_t[0,4]))
print("Bonds sit apart -- SPY-TLT is %.2f, statistically indistinguishable from independence rather" % R_t[0,1])
print("than reliably negative (the Bayesian notebook's posterior for this pair straddles zero). That")
print("near-zero rank correlation is itself the diversification story, and it says nothing yet about")
print("whether TLT stays a diversifier in the tail -- which is the question the copula answers.")
Equities and credit (SPY, HYG, EEM) form a tightly-dependent block: SPY-HYG 0.69, SPY-EEM 0.72. Bonds sit apart -- SPY-TLT is -0.06, statistically indistinguishable from independence rather than reliably negative (the Bayesian notebook's posterior for this pair straddles zero). That near-zero rank correlation is itself the diversification story, and it says nothing yet about whether TLT stays a diversifier in the tail -- which is the question the copula answers.
5. The risk payoff: joint crashes under the wrong copula¶
Does the copula change a risk number? For a single diversified portfolio's headline VaR, only modestly — the fat-tailed marginals dominate the aggregate, and both models share the same ones. Where the copula truly bites is the joint tail: the chance that many assets crash at the same time. That is exactly what tail dependence governs and what a Gaussian copula ($\lambda_L=0$) systematically understates.
Using the same empirical marginals, we draw 200,000 weeks from a Gaussian and from the fitted $t$-copula and compare three things: the equally-weighted portfolio's VaR/ES, the probability of a simultaneous multi-asset crash, and — the sharpest test, anchored to the actual data — given that SPY has a bottom-5% week, how many of the other four assets crash with it.
def simulate(Ucop, X):
return np.column_stack([np.quantile(X[:,k], np.clip(Ucop[:,k], 1e-4, 1-1e-4)) for k in range(X.shape[1])])
w = np.ones(N)/N; m = 200000
Ug = cp.sim_gaussian_copula(R_g, m, rng); Ut = cp.sim_t_copula(R_t, nu, m, rng)
Xg = simulate(Ug, X); Xt = simulate(Ut, X); pg, pt = Xg @ w, Xt @ w
def var_es(r, a=0.01): v = np.quantile(r, a); return v, r[r <= v].mean()
vg, eg = var_es(pg); vt, et = var_es(pt)
thr = np.array([np.quantile(X[:,k], 0.05) for k in range(N)]) # each asset's own 5% crash level
def n_crash(Xd): return (Xd < thr).sum(1)
crash_d, crash_g, crash_t = [np.mean(n_crash(Z) >= 3) for Z in (X, Xg, Xt)]
all_d, all_g, all_t = [np.mean(n_crash(Z) == N) for Z in (X, Xg, Xt)]
sp = ASSETS.index("SPY"); others = [k for k in range(N) if k != sp]
def cocrash(Xd): # E[# of others crashing | SPY crashes]
m_ = Xd[:, sp] < thr[sp]; return (Xd[np.ix_(m_, others)] < thr[others]).sum(1).mean()
co_data, co_g, co_t = cocrash(X), cocrash(Xg), cocrash(Xt)
n_spy = int((X[:, sp] < thr[sp]).sum())
print("Same marginals, different copula. Effect on the DIVERSIFIED portfolio's headline risk:")
print(" 1%% VaR : Gaussian %.2f%% t-copula %.2f%% (marginals dominate -> nearly identical)" % (vg, vt))
print(" 1%% ES : Gaussian %.2f%% t-copula %.2f%%\n" % (eg, et))
print("Where the copula bites -- JOINT tail events, every row anchored to the %d observed weeks:" % T)
print(" P(>=3 of 5 crash together): DATA %.2f%% Gaussian %.2f%% t-copula %.2f%% (t/G %.1fx)"
% (100*crash_d, 100*crash_g, 100*crash_t, crash_t/crash_g))
print(" P(all 5 crash together) : DATA %.2f%% Gaussian %.2f%% t-copula %.2f%% (t/G %.1fx)"
% (100*all_d, 100*all_g, 100*all_t, all_t/all_g))
print(" Given SPY crashes, E[# of the other 4 also crashing] (%d observed SPY-crash weeks):" % n_spy)
print(" DATA %.2f | t-copula %.2f (%+.0f%%) | Gaussian %.2f (%+.0f%%)"
% (co_data, co_t, 100*(co_t-co_data)/co_data, co_g, 100*(co_g-co_data)/co_data))
print("\nTwo things follow. The deeper the joint tail, the more the copula choice matters: the")
print("t/Gaussian ratio grows from the 3-asset to the all-5 event. And BOTH copulas still fall")
print("short of the data -- the t-copula roughly halves the Gaussian's contagion error without")
print("closing it, so the honest claim is that the t is the better of the two, not that it is right.")
Same marginals, different copula. Effect on the DIVERSIFIED portfolio's headline risk:
1% VaR : Gaussian -3.42% t-copula -3.44% (marginals dominate -> nearly identical)
1% ES : Gaussian -4.41% t-copula -4.65%
Where the copula bites -- JOINT tail events, every row anchored to the 521 observed weeks:
P(>=3 of 5 crash together): DATA 2.11% Gaussian 1.49% t-copula 1.86% (t/G 1.2x)
P(all 5 crash together) : DATA 0.19% Gaussian 0.02% t-copula 0.05% (t/G 2.6x)
Given SPY crashes, E[# of the other 4 also crashing] (26 observed SPY-crash weeks):
DATA 1.19 | t-copula 1.05 (-12%) | Gaussian 0.92 (-23%)
Two things follow. The deeper the joint tail, the more the copula choice matters: the
t/Gaussian ratio grows from the 3-asset to the all-5 event. And BOTH copulas still fall
short of the data -- the t-copula roughly halves the Gaussian's contagion error without
closing it, so the honest claim is that the t is the better of the two, not that it is right.
fig, ax = plt.subplots(1, 2, figsize=(13, 4.4))
lo = min(pg.min(), pt.min())
ax[0].hist(pg, bins=np.linspace(lo, 8, 120), density=True, alpha=.5, color=GREY, label="Gaussian copula")
ax[0].hist(pt, bins=np.linspace(lo, 8, 120), density=True, alpha=.5, color=BLUE, label="Student-t copula")
ax[0].axvline(vg, color=GREY, ls="--"); ax[0].axvline(vt, color=BLUE, ls="--")
ax[0].set_xlim(lo, 8); ax[0].set_xlabel("portfolio weekly return (%)"); ax[0].set_ylabel("density")
ax[0].set_title("Portfolio return distribution (1% VaR dashed)"); ax[0].legend()
# number of simultaneous crashes distribution
kk = np.arange(0, N+1)
dg = [np.mean((Xg<thr).sum(1)==k) for k in kk]; dt = [np.mean((Xt<thr).sum(1)==k) for k in kk]
ax[1].bar(kk-0.18, dg, 0.36, color=GREY, label="Gaussian"); ax[1].bar(kk+0.18, dt, 0.36, color=BLUE, label="Student-t")
ax[1].set_yscale("log"); ax[1].set_xlabel("# assets crashing in the same week (of 5)")
ax[1].set_ylabel("probability (log)"); ax[1].set_title("Joint crashes: t-copula has a much fatter cluster"); ax[1].legend()
plt.tight_layout(); plt.show()
6. Summary¶
- Sklar's theorem splits a joint distribution into marginals (each asset alone) and a copula (pure dependence); pseudo-observations expose the copula by erasing the marginals.
- Copulas with the same correlation can have completely different tail behaviour — the Gaussian and Frank have none, the Student-$t$ symmetric, Clayton lower (crashes), Gumbel upper (rallies).
- Tail dependence is the risk-critical feature the Gaussian copula sets to zero; the real cross-asset data (SPY–HYG especially) shows strong lower-tail dependence, and the Student-$t$ copula wins on AIC.
- The payoff: with identical fat-tailed marginals, swapping a Gaussian copula for the fitted $t$-copula barely moves the headline VaR (the marginals dominate it) but raises the probability of a simultaneous multi-asset crash — modestly for three-of-five, and by more as the event gets deeper. Both copulas still understate what the 521 weeks actually delivered, so the $t$ is the better of the two rather than the right one — the exact risk that sinks diversified portfolios, and the reason a single correlation is not enough.
Where this connects. Copulas are the general dependence layer beneath the whole arc: the Gaussian/$t$ covariance of the shrinkage and Bayesian-estimation projects is the elliptical special case, and a $t$-copula is the honest upgrade when joint crashes matter. They also pair naturally with the fat-tailed marginals of the Student-$t$ regression and stochastic-volatility work elsewhere in the collection — and set up a copula-GARCH model for time-varying dependence.