Coherent Risk: VaR, ES & Extreme-Value Theory — the R engine¶

The independent R implementation. VaR and Expected Shortfall are coded from scratch (they are one-liners), and the tail is fit with the classic evir package (Embrechts–McNeil GPD peaks-over-threshold) as the authoritative cross-check. It reproduces the from-scratch Python results: the Gaussian understates the tail, Cornish–Fisher explodes at high kurtosis, EVT extrapolates the far tail, VaR fails sub-additivity while ES is coherent, and a static VaR fails its backtest through clustered violations. Self-contained; no prior reading required.

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

R <- read.csv("crossasset_daily.csv", row.names=1)
A <- colnames(R); X <- as.matrix(R); spy <- R$SPY
cat(sprintf("Daily cross-asset panel: %d days, %s\n", nrow(R), paste(A, collapse=", ")))
cat(sprintf("SPY daily: mean %.3f  sd %.3f  worst %.2f  excess-kurtosis %.1f\n",
            mean(spy), sd(spy), min(spy), mean(((spy-mean(spy))/sd(spy))^4)-3))

Daily cross-asset panel: 3772 days, SPY, TLT, GLD, HYG, EEM

SPY daily: mean 0.051  sd 1.077  worst -11.59  excess-kurtosis 11.5

The dataset¶

Daily log-returns (%) of five cross-asset ETFs, 2010–2024 (3,772 days): SPY, TLT, GLD, HYG, EEM. Daily frequency gives the depth the tail needs; the five assets give a portfolio for coherence and contributions. SPY's worst day is the −11.6% COVID crash.

1. VaR and Expected Shortfall, four ways¶

VaR$_\alpha$ = minus the $\alpha$-quantile of returns; ES$_\alpha$ = the average loss beyond it. Historical (empirical), Gaussian (normal), Cornish–Fisher (skew/kurtosis-adjusted), and Student-$t$ (fat-tailed fit).

In [2]:
skew_ <- function(x) mean(((x-mean(x))/sd(x))^3); exk_ <- function(x) mean(((x-mean(x))/sd(x))^4)-3
hist_var <- function(r,a) -quantile(r,a,names=FALSE)
hist_es  <- function(r,a){ s<-sort(r); k<-max(1,ceiling(a*length(r))); -mean(s[1:k]) }
norm_var <- function(r,a) -(mean(r)+sd(r)*qnorm(a));  norm_es <- function(r,a) -(mean(r)-sd(r)*dnorm(qnorm(a))/a)
cf_z <- function(a,s,k){ z<-qnorm(a); z+(z^2-1)*s/6+(z^3-3*z)*k/24-(2*z^3-5*z)*s^2/36 }
cf_var <- function(r,a) -(mean(r)+sd(r)*cf_z(a,skew_(r),exk_(r)))

tab <- t(sapply(c(0.05,0.01,0.001), function(a) c(alpha=a,
  hist=hist_var(spy,a), norm=norm_var(spy,a), CF=cf_var(spy,a),
  ES.hist=hist_es(spy,a), ES.norm=norm_es(spy,a))))
cat("SPY daily VaR / ES (positive = loss %):\n"); print(round(tab,2))
cat("\nGaussian understates every tail; Cornish-Fisher explodes at 0.1% (excess-kurtosis ~12 breaks it);\n")
cat("historical is the trustworthy benchmark. Same pattern as the Python engine.\n")

SPY daily VaR / ES (positive = loss %):

     alpha hist norm    CF ES.hist ES.norm
[1,]  0.05 1.67 1.72  1.68    2.66    2.17
[2,]  0.01 3.10 2.46  5.72    4.48    2.82
[3,]  0.00 6.12 3.28 14.17    9.13    3.58

Gaussian understates every tail; Cornish-Fisher explodes at 0.1% (excess-kurtosis ~12 breaks it);

historical is the trustworthy benchmark. Same pattern as the Python engine.

2. Extreme-Value Theory with evir¶

The evir package fits the Generalized Pareto Distribution to peaks over a high threshold and returns VaR (quantile) and ES (sfall) at any level. We select the threshold with the mean-excess plot and read off the far-tail risk that history alone cannot.

In [3]:
L <- -spy; u <- quantile(L, 0.90)
fit <- gpd(L, threshold=u)
cat(sprintf("evir GPD fit: xi = %.3f (>0 heavy tail), beta = %.3f, %d exceedances over u=%.2f\n\n",
            fit$par.ests["xi"], fit$par.ests["beta"], length(fit$data), u))
rmv <- riskmeasures(fit, c(0.99, 0.999))          # 99% and 99.9% -> VaR & ES
cat(sprintf("evir riskmeasures:  99%% VaR %.2f  ES %.2f  |  99.9%% VaR %.2f  ES %.2f\n\n",
    rmv[1,2], rmv[1,3], rmv[2,2], rmv[2,3]))
par(mfrow=c(1,2))
meplot(L, main="Mean-excess plot (linear -> GPD tail)"); abline(v=u, col=RED, lty=2)
# VaR term structure: historical vs Gaussian vs EVT
alphas <- 10^seq(-3.3,-1.3,length.out=25)
vh<-sapply(alphas,function(a)hist_var(spy,a)); vn<-sapply(alphas,function(a)norm_var(spy,a))
ve<-riskmeasures(fit, 1-alphas)[,2]
plot(alphas, ve, type="l", col=GREEN, lwd=3, log="x", xlim=rev(range(alphas)),
     xlab="tail probability alpha (rarer ->)", ylab="VaR (loss %)", main="Into the far tail")
lines(alphas, vh, type="b", pch=19, col=BLUE); lines(alphas, vn, type="b", pch=15, col=RED)
legend("topleft", c("EVT (GPD)","historical","Gaussian"), col=c(GREEN,BLUE,RED), lwd=c(3,1,1), pch=c(NA,19,15), bty="n")
par(mfrow=c(1,1))
cat(sprintf("At 0.1%%: historical %.2f  Gaussian %.2f  EVT %.2f -- EVT extrapolates where data runs out.\n",
            hist_var(spy,0.001), norm_var(spy,0.001), rmv[2,2]))

evir GPD fit: xi = 0.155 (>0 heavy tail), beta = 0.763, 378 exceedances over u=1.10


evir riskmeasures:  99% VaR 3.21  ES 4.50  |  99.9% VaR 6.23  ES 8.07


At 0.1%: historical 6.12  Gaussian 3.28  EVT 6.23 -- EVT extrapolates where data runs out.
No description has been provided for this image

3. Coherence: VaR fails sub-additivity, ES does not¶

A coherent risk measure is sub-additive — diversifying never raises risk. VaR breaks this; ES does not. Two independent defaultable bonds (rare large loss, small coupon): each bond's default is rarer than the VaR level so its VaR looks safe, but combined they push a loss into the VaR window — the diversified book reports a larger VaR.

In [4]:
set.seed(0); n <- 200000
bond <- function(p=0.04) ifelse(runif(n)<p, -100, 2)
Aa <- bond(); Bb <- bond(); port <- 0.5*Aa + 0.5*Bb; al <- 0.05
vA<-hist_var(Aa,al); vB<-hist_var(Bb,al); vP<-hist_var(port,al)
eA<-hist_es(Aa,al);  eB<-hist_es(Bb,al);  eP<-hist_es(port,al)
par(mfrow=c(1,2))
barplot(c(vA+vB, vP), names.arg=c("VaR(A)+VaR(B)","VaR(A+B)"), col=c(GREY,RED),
        main="VaR FAILS sub-additivity", ylab="VaR")
barplot(c(eA+eB, eP), names.arg=c("ES(A)+ES(B)","ES(A+B)"), col=c(GREY,GREEN),
        main="ES is coherent", ylab="ES")
par(mfrow=c(1,1))
cat(sprintf("VaR: sum %.1f  vs  diversified %.1f  -> %s\n", vA+vB, vP, ifelse(vP>vA+vB,"VIOLATED (diversification punished!)","ok")))
cat(sprintf("ES : sum %.1f  vs  diversified %.1f  -> coherent (%s)\n", eA+eB, eP, ifelse(eP<=eA+eB,"yes","no")))

VaR: sum -4.0  vs  diversified 49.0  -> VIOLATED (diversification punished!)

ES : sum 161.3  vs  diversified 50.5  -> coherent (yes)
No description has been provided for this image

4. Risk contributions and backtesting¶

Component VaR (Euler) splits portfolio risk into per-asset pieces summing to the total. And a VaR is only trustworthy if it survives a backtest: the Kupiec test checks the violation rate, Christoffersen checks independence.

In [5]:
# component VaR (Gaussian Euler) for the equal-weight portfolio
w <- rep(1/length(A), length(A)); mu <- colMeans(X); Sig <- cov(X)
sdp <- sqrt(as.numeric(t(w)%*%Sig%*%w)); cvar <- w * (-mu - qnorm(0.01)*as.numeric(Sig%*%w)/sdp)
cat("Component VaR (% of total):", paste(sprintf("%s %.0f%%", A, 100*cvar/sum(cvar)), collapse="  "), "\n\n")

# rolling 1% VaR backtest
kupiec <- function(x,n,a){ pi<-x/n; ll0<-(n-x)*log(1-a)+x*log(a); ll1<-if(pi>0&&pi<1)(n-x)*log(1-pi)+x*log(pi) else 0
  LR<--2*(ll0-ll1); c(LR=LR, p=1-pchisq(LR,1)) }
christo <- function(v){ n00<-sum(v[-length(v)]==0&v[-1]==0);n01<-sum(v[-length(v)]==0&v[-1]==1)
  n10<-sum(v[-length(v)]==1&v[-1]==0);n11<-sum(v[-length(v)]==1&v[-1]==1)
  p01<-n01/max(n00+n01,1);p11<-n11/max(n10+n11,1);p<-(n01+n11)/max(length(v)-1,1)
  s<-function(a,b,pp) if(pp>0&&pp<1) a*log(1-pp)+b*log(pp) else 0
  LR<--2*((s(n00,n01,p)+s(n10,n11,p))-(s(n00,n01,p01)+s(n10,n11,p11))); c(LR=LR,p=1-pchisq(LR,1)) }
bt <- function(r,W,method,a=0.01){ v<-c(); for(t in (W+1):length(r)){ past<-r[(t-W):(t-1)]
  vv<-if(method=="normal") norm_var(past,a) else hist_var(past,a); v<-c(v, as.integer(r[t] < -vv)) }; v }
cat("Rolling 1% VaR backtest (window=500):\n")
for(meth in c("normal","historical")){ v<-bt(spy,500,meth); k<-kupiec(sum(v),length(v),0.01); ch<-christo(v)
  cat(sprintf("  %-11s %4d viol (%.2f%%)  Kupiec p=%.3f  Christoffersen p=%.3f\n",
              meth, sum(v), 100*mean(v), k["p"], ch["p"])) }
cat("\nGaussian breached too often (fails Kupiec); historical rate is fine but BOTH cluster violations\n")
cat("(fail Christoffersen) -> volatility clustering, which only a dynamic GARCH-VaR removes.\n")

Component VaR (% of total): SPY 27%  TLT 4%  GLD 17%  HYG 14%  EEM 38% 


Rolling 1% VaR backtest (window=500):

  normal        82 viol (2.51%)  Kupiec p=0.000  Christoffersen p=0.000
  historical    43 viol (1.31%)  Kupiec p=0.085  Christoffersen p=0.000

Gaussian breached too often (fails Kupiec); historical rate is fine but BOTH cluster violations

(fail Christoffersen) -> volatility clustering, which only a dynamic GARCH-VaR removes.

5. Summary¶

  • From-scratch VaR/ES in R reproduce the Python engine: the Gaussian understates the tail, Cornish–Fisher explodes at equity kurtosis, historical is the reliable benchmark.
  • evir confirms the tail: a Generalized-Pareto fit with heavy shape $\xi\approx0.16$ extrapolates the far-tail VaR/ES that history cannot.
  • VaR fails sub-additivity (punishes diversification); ES is coherent — the reason regulation moved to ES.
  • A static VaR fails its backtest through clustered violations (volatility clustering) — a pointer to dynamic GARCH-VaR.

Three engines — from-scratch Python, Bayesian PyMC (posterior of the far-tail VaR), and R (evir) — agree, closing the risk-measurement chapter of the Asset Risk arc.