h5i-db / Cookbook / Risk & production / VaR and Expected Shortfall with an auditable risk table

VaR and Expected Shortfall with an auditable risk table

A risk number nobody can reproduce is a liability. This recipe computes historical and parametric VaR / Expected Shortfall for a fixed equity book on real daily data, backtests the VaR with Kupiec's POF test, and then does the part most notebooks skip: it writes the daily risk metrics into a versioned h5i-db table, one commit per production day, so that "what VaR did we report on date X, and from what inputs?" is a query, not an email thread.

h5i-db features doing the work: SQL window functions for the return pipeline, append-only commits as the audit trail (versions() shows every EOD run with its note and wall-clock commit time), and compact() to keep the table healthy after many small daily appends.

import numpy as np
import pandas as pd
import pyarrow as pa
import h5i_db

import cookbook_utils as cu

db = h5i_db.Database(cu.fresh_db("prod_var"), create=True)

1. Load real prices into a versioned table#

Eight-plus years of daily data for 30 liquid US names (cached Yahoo data). We store adj_close, a total-return-ish series, which is what P&L risk should be computed on. One bulk append = one atomic commit.

One strictness to know about: with sort_key=["ts", "symbol"], append requires the input sorted by the full key. Daily panels have 30 symbols sharing each timestamp, so a ts-only sort is not enough; sort by (ts, symbol) before appending.

real = cu.fetch_daily(cu.SP500_EXAMPLES, start="2018-01-01", end="2026-07-01")

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.sql("SELECT count(*) AS rows, count(DISTINCT symbol) AS names FROM prices").to_pandas()
output
rows names
0 64020 30

2. Portfolio P&L series in one SQL statement#

The book is a fixed 12-name portfolio (weights sum to 1, $10M notional). Daily simple returns come from a lag() window per symbol; the portfolio return is a weighted sum via a join against an inline VALUES weights relation. Everything streams on time-sorted storage, with no pandas reshaping until we actually need rolling quantiles.

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

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")
output
2,133 trading days, ann. vol 20.9%, worst day -1,154,191 USD

3. Rolling historical VaR / ES, plus parametric flavors#

All measures use a trailing 252-day window shifted by one day: the VaR reported for day t only uses information through t−1. That shift is the difference between a risk model and a hindsight model, and it is exactly the kind of detail an auditor will ask about.

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)
output
var95 es95 var95_norm var95_t
ts
2026-06-26 20:00:00+00:00 123777.0 166664.0 114636.0 93226.0
2026-06-29 20:00:00+00:00 123777.0 166664.0 114700.0 93303.0
2026-06-30 20:00:00+00:00 123777.0 166664.0 114576.0 93233.0

4. Kupiec POF backtest: is 95% really 95%?#

Count days where the realized loss exceeded the previous day's VaR and run Kupiec's proportion-of-failures likelihood-ratio test (asymptotically chi-squared with 1 dof). A well-calibrated 95% VaR should be breached on about 5% of days; rejecting the null means the model's coverage is wrong.

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
output
n_days breaches expected breach_rate LR_pof p_value
hist 95% 1881.0 98.0 94.1 0.0521 0.17 0.678
hist 99% 1881.0 33.0 18.8 0.0175 8.83 0.003
normal 95% 1881.0 98.0 94.1 0.0521 0.17 0.678
t 95% 1881.0 136.0 94.1 0.0723 17.41 0.000

Reading the table honestly: the historical and normal models hold coverage at 95%, but the 99% historical VaR under-covers (breaches cluster in vol-regime shifts, which a trailing empirical quantile is slow to catch; Kupiec tests coverage only, not independence, and Christoffersen's test would flag the clustering). And note the trap in the "fat-tail fix": a Student-t rescaled to unit variance has less extreme quantiles than the normal at the 5% level, and the fat tails only cross over around 1%. A t-model naively swapped in at 95% is anti-conservative, and the test catches it.

5. The audit trail: a versioned risk_metrics table#

Production pattern: backfill history in one commit, then one commit per EOD run. Each daily append is atomic and carries a note; versions() then answers "which risk numbers were published on which day, and when exactly", with wall-clock commit timestamps and without any extra logging infrastructure. We simulate the last 60 production days as individual commits.

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:]
]
output
[{'sequence': 58, 'op': 'append', 'rows': 1878, 'note': 'eod risk 2026-06-25'},
 {'sequence': 59, 'op': 'append', 'rows': 1879, 'note': 'eod risk 2026-06-26'},
 {'sequence': 60, 'op': 'append', 'rows': 1880, 'note': 'eod risk 2026-06-29'},
 {'sequence': 61, 'op': 'append', 'rows': 1881, 'note': 'eod risk 2026-06-30'}]

Sixty tiny commits also means sixty tiny Parquet segments. compact() merges them into efficient segments, as its own audited commit, so even maintenance shows up in the history.

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']}")
output
segments: 61 -> 1, op=compact, sequence=62

The risk table is now queryable like any other: here, every breach day of the last two years straight from SQL:

db.sql(
    """
    SELECT ts, round(pnl) AS pnl, round(-var95) AS var95_floor, round(-var99) AS var99_floor
    FROM risk_metrics
    WHERE pnl < -var95 AND ts >= '2024-07-01T00:00:00Z'
    ORDER BY ts DESC
    LIMIT 8
    """
).to_pandas()
output
ts pnl var95_floor var99_floor
0 2026-06-17 20:00:00+00:00 -120366.0 -115987.0 -184759.0
1 2026-06-10 20:00:00+00:00 -148196.0 -104481.0 -184759.0
2 2026-06-05 20:00:00+00:00 -138137.0 -102827.0 -184759.0
3 2026-03-27 20:00:00+00:00 -157587.0 -147858.0 -274463.0
4 2026-03-26 20:00:00+00:00 -147359.0 -147281.0 -274463.0
5 2026-03-12 20:00:00+00:00 -151007.0 -138104.0 -274463.0
6 2026-02-12 20:00:00+00:00 -168461.0 -158913.0 -281868.0
7 2026-02-05 20:00:00+00:00 -171680.0 -147281.0 -281868.0

6. P&L against the VaR band, breaches marked#

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()
output
output figure

Takeaways#

db.close()