21. Copulas (intro)
A joint distribution bundles two separate pieces of information:
- The marginals — how each variable behaves on its own.
- The dependence — how they move together.
The multivariate Normal forces both at once: choose it, and you've chosen Normal marginals and a specific dependence structure. But real data rarely cooperates. Insurance losses are heavy-tailed, not Normal, yet they still move together.
A copula separates these two concerns, letting you pick each independently.
The probability integral transform
Everything rests on one fact. If X is continuous with CDF F, then
U = F(X) \sim \text{Uniform}(0, 1)
Proof (one line, using that F is increasing):
P(U \le u) = P(F(X) \le u) = P\big(X \le F^{-1}(u)\big) = F\big(F^{-1}(u)\big) = u
which is the uniform CDF.
Feeding any continuous variable through its own CDF gives a uniform. And running it backwards, X = F^{-1}(U), turns a uniform into any distribution you like — that's the inverse transform sampling method of §11, and it's the same fact used twice.
So the CDF is a universal translator: every continuous variable can be converted to a uniform, manipulated, and converted back.
The definition
A copula is a joint CDF on [0,1]^k whose marginals are all uniform. It contains only dependence information — the marginals have been standardised away.
Sklar's theorem makes this precise. Any joint CDF H with continuous marginals F_1, F_2 can be written
H(x_1, x_2) = C\big(F_1(x_1),\; F_2(x_2)\big)
for a unique copula C. And conversely, any copula combined with any marginals produces a valid joint distribution.
That is the whole point:
\text{joint distribution} = \text{marginals} + \text{copula}
and the two halves can be chosen independently. Model each variable's own behaviour however you like, then attach a dependence structure separately.
Some standard copulas
Independence copula. C(u,v) = uv — the variables are independent. (It's just the factorisation rule written in copula form.)
Gaussian copula. Take the dependence structure of a bivariate Normal with correlation \rho, but attach whatever marginals you want:
C(u,v) = \Phi_\rho\big(\Phi^{-1}(u),\; \Phi^{-1}(v)\big)
Its defining feature — and its defining weakness — is zero tail dependence. For any \rho < 1, extreme events become asymptotically independent: given one variable hits a once-in-a-century level, the chance the other does too tends to zero.
t copula. Same construction from a multivariate t. Because the t has heavier tails, it has positive tail dependence — extremes genuinely tend to co-occur.
Archimedean copulas (Clayton, Gumbel, Frank) are built from a generator function and allow asymmetric tail dependence. Clayton couples the lower tail (joint crashes); Gumbel couples the upper tail (joint booms).
Tail dependence, and why it mattered
The lower tail dependence coefficient is
\lambda_L = \lim_{q \to 0^+} P\big(U \le q \mid V \le q\big)
"Given one variable is in its worst q%, how likely is the other to be too?"
- Gaussian copula: \lambda_L = 0 for every \rho < 1
- t copula: \lambda_L > 0, increasing as the degrees of freedom fall
This distinction is not academic. Gaussian copulas were used to price mortgage CDOs before 2008. Correlation was calibrated to normal-times data, and the model implied that simultaneous defaults across many mortgages were vanishingly unlikely. Defaults are precisely the tail, and precisely where the Gaussian copula assumes independence. When the housing market turned, they defaulted together.
The lesson generalises: correlation is a summary of the middle of a distribution. It tells you almost nothing about behaviour in the tails, which is usually what you care about in risk.
Worked example
Build a joint distribution with an Exponential(1) margin, a Beta(2,5) margin, and Gaussian dependence with \rho = 0.7.
The construction is three steps:
- Draw (Z_1, Z_2) from a bivariate Normal with correlation 0.7.
- Transform to uniforms: U_1 = \Phi(Z_1), U_2 = \Phi(Z_2). These are uniform marginally but dependent, and their dependence is the Gaussian copula.
- Transform to the targets with inverse CDFs: X = F_{\text{Exp}}^{-1}(U_1), Y = F_{\text{Beta}}^{-1}(U_2).
X is exactly Exponential(1), Y is exactly Beta(2,5), and they're dependent — none of which any standard named bivariate distribution offers.
One caveat: the rank correlation is preserved by steps 2 and 3 (they're monotone), but the ordinary Pearson correlation is not. Feeding in \rho = 0.7 does not give you Pearson 0.7 in the output; it gives you a specific rank correlation. Spearman's \rho is the natural measure here for exactly that reason.
Doing it in Python
The probability integral transform, in both directions:
import numpy as np
from scipy.stats import norm, expon, kstest
rng = np.random.default_rng(0)
# Forward: any continuous variable -> uniform
X = rng.exponential(scale=2.0, size=200_000)
U = expon(scale=2.0).cdf(X)
print("F(X) uniform? KS p-value:", round(kstest(U, "uniform").pvalue, 4))
print(" mean", round(U.mean(), 4), "(theory 0.5) var", round(U.var(), 4), "(theory 0.0833)")
# Backward: uniform -> any distribution
V = rng.uniform(size=200_000)
Y = norm.ppf(V)
print("\nF^-1(U) normal? KS p-value:", round(kstest(Y, "norm").pvalue, 4))
Now the full copula construction:
import numpy as np
from scipy.stats import norm, expon, beta, spearmanr, kstest
rng = np.random.default_rng(1)
n, rho = 200_000, 0.7
# 1. correlated normals
cov = [[1, rho], [rho, 1]]
Z = rng.multivariate_normal([0, 0], cov, size=n)
# 2. -> dependent uniforms (this pair IS the Gaussian copula)
U = norm.cdf(Z)
# 3. -> the marginals we actually want
X = expon(scale=1).ppf(U[:, 0])
Y = beta(2, 5).ppf(U[:, 1])
print("X really Exponential(1)? KS p:", round(kstest(X, expon(scale=1).cdf).pvalue, 4))
print("Y really Beta(2,5)? KS p:", round(kstest(Y, beta(2, 5).cdf).pvalue, 4))
print()
print("but they are dependent:")
print(" Spearman(X, Y) =", round(spearmanr(X, Y).statistic, 4))
print(" Pearson (X, Y) =", round(np.corrcoef(X, Y)[0, 1], 4), "(NOT the 0.7 we fed in)")
And the tail dependence difference — the 2008 lesson in ten lines:
import numpy as np
from scipy.stats import norm, t as student_t
rng = np.random.default_rng(2)
n, rho = 500_000, 0.7
cov = [[1, rho], [rho, 1]]
# Gaussian copula
U_gauss = norm.cdf(rng.multivariate_normal([0, 0], cov, size=n))
# t copula with 3 degrees of freedom (same rho, heavier tails)
df = 3
Z = rng.multivariate_normal([0, 0], cov, size=n)
W = rng.chisquare(df, size=n) / df
U_t = student_t(df).cdf(Z / np.sqrt(W)[:, None])
print(f"{'quantile q':>12} {'gaussian':>12} {'t (df=3)':>12}")
for q in (0.10, 0.05, 0.01, 0.005):
g = ((U_gauss[:, 0] < q) & (U_gauss[:, 1] < q)).sum() / (U_gauss[:, 0] < q).sum()
tt = ((U_t[:, 0] < q) & (U_t[:, 1] < q)).sum() / (U_t[:, 0] < q).sum()
print(f"{q:>12} {g:>12.4f} {tt:>12.4f}")
print("\nP(both extreme | one extreme): Gaussian fades toward 0, t does not.")
Same correlation, radically different behaviour where it matters.
Your turn
1. X \sim \text{Exponential}(2). What's the distribution of F(X)?
2. Why can't the Gaussian copula model "these two stocks crash together"?
3. You have Uniform(0,1) random numbers and want Exponential(\lambda) draws. How?
Solutions
1. Uniform(0,1) — by the probability integral transform, regardless of the rate. The rate 2 is irrelevant: any continuous CDF applied to its own variable gives a standard uniform. That universality is what makes the CDF usable as a translator between distributions.
2. Because the Gaussian copula has zero tail dependence for any \rho < 1. As you look at more and more extreme joint events, it implies
P(\text{both crash} \mid \text{one crashes}) \to 0
So no matter how high you set \rho, the model says simultaneous extreme losses are asymptotically independent. You can raise the correlation to fit everyday co-movement and the model will still understate joint disasters — the parameter has no way to express that behaviour.
A t copula (or a Clayton copula, which loads the lower tail specifically) keeps \lambda_L > 0 and can represent crashing together.
3. Invert the CDF. For an Exponential, F(x) = 1 - e^{-\lambda x}; set u = F(x) and solve:
u = 1 - e^{-\lambda x} \implies x = -\frac{\ln(1 - u)}{\lambda}
And since 1 - U is uniform whenever U is, the usual simplification is
X = -\frac{\ln U}{\lambda}
which is exactly the transformation we derived in lesson 3 of this module. It appeared there as a Jacobian exercise; here it's revealed as one instance of a completely general recipe.
Check yourself in code
Verify the probability integral transform and then build a joint distribution with an Exponential margin and a Beta margin coupled by a Gaussian copula, confirming the marginals survive the construction.
Print exactly this:
U mean 0.5
U var 0.0829
X is exponential: True
Y is beta: True
dependent: True
Round the uniform's mean to 1 decimal place and its variance to 4. For the
marginal checks use a Kolmogorov–Smirnov test and report whether the p-value
exceeds 0.01. For the last line, report whether the Spearman correlation
exceeds 0.5. Use default_rng(1) and 200000 draws.
import numpy as np
from scipy.stats import norm, expon, beta, spearmanr, kstest
rng = np.random.default_rng(1)
n, rho = 200_000, 0.7
Z = rng.multivariate_normal([0, 0], [[1, rho], [rho, 1]], size=n)
U = norm.cdf(Z)
print("U mean", round(U[:, 0].mean(), 1))
print("U var", round(U[:, 0].var(), 4))
# Transform column 0 to Exponential(1) and column 1 to Beta(2, 5),
# check each marginal with kstest (p > 0.01), then check they are dependent
# via spearmanr > 0.5.
import numpy as np
from scipy.stats import norm, expon, beta, spearmanr, kstest
rng = np.random.default_rng(1)
n, rho = 200_000, 0.7
Z = rng.multivariate_normal([0, 0], [[1, rho], [rho, 1]], size=n)
U = norm.cdf(Z)
print("U mean", round(U[:, 0].mean(), 1))
print("U var", round(U[:, 0].var(), 4))
X = expon(scale=1).ppf(U[:, 0])
Y = beta(2, 5).ppf(U[:, 1])
print("X is exponential:", kstest(X, expon(scale=1).cdf).pvalue > 0.01)
print("Y is beta:", kstest(Y, beta(2, 5).cdf).pvalue > 0.01)
print("dependent:", bool(spearmanr(X, Y).statistic > 0.5))
A copula strips the marginals out of a joint distribution and keeps only the dependence. Sklar's theorem guarantees the split is always possible and unique, so you can model each variable's own behaviour separately from how they move together. And the choice of copula decides what happens in the tails — which correlation alone will never tell you.
That closes §2. Next: what happens to averages as samples grow — the limit theorems that make statistics possible at all.