Bayesian Copulas & Tail Dependence — the PyMC engine¶

The from-scratch notebook fit copulas by maximum likelihood and read off a single number for the tail dependence. But that number is estimated — and a risk manager betting on how often assets crash together should want to know how sure we are of it. This notebook does copula estimation in PyMC, delivering the posterior distribution of the dependence structure and, most importantly, of the tail-dependence coefficient itself.

Self-contained; no prior reading required.

Two Bayesian models¶

  1. Gaussian copula — an LKJ prior on the correlation matrix gives the full posterior over the dependence structure of all five assets, with credible intervals.
  2. Tail dependence (Clayton copula) — fit to the equity–credit pair with the strongest joint-crash behaviour, giving the posterior of the lower-tail-dependence coefficient $\lambda_L$. Because $\lambda_L=0$ under a Gaussian, its posterior tells us directly how firmly the data reject the "no joint crashes" assumption. It also supplies a cautionary lesson, which we make explicit rather than gloss: a tight posterior is not the same thing as a well-known quantity, because it is conditional on the copula family being right.
In [1]:
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import pymc as pm, arviz as az, pytensor.tensor as pt
from scipy.stats import norm
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"
RNG = 7
print("PyMC", pm.__version__)

R = pd.read_csv("crossasset_weekly.csv", index_col=0, parse_dates=True)
ASSETS = list(R.columns); X = R.values; T, N = X.shape
U = cp.pseudo_obs(X)                                    # uniform pseudo-observations
print("Cross-asset panel:", T, "weeks,", ASSETS)

g++ not available, if using conda: `conda install gxx`

PyMC 6.0.1
Cross-asset panel: 521 weeks, ['SPY', 'TLT', 'GLD', 'HYG', 'EEM']

The dataset¶

Weekly log-returns (%) of five cross-asset ETFs, 2015–2024 (521 weeks): SPY (US equity), TLT (Treasuries), GLD (gold), HYG (high-yield credit), EEM (EM equity). Chosen for their variety of dependence — equities and credit crash together, bonds and gold are near-independent of equities. We work with the pseudo-observations (ranks in $(0,1)$), which carry only the dependence.

1. Bayesian Gaussian copula: the correlation with error bars¶

The Gaussian copula's dependence is a correlation matrix $R$ acting on the normal scores $z=\Phi^{-1}(u)$. Since the scores are fixed data, placing an LKJ prior on $R$ and modelling $z\sim\mathcal N(0,R)$ yields the copula's posterior directly (the fixed normal marginal terms don't involve $R$). We get every pairwise dependence with a credible interval — the Bayesian upgrade of the point estimate.

In [2]:
Z = norm.ppf(np.clip(U, 1e-4, 1-1e-4))
with pm.Model() as gauss_cop:
    chol, corr, sds = pm.LKJCholeskyCov("L", n=N, eta=2.0,
                                        sd_dist=pm.HalfNormal.dist(1.0), compute_corr=True)
    pm.Deterministic("R", corr)
    pm.MvNormal("z", mu=np.zeros(N), chol=chol, observed=Z)
    idata_g = pm.sample(800, tune=1000, chains=4, cores=1, target_accept=0.9,
                        progressbar=False, random_seed=RNG)
Rpost = idata_g.posterior["R"].mean(("chain","draw")).values
# R-hat over the off-diagonal only: the correlation matrix's diagonal is a constant 1,
# for which between/within chain variance is 0/0 and R-hat is undefined.
iu = np.triu_indices(N, 1)
off = idata_g.posterior["R"].values[:, :, iu[0], iu[1]]
print("max R-hat (off-diagonal):", float(az.rhat(az.convert_to_dataset(off))["x"].max()))

Initializing NUTS using jitter+adapt_diag...

Sequential sampling (4 chains in 1 job)

NUTS: [L]

Sampling 4 chains for 1_000 tune and 800 draw iterations (4_000 + 3_200 draws total) took 3 seconds.

max R-hat (off-diagonal): 1.0032110191731909
In [3]:
# posterior mean correlation + 94% credible interval for each pair
hdi = az.hdi(idata_g, var_names=["R"])["R"].values
print("Gaussian-copula dependence with uncertainty:\n")
print("%-10s %10s %18s" % ("pair", "post. mean", "94% credible"))
for a in range(N):
    for b in range(a+1, N):
        m = Rpost[a,b]; lo, hi = hdi[a,b,0], hdi[a,b,1]
        print("%-10s %10.2f %14s" % (ASSETS[a]+"-"+ASSETS[b], m, "[%.2f, %.2f]" % (lo, hi)))
fig, ax = plt.subplots(figsize=(5.2,4.4))
im = ax.imshow(Rpost, 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"%Rpost[a,b],ha="center",va="center",fontsize=9)
ax.set_title("Posterior mean Gaussian-copula correlation"); fig.colorbar(im, fraction=.046); plt.show()

Gaussian-copula dependence with uncertainty:

pair       post. mean       94% credible
SPY-TLT         -0.02  [-0.09, 0.05]
SPY-GLD          0.12   [0.06, 0.19]
SPY-HYG          0.69   [0.66, 0.73]
SPY-EEM          0.71   [0.68, 0.74]
TLT-GLD          0.36   [0.30, 0.42]
TLT-HYG          0.17   [0.11, 0.24]
TLT-EEM         -0.02  [-0.09, 0.05]
GLD-HYG          0.21   [0.14, 0.28]
GLD-EEM          0.29   [0.23, 0.35]
HYG-EEM          0.59   [0.54, 0.63]
No description has been provided for this image

2. The posterior of tail dependence (Clayton copula)¶

The Gaussian copula cannot express tail dependence at all. To estimate it we fit a Clayton copula — the natural model for joint crashes (lower-tail dependence only) — to the pair with the strongest downside co-movement, SPY–HYG (equity and high-yield credit). Its density is closed-form in a single parameter $\theta>0$, $$c(u,v;\theta)=(1+\theta)(uv)^{-1-\theta}\big(u^{-\theta}+v^{-\theta}-1\big)^{-2-1/\theta},$$ and the lower-tail-dependence coefficient is $\lambda_L=2^{-1/\theta}$. We put a prior on $\theta$, add the Clayton log-likelihood as a pm.Potential, and read off the posterior of $\lambda_L$ — the probability of a joint crash, with honest error bars.

In [4]:
i, j = ASSETS.index("SPY"), ASSETS.index("HYG")
u = np.clip(U[:,i], 1e-4, 1-1e-4); v = np.clip(U[:,j], 1e-4, 1-1e-4)
with pm.Model() as clayton:
    theta = pm.HalfNormal("theta", 3.0)
    logc = (pt.log1p(theta) - (1+theta)*(np.log(u)+np.log(v))
            - (2 + 1/theta)*pt.log(u**(-theta) + v**(-theta) - 1))
    pm.Potential("clayton_ll", logc.sum())
    lamL = pm.Deterministic("lambda_L", 2.0**(-1.0/theta))
    tau  = pm.Deterministic("tau", theta/(theta+2))
    idata_c = pm.sample(1500, tune=1500, chains=4, cores=1, target_accept=0.9,
                        progressbar=False, random_seed=RNG)
print("max R-hat:", float(az.summary(idata_c, var_names=["theta","lambda_L"])["r_hat"].max()))
print(az.summary(idata_c, var_names=["theta","tau","lambda_L"]).to_string())

Initializing NUTS using jitter+adapt_diag...

Sequential sampling (4 chains in 1 job)

NUTS: [theta]

Sampling 4 chains for 1_500 tune and 1_500 draw iterations (6_000 + 6_000 draws total) took 1 seconds.

max R-hat: 1.0
            mean      sd eti89_lb eti89_ub  ess_bulk  ess_tail r_hat mcse_mean  mcse_sd
theta      1.493   0.105      1.3      1.7      1910      2377  1.00    0.0024   0.0016
tau       0.4268  0.0172      0.4     0.45      1910      2377  1.00   0.00039  0.00027
lambda_L  0.6274  0.0206     0.59     0.66      1910      2377  1.00   0.00048  0.00033
In [5]:
lam = idata_c.posterior["lambda_L"].values.flatten()
emp = cp.empirical_tail_dep(U[:,i], U[:,j], 0.10, True)
fig, ax = plt.subplots(1, 2, figsize=(13, 4.3))
ax[0].plot(idata_c.posterior["theta"].values[0][:1200], color=BLUE, lw=.6)
ax[0].set_title("MCMC trace: Clayton $\\theta$"); ax[0].set_xlabel("iteration")
ax[1].hist(lam, bins=50, density=True, color=BLUE, alpha=.75)
ax[1].axvline(lam.mean(), color=RED, lw=2, label="posterior mean %.2f" % lam.mean())
ax[1].axvline(0, color="k", lw=2, ls="--", label="Gaussian ($\\lambda_L=0$)")
ax[1].set_xlabel("lower-tail dependence $\\lambda_L$"); ax[1].set_title("Posterior of SPY–HYG joint-crash probability")
ax[1].legend(); plt.tight_layout(); plt.show()
print("Posterior of lower-tail dependence: mean %.2f, 94%% HDI [%.2f, %.2f]." %
      (lam.mean(), np.percentile(lam,3), np.percentile(lam,97)))
print("P(lambda_L > 0.3 | data) = %.0f%%.  The entire posterior sits far from the Gaussian value of 0," % (100*np.mean(lam>0.3)))
print("so within the Clayton family the data reject 'no joint crashes' overwhelmingly -- and that")
print("qualitative verdict is the robust part. The interval itself is narrow, which is a statement")
print("about theta, not about how well we know the joint-crash probability: the next cell shows the")
print("choice of copula family moves this number several times further than the posterior spread does.")
No description has been provided for this image 
Posterior of lower-tail dependence: mean 0.63, 94% HDI [0.59, 0.66].
P(lambda_L > 0.3 | data) = 100%.  The entire posterior sits far from the Gaussian value of 0,
so within the Clayton family the data reject 'no joint crashes' overwhelmingly -- and that
qualitative verdict is the robust part. The interval itself is narrow, which is a statement
about theta, not about how well we know the joint-crash probability: the next cell shows the
choice of copula family moves this number several times further than the posterior spread does.

How much of that number is the data, and how much is the model?¶

The posterior above is narrow, and it would be easy to read the interval as the uncertainty in the joint-crash probability. It is not — it is the uncertainty given that the copula is Clayton.

Clayton has a single parameter, so it ties the strength of dependence in the body to the strength in the tail. Fitted to 521 weeks that are mostly body, it can only express a heavy tail by inflating $\theta$. The from-scratch notebook's fitted $t$-copula gives $\lambda_L$ for this same pair from two parameters that are free to disagree. Comparing the two — and, more to the point, comparing both against the conditional-crash frequency actually observed at each depth $q$ — settles which number a risk manager should carry.

In [6]:
from scipy.optimize import minimize_scalar
import copulas as cp

# --- the same pair under the t-copula, and both against the observed curve ---
R_t, nu = cp.fit_t_copula(U)
lam_t = cp.tail_dep_t(R_t[i, j], nu)
lam_c = float(np.mean(lam))
th_hat = float(np.mean(idata_c.posterior["theta"].values))

print("Lower-tail dependence for SPY-HYG, lambda_L:")
print("   Clayton posterior mean %.3f   94%% credible [%.3f, %.3f]   (width %.3f)"
      % (lam_c, np.percentile(lam,3), np.percentile(lam,97), np.percentile(lam,97)-np.percentile(lam,3)))
print("   t-copula (from-scratch notebook, nu=%.0f)  %.3f" % (nu, lam_t))
print("   -> the two FAMILIES differ by %.3f; the Clayton posterior spans %.3f."
      % (abs(lam_c-lam_t), np.percentile(lam,97)-np.percentile(lam,3)))
print("      Model choice moves this number about %.0fx more than parameter uncertainty does.\n"
      % (abs(lam_c-lam_t) / (np.percentile(lam,97)-np.percentile(lam,3))))

rng = np.random.default_rng(11); big = 400000
Uc = cp.sim_clayton(th_hat, big, rng)
Ut = cp.sim_t_copula(np.array([[1, R_t[i,j]], [R_t[i,j], 1]]), nu, big, rng)
print("So ask the data which family has the right SHAPE -- P(V<q | U<q) at each depth q:\n")
print("      q      data   Clayton   t-cop")
ec = et = 0.0
for q in (0.20, 0.15, 0.10, 0.05, 0.02):
    d = cp.empirical_tail_dep(u, v, q); c = cp.empirical_tail_dep(Uc[:,0], Uc[:,1], q)
    t = cp.empirical_tail_dep(Ut[:,0], Ut[:,1], q)
    ec += abs(c-d); et += abs(t-d)
    print("   %5.2f    %.2f     %.2f     %.2f" % (q, d, c, t))
print("\n   mean |error|      %.3f     %.3f" % (ec/5, et/5))

nll = lambda th: 1e10 if th <= 0 else -np.sum(
    np.log1p(th) - (1+th)*(np.log(u)+np.log(v)) - (2+1/th)*np.log(u**-th + v**-th - 1))
ll_c = -nll(minimize_scalar(nll, bounds=(0.05, 10), method="bounded").x)
print("\nOn this pair the observed curve DECAYS with depth (%.2f at q=0.20 down to %.2f at q=0.02)."
      % (cp.empirical_tail_dep(u,v,0.20), cp.empirical_tail_dep(u,v,0.02)))
print("Clayton's is flat -- that is its single parameter forcing body and tail to share one number,")
print("and it is why lambda_L comes out high. The t-copula tracks the decay and fits better on this")
print("pair (AIC %.1f vs Clayton %.1f). The %.2f is a Clayton statement, not a data statement."
      % (4-2*cp._t_copula_loglik_given(np.column_stack([u,v]),
                                       np.array([[1,R_t[i,j]],[R_t[i,j],1]]), nu),
         2-2*ll_c, lam_c))

Lower-tail dependence for SPY-HYG, lambda_L:
   Clayton posterior mean 0.627   94% credible [0.586, 0.663]   (width 0.078)
   t-copula (from-scratch notebook, nu=10)  0.185
   -> the two FAMILIES differ by 0.443; the Clayton posterior spans 0.078.
      Model choice moves this number about 6x more than parameter uncertainty does.

So ask the data which family has the right SHAPE -- P(V<q | U<q) at each depth q:

      q      data   Clayton   t-cop
    0.20    0.64     0.65     0.57
    0.15    0.62     0.64     0.53
    0.10    0.52     0.63     0.48
    0.05    0.46     0.63     0.42
    0.02    0.40     0.64     0.36

   mean |error|      0.111     0.057

On this pair the observed curve DECAYS with depth (0.64 at q=0.20 down to 0.40 at q=0.02).
Clayton's is flat -- that is its single parameter forcing body and tail to share one number,
and it is why lambda_L comes out high. The t-copula tracks the decay and fits better on this
pair (AIC -350.2 vs Clayton -321.9). The 0.63 is a Clayton statement, not a data statement.

3. Summary¶

  • A Gaussian copula with an LKJ prior delivers the full posterior over the dependence structure — every pairwise correlation with a credible interval, the Bayesian upgrade of a point estimate.
  • Fitting a Clayton copula to the equity–credit pair yields the posterior of the lower-tail-dependence coefficient $\lambda_L$. Its posterior sits far from the Gaussian value of zero, so the data reject the no-tail-dependence assumption with near-certainty — that verdict is robust.
  • The magnitude is not. The posterior is narrow, but the $t$-copula fitted to the same pair puts $\lambda_L$ several times lower, and it is the $t$ that reproduces the observed conditional-crash curve while Clayton's runs flat. Model uncertainty dwarfs parameter uncertainty here, and a credible interval reports only the second. Quoting the interval as if it bounded the joint-crash probability would be exactly the overconfidence copulas are supposed to cure.
  • Doing copulas Bayesianly turns "how dependent are these assets in a crash?" from a fragile point estimate into a posterior distribution — provided the answer is read as conditional on the model, and the family is checked against the data rather than assumed.

This is the dependence layer, estimated the way the rest of the Risk and Asset Allocation arc estimates everything else — with a prior, a posterior, and error bars that are honest about what they do and do not cover.