""" riskmeasures.py -- Value-at-Risk, Expected Shortfall, Extreme-Value Theory, Cornish-Fisher, coherence, risk contributions and VaR backtesting. Backs the notebooks in "Coherent Risk: VaR, Expected Shortfall & Extreme-Value Theory". Follows Meucci, "Risk and Asset Allocation" (Springer 2005), chapter 5-C/D. CONVENTIONS ----------- r : returns (a bad outcome is a large negative return). A risk measure is reported as a POSITIVE loss. `alpha` is the TAIL PROBABILITY (e.g. 0.01 = the 99% level). Loss L = -r. The measures ------------ VALUE-AT-RISK (VaR_alpha) = the loss the portfolio exceeds with probability alpha; i.e. minus the alpha-quantile of returns. Intuitive, ubiquitous -- and NOT coherent: it can PENALISE diversification (fails sub-additivity). EXPECTED SHORTFALL (ES_alpha, a.k.a. CVaR) = the average loss GIVEN that VaR is breached. It looks past the quantile into the tail, and it IS coherent (sub-additive), so it never punishes diversification. We estimate both four ways -- historical, Gaussian, Cornish-Fisher (moment based), and Extreme-Value-Theory (peaks-over-threshold) -- because the far tail, where the scary numbers live, is exactly where the naive methods fail. """ import numpy as np from scipy.stats import norm, t as student_t, genpareto, chi2, skew, kurtosis # --------------------------------------------------------------------------- # # VaR and ES: the four estimators # # --------------------------------------------------------------------------- # def historical_var(r, alpha=0.01): return -np.quantile(np.asarray(r, float), alpha) def historical_es(r, alpha=0.01): """Average of the worst ceil(alpha*n) returns -- correct even for discrete / point-mass distributions (unlike an `r <= quantile` mask).""" r = np.sort(np.asarray(r, float)) k = max(1, int(np.ceil(alpha * len(r)))) return -r[:k].mean() def normal_var(r, alpha=0.01): mu, sd = np.mean(r), np.std(r, ddof=1) return -(mu + sd * norm.ppf(alpha)) def normal_es(r, alpha=0.01): mu, sd = np.mean(r), np.std(r, ddof=1) return -(mu - sd * norm.pdf(norm.ppf(alpha)) / alpha) def cornish_fisher_z(alpha, s, k): """CF-adjusted normal quantile for skewness s and EXCESS kurtosis k.""" z = norm.ppf(alpha) return (z + (z**2 - 1) * s / 6 + (z**3 - 3*z) * k / 24 - (2*z**3 - 5*z) * s**2 / 36) def cornish_fisher_var(r, alpha=0.01): mu, sd = np.mean(r), np.std(r, ddof=1) z = cornish_fisher_z(alpha, skew(r), kurtosis(r)) return -(mu + sd * z) def t_fit_var_es(r, alpha=0.01): """Parametric Student-t VaR & ES (fit nu, mu, scale by MLE).""" nu, mu, sc = student_t.fit(r) q = student_t.ppf(alpha, nu, mu, sc) var = -q # ES of a location-scale t: xa = student_t.ppf(alpha, nu) es_std = -(nu + xa**2) / (nu - 1) * student_t.pdf(xa, nu) / alpha es = -(mu + sc * es_std) return var, es, dict(nu=nu, mu=mu, scale=sc) # --------------------------------------------------------------------------- # # Extreme-Value Theory: peaks over threshold (GPD) # # --------------------------------------------------------------------------- # def mean_excess(r, thresholds): """Mean-excess function e(u)=E[L-u | L>u] of the losses. Linear-in-u above some point signals a Generalized-Pareto tail (threshold-selection diagnostic).""" L = -np.asarray(r, float) return np.array([L[L > u].mean() - u if (L > u).any() else np.nan for u in thresholds]) def fit_gpd_pot(r, thr_q=0.90): """Fit a Generalized Pareto Distribution to the losses that exceed a high threshold u (the `thr_q` quantile of losses). Returns the tail parameters.""" L = -np.asarray(r, float) u = np.quantile(L, thr_q) exc = L[L > u] - u xi, _, beta = genpareto.fit(exc, floc=0) return dict(u=u, xi=float(xi), beta=float(beta), Nu=len(exc), n=len(L)) def evt_var(fit, alpha=0.001): """POT tail estimator of VaR at (small) tail probability alpha.""" u, xi, beta, Nu, n = fit["u"], fit["xi"], fit["beta"], fit["Nu"], fit["n"] return u + (beta / xi) * ((n * alpha / Nu) ** (-xi) - 1) def evt_es(fit, alpha=0.001): var = evt_var(fit, alpha) xi, beta, u = fit["xi"], fit["beta"], fit["u"] return (var + beta - xi * u) / (1 - xi) # --------------------------------------------------------------------------- # # Risk contributions (Euler allocation) # # --------------------------------------------------------------------------- # def component_var_normal(w, mu, Sigma, alpha=0.01): """Euler decomposition of Gaussian portfolio VaR into per-asset components (they sum to the total VaR).""" w = np.asarray(w, float); z = norm.ppf(alpha) sd = np.sqrt(w @ Sigma @ w) d = -mu - z * (Sigma @ w) / sd # marginal VaR return w * d # component VaR def component_es_historical(R, w, alpha=0.01): """Historical Euler ES contributions: average each asset's loss over the scenarios where the PORTFOLIO breaches its VaR.""" R = np.asarray(R, float); w = np.asarray(w, float) p = R @ w; q = np.quantile(p, alpha); tail = p <= q return -w * R[tail].mean(axis=0) # component ES (sum = total ES) # --------------------------------------------------------------------------- # # Backtesting VaR # # --------------------------------------------------------------------------- # def kupiec_pof(n_viol, n, alpha): """Kupiec proportion-of-failures LR test: is the violation rate = alpha? Returns (LR statistic, p-value) under chi-square(1).""" pi = n_viol / n ll0 = (n - n_viol) * np.log(1 - alpha) + n_viol * np.log(alpha) ll1 = (n - n_viol) * np.log(1 - pi) + n_viol * np.log(pi) if 0 < pi < 1 else 0.0 LR = -2 * (ll0 - ll1) return LR, 1 - chi2.cdf(LR, 1) def christoffersen_independence(viol): """Christoffersen test that VaR violations are independent (not clustered).""" v = np.asarray(viol, int) n00 = np.sum((v[:-1] == 0) & (v[1:] == 0)); n01 = np.sum((v[:-1] == 0) & (v[1:] == 1)) n10 = np.sum((v[:-1] == 1) & (v[1:] == 0)); n11 = np.sum((v[:-1] == 1) & (v[1:] == 1)) p01 = n01 / max(n00 + n01, 1); p11 = n11 / max(n10 + n11, 1); p = (n01 + n11) / max(len(v) - 1, 1) def s(a, b, pp): return a * np.log(1 - pp) + b * np.log(pp) if 0 < pp < 1 else 0.0 LR = -2 * ((s(n00, n01, p) + s(n10, n11, p)) - (s(n00, n01, p01) + s(n10, n11, p11))) return LR, 1 - chi2.cdf(LR, 1)