Value at Risk and Expected Shortfall

Deep dive

How it works

Parts:

• Scenario source — Σ (parametric), a window of past return vectors (historical), or draws L·z (Monte Carlo).
• Revaluation — scenario to loss: a dot product, or a pricer (Black-Scholes) for options.
• Tail estimator — VaR is the ⌈αn⌉-th smallest loss; ES the mean from there up.
• Backtest — counts days the loss exceeded VaR.
• Allocator — Euler components wᵢ·∂VaR/∂wᵢ.

1. Estimate Σ: sample window, or EWMA Σ ← λΣ + (1−λ)rrᵀ (RiskMetrics: daily λ = 0.94; exponential smoothing); shrink it when d approaches the window length.
2. Parametric: VaR = z_α·σ_p·√h, σ_p = √(wᵀΣw); √h assumes independent, identically distributed returns.
3. Historical: apply each past return vector to today's weights. Monte Carlo: factor Σ = LLᵀ (Cholesky), draw r = L·z, revalue.
4. Select the k-th smallest of n losses, k = ⌈αn⌉; ES averages elements k and up.
5. Backtest: x exceptions in T days, p = 1−α. Kupiec's LR = −2 ln[(1−p)^{T−x}p^x] + 2 ln[(1−x/T)^{T−x}(x/T)^x] is χ²(1) if coverage is right (reject above 3.84 at 5 %).
6. Allocate: component_i = z√h·wᵢ(Σw)ᵢ/σ_p; they sum to VaR.

 Σ / history / L·z ──► scenarios (n×d) ──► revalue(w) ──► losses (n) ──► select k ──► VaR, ES
                              ▲  reused for every portfolio                              │
                                                              exceptions per day ◄───────┘ ──► Kupiec

Interpolating between order statistics k and k+1 moves a 250-day, 99 % VaR; state the rule.

Code

pub struct Rng(u64); // xorshift64*
impl Rng {
    fn unif(&mut self) -> f64 { // in (0,1)
        self.0 ^= self.0 >> 12; self.0 ^= self.0 << 25; self.0 ^= self.0 >> 27;
        ((self.0.wrapping_mul(0x2545F4914F6CDD1D) >> 11) as f64 + 0.5) / (1u64 << 53) as f64
    }
    fn normals(&mut self) -> (f64, f64) { // Box-Muller
        let (r, t) = ((-2.0 * self.unif().ln()).sqrt(), 6.283185307179586 * self.unif());
        (r * t.cos(), r * t.sin())
    }
}

/// Lower Cholesky factor of an n*n row-major covariance matrix.
pub fn cholesky(a: &[f64], n: usize) -> Vec<f64> {
    let mut l = vec![0.0; n * n];
    for i in 0..n { for j in 0..=i {
        let s: f64 = (0..j).map(|k| l[i * n + k] * l[j * n + k]).sum();
        l[i * n + j] = if i == j { (a[i * n + i] - s).sqrt() } else { (a[i * n + j] - s) / l[j * n + j] };
    } }
    l
}

/// Delta-normal VaR (z = z_α: 2.3263 at 99 %, 1.9600 at 97.5 %; cov = daily returns; zero mean) and Euler components.
pub fn parametric(w: &[f64], cov: &[f64], z: f64, h: f64) -> (f64, Vec<f64>) {
    let n = w.len();
    let sw: Vec<f64> = (0..n).map(|i| (0..n).map(|j| cov[i * n + j] * w[j]).sum()).collect(); // Σw
    let sigma = (0..n).map(|i| w[i] * sw[i]).sum::<f64>().sqrt(); // sqrt(wᵀΣw)
    let k = z * h.sqrt();
    (k * sigma, (0..n).map(|i| k * w[i] * sw[i] / sigma).collect())
}

/// n scenarios of d correlated daily returns, row-major: r = L z.
pub fn scenarios_mc(chol: &[f64], d: usize, n: usize, rng: &mut Rng) -> Vec<f64> {
    let mut out = Vec::with_capacity(n * d);
    let mut z = vec![0.0; d + 1];
    for _ in 0..n {
        for p in 0..(d + 1) / 2 { (z[2 * p], z[2 * p + 1]) = rng.normals(); }
        for i in 0..d { out.push((0..=i).map(|k| chol[i * d + k] * z[k]).sum()); }
    }
    out
}
/// Loss per scenario row (linear book; a pricer goes here otherwise).
pub fn losses(w: &[f64], scen: &[f64]) -> Vec<f64> {
    scen.chunks(w.len()).map(|r| -r.iter().zip(w).map(|(a, b)| a * b).sum::<f64>()).collect()
}
/// (VaR, ES): the ceil(αn)-th smallest loss, and the mean from there up. O(n) selection; reorders `l`.
pub fn tail(l: &mut [f64], alpha: f64) -> (f64, f64) {
    let k = (alpha * l.len() as f64).ceil() as usize - 1;
    let (_, v, hi) = l.select_nth_unstable_by(k, |a, b| a.total_cmp(b));
    let v = *v;
    (v, (v + hi.iter().sum::<f64>()) / (1 + hi.len()) as f64)
}
/// Kupiec POF likelihood ratio: x exceptions in t days at level p = 1-α; > 3.84 rejects at 5 %.
pub fn kupiec(x: f64, t: f64, p: f64) -> f64 {
    let ll = |q: f64| (t - x) * (1.0 - q).ln() + if x > 0.0 { x * q.ln() } else { 0.0 };
    2.0 * (ll(x / t) - ll(p))
}
fn main() {
    let vol = [0.01, 0.02, 0.015];
    let rho = [[1.0, 0.3, 0.2], [0.3, 1.0, 0.5], [0.2, 0.5, 1.0]];
    let cov: Vec<f64> = (0..9).map(|k| vol[k / 3] * vol[k % 3] * rho[k / 3][k % 3]).collect();
    let (w, a, z) = ([1e6, 2e6, 5e5], 0.99, 2.3263478740408408);
    let (pv, comp) = parametric(&w, &cov, z, 1.0);
    let mut rng = Rng(0x9E3779B97F4A7C15);
    let sc = scenarios_mc(&cholesky(&cov, 3), 3, 400_000, &mut rng);
    let (mv, es) = tail(&mut losses(&w, &sc), a);
    assert!((mv / pv - 1.0).abs() < 0.01 && es >= mv);   // MC ~ parametric
    let mut s: Vec<f64> = (1..=1000).rev().map(f64::from).collect();
    assert_eq!(tail(&mut s, 0.99), (990.0, 995.0));              // order statistic; mean of 990..=1000
    let w2: Vec<f64> = w.iter().map(|x| 2.0 * x).collect();
    assert!((tail(&mut losses(&w2, &sc), a).0 / mv - 2.0).abs() < 1e-12);  // linear
    let alone: f64 = (0..3).map(|i| { let mut x = [0.0; 3]; x[i] = w[i]; parametric(&x, &cov, z, 1.0).0 }).sum();
    assert!(pv < alone);                                // diversified
    assert!((comp.iter().sum::<f64>() / pv - 1.0).abs() < 1e-12);           // Euler
    assert!((parametric(&w, &cov, z, 10.0).0 / pv - 10f64.sqrt()).abs() < 1e-12); // √h
    assert!(kupiec(2.5, 250.0, 0.01).abs() < 1e-9 && kupiec(8.0, 250.0, 0.01) > 3.84);
}

Compiled and run (rustc 1.99.1), asserts pass: Monte Carlo within 1 % of parametric; ES ≥ VaR; 1..1000 at 99 % gives 990 and 995; VaR linear in size; diversified < standalone sum; Euler parts sum to VaR; √h; Kupiec.

Complexity and constants

Parametric is O(d²). Scenario methods cost O(n·d) per portfolio plus O(n) selection; Monte Carlo generation, O(n·d²), dominates. The scenario matrix is n·d·8 bytes (1M × 50 = 400 MB), so revaluation streams memory. A Monte Carlo quantile has standard error √(α(1−α)/n)/f(q): order of 0.25 % of VaR at n = 400k, α = 0.99.

Efficiency

• Cost: Apple M4 Pro, one core, measured, 1M scenarios: generation 39 ns per scenario at d = 3, 694 ns at d = 50; revaluing a portfolio 3.5 ns (d = 3), 43 ns (d = 50, includes allocating the loss vector); selection plus ES 0.9-1.2 ns. Generation is compute-bound.
• Levers:
  • Quickselect instead of sort — 12.5 → 0.9 ns/scenario at 1M, 7.1 → 0.8 at 10k (measured) — any tail quantile.
  • One scenario matrix for all portfolios — d = 50: 43 → 9.6 ns per scenario-portfolio over 100 portfolios (measured) — many books on shared factors; a matrix product, so BLAS.
  • Delta or delta-gamma instead of full revaluation — a pricer call becomes a few flops (not measured) — smooth payoffs.
  • Low-discrepancy sequences or importance sampling — fewer scenarios for the same tail error (not measured) — Monte Carlo of rare tails.
  • Incremental window or EWMA Σ — O(d²) per day, not O(window·d²) (not measured) — daily runs.
• Naive → standard → state of the art:
  • Full sort, full revaluation · 12.5 ns/scenario for the sort alone at 1M (measured)
  • Selection, shared scenario matrix · ~1 ns per scenario for the tail, ~10 ns per scenario-portfolio at d = 50 (measured)
  • Delta-gamma proxies, GEMM, quasi-random or importance-sampled tails · not measured

Implementation notes

• Report VaR as a positive loss; sign slips are the commonest bug.
• Selection reorders the buffer; copy first if needed again.
• Fix the seed and report n: sampling error is part of the number. Use a real generator.
• Store scenarios as n contiguous rows of d returns.

Pitfalls

• √h scaling fails with autocorrelation, volatility clustering or intra-horizon trading; Basel uses per-risk-factor liquidity horizons instead.
• Normal tails understate fat-tailed returns.
• VaR ignores losses past the threshold and can rise when books merge; ES does not.
• Exceptions cluster; Kupiec does not test independence.

Variants and alternatives

In production

• Basel market-risk standard: ES at the 97.5th percentile, one-tailed; backtests at 97.5 % and 99 % over 12 months (MAR33).
• Solvency II: solvency capital targets 99.5 % over one year.
• Cloud-cost, energy-procurement and project-contingency tools report P80/P95 from Monte Carlo.

Implementations

Sources

• Artzner, Delbaen, Eber, Heath, Coherent Measures of Risk, Mathematical Finance 9(3), 1999, 203-228. paper
• Kupiec, Techniques for Verifying the Accuracy of Risk Measurement Models, Journal of Derivatives 3(2), 1995, 73-84. paper
• Basel Committee, Minimum capital requirements for market risk, January 2019. BIS
• J.P. Morgan/Reuters, RiskMetrics Technical Document, 4th ed., 1996. MSCI
• Directive 2009/138/EC (Solvency II), Art. 101.