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

Python · PyMC · R · Download risk-measures module

Two Measures

The arc so far has been about building portfolios. This project is about measuring what they can lose — and about the ways the standard measure misleads. Value-at-Risk is a quantile of the loss distribution: the threshold breached with probability $\alpha$ . Expected Shortfall is the average loss given that the threshold is breached. The difference sounds technical and turns out to be decisive, both statistically and mathematically.

\text{VaR}_\alpha=\inf\{\ell:\Pr(L>\ell)\le\alpha\},\qquad \text{ES}_\alpha=\mathbb{E}\big[L\mid L>\text{VaR}_\alpha\big]

Four estimators, and where they break

The first section estimates both four ways on the same daily S&P returns, and the disagreement is the point. The Gaussian understates every tail — 2.46% against a historical 3.10% at the 1% level, and 3.28% against 6.12% at 0.1%. Cornish–Fisher, the standard moment-based correction, is reasonable at 5% and then explodes to 14.18% at 0.1%: the expansion is only valid for mild non-normality, and this series has excess kurtosis around 12. The Student- $t$ and historical estimates are the trustworthy pair. Python and R reproduce every one of these figures identically.

S&P daily loss, VaR (%)	historical	Gaussian	Cornish–Fisher	Student- $t$
$\alpha=5\%$	1.67	1.72	1.68	1.46
$\alpha=1\%$	3.10	2.46	5.72	3.03
$\alpha=0.1\%$	6.12	3.28	14.18	7.40

Extreme-value theory, past the edge of the data

Where the data runs out, extreme-value theory takes over. Fitting a Generalized Pareto distribution to the 378 exceedances over a 90th-percentile threshold gives a shape parameter $\xi=0.155$ — decisively positive, meaning a genuine power-law tail rather than a thin one. The value of this is extrapolation: at the 0.1% level EVT gives 6.23% where the Gaussian says 3.28%, and unlike the historical estimate it keeps working past the point where the sample simply contains no further observations. R's evir package reproduces the fit to three decimals.

Why VaR Is Not Coherent

Then the mathematical objection, demonstrated rather than described. A risk measure ought to be sub-additive — combining two positions should never increase measured risk, because diversification is supposed to help. VaR fails this. The notebooks construct a case where $\text{VaR}(A)+\text{VaR}(B)=-4.0$ while the diversified $\text{VaR}(A+B)=49.0$ : on VaR's own arithmetic, diversifying made things worse. Expected Shortfall on the same portfolios gives 157.0 against 50.7 — coherent, and rewarding diversification as it should. This is the formal reason regulators moved from VaR to ES.

\rho(A+B)\le\rho(A)+\rho(B)\qquad\text{(sub-additivity — ES satisfies it, VaR does not)}

Does the number hold up?

The closing section asks the question a risk manager actually cares about: does the number hold up? A rolling 1% VaR backtest on 500-day windows is scored with the Kupiec test, which asks whether the number of breaches matches the 1% promised, and the Christoffersen test, which asks whether they arrived independently or in clusters — a model can take exactly the right number of hits and still be useless if they all land in one bad fortnight. Gaussian VaR is breached 2.51% of the time against a 1% target and fails both. Historical VaR gets the rate about right at 1.31% and passes Kupiec — but both fail Christoffersen, because their breaches arrive in clusters. That clustering is volatility clustering, and no static VaR can remove it; it needs a dynamic model, which is exactly what the GARCH-VaR work elsewhere in the collection provides. Ending on a diagnosis that points outside itself is the right way to close.

Honest Error Bars on the Far Tail

A PyMC companion fits the same GPD tail Bayesianly, and delivers the most quietly devastating number in the project. The posterior for the tail shape is $\xi=0.164$ with a 94% interval of $[0.069,\,0.275]$ — heavy for certain, but how heavy is genuinely uncertain. Propagate that into the far tail and the 1-in-10,000-day VaR has a posterior mean of 11.1% with an interval of [8.3%, 15.6%] — a factor of 1.9 between the plausible extremes. A single plug-in figure of 10.5% conceals all of it. The tail is precisely where point estimates pretend to precision they do not have.

Notebooks

Python VaR and ES four ways, GPD extreme-value fitting, the sub-additivity violation, Euler risk contributions, and Kupiec / Christoffersen backtesting (Python) Open → Python Bayesian extreme-value risk — the GPD tail with a posterior, and far-tail VaR and ES reported with the credible intervals a point estimate hides (PyMC) Open → R Independent R implementation — the same four estimators, EVT via the evir package, the coherence demonstration and the same backtests (R) Open →

Downloads

riskmeasures.py Historical, Gaussian, Cornish–Fisher and Student-t VaR/ES; mean-excess and peaks-over-threshold GPD fitting; Euler risk contributions; Kupiec and Christoffersen backtests (NumPy / SciPy) crossasset_daily.csv Daily returns for SPY, TLT, GLD, HYG and EEM — 3,772 trading days

Risk-Measures Module — Source Code

"""
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)

References

Artzner, P., Delbaen, F., Eber, JM. & Heath, D. (1999). Coherent measures of risk. Mathematical Finance 9(3), 203–228. — the coherence axioms, and the sub-additivity failure demonstrated here
McNeil, A. J., Frey, R. & Embrechts, P. (2015). Quantitative Risk Management (2nd ed.). Princeton University Press. — the standard treatment of VaR, ES and peaks-over-threshold EVT
Pickands, J. (1975). Statistical inference using extreme order statistics. Annals of Statistics 3(1), 119–131. — the Generalized Pareto limit for threshold exceedances
Kupiec, P. H. (1995). Techniques for verifying the accuracy of risk measurement models. Journal of Derivatives 3(2), 73–84. — the proportion-of-failures test used in the backtest
Christoffersen, P. F. (1998). Evaluating interval forecasts. International Economic Review 39(4), 841–862. — the independence test that both static VaR methods fail here
Basel Committee on Banking Supervision (2016). Minimum Capital Requirements for Market Risk. — the regulatory move from VaR to Expected Shortfall
Meucci, A. (2005). Risk and Asset Allocation. Springer. — risk measures and their role in the allocation framework of this arc