Add VaR model validation: rolling forecasts and backtests

.github/workflows/ci.yml

@@ -54,9 +54,13 @@ jobs:
           uv run quantkit simulate --output "$RUNNER_TEMP/prices.csv"
           uv run quantkit backtest "$RUNNER_TEMP/prices.csv" --strategy momentum
           uv run quantkit frontier "$RUNNER_TEMP/prices.csv" --plot "$RUNNER_TEMP/frontier.png"
+          uv run quantkit var-backtest "$RUNNER_TEMP/prices.csv" --method ewma --level 0.99
+          uv run quantkit var-backtest "$RUNNER_TEMP/prices.csv" --es-simulations 999
 
       - name: Run the example scripts
         run: |
           uv run python examples/binomial_convergence.py
           uv run python examples/efficient_frontier.py
           uv run python examples/backtest_trend.py
+          uv run python examples/var_backtest.py
+          uv run python examples/var_size_power.py

CHANGELOG.md

@@ -6,17 +6,40 @@ All notable changes to this project are documented here. The format follows
 
 ## Unreleased
 
-The `quantkit` package itself is unchanged since 0.1.0.
-
 ### Added
 
+- `quantkit.var`, VaR model validation:
+  - `rolling_var_es`: rolling one-day VaR and expected shortfall forecasts by
+    historical simulation, a Gaussian model, EWMA (RiskMetrics, decay 0.94) and
+    filtered historical simulation. Each forecast uses only the window before
+    the day it forecasts; `RiskForecast.simulate` draws from the forecast
+    distributions.
+  - `backtest_var`, `kupiec_test`, `christoffersen_test` and
+    `conditional_coverage_test`: exception counts and likelihood-ratio tests
+    with asymptotic p-values, finite with no exceptions, only exceptions or no
+    consecutive pairs. `null_rejection_rates` computes the exact probability
+    that each test rejects a correct model at a given sample length.
+  - `basel_traffic_light`: the Basel Committee's 1996 zones and capital plus
+    factors for 250 days of 99% VaR.
+  - `acerbi_szekely_z2` and `es_backtest`: Acerbi and Székely's Z2 test of
+    expected shortfall with a simulated null distribution.
+- `quantkit var-backtest` command, and the `examples/var_backtest.py` (figure)
+  and `examples/var_size_power.py` (size and power study) scripts, which CI
+  now runs.
+
 - `examples/walk_forward_risk_parity.py`: a walk-forward example that rebuilds a
   risk-parity portfolio each month from trailing returns, compares it with an
   equal-weight benchmark and exports the results as CSVs. The walkthrough is in
   `docs/WALK_FORWARD.md`.
 - A release workflow that builds the wheel and sdist from a version tag and
   attaches them, with SHA-256 checksums, to the GitHub release.
 
+### Changed
+
+- The tail size `ceil(n * (1 - level))` of historical VaR moved to a shared
+  helper, so single-sample and rolling historical VaR always use the same rule.
+  Results are unchanged.
+
 ## 0.1.0 - 2026-10-23
 
 First release.

examples/var_backtest.py

@@ -0,0 +1,128 @@
+"""Four VaR models on returns with volatility clustering, and the backtests that judge them.
+
+Run from the repository root:
+
+    uv run python examples/var_backtest.py
+
+The returns are SYNTHETIC: a GARCH(1,1) process (volatility that rises after
+large moves and decays slowly) with fat-tailed Student-t shocks, which mimics
+two well-known features of daily equity returns. Historical simulation and the
+Gaussian model weight a year of history equally, so they react slowly when
+volatility jumps and their exceptions cluster. EWMA and filtered historical
+simulation follow the volatility. The script prints the backtest of each model
+and writes docs/images/var_backtest.png.
+"""
+
+import math
+from pathlib import Path
+
+import matplotlib
+
+matplotlib.use("Agg")  # write files only; never open a window
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+
+from quantkit.var import (
+    METHODS,
+    RiskForecast,
+    backtest_var,
+    basel_traffic_light,
+    es_backtest,
+    rolling_var_es,
+)
+
+ROOT = Path(__file__).resolve().parents[1]
+OUTPUT = ROOT / "docs" / "images" / "var_backtest.png"
+LEVEL, WINDOW = 0.99, 250
+LABELS = {"historical": "historical", "gaussian": "Gaussian", "ewma": "EWMA", "fhs": "FHS"}
+COLORS = {
+    "historical": "tab:red",
+    "gaussian": "tab:orange",
+    "ewma": "tab:green",
+    "fhs": "tab:blue",
+}
+
+
+def simulate_garch_returns(n_days: int, seed: int) -> pd.Series:
+    """GARCH(1,1) returns: sigma2[t] = omega + a r[t-1]^2 + b sigma2[t-1], t(5) shocks."""
+    a, b, long_run_vol = 0.08, 0.90, 0.01
+    omega = long_run_vol**2 * (1 - a - b)
+    rng = np.random.default_rng(seed)
+    shocks = rng.standard_t(5, n_days) * math.sqrt(3 / 5)  # unit variance
+    returns = np.empty(n_days)
+    variance = long_run_vol**2
+    for t in range(n_days):
+        returns[t] = math.sqrt(variance) * shocks[t]
+        variance = omega + a * returns[t] ** 2 + b * variance
+    return pd.Series(returns, index=pd.bdate_range("2018-01-01", periods=n_days, name="date"))
+
+
+def summarize(forecast: RiskForecast) -> dict[str, object]:
+    result = backtest_var(forecast.exceptions, level=forecast.level)
+    light = basel_traffic_light(int(forecast.last(250).exceptions.sum()))
+    z2 = es_backtest(forecast, n_simulations=9_999, seed=0)
+    return {
+        "exceptions": result.n_exceptions,
+        "expected": round(result.expected_exceptions, 1),
+        "Kupiec p": round(result.unconditional.p_value, 4),
+        "independence p": round(result.independence.p_value, 4),
+        "cond. coverage p": round(result.conditional.p_value, 4),
+        "Basel zone (last 250 days)": f"{light.zone} ({light.n_exceptions})",
+        "Z2": round(z2.statistic, 2),
+        "Z2 p": round(z2.p_value, 4),
+    }
+
+
+def main() -> None:
+    """Forecast with all four methods, backtest them on common days, and plot."""
+    returns = simulate_garch_returns(1_750, seed=3)
+    forecasts = {m: rolling_var_es(returns, method=m, level=LEVEL, window=WINDOW) for m in METHODS}
+    n_days = min(len(f) for f in forecasts.values())  # FHS needs two windows of history
+    forecasts = {m: f.last(n_days) for m, f in forecasts.items()}
+    days = forecasts["fhs"].var.index
+    print(
+        f"SYNTHETIC GARCH(1,1) returns with t(5) shocks; {LEVEL:.0%} one-day VaR and ES, "
+        f"{WINDOW}-day window; backtest over {n_days} days ({days[0].date()} to "
+        f"{days[-1].date()})\n"
+    )
+    table = pd.DataFrame({LABELS[m]: summarize(f) for m, f in forecasts.items()})
+    print(table.to_string())
+
+    fig, (top, bottom) = plt.subplots(
+        2, 1, figsize=(9, 6.5), sharex=True, gridspec_kw={"height_ratios": [3, 2]}
+    )
+    realized = forecasts["fhs"].returns
+    top.plot(days, realized, color="0.65", lw=0.6, label="daily return")
+    for method in ("historical", "fhs"):
+        forecast, color = forecasts[method], COLORS[method]
+        top.plot(days, -forecast.var, color=color, lw=1.1, label=f"-VaR, {LABELS[method]}")
+        hit = forecast.exceptions
+        top.scatter(days[hit], realized[hit], color=color, s=14, zorder=3)
+    top.set(ylabel="return", title="SYNTHETIC data: 99% one-day VaR and its exceptions (dots)")
+    top.yaxis.set_major_formatter("{x:.0%}")
+    top.legend(loc="lower left", fontsize=8)
+    top.grid(alpha=0.3)
+
+    for method, forecast in forecasts.items():
+        bottom.step(
+            days,
+            forecast.exceptions.cumsum(),
+            where="post",
+            color=COLORS[method],
+            label=LABELS[method],
+        )
+    expected = (1 - LEVEL) * np.arange(1, n_days + 1)
+    bottom.plot(days, expected, color="black", ls="--", lw=1, label="expected (1% of days)")
+    bottom.set(ylabel="cumulative exceptions")
+    bottom.legend(loc="upper left", fontsize=8, ncol=3)
+    bottom.grid(alpha=0.3)
+    fig.tight_layout()
+    OUTPUT.parent.mkdir(parents=True, exist_ok=True)
+    fig.savefig(OUTPUT, dpi=110)
+    print(f"\nsaved {OUTPUT.relative_to(ROOT)}")
+
+
+if __name__ == "__main__":
+    main()

examples/var_size_power.py

@@ -0,0 +1,115 @@
+"""How often do the VaR and ES backtests reject a correct model, and a wrong one?
+
+Run from the repository root:
+
+    uv run python examples/var_size_power.py
+
+All data are SYNTHETIC. The script simulates many independent return paths,
+forecasts VaR and ES on each with a 250-day window, runs the 5% tests on the
+last 250 or 1,000 days, and prints the share of paths each test rejects.
+
+* Size: returns drawn from the EWMA (RiskMetrics) model itself, so the EWMA
+  forecast is the true conditional distribution and every rejection is a false
+  alarm. The exact rejection probability of each likelihood-ratio test comes
+  from quantkit.var.null_rejection_rates. The Z2 Monte Carlo test with 199
+  simulations has an exact 5% size, because 5% of 199 + 1 is a whole number.
+* Power: Gaussian forecasts on Student-t returns with 4 degrees of freedom,
+  whose tails are fatter than the model assumes, next to the same forecasts on
+  normal returns (the Gaussian model's own assumption, with estimation error).
+"""
+
+import math
+from collections.abc import Callable
+
+import numpy as np
+
+from quantkit.var import (
+    Method,
+    backtest_var,
+    es_backtest,
+    ewma_weights,
+    null_rejection_rates,
+    rolling_var_es,
+)
+
+N_PATHS = 1_000
+WINDOW = 250
+Z2_SIMULATIONS = 199  # 0.05 * (199 + 1) = 10, so 5% is an exact Monte Carlo cutoff
+Sampler = Callable[[np.random.Generator, int], np.ndarray]
+
+
+def ewma_model(rng: np.random.Generator, n_days: int) -> np.ndarray:
+    """Returns from the EWMA model: r[t] = sigma[t] z[t], sigma from the past 250 returns."""
+    weights = ewma_weights(WINDOW)
+    paths = np.empty((N_PATHS, n_days))
+    paths[:, :WINDOW] = 0.01 * rng.standard_normal((N_PATHS, WINDOW))
+    for t in range(WINDOW, n_days):
+        sigma = np.sqrt(np.square(paths[:, t - WINDOW : t]) @ weights)
+        paths[:, t] = sigma * rng.standard_normal(N_PATHS)
+    return paths
+
+
+def normal(rng: np.random.Generator, n_days: int) -> np.ndarray:
+    return 0.01 * rng.standard_normal((N_PATHS, n_days))
+
+
+def student_t4(rng: np.random.Generator, n_days: int) -> np.ndarray:
+    return 0.01 * rng.standard_t(4, (N_PATHS, n_days)) * math.sqrt(2 / 4)  # 1% daily vol
+
+
+def rejection_rates(
+    sampler: Sampler, method: Method, level: float, n_test: int, rng: np.random.Generator
+) -> list[float]:
+    """Share of paths rejected by LR_uc, LR_ind, LR_cc and Z2 at 5%, and the exception rate."""
+    rejected = np.zeros(4)
+    exceptions = 0
+    for path in sampler(rng, WINDOW + n_test):
+        forecast = rolling_var_es(path, method=method, level=level, window=WINDOW)
+        result = backtest_var(forecast.exceptions, level=level)
+        z2 = es_backtest(forecast, n_simulations=Z2_SIMULATIONS, seed=rng)
+        rejected += [
+            result.unconditional.p_value < 0.05,
+            result.independence.p_value < 0.05,
+            result.conditional.p_value < 0.05,
+            z2.p_value <= 0.05,
+        ]
+        exceptions += result.n_exceptions
+    return [*(rejected / N_PATHS), exceptions / (N_PATHS * n_test)]
+
+
+def row(label: str, values: list[float]) -> str:
+    cells = "".join(f"{v:>9.1%}" for v in values[:4])
+    rate = f"{values[4]:>11.2%}" if len(values) > 4 else ""
+    return f"{label:<46}{cells}{rate}"
+
+
+def main() -> None:
+    """Run the size and power simulations and print one table."""
+    rng = np.random.default_rng(2026)
+    print(
+        f"Share of {N_PATHS:,} SYNTHETIC paths rejected by each 5% test "
+        f"(Monte Carlo standard error at most {0.5 / math.sqrt(N_PATHS):.1%})\n"
+    )
+    print(f"{'':<46}{'LR_uc':>9}{'LR_ind':>9}{'LR_cc':>9}{'Z2':>9}{'exceptions':>11}")
+    print("size: EWMA forecasts on returns from the EWMA model, 99%")
+    for n_test in (250, 1_000):
+        measured = rejection_rates(ewma_model, "ewma", 0.99, n_test, rng)
+        exact = null_rejection_rates(n_test, level=0.99)
+        print(row(f"  {n_test:,} days, measured", measured))
+        print(row(f"  {n_test:,} days, exact", [*vars_of(exact), 0.05]))
+    for level in (0.99, 0.975):
+        print(f"power: Gaussian forecasts, {100 * level:g}%")
+        tests = (250, 1_000) if level == 0.99 else (1_000,)
+        for n_test in tests:
+            for name, sampler in (("normal", normal), ("Student-t(4)", student_t4)):
+                values = rejection_rates(sampler, "gaussian", level, n_test, rng)
+                print(row(f"  {name} returns, {n_test:,} days", values))
+
+
+def vars_of(rates: object) -> list[float]:
+    """The three exact rejection probabilities, in table order."""
+    return [getattr(rates, name) for name in ("unconditional", "independence", "conditional")]
+
+
+if __name__ == "__main__":
+    main()

pyproject.toml

@@ -1,7 +1,7 @@
 [project]
 name = "quantkit"
 version = "0.1.0"
-description = "Options pricing, portfolio optimization and backtesting in Python, with every model checked against closed-form results, textbook values or an independent method."
+description = "Options pricing, portfolio optimization, backtesting and VaR model validation in Python, with every model checked against closed-form results, textbook values or an independent method."
 readme = "README.md"
 license = "MIT"
 license-files = ["LICENSE"]
@@ -14,6 +14,7 @@ keywords = [
     "portfolio-optimization",
     "risk-parity",
     "backtesting",
+    "value-at-risk",
 ]
 classifiers = [
     "Development Status :: 3 - Alpha",

src/quantkit/__init__.py

@@ -1,4 +1,4 @@
-"""quantkit: options pricing, portfolio optimization and backtesting.
+"""quantkit: options pricing, portfolio optimization, backtesting and VaR validation.
 
 Subpackages:
     quantkit.pricing: Black-Scholes-Merton, implied volatility, binomial trees and
@@ -7,6 +7,8 @@
         optimization, risk parity and risk metrics.
     quantkit.backtest: a daily vectorized backtester, example strategies and
         performance metrics.
+    quantkit.var: rolling VaR and expected shortfall forecasts, exception
+        backtests, the Basel traffic light and an expected shortfall backtest.
     quantkit.data: synthetic price generation and CSV loading.
 
 This is educational and research software, not investment advice.

src/quantkit/_validation.py

@@ -20,6 +20,20 @@
 SeedLike = int | np.random.Generator | None
 """Anything accepted by :func:`numpy.random.default_rng`: a seed, a generator or None."""
 
+# (1 - 0.99) * 500 evaluates to 5.000000000000004 in binary floating point; this
+# slack stops math.ceil from turning that into 6.
+_FLOAT_SLACK = 1e-9
+
+
+def tail_count(n_obs: int, level: float) -> int:
+    """Number of observations in the empirical ``1 - level`` tail: ``ceil(n * (1 - level))``.
+
+    At least one observation is always in the tail. Historical VaR and CVaR use
+    this count both for a single sample and for every rolling window, so the two
+    always agree.
+    """
+    return max(1, math.ceil(n_obs * (1.0 - level) - _FLOAT_SLACK))
+
 
 def require_finite(name: str, value: float) -> float:
     """Return ``value`` as a float, or raise if it is NaN or infinite."""