Bayesian Copula-GARCH — the PyMC engine

The from-scratch notebook fit the residual dependence by maximum likelihood and observed that GARCH reduces the apparent tail dependence — some of the raw joint-crash signal was merely shared volatility. But is that reduction real, or could it be noise? This notebook answers Bayesianly. We keep the GARCH marginals from maximum likelihood (GARCH in PyMC is impractical on this platform) and fit the copula in PyMC, delivering the posterior of the tail-dependence coefficient for the raw returns and for the GARCH residuals — and hence the posterior of the difference.

Self-contained; no prior reading required.

The question, made Bayesian

The tail-dependence coefficient $\lambda_L$ is the probability two assets crash together. We model the equity–credit pair SPY–HYG as bivariate Student-$t$ (whose dependence is the $t$-copula from the from-scratch notebook), fitting the correlation $\rho$ and the degrees of freedom $\nu$ — and the tail dependence follows, $\lambda_L=2\,t_{\nu+1}\!\big(-\sqrt{(\nu+1)(1-\rho)/(1+\rho)}\big)$. We fit it once to the raw returns and once to the GARCH residuals, and compare the posteriors of $\lambda_L$. If $\lambda_L^{\text{resid}}$ sits credibly below $\lambda_L^{\text{raw}}$, GARCH has genuinely reallocated joint-crash risk from "shared volatility" to "true dependence".

import numpy as np, pandas as pd, matplotlib.pyplot as plt
import pymc as pm, arviz as az, pytensor.tensor as pt
import copulagarch as cg, 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 = 4

R = pd.read_csv("crossasset_daily.csv", index_col=0, parse_dates=True)
ASSETS = list(R.columns); X = R.values; N = X.shape[1]
marg = cg.fit_all(X)                                    # GARCH(1,1)-t marginals (MLE)
Z = cg.std_resid_matrix(marg)                           # standardized residuals (i.i.d.)
print("GARCH marginals fitted; persistence:", {a: round(m["persistence"],3) for a,m in zip(ASSETS,marg)})

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

GARCH marginals fitted; persistence: {'SPY': 0.99, 'TLT': 0.986, 'GLD': 0.989, 'HYG': 1.0, 'EEM': 0.97}

The dataset

Daily log-returns (%) of five cross-asset ETFs, 2010–2024 (SPY, TLT, GLD, HYG, EEM). We reduce each to its GARCH(1,1)-$t$ standardized residuals (volatility removed) and work with the pseudo-observations of both the raw returns and the residuals for the pair SPY–HYG.

1. Bayesian bivariate-$t$: raw returns vs GARCH residuals

We fit a bivariate Student-$t$ to the standardized pair, with priors $\rho\sim\text{Uniform}(-0.95,0.95)$ and $\nu\sim\text{Gamma}$ (mean $\approx10$). The degrees of freedom $\nu$ control the joint tail fatness: a small $\nu$ means the two assets have violent joint extremes. From each posterior draw of $(\rho,\nu)$ we compute the $t$-copula tail dependence $\lambda_L$ in closed form.

from scipy.stats import t as st
def fit_bivariate_t(x, y, seed):
    d = np.column_stack([(x-x.mean())/x.std(), (y-y.mean())/y.std()])
    with pm.Model():
        rho = pm.Uniform("rho", -0.95, 0.95)
        nu  = pm.Gamma("nu", 3.0, 0.3)                          # dof, mean ~10
        chol = pt.stack([pt.stack([pt.ones(()), pt.zeros(())]),
                         pt.stack([rho, pt.sqrt(1 - rho**2)])])  # correlation Cholesky
        pm.MvStudentT("obs", nu=nu, mu=np.zeros(2), chol=chol, observed=d)
        idata = pm.sample(1500, tune=1500, chains=4, cores=1, target_accept=0.9,
                          progressbar=False, random_seed=seed)
    rd = idata.posterior["rho"].values.flatten(); nd = idata.posterior["nu"].values.flatten()
    lam = 2 * st.cdf(-np.sqrt((nd+1)*(1-rd)/(1+rd)), nd+1)       # t-copula lower-tail dependence
    return rd, nd, lam

i, j = ASSETS.index("SPY"), ASSETS.index("HYG")
rho_raw, nu_raw, lam_raw = fit_bivariate_t(X[:,i], X[:,j], RNG)
rho_res, nu_res, lam_res = fit_bivariate_t(Z[:,i], Z[:,j], RNG+1)
print("SPY-HYG bivariate-t posteriors (mean [94%% HDI]):")
print("  raw returns    : rho %.2f, nu %.1f  ->  tail dep lambda_L %.2f [%.2f, %.2f]"
      % (rho_raw.mean(), nu_raw.mean(), lam_raw.mean(), np.percentile(lam_raw,3), np.percentile(lam_raw,97)))
print("  GARCH residuals: rho %.2f, nu %.1f  ->  tail dep lambda_L %.2f [%.2f, %.2f]"
      % (rho_res.mean(), nu_res.mean(), lam_res.mean(), np.percentile(lam_res,3), np.percentile(lam_res,97)))
print("\nThe residual nu is much larger (thinner joint tails) -> lambda_L drops, matching the from-scratch fit.")

Initializing NUTS using jitter+adapt_diag...

Sequential sampling (4 chains in 1 job)

NUTS: [rho, nu]

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

Initializing NUTS using jitter+adapt_diag...

Sequential sampling (4 chains in 1 job)

NUTS: [rho, nu]

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

SPY-HYG bivariate-t posteriors (mean [94%% HDI]):
  raw returns    : rho 0.85, nu 5.7  ->  tail dep lambda_L 0.49 [0.46, 0.52]
  GARCH residuals: rho 0.78, nu 10.9  ->  tail dep lambda_L 0.25 [0.20, 0.29]

The residual nu is much larger (thinner joint tails) -> lambda_L drops, matching the from-scratch fit.

# posterior of the REDUCTION (pair independent draws)
n = min(len(lam_raw), len(lam_res))
diff = lam_raw[:n] - lam_res[:n]
fig, ax = plt.subplots(1, 2, figsize=(13, 4.3))
ax[0].hist(lam_raw, bins=50, density=True, color=GREY, alpha=.6, label="raw returns")
ax[0].hist(lam_res, bins=50, density=True, color=BLUE, alpha=.6, label="GARCH residuals")
ax[0].set_xlabel("lower-tail dependence $\\lambda_L$"); ax[0].set_title("Posterior of joint-crash probability"); ax[0].legend()
ax[1].hist(diff, bins=50, density=True, color=GREEN, alpha=.8)
ax[1].axvline(0, color="k", lw=2, ls="--"); ax[1].axvline(diff.mean(), color=RED, lw=2, label="mean %.2f" % diff.mean())
ax[1].set_xlabel("$\\lambda_L^{raw} - \\lambda_L^{resid}$"); ax[1].set_title("Posterior of the reduction from GARCH"); ax[1].legend()
plt.tight_layout(); plt.show()
print("P(lambda_L is lower for residuals than raw) = %.0f%%." % (100*np.mean(diff > 0)))
print("The reduction is credible, not noise: removing common volatility genuinely lowers the estimated")
print("joint-crash probability -- with the whole posterior on the positive side. A single ML number could")
print("not tell you whether that drop was real; the Bayesian comparison can.")

P(lambda_L is lower for residuals than raw) = 100%.
The reduction is credible, not noise: removing common volatility genuinely lowers the estimated
joint-crash probability -- with the whole posterior on the positive side. A single ML number could
not tell you whether that drop was real; the Bayesian comparison can.

2. The residual dependence structure, with uncertainty

Beyond the single pair, the whole residual dependence can be given error bars: a Gaussian copula with an LKJ prior on the residuals returns the posterior correlation of all five assets — the clean, volatility-adjusted co-movement that feeds a risk model.

from scipy.stats import norm
Uz = cp.pseudo_obs(Z)                       # pseudo-observations of the GARCH residuals
Zc = norm.ppf(np.clip(Uz, 1e-4, 1-1e-4))
with pm.Model() as gc:
    chol, corr, sds = pm.LKJCholeskyCov("L", n=N, eta=2.0, sd_dist=pm.HalfNormal.dist(1.0), compute_corr=True)
    pm.Deterministic("Rc", corr)
    pm.MvNormal("z", mu=np.zeros(N), chol=chol, observed=Zc)
    idata_g = pm.sample(800, tune=1000, chains=4, cores=1, target_accept=0.9, progressbar=False, random_seed=RNG)
Rc = idata_g.posterior["Rc"].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["Rc"].values[:, :, iu[0], iu[1]]
print("max R-hat (off-diagonal):", float(az.rhat(az.convert_to_dataset(off))["x"].max()))

fig, ax = plt.subplots(figsize=(5.2,4.4))
im = ax.imshow(Rc, 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"%Rc[a,b],ha="center",va="center",fontsize=9)
ax.set_title("Posterior residual-copula correlation"); fig.colorbar(im, fraction=.046); plt.show()

print("Residual-copula dependence, posterior mean and 94% credible interval for every pair:\n")
print("  pair        post. mean        94% credible     width")
widths = {}
for k, (a, b) in enumerate(zip(iu[0], iu[1])):
    d = off[:, :, k].flatten(); lo, hi = np.percentile(d, 3), np.percentile(d, 97)
    widths["%s-%s" % (ASSETS[a], ASSETS[b])] = hi - lo
    print("  %-4s-%-4s   %11.2f   [%5.2f, %5.2f]     %.3f" % (ASSETS[a], ASSETS[b], d.mean(), lo, hi, hi-lo))
tight = min(widths, key=widths.get); wide = max(widths, key=widths.get)
print("\nThis is the dependence AFTER stripping each asset's own volatility dynamics -- the honest")
print("co-movement structure, delivered as a posterior with a credible interval on every entry.")
print("The intervals are not uniform: %s is pinned to within %.3f while %s spans %.3f, so the"
      % (tight, widths[tight], wide, widths[wide]))
print("residual dependence structure is known far more precisely for some pairs than for others --")
print("a distinction a single point-estimate correlation matrix cannot express at all.\n")
k_st = [k for k,(a,b) in enumerate(zip(iu[0],iu[1]))
        if {ASSETS[a],ASSETS[b]} == {"SPY","TLT"}][0]
d_st = off[:,:,k_st].flatten()
print("Note SPY-TLT: posterior mean %+.2f, entirely below zero [%+.2f, %+.2f]. Bonds hedge equities"
      % (d_st.mean(), np.percentile(d_st,3), np.percentile(d_st,97)))
print("across this sample and the posterior is emphatic about it. But that is an average over")
print("2010-2024, and the from-scratch notebook shows the relationship CHANGED SIGN in 2023.")
print("A posterior this tight around a parameter that is not constant is confident about the")
print("wrong thing -- precision on the average, silence on the drift. It is the same lesson the")
print("copulas notebook drew about model choice, arriving here from the other direction.")

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 10 seconds.

max R-hat (off-diagonal): 1.0018135444613485

Residual-copula dependence, posterior mean and 94% credible interval for every pair:

  pair        post. mean        94% credible     width
  SPY -TLT          -0.28   [-0.31, -0.25]     0.057
  SPY -GLD           0.05   [ 0.02,  0.08]     0.063
  SPY -HYG           0.73   [ 0.71,  0.74]     0.029
  SPY -EEM           0.74   [ 0.73,  0.75]     0.029
  TLT -GLD           0.26   [ 0.23,  0.29]     0.058
  TLT -HYG          -0.06   [-0.09, -0.03]     0.061
  TLT -EEM          -0.21   [-0.24, -0.18]     0.059
  GLD -HYG           0.14   [ 0.11,  0.17]     0.060
  GLD -EEM           0.19   [ 0.16,  0.22]     0.060
  HYG -EEM           0.64   [ 0.62,  0.66]     0.036

This is the dependence AFTER stripping each asset's own volatility dynamics -- the honest
co-movement structure, delivered as a posterior with a credible interval on every entry.
The intervals are not uniform: SPY-EEM is pinned to within 0.029 while SPY-GLD spans 0.063, so the
residual dependence structure is known far more precisely for some pairs than for others --
a distinction a single point-estimate correlation matrix cannot express at all.

Note SPY-TLT: posterior mean -0.28, entirely below zero [-0.31, -0.25]. Bonds hedge equities
across this sample and the posterior is emphatic about it. But that is an average over
2010-2024, and the from-scratch notebook shows the relationship CHANGED SIGN in 2023.
A posterior this tight around a parameter that is not constant is confident about the
wrong thing -- precision on the average, silence on the drift. It is the same lesson the
copulas notebook drew about model choice, arriving here from the other direction.

3. Summary

• Keeping GARCH marginals from maximum likelihood, a Bayesian copula on the residuals turns the from-scratch notebook's observation into an inference: the posterior of the tail-dependence reduction from GARCH sits credibly above zero — the drop in apparent joint-crash risk is real, not sampling noise.
• An LKJ Gaussian copula gives the full residual dependence structure with error bars — the volatility-adjusted co-movement, ready to feed a risk model.

Copula-GARCH, done Bayesianly, separates and quantifies the two moods of the market — common turbulence and true dependence — with the honesty about uncertainty that runs through the whole Risk and Asset Allocation arc.