Bayesian Extreme-Value Risk — the PyMC engine¶
The from-scratch notebook fit the tail with a single Generalized-Pareto estimate and reported one number for the 1-in-1000-day loss. But that is the most uncertain number in all of risk management: it is extrapolated from a few dozen extreme observations. A point estimate here is almost dishonest. This notebook fits the tail in PyMC, delivering the posterior distribution of the far-tail VaR and Expected Shortfall — error bars that widen, exactly as they should, the deeper into the tail we look.
Self-contained; no prior reading required.
The model: Bayesian peaks-over-threshold¶
Extreme-Value Theory says the exceedances of the losses over a high threshold $u$ follow a Generalized Pareto Distribution with shape $\xi$ (tail heaviness) and scale $\beta$: $$f(y)=\tfrac1\beta\Big(1+\xi\tfrac{y}{\beta}\Big)^{-1/\xi-1},\qquad y=\text{loss}-u>0.$$ We place priors on $(\xi,\beta)$, fit by MCMC, and propagate the posterior into every tail risk number. The shape $\xi$ is the crux — it governs how fast the tail decays — and, being estimated from the scarcest data, it is genuinely uncertain.
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import pymc as pm, arviz as az, pytensor.tensor as pt
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"
RNG = 5
R = pd.read_csv("crossasset_daily.csv", index_col=0, parse_dates=True)
spy = R["SPY"].values
L = -spy # losses
u = np.quantile(L, 0.90) # threshold (90th percentile of losses)
y = L[L > u] - u # exceedances
n, Nu = len(L), len(y)
print("Losses: %d days, threshold u=%.2f%%, %d exceedances above it" % (n, u, Nu))
g++ not available, if using conda: `conda install gxx`
Losses: 3772 days, threshold u=1.10%, 378 exceedances above it
The dataset¶
Daily SPY losses (%), 2010–2024 (3,772 days). We keep only the exceedances over the 90th-percentile loss threshold ($u\approx$ 1.1%) — the $\sim$377 observations that carry the tail information. Everything downstream depends on how those few extremes are modelled.
1. Bayesian fit of the tail¶
Weakly-informative priors: $\xi\sim\mathcal N(0.1,0.3)$ (centred on a mildly heavy tail but open to either sign) and $\beta\sim\text{HalfNormal}$. The GPD log-likelihood enters as a pm.Potential, with the support constraint $1+\xi y/\beta>0$ enforced.
with pm.Model() as gpd:
xi = pm.Normal("xi", 0.1, 0.3)
beta = pm.HalfNormal("beta", 2.0)
z = 1.0 + xi * y / beta
logp = pt.switch(z > 1e-12, -pt.log(beta) - (1.0 + 1.0/xi) * pt.log(pt.maximum(z, 1e-12)), -1e12)
pm.Potential("gpd_ll", logp.sum())
idata = pm.sample(2000, tune=2000, chains=4, cores=1, target_accept=0.95,
progressbar=False, random_seed=RNG)
print("max R-hat:", float(az.summary(idata, var_names=["xi","beta"])["r_hat"].max()))
print(az.summary(idata, var_names=["xi","beta"]).iloc[:, :5].to_string())
xi_d = idata.posterior["xi"].values.flatten(); beta_d = idata.posterior["beta"].values.flatten()
xi_mle = rm.fit_gpd_pot(spy, 0.90)["xi"]
print("\nPosterior tail shape xi: mean %.3f, 94%% HDI [%.3f, %.3f] (MLE point estimate %.3f)"
% (xi_d.mean(), np.percentile(xi_d,3), np.percentile(xi_d,97), xi_mle))
print("xi > 0 across essentially the whole posterior -> the tail is genuinely heavy (power-law), not")
print("thin -- but HOW heavy is uncertain, and that uncertainty must flow into the risk numbers.")
Initializing NUTS using jitter+adapt_diag...
Sequential sampling (4 chains in 1 job)
NUTS: [xi, beta]
Sampling 4 chains for 2_000 tune and 2_000 draw iterations (8_000 + 8_000 draws total) took 3 seconds.
max R-hat: 1.0
mean sd eti89_lb eti89_ub ess_bulk
xi 0.164 0.055 0.081 0.25 2953
beta 0.765 0.057 0.68 0.86 3419
Posterior tail shape xi: mean 0.164, 94% HDI [0.069, 0.275] (MLE point estimate 0.155)
xi > 0 across essentially the whole posterior -> the tail is genuinely heavy (power-law), not
thin -- but HOW heavy is uncertain, and that uncertainty must flow into the risk numbers.
2. The posterior of far-tail VaR and Expected Shortfall¶
For each posterior draw $(\xi,\beta)$ the peaks-over-threshold formulas give a VaR and ES at any tail probability $\alpha$: $$\text{VaR}_\alpha=u+\tfrac{\beta}{\xi}\Big[\big(\tfrac{n\alpha}{N_u}\big)^{-\xi}-1\Big],\qquad \text{ES}_\alpha=\tfrac{\text{VaR}_\alpha+\beta-\xi u}{1-\xi}.$$ Evaluating them across the posterior turns each risk number into a distribution. We look at the 1-in-1000 and 1-in-10000-day losses — and watch the uncertainty grow as the event gets rarer.
def evt_var_vec(xi, beta, alpha): return u + (beta/xi)*((n*alpha/Nu)**(-xi) - 1)
def evt_es_vec(xi, beta, alpha):
v = evt_var_vec(xi, beta, alpha); return (v + beta - xi*u)/(1 - xi)
fig, ax = plt.subplots(1, 2, figsize=(13, 4.4))
for a_, c in [(1e-3, BLUE), (1e-4, RED)]:
vv = evt_var_vec(xi_d, beta_d, a_)
ax[0].hist(vv, bins=60, density=True, alpha=.6, color=c, label="VaR at alpha=%.0e (mean %.1f%%)" % (a_, vv.mean()))
ax[0].set_xlabel("VaR (loss %)"); ax[0].set_title("Posterior of far-tail VaR — rarer = far more uncertain"); ax[0].legend()
alphas = np.logspace(-4, -1.5, 25)
V = np.array([evt_var_vec(xi_d, beta_d, a_) for a_ in alphas])
lo, med, hi = np.percentile(V, [3, 50, 97], axis=1)
ax[1].fill_between(alphas, lo, hi, color="#bee3f8", alpha=.8, label="94% credible band")
ax[1].plot(alphas, med, color=BLUE, lw=2, label="posterior median")
ax[1].plot(alphas, [rm.evt_var(rm.fit_gpd_pot(spy,0.90), a_) for a_ in alphas], "--", color=RED, lw=1.5, label="MLE plug-in")
ax[1].set_xscale("log"); ax[1].invert_xaxis(); ax[1].set_xlabel("tail probability alpha (rarer ->)")
ax[1].set_ylabel("VaR (loss %)"); ax[1].set_title("VaR term structure with credible band"); ax[1].legend()
plt.tight_layout(); plt.show()
for a_ in (1e-3, 1e-4):
vv = evt_var_vec(xi_d, beta_d, a_); ee = evt_es_vec(xi_d, beta_d, a_)
print(" alpha=%.0e: VaR %.1f%% [%.1f, %.1f] ES %.1f%% [%.1f, %.1f]"
% (a_, vv.mean(), np.percentile(vv,3), np.percentile(vv,97), ee.mean(), np.percentile(ee,3), np.percentile(ee,97)))
alpha=1e-03: VaR 6.4% [5.5, 7.7] ES 8.4% [6.7, 11.1] alpha=1e-04: VaR 11.1% [8.3, 15.6] ES 14.2% [9.7, 22.0]
# The honest headline: how wide is the 1-in-10000 loss estimate?
v1e4 = evt_var_vec(xi_d, beta_d, 1e-4)
print("The 1-in-10,000-day VaR has a posterior mean of %.1f%% but a 94%% credible interval of"
% v1e4.mean())
print("[%.1f%%, %.1f%%] -- a factor-of-%.1f range. A single plug-in number (%.1f%%) hides that the data"
% (np.percentile(v1e4,3), np.percentile(v1e4,97), np.percentile(v1e4,97)/np.percentile(v1e4,3),
rm.evt_var(rm.fit_gpd_pot(spy,0.90), 1e-4)))
print("simply cannot pin down the most extreme risks. Bayesian EVT reports that honestly; a point")
print("estimate pretends to a precision the tail does not possess.")
The 1-in-10,000-day VaR has a posterior mean of 11.1% but a 94% credible interval of [8.3%, 15.6%] -- a factor-of-1.9 range. A single plug-in number (10.5%) hides that the data simply cannot pin down the most extreme risks. Bayesian EVT reports that honestly; a point estimate pretends to a precision the tail does not possess.
3. Summary¶
- Extreme risk is extreme uncertainty. The far-tail VaR/ES is extrapolated from a few dozen exceedances; a point estimate misrepresents how little we know.
- A Bayesian peaks-over-threshold model returns the posterior of the tail-shape $\xi$ (heavy, $>0$, but with real spread) and, through it, the posterior of VaR and ES at any tail probability.
- The credible band widens as the event gets rarer — the 1-in-10,000-day loss spans a large multiplicative range. That is not a defect of the method; it is the truth about the tail, made visible.
Bayesian EVT is the natural meeting point of this arc's two threads — tail risk (this project) and estimation uncertainty (Projects 2 and the horizon notebook): the scariest number in finance, reported with the humility it deserves.