"""value-at-risk-three-ways.py — one-day Value at Risk computed three ways on FRED's S&P 500 series, then backtested.

Input : datasets/value-at-risk-three-ways.csv - the derived daily simple return series computed from
        FRED's SP500 daily close (a price index without dividends, on FRED's ten-year rolling window) as
        it stood at the 2026-09-08 retrieval. The index levels themselves are not shipped. S&P Dow Jones
        Indices LLC states in the FRED series notes that "Reproduction of S&P 500 in any form is
        prohibited except with the prior written permission of S&P Dow Jones Indices LLC ("S&P")", so
        this file carries only the returns computed from the series, which are this site's own output.
        Source: S&P Dow Jones Indices LLC, https://fred.stlouisfed.org/series/SP500.
        Run with --download to pull SP500 from FRED and rebuild the return file end to end. FRED's
        ten-year window rolls forward, so a pull today covers a later sample than the article reports.
Method: simple daily returns r_t = P_t/P_{t-1} - 1, expressed as a loss L_t = -r_t.
        (1) Historical simulation: VaR_alpha = the empirical alpha-quantile of the loss sample
            (numpy's linear-interpolation quantile).
        (2) Parametric normal: VaR_alpha = mu_L + z_alpha * sigma_L with z from the inverse normal CDF
            (bisection on math.erfc, so no scipy) and sigma using ddof=1.
        (3) Cornish-Fisher: the same formula with z replaced by the Cornish-Fisher expansion
            z + (z^2-1)S/6 + (z^3-3z)K/24 - (2z^3-5z)S^2/36, S = sample skewness, K = sample excess
            kurtosis (both Fisher-Pearson, as pandas reports them).
        Expected shortfall is the mean loss beyond the historical VaR, and the closed-form normal ES
        phi(z)/(1-alpha) * sigma + mu for the parametric case.
        Backtest: a rolling estimation window of WINDOW trading days produces a one-day-ahead VaR for the
        next day; an exception is a realised loss above that forecast. The Kupiec proportion-of-failures
        likelihood-ratio statistic is compared with a chi-square, 1 degree of freedom, whose upper tail is
        erfc(sqrt(LR/2)) exactly.
        Nothing here is a forecast of any future loss and nothing here is investment advice; the numbers
        describe one sample of one price index.
Output: prints the tables used in the article and writes datasets/value-at-risk-three-ways-backtest.csv
        (date, return, loss, the three 99% forecasts, the three exception flags).
Run   : python code/value-at-risk-three-ways.py [--download]
Needs : Python 3.13, pandas 3.0.2, numpy 2.3. Standard library only for the statistics.
"""
import io
import math
import sys
import urllib.request

import numpy as np
import pandas as pd

CSV = "datasets/value-at-risk-three-ways.csv"
OUT = "datasets/value-at-risk-three-ways-backtest.csv"
WINDOW = 500
LEVELS = (0.95, 0.99)


def download():
    """Pull SP500 from FRED and write only the derived returns; the index levels are never stored."""
    req = urllib.request.Request("https://fred.stlouisfed.org/graph/fredgraph.csv?id=SP500",
                                 headers={"User-Agent": "prism-data-lab/1.0", "Accept": "*/*"})
    txt = urllib.request.urlopen(req, timeout=60).read().decode()
    s = pd.read_csv(io.StringIO(txt), na_values=".", index_col=0).iloc[:, 0]
    s.index.name = "date"
    r = s.dropna().pct_change().rename("return_simple")
    r.reindex(s.index).to_frame().to_csv(CSV)


def norm_cdf(x: float) -> float:
    return 0.5 * math.erfc(-x / math.sqrt(2.0))


def norm_ppf(p: float) -> float:
    """Inverse standard normal CDF by bisection on erfc. Accurate to ~1e-12 for p in (1e-9, 1-1e-9)."""
    lo, hi = -12.0, 12.0
    for _ in range(200):
        mid = 0.5 * (lo + hi)
        if norm_cdf(mid) < p:
            lo = mid
        else:
            hi = mid
    return 0.5 * (lo + hi)


def norm_pdf(x: float) -> float:
    return math.exp(-0.5 * x * x) / math.sqrt(2.0 * math.pi)


def cornish_fisher_z(z: float, skew: float, kurt: float) -> float:
    return (z
            + (z * z - 1.0) * skew / 6.0
            + (z ** 3 - 3.0 * z) * kurt / 24.0
            - (2.0 * z ** 3 - 5.0 * z) * skew * skew / 36.0)


def var_three_ways(loss: np.ndarray, alpha: float) -> dict:
    mu = float(np.mean(loss))
    sd = float(np.std(loss, ddof=1))
    s = float(pd.Series(loss).skew())
    k = float(pd.Series(loss).kurt())
    z = norm_ppf(alpha)
    zcf = cornish_fisher_z(z, s, k)
    hist = float(np.quantile(loss, alpha))
    return {
        "historical": hist,
        "normal": mu + z * sd,
        "cornish_fisher": mu + zcf * sd,
        "es_historical": float(np.mean(loss[loss >= hist])),
        "es_normal": mu + sd * norm_pdf(z) / (1.0 - alpha),
        "mu": mu, "sd": sd, "skew": s, "kurt": k, "z": z, "z_cf": zcf,
    }


def kupiec(x: int, T: int, p: float) -> tuple[float, float]:
    """Proportion-of-failures LR statistic and its chi-square(1) upper-tail probability."""
    if x == 0:
        lr = -2.0 * (T * math.log(1.0 - p))
    else:
        pi = x / T
        lr = -2.0 * ((T - x) * math.log(1.0 - p) + x * math.log(p)
                     - (T - x) * math.log(1.0 - pi) - x * math.log(pi))
    return lr, math.erfc(math.sqrt(max(lr, 0.0) / 2.0))


def main():
    if "--download" in sys.argv:
        download()
    r = pd.read_csv(CSV, parse_dates=["date"], index_col="date")["return_simple"].dropna()
    loss = -r
    print(f"sample {r.index[0]:%Y-%m-%d} to {r.index[-1]:%Y-%m-%d}: {len(r)} daily returns")
    print(f"loss mean {100 * loss.mean():.4f}%, sd {100 * loss.std(ddof=1):.4f}%, "
          f"skew {loss.skew():.2f}, excess kurtosis {loss.kurt():.1f}")

    print("\n| level | historical | normal | Cornish-Fisher | z | z (CF) |")
    print("|---|---|---|---|---|---|")
    full = {}
    for a in LEVELS:
        v = var_three_ways(loss.to_numpy(), a)
        full[a] = v
        print(f"| {100 * a:.0f}% | {100 * v['historical']:.2f}% | {100 * v['normal']:.2f}% | "
              f"{100 * v['cornish_fisher']:.2f}% | {v['z']:.3f} | {v['z_cf']:.3f} |")

    print("\n| level | ES historical | ES normal | worst loss in sample |")
    print("|---|---|---|---|")
    for a in LEVELS:
        v = full[a]
        print(f"| {100 * a:.0f}% | {100 * v['es_historical']:.2f}% | {100 * v['es_normal']:.2f}% | "
              f"{100 * loss.max():.2f}% on {loss.idxmax():%Y-%m-%d} |")

    # rolling one-day-ahead backtest
    lv = loss.to_numpy()
    rows = []
    for i in range(WINDOW, len(lv)):
        w = lv[i - WINDOW:i]
        rec = {"date": loss.index[i].date().isoformat(),
               "return_pct": round(100 * float(r.iloc[i]), 4),
               "loss_pct": round(100 * float(lv[i]), 4)}
        for a in LEVELS:
            v = var_three_ways(w, a)
            tag = f"{int(100 * a)}"
            for m in ("historical", "normal", "cornish_fisher"):
                rec[f"var{tag}_{m}_pct"] = round(100 * v[m], 4)
                rec[f"exc{tag}_{m}"] = int(lv[i] > v[m])
        rows.append(rec)
    bt = pd.DataFrame(rows)
    bt.to_csv(OUT, index=False)
    T = len(bt)
    print(f"\nbacktest: {WINDOW}-day rolling window, {T} one-day-ahead forecasts "
          f"({bt['date'].iloc[0]} to {bt['date'].iloc[-1]})")
    print("\n| level | method | expected exceptions | observed | rate | Kupiec LR | p |")
    print("|---|---|---|---|---|---|---|")
    for a in LEVELS:
        tag = f"{int(100 * a)}"
        for m, label in (("historical", "historical"), ("normal", "normal"), ("cornish_fisher", "Cornish-Fisher")):
            x = int(bt[f"exc{tag}_{m}"].sum())
            lr, p = kupiec(x, T, 1.0 - a)
            print(f"| {100 * a:.0f}% | {label} | {T * (1 - a):.1f} | {x} | {100 * x / T:.2f}% | "
                  f"{lr:.2f} | {p:.3f} |")

    # clustering: how many 99% historical exceptions fall in the same calendar year
    bt["year"] = bt["date"].str.slice(0, 4)
    print("\n| year | forecasts | 99% exceptions (historical) | 99% exceptions (normal) |")
    print("|---|---|---|---|")
    for y, g in bt.groupby("year"):
        print(f"| {y} | {len(g)} | {int(g['exc99_historical'].sum())} | {int(g['exc99_normal'].sum())} |")
    # independence: an unconditional coverage test cannot see clustering, so count back-to-back exceptions
    for m, label in (("historical", "historical"), ("normal", "normal"), ("cornish_fisher", "Cornish-Fisher")):
        e = bt[f"exc99_{m}"].to_numpy()
        pairs = int(((e[1:] == 1) & (e[:-1] == 1)).sum())
        print(f"99% {label}: {int(e.sum())} exceptions, {pairs} of them on the day after another exception")
    cf = bt["var99_cornish_fisher_pct"]
    print(f"99% Cornish-Fisher forecast range: {cf.min():.2f}% to {cf.max():.2f}% "
          f"(historical {bt['var99_historical_pct'].min():.2f}% to {bt['var99_historical_pct'].max():.2f}%)")
    worst = bt.loc[bt["loss_pct"].idxmax()]
    print(f"\nworst backtest day: {worst['date']} loss {worst['loss_pct']:.2f}% vs 99% forecasts "
          f"hist {worst['var99_historical_pct']:.2f}%, normal {worst['var99_normal_pct']:.2f}%, "
          f"CF {worst['var99_cornish_fisher_pct']:.2f}%")
    print(f"wrote {OUT} ({T} rows)")


if __name__ == "__main__":
    main()