VaR and Expected Shortfall with an auditable risk table

real = cu.fetch_daily(cu.SP500_EXAMPLES, start="2018-01-01", end="2026-07-01")
print(f"{real.num_rows:,} rows x {real.num_columns} columns, "
      f"{len(set(real['symbol'].to_pylist()))} symbols")
real.to_pandas().head()

schema = pa.schema(
    [
        pa.field("ts", pa.timestamp("us", tz="UTC"), nullable=False),
        pa.field("symbol", pa.string()),
        pa.field("adj_close", pa.float64()),
    ]
)
db.create_table("prices", schema, time_column="ts", sort_key=["ts", "symbol"])
db.append(
    "prices",
    real.select(["ts", "symbol", "adj_close"]).sort_by([("ts", "ascending"), ("symbol", "ascending")]),
    note="vendor backfill 2018-2026",
)

db.table("prices").select(rows=count_star(), names=col("symbol").n_unique()).to_pandas()

NOTIONAL = 10_000_000.0

WEIGHTS_SQL = """
    (VALUES ('AAPL', 0.10), ('MSFT', 0.10), ('NVDA', 0.08), ('AMZN', 0.08),
            ('GOOGL', 0.08), ('JPM', 0.10), ('XOM', 0.08), ('UNH', 0.08),
            ('PG', 0.08), ('KO', 0.06), ('CAT', 0.08), ('GS', 0.08)
    ) AS w(symbol, weight)
"""

# An inline VALUES relation has no builder verb, so this one stays SQL.
port = db.sql(
    f"""
    WITH rets AS (
        SELECT ts, symbol,
               adj_close / lag(adj_close) OVER (PARTITION BY symbol ORDER BY ts) - 1 AS ret
        FROM prices
    )
    SELECT r.ts, sum(r.ret * w.weight) AS port_ret, count(*) AS names
    FROM rets r
    JOIN {WEIGHTS_SQL} ON r.symbol = w.symbol
    WHERE r.ret IS NOT NULL
    GROUP BY r.ts
    ORDER BY r.ts
    """
).to_pandas()

# Keep only days where all 12 names have a return (guards against partial days).
port = port[port["names"] == 12].set_index("ts")
port["pnl"] = NOTIONAL * port["port_ret"]
print(f"{len(port):,} trading days, "
      f"ann. vol {port['port_ret'].std() * np.sqrt(252):.1%}, "
      f"worst day {port['pnl'].min():,.0f} USD")

from scipy import stats

r = port["port_ret"]
W = 252

def hist_var(x, q):
    return -np.quantile(x, q)

def hist_es(x, q):
    cut = np.quantile(x, q)
    return -x[x <= cut].mean()

roll = r.rolling(W)
port["var95"] = roll.apply(hist_var, args=(0.05,), raw=True).shift(1)
port["var99"] = roll.apply(hist_var, args=(0.01,), raw=True).shift(1)
port["es95"] = roll.apply(hist_es, args=(0.05,), raw=True).shift(1)
port["es99"] = roll.apply(hist_es, args=(0.01,), raw=True).shift(1)

mu, sd = roll.mean().shift(1), roll.std().shift(1)
port["var95_norm"] = -(mu + sd * stats.norm.ppf(0.05))

# Fit t df on the full standardized history, rescale to unit variance.
z = ((r - r.mean()) / r.std()).to_numpy()
nu = stats.t.fit(z, floc=0)[0]
t_scale = np.sqrt((nu - 2) / nu)
port["var95_t"] = -(mu + sd * stats.t.ppf(0.05, nu) * t_scale)

port[["var95", "es95", "var95_norm", "var95_t"]].dropna().tail(3).mul(NOTIONAL).round(0)

def kupiec_pof(returns: pd.Series, var: pd.Series, p: float) -> dict:
    mask = var.notna()
    ret, v = returns[mask], var[mask]
    n = len(ret)
    x = int((ret < -v).sum())
    pi = x / n
    ll_null = (n - x) * np.log(1 - p) + x * np.log(p)
    ll_alt = ((n - x) * np.log(1 - pi) + x * np.log(pi)) if 0 < x < n else 0.0
    lr = -2 * (ll_null - ll_alt)
    return {"n_days": n, "breaches": x, "expected": round(n * p, 1),
            "breach_rate": round(pi, 4), "LR_pof": round(lr, 2),
            "p_value": round(float(stats.chi2.sf(lr, 1)), 4)}

backtest = pd.DataFrame(
    {
        "hist 95%": kupiec_pof(r, port["var95"], 0.05),
        "hist 99%": kupiec_pof(r, port["var99"], 0.01),
        "normal 95%": kupiec_pof(r, port["var95_norm"], 0.05),
        "t 95%": kupiec_pof(r, port["var95_t"], 0.05),
    }
).T
backtest

metrics = port.dropna(subset=["var95", "var99", "es95", "es99"])[
    ["pnl", "var95", "es95", "var99", "es99"]
].copy()
for c in ["var95", "es95", "var99", "es99"]:
    metrics[c] *= NOTIONAL

mschema = pa.schema(
    [
        pa.field("ts", pa.timestamp("us", tz="UTC"), nullable=False),
        pa.field("pnl", pa.float64()),
        pa.field("var95", pa.float64()),
        pa.field("es95", pa.float64()),
        pa.field("var99", pa.float64()),
        pa.field("es99", pa.float64()),
    ]
)
db.create_table("risk_metrics", mschema, time_column="ts")

hist_part = metrics.iloc[:-60].reset_index()
db.append("risk_metrics", pa.Table.from_pandas(hist_part, schema=mschema, preserve_index=False),
          note="backfill through model go-live")

for day, row in metrics.iloc[-60:].iterrows():
    one = pa.Table.from_pandas(
        pd.DataFrame([row]).rename_axis("ts").reset_index(), schema=mschema, preserve_index=False
    )
    db.append("risk_metrics", one, note=f"eod risk {day.date()}")

[
    {k: v.get(k) for k in ("sequence", "op", "rows", "note")}
    for v in db.versions("risk_metrics")[-4:]
]

before = db.versions("risk_metrics")[-1]
compacted = db.compact("risk_metrics", note="post-EOD-loop compaction")
print(f"segments: {before['segments']} -> {compacted['segments_total']}, "
      f"op={compacted['op']}, sequence={compacted['sequence']}")

(
    db.table("risk_metrics")
    .filter(col("pnl") < -col("var95"), col("ts") >= "2024-07-01T00:00:00Z")
    .select(
        "ts",
        pnl=col("pnl").round(),
        var95_floor=(-col("var95")).round(),
        var99_floor=(-col("var99")).round(),
    )
    .sort("ts", descending=True)
    .limit(8)
    .to_pandas()
)

import matplotlib.pyplot as plt

plot_df = metrics.loc["2024-07-01":]
breach95 = plot_df[plot_df["pnl"] < -plot_df["var95"]]
breach99 = plot_df[plot_df["pnl"] < -plot_df["var99"]]

fig, ax = plt.subplots(figsize=(11, 4.5))
ax.plot(plot_df.index, plot_df["pnl"], lw=0.6, color="0.4", label="daily P&L")
ax.plot(plot_df.index, -plot_df["var95"], lw=1.2, color="tab:orange", label="-VaR 95%")
ax.plot(plot_df.index, -plot_df["var99"], lw=1.2, color="tab:red", label="-VaR 99%")
ax.scatter(breach95.index, breach95["pnl"], s=28, color="tab:orange", zorder=3,
           label=f"95% breach ({len(breach95)})")
ax.scatter(breach99.index, breach99["pnl"], s=36, color="tab:red", marker="x", zorder=4,
           label=f"99% breach ({len(breach99)})")
ax.set_title("Daily P&L vs rolling 252d historical VaR ($10M book)")
ax.set_xlabel("date")
ax.set_ylabel("USD")
ax.legend(loc="lower left", ncols=3, fontsize=8)
fig.tight_layout()