Copulas and Tail Dependence — the R engine¶

The independent R implementation, built on the standard copula package (Hofert, Kojadinovic, Mächler, Yan). It fits and simulates the same copulas as the from-scratch Python engine and confirms the results: the Student-$t$ copula beats the Gaussian, and equity–credit pairs carry strong lower-tail dependence that the Gaussian sets to zero. Self-contained; no prior reading required.

Sklar's theorem. Any joint distribution splits into marginals (each asset alone) and a copula (pure dependence): $F(x)=C(F_1(x_1),\dots,F_N(x_N))$. The copula's tail behaviour — do assets crash together? — is the risk-critical feature that a correlation cannot express and a Gaussian copula sets to zero.

In [1]:
suppressMessages(library(copula))
options(repr.plot.width=10, repr.plot.height=4.6)
BLUE<-"#2b6cb0"; ORANGE<-"#dd6b20"; GREEN<-"#2f855a"; RED<-"#c53030"; GREY<-"#718096"; PURP<-"#6b46c1"

X <- as.matrix(read.csv("crossasset_weekly.csv", row.names=1))
A <- colnames(X); N <- ncol(X); U <- pobs(X)          # pseudo-observations (ranks in (0,1))
cat(sprintf("Cross-asset panel: %d weeks, %s\n", nrow(X), paste(A, collapse=", ")))

Warning message:
"package 'copula' was built under R version 4.6.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).

1. The copula zoo¶

Five classic copulas, each calibrated to the same Kendall's $\tau=0.6$, simulated with rCopula. Same rank correlation, completely different tails: Gaussian and Frank empty their corners (no tail dependence), the $t$ fills both, Clayton the lower (crashes), Gumbel the upper (rallies).

In [2]:
tau <- 0.6
zoo <- list(
  "Gaussian"        = normalCopula(iTau(normalCopula(), tau)),
  "Student-t (df=3)"= tCopula(iTau(tCopula(df=3), tau), df=3),
  "Clayton (lower)" = claytonCopula(iTau(claytonCopula(), tau)),
  "Gumbel (upper)"  = gumbelCopula(iTau(gumbelCopula(), tau)),
  "Frank (no tails)"= frankCopula(iTau(frankCopula(), tau)))
set.seed(1)
par(mfrow=c(1,5), mar=c(2,2,2,1))
for (nm in names(zoo)) {
  S <- rCopula(3000, zoo[[nm]])
  plot(S, pch=19, cex=.15, col=BLUE, xlab="", ylab="", main=nm, xaxt="n", yaxt="n")
}
par(mfrow=c(1,1))
cat("Same Kendall's tau, different corners: the tail behaviour is what a single correlation misses.\n")

Same Kendall's tau, different corners: the tail behaviour is what a single correlation misses.
No description has been provided for this image

2. Fitting Gaussian vs Student-t to the data¶

We fit both elliptical copulas to all five assets by maximum pseudo-likelihood and compare by AIC (lower is better). The $t$-copula's extra parameter is the degrees of freedom $\nu$ (tail fatness); if it wins, the data demand tail dependence.

In [3]:
fit_g <- fitCopula(normalCopula(dim=N, dispstr="un"), U, method="mpl")
fit_t <- fitCopula(tCopula(dim=N, dispstr="un"),   U, method="mpl")
nu <- tail(coef(fit_t), 1)
cat(sprintf("Student-t copula dof nu = %.1f  (small nu = fat joint tails)\n\n", nu))
cat(sprintf("  Gaussian copula : logLik %7.1f   AIC %8.1f\n", logLik(fit_g), AIC(fit_g)))
cat(sprintf("  Student-t copula: logLik %7.1f   AIC %8.1f   <- %s\n",
            logLik(fit_t), AIC(fit_t), ifelse(AIC(fit_t) < AIC(fit_g), "WINS", "loses")))

Warning message in var.mpl(copula, u):
"the covariance matrix of the parameter estimates is computed as if 'df.fixed = TRUE' with df = 8.74148609684518"

Student-t copula dof nu = 8.7  (small nu = fat joint tails)


  Gaussian copula : logLik   461.0   AIC   -902.0

  Student-t copula: logLik   478.9   AIC   -935.8   <- WINS

3. Tail dependence: where the Gaussian fails¶

The lower tail-dependence coefficient $\lambda_L=\lim_{q\to0}\Pr(V<q\mid U<q)$ is $0$ for the Gaussian but positive for the $t$. Note that it is a limit: at any measurable threshold both copulas predict far more joint crashes than their own $\lambda_L$, so we report the empirical frequency at $q=0.10$ beside each model's own frequency at $q=0.10$, and keep the $q\to0$ limits in a separate block.

In [4]:
Rt <- p2P(coef(fit_t)[1:(N*(N-1)/2)])              # fitted t-copula correlation matrix
Rg <- p2P(coef(fit_g))
tdep_t <- function(rho, nu) 2*pt(-sqrt((nu+1)*(1-rho)/(1+rho)), nu+1)
emp_tail <- function(u, v, q) { c <- u < q; if (sum(c)==0) NA else mean(v[c] < q) }

# lambda_L is a q->0 LIMIT and cannot be set against an empirical frequency measured at
# q=0.10: 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.
set.seed(3); big <- 200000
cat("Lower-tail dependence by pair:\n")
cat("                  P(V<q | U<q)  at q = 0.10     |   lambda_L  (q -> 0)\n")
cat("  pair            data    t-cop   Gaussian      |   t-cop   Gaussian\n")
for (a in 1:(N-1)) for (b in (a+1):N) {
  Ut_ <- rCopula(big, tCopula(Rt[a,b], dim=2, df=nu))
  Ug_ <- rCopula(big, normalCopula(Rg[a,b], dim=2))
  cat(sprintf("  %-4s-%-4s      %.2f    %.2f     %.2f       |    %.2f     0.00\n",
              A[a], A[b], emp_tail(U[,a], U[,b], 0.10),
              emp_tail(Ut_[,1], Ut_[,2], 0.10), emp_tail(Ug_[,1], Ug_[,2], 0.10),
              tdep_t(Rt[a,b], nu)))
}
sp <- which(A=="SPY"); hy <- which(A=="HYG")
cat(sprintf("\nSPY-HYG: fitted correlation %.2f, t-copula lambda_L %.2f (the Gaussian's limit is 0).\n",
            Rt[sp,hy], tdep_t(Rt[sp,hy], nu)))
cat("Read the blocks against each other: at q=0.10 the two copulas are nearly indistinguishable\n")
cat("and both sit close to the data. The Gaussian's failure is ASYMPTOTIC -- it opens up only as q\n")
cat("shrinks, which is why the limit and the measurable frequency must be reported separately.\n")

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.12     0.10       |    0.01     0.00
  SPY -GLD       0.25    0.18     0.14       |    0.02     0.00
  SPY -HYG       0.52    0.48     0.47       |    0.21     0.00
  SPY -EEM       0.60    0.50     0.49       |    0.23     0.00
  TLT -GLD       0.31    0.27     0.26       |    0.06     0.00
  TLT -HYG       0.29    0.19     0.16       |    0.03     0.00
  TLT -EEM       0.17    0.12     0.10       |    0.01     0.00
  GLD -HYG       0.23    0.21     0.18       |    0.03     0.00
  GLD -EEM       0.25    0.24     0.21       |    0.04     0.00
  HYG -EEM       0.44    0.41     0.39       |    0.15     0.00

SPY-HYG: fitted correlation 0.69, t-copula lambda_L 0.21 (the Gaussian's limit is 0).

Read the blocks against each other: at q=0.10 the two copulas are nearly indistinguishable

and both sit close to the data. The Gaussian's failure is ASYMPTOTIC -- it opens up only as q

shrinks, which is why the limit and the measurable frequency must be reported separately.

4. The risk payoff: joint crashes under the wrong copula¶

Using the same empirical marginals, we draw from the fitted Gaussian and $t$ copulas (rCopula), invert through each asset's empirical quantiles, and compare the equally-weighted portfolio's Expected Shortfall and — the sharp test — given SPY crashes, how many of the other four crash with it.

In [5]:
set.seed(7); m <- 200000
sim_port <- function(fit) {
  Uc <- rCopula(m, fit@copula)
  Xs <- sapply(1:N, function(k) quantile(X[,k], pmin(pmax(Uc[,k],1e-4),1-1e-4), names=FALSE))
  list(port=rowMeans(Xs), Xs=Xs)
}
g <- sim_port(fit_g); t_ <- sim_port(fit_t)
es <- function(r, a=0.01){ v<-quantile(r,a); mean(r[r<=v]) }
thr <- sapply(1:N, function(k) quantile(X[,k], 0.05))
sp <- which(A=="SPY"); oth <- setdiff(1:N, sp)
ncrash  <- function(Xs) rowSums(sweep(Xs, 2, thr, `<`))
cocrash <- function(Xs){ cr<-Xs[,sp]<thr[sp]; mean(rowSums(sweep(Xs[cr,oth],2,thr[oth],`<`))) }
cd <- cocrash(X); ct <- cocrash(t_$Xs); cg <- cocrash(g$Xs)

cat(sprintf("1%% Expected Shortfall:  Gaussian %.2f%%   t-copula %.2f%%\n", es(g$port), es(t_$port)))
cat(sprintf("\nJoint tail events, every row anchored to the %d observed weeks:\n", nrow(X)))
cat(sprintf("  P(>=3 of 5 crash together): DATA %.2f%%   Gaussian %.2f%%   t-copula %.2f%%\n",
    100*mean(ncrash(X)>=3), 100*mean(ncrash(g$Xs)>=3), 100*mean(ncrash(t_$Xs)>=3)))
cat(sprintf("  P(all 5 crash together)   : DATA %.2f%%   Gaussian %.2f%%   t-copula %.2f%%\n",
    100*mean(ncrash(X)==N), 100*mean(ncrash(g$Xs)==N), 100*mean(ncrash(t_$Xs)==N)))
cat(sprintf("\nGiven SPY crashes, E[# of the other 4 also crashing]  (%d observed SPY-crash weeks):\n",
    sum(X[,sp]<thr[sp])))
cat(sprintf("   DATA %.2f  |  t-copula %.2f (%+.0f%%)  |  Gaussian %.2f (%+.0f%%)\n",
    cd, ct, 100*(ct-cd)/cd, cg, 100*(cg-cd)/cd))
cat("Both copulas understate the observed contagion; the t-copula closes roughly half the\n")
cat("Gaussian's shortfall, so it is the better of the two rather than a match for the data.\n")
# joint-crash count distribution
kk <- 0:N
dg <- sapply(kk, function(k) mean(ncrash(g$Xs)==k))
dt <- sapply(kk, function(k) mean(ncrash(t_$Xs)==k))
barplot(rbind(dg,dt), beside=TRUE, names.arg=kk, col=c(GREY,BLUE), log="y",
        xlab="# assets crashing in the same week (of 5)", ylab="probability (log)",
        main="Joint crashes: t-copula has the fatter cluster")
legend("topright", c("Gaussian","Student-t"), fill=c(GREY,BLUE), bty="n")

1% Expected Shortfall:  Gaussian -4.52%   t-copula -4.72%

Joint tail events, every row anchored to the 521 observed weeks:

  P(>=3 of 5 crash together): DATA 2.11%   Gaussian 1.55%   t-copula 1.92%

  P(all 5 crash together)   : DATA 0.19%   Gaussian 0.03%   t-copula 0.08%

Given SPY crashes, E[# of the other 4 also crashing]  (26 observed SPY-crash weeks):

   DATA 1.19  |  t-copula 1.08 (-9%)  |  Gaussian 0.95 (-21%)

Both copulas understate the observed contagion; the t-copula closes roughly half the

Gaussian's shortfall, so it is the better of the two rather than a match for the data.
No description has been provided for this image

5. Summary¶

  • The R copula package reproduces the from-scratch Python results: the Student-$t$ copula beats the Gaussian on AIC (dof $\nu\approx9$), and equity–credit pairs (SPY–HYG, SPY–EEM) carry strong lower-tail dependence that the Gaussian sets to zero.
  • With identical marginals, swapping in the $t$-copula fattens the portfolio's joint-crash tail. Measured against the observed weeks the Gaussian understates contagion — and so, by less, does the $t$.
  • Three engines — from-scratch Python, Bayesian PyMC, and R — agree: dependence in the tails is real, and modelling it requires a copula, not a correlation.

This adds the missing dependence-modelling layer to the collection and connects to the fat-tailed marginals of the Student-$t$ and stochastic-volatility work elsewhere, pointing toward copula-GARCH for time-varying dependence.