Coherent Risk: VaR, Expected Shortfall & Extreme-Value Theory

The risk-measurement engine of the Asset Risk arc — from scratch

"Value-at-Risk is the most widely used risk number in finance, and it is not even a coherent one: it can punish you for diversifying."

Self-contained; no prior reading of Meucci's Risk and Asset Allocation (Springer 2005, ch. 5) required.

What makes a risk measure good?

Everything in this arc — shrinkage, Bayesian estimation, allocation, copulas, horizon projection — ultimately feeds a single question: how bad can it get? The answer is a risk measure, and the two that dominate practice are:

• Value-at-Risk (VaR) — the loss you exceed with probability $\alpha$ (e.g. the 1-in-100 daily loss). Simple, intuitive, regulator-mandated — and quietly flawed: it ignores how far past the threshold losses go, and it is not coherent — it can report that a diversified portfolio is riskier than its parts, violating the one property any sane risk measure should have.
• Expected Shortfall (ES / CVaR) — the average loss given that VaR is breached. It looks into the tail, and it is coherent (sub-additive), so it never penalises diversification. Basel III moved bank capital from VaR to ES for exactly this reason.

And a third problem cuts across both: the far tail (1-in-1000, 1-in-10000 days) is where the scary numbers live, but it is exactly where we have almost no data. Extreme-Value Theory is the mathematics of extrapolating there safely.

Roadmap

1. VaR & ES and four ways to estimate them — historical, Gaussian, Cornish–Fisher, and where each fails.
2. Extreme-Value Theory (peaks-over-threshold) for the far tail.
3. Coherence — VaR's sub-additivity failure, ES's coherence.
4. Risk contributions — which assets actually drive portfolio risk.
5. Backtesting — does the VaR hold up out of sample?

import numpy as np, pandas as pd, matplotlib.pyplot as plt
import riskmeasures as rm

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"
np.set_printoptions(precision=3, suppress=True)

R = pd.read_csv("crossasset_daily.csv", index_col=0, parse_dates=True)
ASSETS = list(R.columns); X = R.values; spy = R["SPY"].values
print("Daily cross-asset panel:", X.shape, "->", ASSETS)
print("SPY daily: mean %.3f%%  std %.3f%%  worst day %.2f%%  skew %.2f  ex-kurt %.1f"
      % (spy.mean(), spy.std(), spy.min(), pd.Series(spy).skew(), pd.Series(spy).kurt()))

Daily cross-asset panel: (3772, 5) -> ['SPY', 'TLT', 'GLD', 'HYG', 'EEM']
SPY daily: mean 0.051%  std 1.077%  worst day -11.59%  skew -0.72  ex-kurt 11.5

The dataset

Daily log-returns (%) of five cross-asset ETFs, 2010–2024 (3,772 days): SPY (equity), TLT (Treasuries), GLD (gold), HYG (credit), EEM (EM equity). Daily frequency gives us the many observations the far tail needs, and the five assets give us a portfolio for coherence and risk-contribution analysis. SPY's worst day here is the −11.6% COVID crash of March 2020 — a reminder of what the tail contains.

1. VaR and ES, and four ways to estimate them

VaR$_\alpha$ is minus the $\alpha$-quantile of returns (a positive loss); ES$_\alpha$ is the average loss beyond it, $\text{ES}_\alpha=\mathbb E[-r\mid -r>\text{VaR}_\alpha]$. Four estimators, in increasing sophistication:

• Historical — just read the empirical quantile / tail average off the data. Assumption-free, but only as reliable as the data is deep.
• Gaussian — assume normal returns: $\text{VaR}=-(\mu+\sigma z_\alpha)$. Fast, and wrong in the tail — it ignores fat tails and understates risk.
• Cornish–Fisher — correct the Gaussian quantile using the sample skewness and kurtosis. Better for mild non-normality — but it diverges for very fat tails, as we will see.
• Student-$t$ — fit a fat-tailed $t$ by maximum likelihood; a parametric tail that actually has fat tails.

rows = []
for a in (0.05, 0.01, 0.001):
    tv, te, _ = rm.t_fit_var_es(spy, a)
    rows.append([a, rm.historical_var(spy,a), rm.normal_var(spy,a), rm.cornish_fisher_var(spy,a), tv,
                 rm.historical_es(spy,a), rm.normal_es(spy,a), te])
df = pd.DataFrame(rows, columns=["alpha","VaR hist","VaR norm","VaR CF","VaR t","ES hist","ES norm","ES t"]).set_index("alpha")
print("SPY daily VaR / ES (positive = loss %):\n"); print(df.round(2).to_string())
print("\nThe Gaussian understates every tail; Cornish-Fisher is reasonable at 5% but EXPLODES at 0.1%")
print("(CF is only valid for mild non-normality -- SPY's excess kurtosis ~12 breaks it). The Student-t")
print("and historical estimates are the trustworthy ones -- and they diverge from the Gaussian fast.")

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

       VaR hist  VaR norm  VaR CF  VaR t  ES hist  ES norm   ES t
alpha                                                            
0.050      1.67      1.72    1.68   1.46     2.66     2.17   2.58
0.010      3.10      2.46    5.72   3.03     4.48     2.82   4.94
0.001      6.12      3.28   14.18   7.40     9.13     3.58  11.77

The Gaussian understates every tail; Cornish-Fisher is reasonable at 5% but EXPLODES at 0.1%
(CF is only valid for mild non-normality -- SPY's excess kurtosis ~12 breaks it). The Student-t
and historical estimates are the trustworthy ones -- and they diverge from the Gaussian fast.

a = 0.01
fig, ax = plt.subplots(figsize=(11,4.8))
ax.hist(spy, bins=200, density=True, color=GREY, alpha=.6)
for name, v, c in [("historical", rm.historical_var(spy,a), BLUE),
                   ("Gaussian",   rm.normal_var(spy,a),     RED),
                   ("Student-t",  rm.t_fit_var_es(spy,a)[0], GREEN)]:
    ax.axvline(-v, color=c, lw=2, label="VaR$_{1\\%%}$ %s = %.2f%%" % (name, v))
ax.axvline(-rm.historical_es(spy,a), color=PURP, lw=2, ls="--", label="ES$_{1\\%%}$ hist = %.2f%%" % rm.historical_es(spy,a))
ax.set_xlim(-8, 4); ax.set_xlabel("SPY daily return (%)"); ax.set_ylabel("density")
ax.set_title("1% VaR and ES on SPY: the Gaussian sits too close to the centre"); ax.legend()
plt.show()
print("ES (dashed) sits deeper than any VaR: it is the average of the losses in the shaded left tail,")
print("not merely the threshold. That is the information VaR throws away -- and why ES is more prudent.")

ES (dashed) sits deeper than any VaR: it is the average of the losses in the shaded left tail,
not merely the threshold. That is the information VaR throws away -- and why ES is more prudent.

2. Extreme-Value Theory: the far tail done right

At $\alpha=0.001$ (a 1-in-1000-day loss) the historical estimate rests on only a handful of observations, and the Gaussian is hopeless. Extreme-Value Theory provides the principled extrapolation. The peaks-over-threshold approach uses a deep theorem: the exceedances of any distribution over a high threshold converge to a Generalized Pareto Distribution (GPD), with a shape parameter $\xi$ that measures tail heaviness ($\xi>0$ = heavy/power-law tail). We pick a threshold with the mean-excess plot (linear-in-threshold signals the GPD regime), fit the GPD to the exceedances, and read off VaR and ES at any tail probability.

L = -spy                                    # losses
us = np.linspace(np.quantile(L,0.80), np.quantile(L,0.995), 40)
me = rm.mean_excess(spy, us)
fit = rm.fit_gpd_pot(spy, thr_q=0.90)
fig, ax = plt.subplots(1, 2, figsize=(13, 4.4))
ax[0].plot(us, me, "o-", color=BLUE); ax[0].axvline(fit["u"], color=RED, ls="--", label="chosen threshold u=%.2f"%fit["u"])
ax[0].set_xlabel("threshold u (loss %)"); ax[0].set_ylabel("mean excess E[L-u | L>u]")
ax[0].set_title("Mean-excess plot (roughly linear -> GPD tail)"); ax[0].legend()
# tail fit: empirical vs GPD survival above u
from scipy.stats import genpareto
exc = np.sort(L[L > fit["u"]] - fit["u"])
emp_surv = 1 - np.arange(1, len(exc)+1)/(len(exc)+1)
ax[1].plot(exc + fit["u"], emp_surv, "o", color=GREY, ms=3, label="empirical exceedances")
xx = np.linspace(0, exc.max(), 200)
ax[1].plot(xx + fit["u"], genpareto.sf(xx, fit["xi"], 0, fit["beta"]), color=RED, lw=2, label="fitted GPD")
ax[1].set_yscale("log"); ax[1].set_xlabel("loss (%)"); ax[1].set_ylabel("P(loss > x | exceedance)")
ax[1].set_title("GPD fits the tail of exceedances"); ax[1].legend()
plt.tight_layout(); plt.show()
print("Fitted GPD shape xi = %.3f (> 0 -> genuinely heavy, power-law tail); scale beta = %.2f; %d exceedances."
      % (fit["xi"], fit["beta"], fit["Nu"]))

Fitted GPD shape xi = 0.155 (> 0 -> genuinely heavy, power-law tail); scale beta = 0.76; 378 exceedances.

# VaR across the whole tail: historical runs out of data, Gaussian understates, EVT extrapolates
alphas = np.logspace(-3.3, -1.3, 25)
v_hist = [rm.historical_var(spy,a) for a in alphas]
v_norm = [rm.normal_var(spy,a) for a in alphas]
v_evt  = [rm.evt_var(fit,a) for a in alphas]
e_evt  = [rm.evt_es(fit,a) for a in alphas]
plt.figure(figsize=(9.5,5))
plt.plot(alphas, v_hist, "o-", color=BLUE, label="historical VaR")
plt.plot(alphas, v_norm, "s-", color=RED,  label="Gaussian VaR")
plt.plot(alphas, v_evt,  "-",  color=GREEN, lw=2.5, label="EVT VaR (GPD)")
plt.plot(alphas, e_evt,  "--", color=PURP, lw=2, label="EVT Expected Shortfall")
plt.xscale("log"); plt.gca().invert_xaxis(); plt.xlabel("tail probability alpha (rarer ->)"); plt.ylabel("loss (%)")
plt.title("Into the far tail: EVT extrapolates where history is empty and the Gaussian gives up")
plt.legend(); plt.show()
for a in (0.01, 0.001):
    print("  alpha=%.3f:  historical VaR %.2f%%   Gaussian %.2f%%   EVT %.2f%%   EVT-ES %.2f%%"
          % (a, rm.historical_var(spy,a), rm.normal_var(spy,a), rm.evt_var(fit,a), rm.evt_es(fit,a)))

  alpha=0.010:  historical VaR 3.10%   Gaussian 2.46%   EVT 3.21%   EVT-ES 4.50%
  alpha=0.001:  historical VaR 6.12%   Gaussian 3.28%   EVT 6.23%   EVT-ES 8.07%

3. Coherence: why VaR can punish diversification

A risk measure $\rho$ is coherent if, among other properties, it is sub-additive: $\rho(A+B)\le\rho(A)+\rho(B)$ — merging two portfolios never increases risk, because diversification helps. VaR violates this. The cleanest example is two independent defaultable bonds: each pays a small coupon but, with small probability, defaults and loses everything. Individually, the default is rarer than the VaR level, so each bond's VaR looks safe; combine them and the chance that at least one defaults pushes a loss into the VaR window — so the diversified book reports a larger VaR than the sum of the parts. ES never does this.

rng = np.random.default_rng(0); n = 200000
def bond(p=0.04): return np.where(rng.random(n) < p, -100.0, 2.0)   # 4% default (-100), else +2 coupon
A, B = bond(), bond(); port = 0.5*A + 0.5*B
al = 0.05
vA, vB, vP = rm.historical_var(A,al), rm.historical_var(B,al), rm.historical_var(port,al)
eA, eB, eP = rm.historical_es(A,al),  rm.historical_es(B,al),  rm.historical_es(port,al)
fig, ax = plt.subplots(1, 2, figsize=(12, 4.2))
ax[0].bar(["VaR(A)+VaR(B)","VaR(A+B)"], [vA+vB, vP], color=[GREY, RED])
ax[0].set_title("VaR FAILS sub-additivity (5% level)"); ax[0].set_ylabel("VaR (loss)")
ax[0].text(1, vP, "diversified\nis WORSE!", ha="center", va="bottom", color=RED, fontweight="bold")
ax[1].bar(["ES(A)+ES(B)","ES(A+B)"], [eA+eB, eP], color=[GREY, GREEN])
ax[1].set_title("ES is coherent (diversified <= sum)"); ax[1].set_ylabel("ES (loss)")
plt.tight_layout(); plt.show()
print("VaR:  VaR(A)+VaR(B) = %.1f   but   VaR(A+B) = %.1f   -> sub-additivity VIOLATED" % (vA+vB, vP))
print("ES :  ES(A)+ES(B)   = %.1f   and   ES(A+B)  = %.1f   -> coherent (diversification rewarded)" % (eA+eB, eP))

VaR:  VaR(A)+VaR(B) = -4.0   but   VaR(A+B) = 49.0   -> sub-additivity VIOLATED
ES :  ES(A)+ES(B)   = 157.0   and   ES(A+B)  = 50.7   -> coherent (diversification rewarded)

4. Risk contributions: who drives the risk?

For a portfolio, the total risk can be split into per-asset contributions that sum exactly to it (the Euler decomposition, valid because VaR and ES are homogeneous of degree one). The component for asset $i$ answers "how much of the portfolio's tail risk is because of holding $i$?" — the number that actually guides hedging and limits, far more useful than a standalone per-asset VaR.

w = np.ones(len(ASSETS))/len(ASSETS)
mu, Sig = X.mean(0), np.cov(X.T)
cvar = rm.component_var_normal(w, mu, Sig, 0.01)
ces  = rm.component_es_historical(X, w, 0.01)
xx = np.arange(len(ASSETS))
fig, ax = plt.subplots(1, 2, figsize=(13, 4.4))
ax[0].bar(xx, 100*cvar/cvar.sum(), color=BLUE); ax[0].set_xticks(xx); ax[0].set_xticklabels(ASSETS)
ax[0].set_ylabel("% of portfolio VaR"); ax[0].set_title("Component VaR (Gaussian, Euler)")
ax[1].bar(xx, 100*ces/ces.sum(), color=PURP); ax[1].set_xticks(xx); ax[1].set_xticklabels(ASSETS)
ax[1].set_ylabel("% of portfolio ES"); ax[1].set_title("Component ES (historical, Euler)")
plt.tight_layout(); plt.show()
print("Equal-weighted, but NOT equal-risk: the equity/EM/credit names (SPY, EEM, HYG) dominate the tail,")
print("while TLT (Treasuries) contributes little -- often negative in a crisis, the genuine diversifier.")
print("Component VaR sums to %.3f = total VaR %.3f (Euler)." % (cvar.sum(), rm.normal_var(X@w,0.01)))

Equal-weighted, but NOT equal-risk: the equity/EM/credit names (SPY, EEM, HYG) dominate the tail,
while TLT (Treasuries) contributes little -- often negative in a crisis, the genuine diversifier.
Component VaR sums to 1.375 = total VaR 1.375 (Euler).

5. Backtesting: does the VaR actually hold?

A VaR model is only as good as its out-of-sample track record. We roll a 1% one-day VaR through the sample (re-estimating from a trailing window each day) and count violations — days the loss exceeded the VaR. Two formal tests: Kupiec (is the violation rate right — should be ~1%?) and Christoffersen (are violations independent, not clustered in crises?). We will see the Gaussian fail the rate test outright, the historical get the rate about right — and both fail the independence test, because losses cluster in crises. That clustering is the fingerprint of volatility clustering, which a static VaR cannot capture and a dynamic (GARCH) VaR can — a direct pointer to the volatility-modelling arc elsewhere in the collection.

def backtest_var(r, window, method, alpha=0.01):
    viol = []
    for t in range(window, len(r)):
        past = r[t-window:t]
        v = rm.normal_var(past, alpha) if method == "normal" else rm.historical_var(past, alpha)
        viol.append(1 if r[t] < -v else 0)
    return np.array(viol)

W = 500
print("Rolling 1%% VaR backtest on SPY (window=%d days), target violation rate 1.0%%:\n" % W)
print("%-12s %10s %14s %18s" % ("method", "violations", "Kupiec p", "Christoffersen p"))
for method in ("normal", "historical"):
    v = backtest_var(spy, W, method)
    rate = 100*v.mean(); lr_k, p_k = rm.kupiec_pof(v.sum(), len(v), 0.01); lr_c, p_c = rm.christoffersen_independence(v)
    verdict = "PASS" if (p_k > 0.05 and p_c > 0.05) else "FAIL"
    print("%-12s %5d (%.2f%%) %14.3f %14.3f    %s" % (method, v.sum(), rate, p_k, p_c, verdict))
print("\nThe Gaussian VaR is breached far more than 1%% of the time (it understates the tail) -> fails")
print("Kupiec. The historical VaR gets the RATE about right (Kupiec ok) but BOTH cluster their")
print("violations in crises -> both fail Christoffersen independence. That clustering is volatility")
print("clustering: a static VaR cannot remove it -- a dynamic GARCH-VaR (see the volatility arc) can.")

Rolling 1% VaR backtest on SPY (window=500 days), target violation rate 1.0%:

method       violations       Kupiec p   Christoffersen p

normal          82 (2.51%)          0.000          0.000    FAIL

historical      43 (1.31%)          0.085          0.000    FAIL

The Gaussian VaR is breached far more than 1%% of the time (it understates the tail) -> fails
Kupiec. The historical VaR gets the RATE about right (Kupiec ok) but BOTH cluster their
violations in crises -> both fail Christoffersen independence. That clustering is volatility
clustering: a static VaR cannot remove it -- a dynamic GARCH-VaR (see the volatility arc) can.

6. Summary

• VaR is the tail quantile; ES is the average loss beyond it. ES sees how bad the tail is; VaR does not.
• Four estimators, and where they fail: the Gaussian understates every tail; Cornish–Fisher helps for mild non-normality but explodes at equity-like kurtosis; historical and Student-$t$ are trustworthy where data is deep; EVT is the only principled tool for the far tail, extrapolating via a Generalized-Pareto fit with heavy shape $\xi>0$.
• Coherence: VaR fails sub-additivity — it can report a diversified book as riskier than its parts — while ES is coherent. This is why regulation moved to ES.
• Risk contributions (Euler) split portfolio risk into per-asset pieces that sum to the whole, revealing the true risk drivers (equities/credit) versus the diversifiers (Treasuries).
• Backtesting (Kupiec, Christoffersen) is non-negotiable: the Gaussian VaR fails the rate test; the historical gets the rate right; and both fail independence because violations cluster — the signature of volatility clustering that only a dynamic (GARCH) VaR removes.

This closes the measurement side of the arc: having estimated distributions (Projects 1–2), built dependence (copulas) and projected to the horizon, we can now put an honest, coherent number on the downside — and check that it holds.