diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 9f581e8..90617a1 100644 --- a/.github/workflows/ci.yml +++ b/.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 diff --git a/CHANGELOG.md b/CHANGELOG.md index 302ca50..becd821 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,10 +6,27 @@ 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 @@ -17,6 +34,12 @@ The `quantkit` package itself is unchanged since 0.1.0. - 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. diff --git a/README.md b/README.md index 54717d9..120cee7 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # quantkit -Options pricing, portfolio optimization and backtesting in Python. Every model is checked against closed-form results, textbook values or an independent method. +Options pricing, portfolio optimization, backtesting and VaR model validation in Python. Every model is checked against closed-form results, textbook values or an independent method. [![CI](https://github.com/ult666666/quantkit/actions/workflows/ci.yml/badge.svg)](https://github.com/ult666666/quantkit/actions/workflows/ci.yml) [![License: MIT](https://img.shields.io/badge/license-MIT-blue.svg)](LICENSE) @@ -29,8 +29,13 @@ Quantitative finance code fails quietly: a sign slip in a Greek, an off-by-one d - Daily vectorized engine: weights decided at the close of day *t* are held from day *t*+1, with drift-aware turnover and proportional costs in basis points. - Strategies: buy-and-hold, periodically rebalanced equal weight, moving-average crossover, time-series momentum and inverse volatility. - Metrics: CAGR, volatility, Sharpe, Sortino, maximum drawdown and its duration, Calmar, hit rate and turnover, plus a look-ahead checker for custom strategies. +- **VaR model validation** (`quantkit.var`) + - Rolling one-day Value at Risk (VaR) and expected shortfall (ES, the mean loss beyond VaR) forecasts by historical simulation, a Gaussian model, EWMA (RiskMetrics, decay 0.94) and filtered historical simulation, each from a trailing window that never contains the return it is scored against. + - Exception backtests: Kupiec's proportion-of-failures test and Christoffersen's independence and conditional coverage tests, with p-values, plus the exact probability that each test rejects a correct model at a given sample length. + - The Basel traffic light: green, yellow and red zones and the capital plus factors for 250 days of 99% VaR. + - Acerbi and Székely's Z2 test of expected shortfall, with p-values from a null distribution simulated from the forecasts themselves. - **Data** (`quantkit.data`): a seeded correlated-GBM price simulator and a validating CSV loader (long or wide format). -- **CLI**: `quantkit price | iv | simulate | backtest | frontier`. +- **CLI**: `quantkit price | iv | simulate | backtest | frontier | var-backtest`. - Fully typed (`mypy --strict`, ships `py.typed`) and linted with ruff; the CI workflow targets Python 3.11–3.13. ## Quick start @@ -106,6 +111,34 @@ CAGR 0.36%, Sharpe 0.08, max drawdown 14.9% The maximum-Sharpe portfolio leaves out `CCC` because its *estimated* mean return over this sample is negative (−10.9% a year), even though it was simulated with a +8% drift: five years of data pin down expected returns very poorly. That estimation error, not the optimizer, is the main weakness of mean–variance optimization in practice. +Forecast the equal-weight portfolio's one-day 99% VaR and ES each day from the previous 250 days, then backtest the forecasts: + +```python +>>> from quantkit.var import ( +... backtest_var, basel_traffic_light, es_backtest, null_rejection_rates, rolling_var_es, +... ) +>>> portfolio = returns.mean(axis=1) # daily-rebalanced equal weights +>>> forecast = rolling_var_es(portfolio, method="fhs", level=0.99, window=250) +>>> print(f"{len(forecast)} forecasts, {forecast.exceptions.sum()} exceptions; " +... f"latest VaR {forecast.var.iloc[-1]:.2%}, ES {forecast.es.iloc[-1]:.2%}") +759 forecasts, 8 exceptions; latest VaR 1.50%, ES 1.71% +>>> result = backtest_var(forecast.exceptions, level=0.99) +>>> print(f"LR_uc p={result.unconditional.p_value:.3f}, LR_ind p={result.independence.p_value:.3f}, " +... f"LR_cc p={result.conditional.p_value:.3f}") +LR_uc p=0.882, LR_ind p=0.680, LR_cc p=0.908 +>>> basel_traffic_light(int(forecast.last(250).exceptions.sum())) +TrafficLight(zone='green', n_exceptions=3, plus_factor=0.0) +>>> z2 = es_backtest(forecast, n_simulations=999, seed=0) +>>> print(f"Z2 {z2.statistic:.2f}, simulated p-value {z2.p_value:.3f}") +Z2 -0.13, simulated p-value 0.164 +>>> rates = null_rejection_rates(250, level=0.99) # how often a correct model fails at 5% +>>> print(f"{rates.unconditional:.3f} {rates.independence:.3f} {rates.conditional:.3f}") +0.095 0.014 0.011 + +``` + +FHS (filtered historical simulation) needs two 250-day windows of history, so the 1,259 returns give 759 forecasts. The last line is a warning about the tests themselves: with 2.5 exceptions expected in 250 days, the nominal 5% Kupiec test rejects a correct model 9.5% of the time, and the two Christoffersen tests about 1%. + ## Usage ### Command line @@ -195,6 +228,32 @@ max Sharpe 7.2% 16.2% 0.32 0.0% 0.0% 81.7% 18.3% 0.0% Add `--plot frontier.png` to also save a chart (needs the `plot` extra), and, for example, `--min-weight -0.2 --max-weight 0.6` to change the weight bounds. +Backtest rolling VaR and ES forecasts of the equal-weight portfolio (or one column, with `--ticker`): + +```console +$ uv run quantkit var-backtest prices.csv --method ewma --level 0.99 +equal-weight portfolio of 5 assets, rebalanced daily: 1259 daily returns, 2020-01-02 to 2024-10-29 +99% one-day VaR and ES, 250-day window, EWMA decay 0.94 +backtest: 1009 days, 2020-12-17 to 2024-10-29 + + EWMA +exceptions 14 +expected exceptions 10.1 +exception rate 1.39% +Kupiec LR_uc p-value 0.243 +Christoffersen LR_ind p-value 0.185 +conditional coverage LR_cc p-value 0.209 +Basel zone, last 250 days green (3) +Acerbi-Szekely Z2 -0.42 +Z2 simulated p-value 0.102 + +LR p-values use chi-squared limits. Over 1009 days at 99%, a correct model fails the 5% tests +with probability 5.5% (LR_uc), 1.7% (LR_ind) and 2.7% (LR_cc). +Z2 p-values compare with 9,999 paths simulated from each method's own forecasts (seed 0). +``` + +Without `--method` it compares all four methods on the days they all cover. An exception is a day whose loss exceeded that day's VaR forecast; at 99% a correct model expects one in a hundred days, independent of each other. The Basel zone counts the exceptions of the last 250 days, as bank supervisors do. Z2 is negative when losses beyond VaR are larger than the ES forecasts promised. + ### Examples and figures The scripts in [`examples/`](examples) regenerate the figures below (`uv run python examples/.py`). The walk-forward example writes CSVs instead of a figure. @@ -217,6 +276,48 @@ Seed 42 happens to be one of the 14. With no trends in the data, the crossover m ![Equity curves and drawdowns on synthetic data](docs/images/backtest_equity.png) +**VaR backtest** ([`var_backtest.py`](examples/var_backtest.py)), on **synthetic** data: GARCH(1,1) returns, whose volatility jumps after large moves, with fat-tailed Student-t shocks. Historical simulation weights a year of history equally, so its VaR (red) reacts late to the volatility burst and its exceptions cluster; filtered historical simulation (blue) rescales history by the current EWMA volatility. Over 1,250 days at 99%: + +```text + historical Gaussian EWMA FHS +exceptions 20 24 24 16 +expected 12.5 12.5 12.5 12.5 +Kupiec p 0.0499 0.0037 0.0037 0.3403 +independence p 0.0389 0.0833 0.4815 0.5193 +cond. coverage p 0.0172 0.0033 0.0115 0.5142 +Basel zone (last 250 days) yellow (5) yellow (7) green (2) green (2) +Z2 -0.79 -1.25 -1.14 -0.33 +Z2 p 0.0014 0.0002 0.0005 0.0484 +``` + +All three exception tests reject historical simulation, and the independence test shows why: its exceptions bunch up in the volatile period. EWMA follows the volatility (independence p = 0.48) but keeps normal tails, so it has as many exceptions as the Gaussian model, 24 against 12.5 expected. FHS passes the three exception tests; Z2 still rejects it, narrowly, at 5%. + +![VaR forecasts, exceptions and cumulative exception counts on synthetic data](docs/images/var_backtest.png) + +**Size and power of the backtests** ([`var_size_power.py`](examples/var_size_power.py)), on **synthetic** data. Each row simulates 1,000 independent paths, forecasts with a 250-day window and runs the 5% tests on the last 250 or 1,000 days (`uv run python examples/var_size_power.py`, about 40 seconds): + +```text +Share of 1,000 SYNTHETIC paths rejected by each 5% test (Monte Carlo standard error at most 1.6%) + + LR_uc LR_ind LR_cc Z2 exceptions +size: EWMA forecasts on returns from the EWMA model, 99% + 250 days, measured 10.5% 1.8% 1.0% 5.1% 1.00% + 250 days, exact 9.5% 1.4% 1.1% 5.0% + 1,000 days, measured 4.2% 2.2% 2.2% 4.4% 1.00% + 1,000 days, exact 5.5% 1.8% 2.7% 5.0% +power: Gaussian forecasts, 99% + normal returns, 250 days 8.1% 1.3% 0.8% 5.3% 1.02% + Student-t(4) returns, 250 days 13.0% 1.4% 6.8% 42.2% 1.65% + normal returns, 1,000 days 2.2% 1.2% 0.9% 3.9% 1.05% + Student-t(4) returns, 1,000 days 46.9% 1.9% 31.3% 87.3% 1.63% +power: Gaussian forecasts, 97.5% + normal returns, 1,000 days 1.2% 2.3% 2.1% 4.0% 2.57% + Student-t(4) returns, 1,000 days 3.4% 2.0% 3.2% 43.3% 2.70% +``` + +- **Size.** The paths are drawn from the EWMA model itself, so the forecasts are exactly right and every rejection is a false alarm. The measured rates agree with the exact ones from `null_rejection_rates` within the Monte Carlo error, and they are not 5%: with 2.5 expected exceptions in 250 days, zero exceptions (8.1% likely) already fails Kupiec's test, and the independence test rarely has a pair of exceptions to work with. +- **Power.** Student-t returns with 4 degrees of freedom have the same variance as the normal ones but fatter tails, so Gaussian VaR is too low. Exception counts catch this in only 13% of 250-day samples. Z2, which also weighs how large the tail losses are, catches it in 42%, and in 87% of 1,000-day samples. At 97.5% (the level of the Basel FRTB expected shortfall) the Gaussian VaR happens to have about the right number of exceptions, so the count-based tests see nothing (3.4%) while Z2 rejects 43% of samples. + **Walk-forward portfolio construction** ([`walk_forward_risk_parity.py`](examples/walk_forward_risk_parity.py)), on **synthetic** data. This example connects trailing-return estimation, risk-parity allocation and backtesting. Each monthly allocation uses only closes before the rebalance date. A monthly equal-weight benchmark shares the same cash warmup, dates and costs. ```bash @@ -231,6 +332,7 @@ The command exports prices, both target tables, estimation windows, net wealth a flowchart LR CSV["CSV file or synthetic GBM
quantkit.data"] --> PX["wide price table
(dates x tickers)"] PX --> PF["quantkit.portfolio
returns, covariance, optimizers, VaR"] + PF --> VR["quantkit.var
rolling VaR/ES forecasts, exception tests,
Basel traffic light, ES backtest"] PX --> ST["strategy
target weights per close"] ST --> EN["backtest engine
weights held from t+1, costs"] EN --> MT["metrics
CAGR, Sharpe, drawdown, ..."] @@ -247,6 +349,8 @@ flowchart LR The engine is vectorized: between trades the share counts are constant, so each day's weights follow from the weights at the last trade and the price growth since then. +- **VaR model validation.** Every forecaster is a location-scale model `r[t] = mu[t] + sigma[t] e[t]` fitted to the window before day `t`, with normal or resampled innovations `e`; VaR, ES and the simulated null of the ES backtest all come from that one description. The likelihood-ratio tests use `0 log 0 = 0`, so they stay finite with no exceptions or no consecutive pairs, and a dynamic program over the exception sequence gives their exact size. + [docs/DESIGN.md](docs/DESIGN.md) explains the architecture, the main design decisions and their trade-offs in more detail. ## Correctness @@ -263,11 +367,15 @@ flowchart LR | Mean–variance | Closed-form minimum-variance, tangency and target-return portfolios (unbounded case); first-order conditions; no random feasible portfolio beats the optimizer; budget and bounds respected | | Risk parity | Equal risk contributions; inverse-volatility closed forms (uncorrelated or two assets); volatility between minimum-variance and 1/N | | VaR / CVaR | Hand-computed tails; Gaussian closed forms; historical converges to parametric on normal data | +| Rolling VaR / ES forecasts | Every window agrees with the single-sample `historical_var` and `parametric_var`; the RiskMetrics recursion; a Hull–White volatility-rescaling loop for FHS; **no look-ahead**: changing the return a forecast is scored against (and all later ones) leaves the forecast unchanged, and a peeking forecaster fails the same check; exception probabilities on i.i.d. data match closed forms (`m / (window + 1)` for historical simulation, a Student-t prediction interval for the Gaussian model) | +| Exception tests | Hand-computed Kupiec and Christoffersen statistics; SciPy's binomial log-likelihoods and `chi2` p-values; `LR_cc = LR_uc + LR_ind`; exact sizes against brute-force enumeration of all 2^T exception sequences; simulated sizes on returns drawn from the forecast model itself | +| Basel traffic light | The zones, plus factors and cumulative probabilities of the 1996 Basel table, and the 95% / 99.99% binomial rule behind them | +| ES backtest (Z2) | A hand calculation; mean zero when returns follow the forecasts; a null distribution that does not depend on the volatility path; power against fat tails | | Backtester | An independent day-by-day reference loop; buy-and-hold equals the asset return; costs reduce wealth by exactly the computed amount; a hand-computed two-day example; **no look-ahead**: perturbing prices after day *t* leaves every weight in effect through day *t*+1 unchanged, and a deliberately cheating strategy fails the same check | | Metrics | Hand-computed series (drawdown and its duration, CAGR, Sharpe, Sortino, Calmar, hit rate) | | Data | Simulated means, volatilities and correlations match their targets; CSV round trips and specific error messages | -Monte Carlo and other statistical tests use fixed seeds, so they cannot fail by chance; the Hypothesis property tests draw fresh random inputs on every run (200 per test, or 5,000 with `pytest --hypothesis-profile=thorough`). The CI workflow runs ruff, `mypy --strict`, the tests with a 90% coverage floor, the README doctests and CLI smoke tests on Python 3.11, 3.12 and 3.13. +Monte Carlo and other statistical tests use fixed seeds, so they cannot fail by chance; the Hypothesis property tests draw fresh random inputs on every run (200 per test, or 5,000 with `pytest --hypothesis-profile=thorough`). The CI workflow runs ruff, `mypy --strict`, the tests with a 90% coverage floor, the README doctests, CLI smoke tests and the example scripts on Python 3.11, 3.12 and 3.13. ## Development @@ -286,11 +394,13 @@ uv run python examples/backtest_trend.py # regenerate a figure in docs/images **Portfolio.** Expected returns are plain historical means, which are very noisy. The same weight bounds apply to every asset, with no per-asset, sector or turnover constraints. SLSQP is a general-purpose solver that works well for tens of assets; a dedicated QP solver would scale better. Risk parity is long-only. VaR is one-period and in return units. +**VaR models and backtests.** Forecasts cover one day and one return series, in return units: there is no multi-day scaling, no P&L from positions and no revaluation of options. The EWMA decay is fixed rather than estimated (no GARCH fit), and FHS filters with that EWMA volatility. The likelihood-ratio p-values are asymptotic; `null_rejection_rates` reports how far the true size is from the nominal one, but the tests do not yet offer exact or Monte Carlo p-values. The independence test only looks at consecutive days. The traffic light implements the 1996 table for 250 days of 99% VaR, not the 97.5% zones or the P&L attribution test of the Basel FRTB standards. Z2's null distribution is simulated from each method's own forecasts, so it tests the forecast distribution's tail as a whole and needs positive ES forecasts. + **Backtesting.** Daily close-to-close only: trades fill at the same close that produced the signal, and the only friction is a proportional cost (no slippage, market impact, borrowing or financing costs). Idle cash earns nothing. Prices must be complete (no missing values and no assets entering mid-sample) and should already be adjusted for splits and dividends. The no-look-ahead guarantee covers the engine and the built-in strategies; custom strategies can be checked with `check_no_lookahead`. **Data and verification.** There is deliberately no data downloader: bring your own CSV with ISO dates. CI runs lint, strict type checks and the full suite with README doctests on Python 3.11, 3.12 and 3.13. -**Roadmap ideas:** per-asset and linear constraints, a vectorized BSM interface, Longstaff–Schwartz for American options, Black–Litterman expected returns, and an engine option to fill trades after the signal close. The walk-forward example lags its estimation window by one close, but its trades still fill at the signal close. +**Roadmap ideas:** per-asset and linear constraints, a vectorized BSM interface, Longstaff–Schwartz for American options, Black–Litterman expected returns, an engine option to fill trades after the signal close, Monte Carlo p-values and duration-based tests (Christoffersen & Pelletier 2004) for VaR exceptions, and GARCH-filtered historical simulation. The walk-forward example lags its estimation window by one close, but its trades still fill at the signal close. ## References @@ -307,6 +417,17 @@ uv run python examples/backtest_trend.py # regenerate a figure in docs/images - Maillard, S., Roncalli, T. & Teïletche, J. (2010). The properties of equally weighted risk contribution portfolios. *Journal of Portfolio Management*, 36(4), 60–70. - Griveau-Billion, T., Richard, J.-C. & Roncalli, T. (2013). A fast algorithm for computing high-dimensional risk parity portfolios. arXiv:1311.4057. - Moskowitz, T. J., Ooi, Y. H. & Pedersen, L. H. (2012). Time series momentum. *Journal of Financial Economics*, 104(2), 228–250. +- J.P. Morgan/Reuters (1996). *RiskMetrics Technical Document*, 4th edition. +- Hull, J. & White, A. (1998). Incorporating volatility updating into the historical simulation method for value-at-risk. *Journal of Risk*, 1(1), 5–19. +- Barone-Adesi, G., Giannopoulos, K. & Vosper, L. (1999). VaR without correlations for portfolios of derivative securities. *Journal of Futures Markets*, 19(5), 583–602. +- Kupiec, P. H. (1995). Techniques for verifying the accuracy of risk measurement models. *Journal of Derivatives*, 3(2), 73–84. +- Christoffersen, P. F. (1998). Evaluating interval forecasts. *International Economic Review*, 39(4), 841–862. +- Christoffersen, P. & Pelletier, D. (2004). Backtesting value-at-risk: A duration-based approach. *Journal of Financial Econometrics*, 2(1), 84–108. +- Basel Committee on Banking Supervision (1996). *Supervisory framework for the use of "backtesting" in conjunction with the internal models approach to market risk capital requirements*. [bis.org/publ/bcbs22.htm](https://www.bis.org/publ/bcbs22.htm) +- Acerbi, C. & Székely, B. (2014). Backtesting expected shortfall. *Risk*, December 2014. +- Gneiting, T. (2011). Making and evaluating point forecasts. *Journal of the American Statistical Association*, 106(494), 746–762. +- Jorion, P. (2007). *Value at Risk: The New Benchmark for Managing Financial Risk*, 3rd edition. McGraw-Hill (chapter on backtesting VaR). +- Davison, A. C. & Hinkley, D. V. (1997). *Bootstrap Methods and their Application*. Cambridge University Press (Monte Carlo tests). ## License diff --git a/docs/DESIGN.md b/docs/DESIGN.md index 0cb1102..f7e53f0 100644 --- a/docs/DESIGN.md +++ b/docs/DESIGN.md @@ -6,7 +6,7 @@ This document explains how quantkit is put together, why it is built that way, a **Goals** -- Implement the standard models of an introductory quantitative finance course (option pricing, mean–variance and risk-parity portfolios, backtesting) in code that a reader can follow line by line. +- Implement the standard models of an introductory quantitative finance course (option pricing, mean–variance and risk-parity portfolios, backtesting) and the standard way risk teams validate a VaR model, in code that a reader can follow line by line. - Verify every model against something it did not produce itself: a textbook value, a closed form, an independent implementation, or an identity. - Fail loudly: validate inputs at the public boundary and raise typed errors with messages that name the problem. @@ -27,11 +27,12 @@ src/quantkit/ │ monte_carlo.py · asian.py ├── portfolio/ estimators.py · mean_variance.py · risk_parity.py · risk.py ├── backtest/ engine.py · strategies.py · metrics.py +├── var/ forecasts.py · backtests.py · traffic_light.py · es_backtest.py ├── data.py synthetic GBM prices, CSV loader └── cli.py, __main__.py argparse front end (`quantkit ...`, `python -m quantkit`) ``` -Dependencies point one way: `pricing`, `portfolio` and `backtest` depend only on `_validation` and `exceptions` (and `backtest.engine` on `backtest.metrics`); `data` depends on `_validation`; only `cli` depends on everything, and only `__main__` imports `cli`. This keeps each subpackage testable on its own. +Dependencies point one way: `pricing`, `portfolio`, `backtest` and `var` depend only on `_validation` and `exceptions` (and `backtest.engine` on `backtest.metrics`; inside `var`, `forecasts` uses `backtests` and `es_backtest` uses `forecasts`); `data` depends on `_validation`; only `cli` depends on everything, and only `__main__` imports `cli`. This keeps each subpackage testable on its own. The tail size shared by historical VaR and rolling historical simulation, `ceil(n (1 - level))`, lives in `_validation.tail_count` so that the two cannot drift apart. **Conventions used everywhere** @@ -182,30 +183,87 @@ for the last trade day `s ≤ t`, computed for all days with a forward fill. The Sample standard deviations use `ddof=1`. The annual risk-free rate is converted to a per-period rate that compounds to it. Drawdowns are measured from a running peak that starts at the initial capital, so a loss on the first day counts. Maximum drawdown duration is the longest run of periods below a previous peak, until recovery or the end of the data. Hit rate ignores zero-return periods (for example, days spent in cash). Ratios with a zero denominator, or one that is zero up to round-off, return NaN rather than infinity. -## 6. How correctness is verified +## 6. VaR model validation (`quantkit.var`) + +A VaR model is only as good as its backtest. Banks forecast tomorrow's 99% one-day VaR every evening, count the days on which the loss exceeded it (*exceptions*), and test whether those exceptions behave like a correct model's: about 1% of days, independent of each other. Supervisors map the count of the last 250 days to a capital penalty. The subpackage implements that loop end to end. + +### Rolling forecasts (`forecasts.py`) + +**Timing contract.** The forecast for day `t` is made at the close of day `t - 1` from the returns `r[t-h], ..., r[t-1]`, where `h` is the method's history. It never sees `r[t]`, the return it is scored against, and nothing older than `h` days. Forecasts are indexed by the day they forecast, next to the realized return of that day, so exceptions are an elementwise comparison with no shifting left to the caller. + +**One representation for every method.** Each method describes tomorrow's return as `r[t] = mu[t] + sigma[t] e[t]`: + +| method | `mu[t]` | `sigma[t]` | innovation `e[t]` | history `h` | +| --- | --- | --- | --- | --- | +| historical | 0 | 1 | resampled from the window's returns | window | +| gaussian | window mean | window std (`ddof=1`) | N(0, 1) | window | +| ewma | 0 | EWMA volatility | N(0, 1) | window | +| fhs | 0 | EWMA volatility | resampled standardized returns `r[s] / sigma[s]` | 2 × window | + +VaR and ES follow from the representation: `z sigma - mu` and `sigma phi(z) / (1 - c) - mu` for normal innovations, and `-(mu + sigma e_(m))` and `-(mu + sigma mean(e_(1..m)))` for an empirical pool, with `m = ceil(window (1 - c))`. For historical simulation that is exactly `historical_var` and `historical_cvar` of the window, and for the Gaussian method exactly `parametric_var` and `parametric_cvar`; the tests check every window against those functions. The same representation defines `RiskForecast.simulate`, which draws returns from each day's forecast distribution for the ES backtest. Keeping VaR, ES and simulation on one description means they cannot disagree about what the model says. + +**Vectorization.** `numpy.lib.stride_tricks.sliding_window_view` gives every window as a read-only view without copying. Historical VaR is `np.partition` along each row (the m-th smallest element and the m - 1 below it, in linear time per row); the Gaussian method uses row means and standard deviations; EWMA is one matrix-vector product of the squared windows with the weight vector. Each row depends only on its own window, so changing a later return cannot change an earlier forecast even by round-off, and the no-look-ahead test can demand bit-identical forecasts. + +**EWMA (RiskMetrics).** RiskMetrics forecasts variance as `sigma²[t] = λ sigma²[t-1] + (1 - λ) r²[t-1]` with `λ = 0.94` for daily data, which weights the return `k` days back by `(1 - λ) λ^k` and assumes a zero mean. quantkit applies the same weights inside the window and normalizes them to sum to one. *Trade-off:* the infinite recursion depends on how it was started and on every return since; the windowed form gives every method the same finite information set and memory, and with a 250-day window the omitted weight `0.94^250` is about 2e-7, so the two agree to that precision (the tests run the recursion over each window as an independent check). + +**Filtered historical simulation.** Hull & White (1998) rescale each historical return by the ratio of today's volatility to the volatility on its own day, `r[s] sigma[t] / sigma[s]`; Barone-Adesi, Giannopoulos & Vosper (1999) call the resampling of such standardized returns filtered historical simulation. It keeps the empirical shape of the tails (unlike EWMA) and reacts to volatility (unlike plain historical simulation). Each `sigma[s]` is itself an EWMA forecast that needs a window, so FHS needs two windows of history. *Alternative not taken:* a GARCH filter fitted by maximum likelihood, which would add an optimizer and its convergence failures to every forecast. + +**Exceptions** use a strict inequality, `r[t] < -VaR[t]`: a loss exactly equal to VaR is not an exception. + +### Exception tests (`backtests.py`) + +Under a correct model the exception indicators are i.i.d. Bernoulli(`p = 1 - c`). With `x` exceptions in `T` days and `n_ij` the number of days in state `i` followed by state `j`: + +- **Kupiec (1995), unconditional coverage:** `LR_uc = -2 ln[(1-p)^(T-x) p^x / ((1-π̂)^(T-x) π̂^x)]`, `π̂ = x/T`, asymptotically `χ²(1)`. The test is two-sided: too few exceptions mean the VaR is too high, which wastes capital. +- **Christoffersen (1998), independence:** a first-order Markov chain with `π01 = n01 / (n00 + n01)` and `π11 = n11 / (n10 + n11)` against one pooled `π`; `LR_ind ~ χ²(1)`. +- **Christoffersen (1998), conditional coverage:** the Markov chain against i.i.d. Bernoulli(`p`); `LR_cc ~ χ²(2)`. Both Christoffersen tests use the `T - 1` transitions, so `LR_cc` equals `LR_ind` plus Kupiec's statistic on days 2..T, an identity the tests check on random sequences. `kupiec_test` itself uses all `T` days, as Kupiec defines it. + +**Edge cases.** Every likelihood term is written as `count × ln(probability)` with `0 × ln 0 = 0`, and an empty set of transitions contributes an empty product. With no exceptions, Kupiec's statistic is `-2T ln(c)`, finite, and `LR_ind = 0` because there is nothing to compare. With only exceptions the same holds the other way round. With isolated exceptions (`n11 = 0`) the alternative fits `π11 = 0` exactly and `LR_ind` stays finite. A single observation is enough for Kupiec; the Christoffersen tests need two. + +**P-values.** `χ²(1)` is the square of a standard normal and `χ²(2)` is exponential with mean 2, so the p-values are `erfc(sqrt(x/2))` and `exp(-x/2)`; the tests compare them with SciPy's `chi2.sf`. The chi-squared limits are poor when few exceptions are expected (2.5 in 250 days at 99%). `null_rejection_rates` measures how poor: a dynamic program walks the days and keeps the joint distribution of (first day's state, last day's state, number of exceptions, number of consecutive exception pairs). Those four numbers determine all transition counts, so every statistic, and summing the probability of the states each test rejects gives its exact size. The work is `T × (x_max + 1)²` for `x_max` the largest exception count with non-negligible probability (`binom.isf(1e-15)`), a few milliseconds for 250 or 1,000 days. The tests check it against brute-force enumeration of all `2^T` sequences for small `T`, and against the binomial distribution for Kupiec. *Alternative not taken:* exact or Monte Carlo p-values for the tests themselves (Dufour-style Monte Carlo tests, as in Christoffersen & Pelletier 2004); the dynamic program would support exact p-values, which is on the roadmap. + +### Basel traffic light (`traffic_light.py`) + +Table 2 of the Basel Committee's 1996 backtesting framework maps the exceptions of 250 daily 99% VaR forecasts to zones: green for 0 to 4, yellow for 5 to 9 with plus factors 0.40, 0.50, 0.65, 0.75 and 0.85 added to the capital multiplier of 3, and red (plus factor 1) from 10. The boundaries come from the binomial(250, 1%) distribution: green while the probability of that many exceptions or fewer is below 95%, red from 99.99%. The tests reproduce the published cumulative probabilities from `scipy.stats.binom` and re-derive the zones from that rule. The function deliberately takes only the count for 250 days at 99%, the case the table defines. + +### Expected shortfall backtest (`es_backtest.py`) + +VaR says nothing about how large the losses beyond it are; expected shortfall does, and the Basel FRTB standards moved market risk capital to 97.5% ES. ES is harder to backtest because it is not *elicitable* (no scoring function is minimized by it, Gneiting 2011), so there is no simple exception count. Acerbi & Székely (2014) propose + +```text +Z2 = Σ_t r[t] I[t] / (T α ES[t]) + 1, I[t] = 1{r[t] < -VaR[t]}, α = 1 - c +``` + +If each day's forecast distribution is right and continuous, `E[r I] = -α ES` every day, so `E[Z2] = 0`; underestimated ES pushes `Z2` below zero. Of their three statistics, Z2 needs only the VaR and ES forecasts and the realized returns. + +`Z2` has no standard distribution, so `es_backtest` simulates it under the null, as the paper prescribes: it draws `N` return paths from the forecast distributions with `RiskForecast.simulate`, computes `Z2` on each and reports `(1 + #{Z_sim ≤ Z_obs}) / (N + 1)`. Counting the observed path among the draws makes this a Monte Carlo test with exact size whenever `α_test (N + 1)` is a whole number (Davison & Hinkley 1997), whatever the sample length; ties, such as `Z2 = 1` on paths without exceptions, count against the model. Paths are simulated in batches of at most two million path-days to bound memory, consuming the random stream exactly as one large draw would (a test checks this). With zero-mean normal forecasts `r / sigma` is standard normal and `Z2` depends only on it, so the null distribution does not depend on the volatility path; with the same seed, two different volatility paths give identical critical values, which another test checks. + +*Trade-off:* simulating the null from each method's own forecasts tests the whole tail of the forecast distribution, not ES alone, and for historical simulation and FHS the forecast distribution is discrete (resampled returns), so `E[Z2]` is not exactly zero there; the simulated critical values account for that. The README's size and power study measures how the tests behave. + +## 7. How correctness is verified | Kind of check | Examples | | --- | --- | -| Textbook values | Hull's BSM call and put; Hull's Greeks example; Hull's American put on 5–500-step trees | -| Closed forms | Unbounded mean–variance portfolios; inverse-volatility risk parity for uncorrelated or two assets; geometric Asian with one fixing = BSM; Gaussian VaR | -| Independent implementations | scikit-learn's Ledoit–Wolf; a day-by-day reference backtest; a separate simulation of the geometric average; finite-difference Greeks | +| Textbook values | Hull's BSM call and put; Hull's Greeks example; Hull's American put on 5–500-step trees; the Basel traffic-light table | +| Closed forms | Unbounded mean–variance portfolios; inverse-volatility risk parity for uncorrelated or two assets; geometric Asian with one fixing = BSM; Gaussian VaR; Kupiec and Christoffersen statistics by hand; exception probabilities of rolling VaR on i.i.d. data | +| Independent implementations | scikit-learn's Ledoit–Wolf; a day-by-day reference backtest; a separate simulation of the geometric average; finite-difference Greeks; single-window VaR for every rolling window; the RiskMetrics recursion; a Hull–White rescaling loop; SciPy's binomial and chi-squared distributions; brute-force enumeration of exception sequences | | Identities and properties (Hypothesis) | Put–call parity; no-arbitrage bounds; IV round trips; American ≥ European; CVaR ≥ VaR | -| Statistical checks with fixed seeds | Monte Carlo within 3 standard errors; confidence-interval coverage; simulated moments and correlations | +| Statistical checks with fixed seeds | Monte Carlo within 3 standard errors; confidence-interval coverage; simulated moments and correlations; backtest sizes on returns drawn from the forecast model itself, against their exact values; power against fat tails; `E[Z2] = 0` under the null | | Optimality | No random feasible portfolio beats the minimum-variance or maximum-Sharpe result; first-order conditions of the closed forms | -| Negative controls | A peeking strategy must fail the look-ahead checks; prices violating no-arbitrage must raise | +| Negative controls | A peeking strategy or VaR forecaster must fail the look-ahead checks; prices violating no-arbitrage must raise | | Documentation | The README's Python examples run as doctests | Monte Carlo and simulation tests use fixed seeds, so each either always passes or always fails; they cannot fail by chance in CI. Hypothesis property tests draw new random inputs on every run (200 per test, or 5,000 with `pytest --hypothesis-profile=thorough`). Tolerances are tied to the method's error: for example, the binomial test allows an error below `0.05 · S / steps` because CRR converges at rate 1/steps, and the gamma finite difference allows a relative error of 1e-4 because a second difference loses roughly half the digits to round-off. -## 7. Tooling +## 8. Tooling - **uv** manages the environment; `uv.lock` pins every dependency for every supported Python. The minimum versions in `pyproject.toml` are the oldest versions the lockfile actually exercises (the Python 3.11 job resolves NumPy 2.4 and SciPy 1.17). - **mypy --strict** covers `src`, `tests` and `examples`. SciPy and scikit-learn ship no type information, so a single override ignores their missing stubs, and their results are wrapped in `float()` or `np.asarray()`. Tests that deliberately pass wrong types use `typing.cast` instead of `type: ignore`. - **ruff** handles linting (including pydocstyle with Google-style docstrings) and formatting. - **pytest** runs unit tests, Hypothesis properties and the README doctests, with warnings turned into errors so deprecations surface immediately. Coverage must stay at or above 90% in CI. -- **CI** (GitHub Actions) runs lint, format check, type check, tests with coverage, and CLI smoke tests on Python 3.11, 3.12 and 3.13, with read-only permissions and cancellation of superseded runs. +- **CI** (GitHub Actions) runs lint, format check, type check, tests with coverage, CLI smoke tests and the example scripts on Python 3.11, 3.12 and 3.13, with read-only permissions and cancellation of superseded runs. -## 8. Questions worth being able to answer +## 9. Questions worth being able to answer - *Why does the binomial error oscillate instead of falling smoothly?* The strike falls at different positions between terminal nodes as the step count changes; the payoff kink is not resolved consistently. - *Why Brent rather than Newton for implied volatility?* Guaranteed convergence on a bracket; Newton divides by vega, which can be nearly zero. @@ -214,6 +272,10 @@ Monte Carlo and simulation tests use fixed seeds, so each either always passes o - *Why does the maximum-Sharpe portfolio often concentrate in one or two assets?* It leans on the noisiest input, the expected returns; small estimation errors in `μ` move the weights a lot (see the README examples). - *How do you know the backtest has no look-ahead?* The contract is explicit, the engine shifts weights by one day, the schedule uses first trading days, and a perturbation test shows that future prices cannot change past weights, with a cheating strategy as a negative control. - *Why is turnover compared with drifted weights?* So that buy-and-hold is not charged for price moves. -- *What would you change for production use?* Execution lag and slippage, per-asset constraints, a QP solver, vectorized pricing, and robust estimates of expected returns. +- *Why does a 99% VaR model with no exceptions in 250 days fail Kupiec's test?* The test is two-sided, and zero exceptions is unlikely (0.99^250 = 8.1%) for a correct model: its VaR is probably too high. This is also why the nominal 5% test rejects a correct model 9.5% of the time at that sample size. +- *Why does 250-day historical simulation have more than 1% exceptions at 99%?* Its VaR is the 3rd-worst return of 250, and on i.i.d. data a new return falls below the 3rd-smallest of 250 others with probability 3/251, about 1.2%. +- *How do you know the VaR forecasts have no look-ahead?* Forecasts are a function of the window before the forecast day only; a perturbation test changes the scored return and every later one and demands bit-identical forecasts, and a forecaster that peeks fails it. +- *Why simulate Z2's critical values from the forecasts themselves?* Under the null the realized returns are one draw from the forecast distributions, so ranking the observed statistic among simulated ones gives a test of exact size without asymptotics. +- *What would you change for production use?* Execution lag and slippage, per-asset constraints, a QP solver, vectorized pricing, robust estimates of expected returns, and VaR from position-level P&L with full revaluation. See the README's *Limitations and roadmap* section for what is out of scope. diff --git a/docs/images/var_backtest.png b/docs/images/var_backtest.png new file mode 100644 index 0000000..3b691d1 Binary files /dev/null and b/docs/images/var_backtest.png differ diff --git a/examples/var_backtest.py b/examples/var_backtest.py new file mode 100644 index 0000000..af3e4f3 --- /dev/null +++ b/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() diff --git a/examples/var_size_power.py b/examples/var_size_power.py new file mode 100644 index 0000000..49e83d3 --- /dev/null +++ b/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() diff --git a/pyproject.toml b/pyproject.toml index c2a5f1c..7f43185 100644 --- a/pyproject.toml +++ b/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", diff --git a/src/quantkit/__init__.py b/src/quantkit/__init__.py index 71525a0..a502904 100644 --- a/src/quantkit/__init__.py +++ b/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. diff --git a/src/quantkit/_validation.py b/src/quantkit/_validation.py index f992a0f..352d667 100644 --- a/src/quantkit/_validation.py +++ b/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.""" diff --git a/src/quantkit/cli.py b/src/quantkit/cli.py index 17dd99e..aecc0bb 100644 --- a/src/quantkit/cli.py +++ b/src/quantkit/cli.py @@ -1,4 +1,4 @@ -"""Command-line interface: ``quantkit price | iv | simulate | backtest | frontier``. +"""Command-line interface: ``quantkit price | iv | simulate | backtest | frontier | var-backtest``. Run ``quantkit --help`` for the options of each command. """ @@ -42,6 +42,18 @@ implied_vol, monte_carlo_price, ) +from quantkit.var import ( + BASEL_DAYS, + BASEL_LEVEL, + METHODS, + RISKMETRICS_DECAY, + RiskForecast, + backtest_var, + basel_traffic_light, + es_backtest, + null_rejection_rates, + rolling_var_es, +) SIMULATED_TICKERS = ("AAA", "BBB", "CCC", "DDD", "EEE") SIMULATED_DRIFTS = (0.03, 0.05, 0.07, 0.09, 0.11) @@ -50,6 +62,13 @@ STRATEGY_NAMES = ("buy-and-hold", "equal-weight", "ma-crossover", "momentum", "inverse-vol") +VAR_METHOD_LABELS = { + "historical": "historical", + "gaussian": "Gaussian", + "ewma": "EWMA", + "fhs": "FHS", +} + def main(argv: Sequence[str] | None = None) -> int: """Run the command line interface and return the process exit code.""" @@ -282,6 +301,85 @@ def _save_frontier_plot(table: pd.DataFrame, mu: pd.Series, cov: pd.DataFrame, p plt.close(fig) +def _cmd_var_backtest(args: argparse.Namespace) -> None: + prices = load_prices_csv(args.prices) + asset_returns = compute_returns(prices) + if args.ticker is None: + returns = asset_returns.mean(axis=1) + subject = f"equal-weight portfolio of {prices.shape[1]} assets, rebalanced daily" + elif args.ticker in asset_returns.columns: + returns = asset_returns[args.ticker] + subject = args.ticker + else: + raise QuantkitError( + f"ticker {args.ticker!r} is not in {args.prices}; " + f"columns: {', '.join(map(str, asset_returns.columns))}" + ) + methods = list(dict.fromkeys(args.method or METHODS)) + forecasts = { + name: rolling_var_es( + returns, method=name, level=args.level, window=args.window, decay=args.decay + ) + for name in methods + } + # Score every method on the same days: the period that the longest history covers. + n_days = min(len(forecast) for forecast in forecasts.values()) + forecasts = {name: forecast.last(n_days) for name, forecast in forecasts.items()} + days = next(iter(forecasts.values())).var.index + decay = f", EWMA decay {args.decay:g}" if {"ewma", "fhs"} & set(methods) else "" + level = f"{100 * args.level:g}%" + print( + f"{subject}: {len(returns)} daily returns, {_day(returns.index[0])} to " + f"{_day(returns.index[-1])}\n" + f"{level} one-day VaR and ES, {args.window}-day window{decay}\n" + f"backtest: {n_days} days, {_day(days[0])} to {_day(days[-1])}" + f"{', the period every method covers' if len(methods) > 1 else ''}\n" + ) + table = pd.DataFrame( + { + VAR_METHOD_LABELS[name]: _var_backtest_column(forecast, args) + for name, forecast in forecasts.items() + } + ) + print(table.to_string()) + sizes = null_rejection_rates(n_days, level=args.level) + print( + f"\nLR p-values use chi-squared limits. Over {n_days} days at {level}, a correct model " + f"fails the 5% tests\nwith probability {sizes.unconditional:.1%} (LR_uc), " + f"{sizes.independence:.1%} (LR_ind) and {sizes.conditional:.1%} (LR_cc).\n" + f"Z2 p-values compare with {args.es_simulations:,} paths simulated from each method's " + f"own forecasts (seed {args.seed})." + ) + + +def _var_backtest_column(forecast: RiskForecast, args: argparse.Namespace) -> dict[str, str]: + result = backtest_var(forecast.exceptions, level=forecast.level) + if forecast.level != BASEL_LEVEL: + zone = "n/a (99% only)" + elif len(forecast) < BASEL_DAYS: + zone = f"n/a (< {BASEL_DAYS} days)" + else: + light = basel_traffic_light(int(forecast.last(BASEL_DAYS).exceptions.sum())) + zone = f"{light.zone} ({light.n_exceptions})" + es = es_backtest(forecast, n_simulations=args.es_simulations, seed=args.seed) + return { + "exceptions": f"{result.n_exceptions}", + "expected exceptions": f"{result.expected_exceptions:.1f}", + "exception rate": f"{result.exception_rate:.2%}", + "Kupiec LR_uc p-value": f"{result.unconditional.p_value:.3f}", + "Christoffersen LR_ind p-value": f"{result.independence.p_value:.3f}", + "conditional coverage LR_cc p-value": f"{result.conditional.p_value:.3f}", + f"Basel zone, last {BASEL_DAYS} days": zone, + "Acerbi-Szekely Z2": f"{es.statistic:.2f}", + "Z2 simulated p-value": f"{es.p_value:.3f}", + } + + +def _day(label: object) -> str: + """A date label for messages (the index of a CSV price table is a DatetimeIndex).""" + return str(pd.Timestamp(str(label)).date()) + + # ------------------------------------------------------------------------ parser @@ -356,4 +454,36 @@ def _build_parser() -> argparse.ArgumentParser: ) frontier.add_argument("--plot", type=Path, help="also save a PNG chart (needs matplotlib)") frontier.set_defaults(handler=_cmd_frontier) + + var_backtest = commands.add_parser( + "var-backtest", + help="backtest rolling VaR and ES forecasts on a CSV of daily closing prices", + description="Forecast one-day VaR and ES from a trailing window, then test the " + "forecasts: exception counts (Kupiec, Christoffersen), the Basel traffic light " + "and the Acerbi-Szekely Z2 test of expected shortfall.", + ) + var_backtest.add_argument( + "prices", type=Path, help="CSV in long (date,ticker,close) or wide form" + ) + var_backtest.add_argument( + "--method", + choices=METHODS, + action="append", + help="forecasting method; repeat to compare several (default: all four)", + ) + var_backtest.add_argument( + "--level", type=float, default=0.99, help="VaR and ES confidence level" + ) + var_backtest.add_argument("--window", type=int, default=250, help="estimation window (days)") + var_backtest.add_argument( + "--decay", type=float, default=RISKMETRICS_DECAY, help="EWMA decay for ewma and fhs" + ) + var_backtest.add_argument( + "--ticker", help="backtest one column instead of the equal-weight portfolio" + ) + var_backtest.add_argument( + "--es-simulations", type=int, default=9_999, help="null paths for the Z2 p-value" + ) + var_backtest.add_argument("--seed", type=int, default=0, help="seed for the Z2 simulation") + var_backtest.set_defaults(handler=_cmd_var_backtest) return parser diff --git a/src/quantkit/portfolio/risk.py b/src/quantkit/portfolio/risk.py index fa8d3ce..03472f5 100644 --- a/src/quantkit/portfolio/risk.py +++ b/src/quantkit/portfolio/risk.py @@ -15,15 +15,17 @@ import numpy.typing as npt import pandas as pd -from quantkit._validation import FloatArray, as_1d_array, require_level, validate_covariance +from quantkit._validation import ( + FloatArray, + as_1d_array, + require_level, + tail_count, + validate_covariance, +) from quantkit.exceptions import InvalidParameterError _STANDARD_NORMAL = NormalDist() -# (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 portfolio_volatility(weights: pd.Series, cov: pd.DataFrame) -> float: """Portfolio volatility ``sqrt(w' Sigma w)``, in the units of ``cov``.""" @@ -37,9 +39,7 @@ def portfolio_volatility(weights: pd.Series, cov: pd.DataFrame) -> float: def _tail(returns: npt.ArrayLike, level: float) -> FloatArray: """The ``ceil(n * (1 - level))`` worst returns, worst first.""" r = np.sort(as_1d_array("returns", returns, min_length=2)) - level = require_level("level", level) - n_tail = max(1, math.ceil(len(r) * (1.0 - level) - _FLOAT_SLACK)) - return r[:n_tail] + return r[: tail_count(len(r), require_level("level", level))] def historical_var(returns: npt.ArrayLike, *, level: float = 0.95) -> float: diff --git a/src/quantkit/var/__init__.py b/src/quantkit/var/__init__.py new file mode 100644 index 0000000..4d00c45 --- /dev/null +++ b/src/quantkit/var/__init__.py @@ -0,0 +1,74 @@ +"""VaR model validation: rolling VaR and ES forecasts and the backtests that judge them. + +* :func:`rolling_var_es` produces one-day VaR and expected shortfall forecasts by + historical simulation, a Gaussian model, EWMA (RiskMetrics) or filtered + historical simulation, each from a trailing window that never includes the + return being forecast. +* :func:`backtest_var` counts the exceptions (days whose loss exceeded VaR) and + runs Kupiec's and Christoffersen's likelihood-ratio tests; + :func:`null_rejection_rates` gives their exact size at a given sample length. +* :func:`basel_traffic_light` maps 250 days of 99% exceptions to the Basel + green, yellow and red zones and their capital plus factors. +* :func:`es_backtest` tests expected shortfall with Acerbi and Székely's ``Z2`` + statistic against its simulated null distribution. + +Returns, VaR and ES are decimals; VaR and ES are positive loss fractions. +""" + +from quantkit.var.backtests import ( + LikelihoodRatioTest, + RejectionRates, + TransitionCounts, + VaRBacktest, + backtest_var, + christoffersen_test, + conditional_coverage_test, + kupiec_test, + null_rejection_rates, + transition_counts, + var_exceptions, +) +from quantkit.var.es_backtest import ESBacktest, acerbi_szekely_z2, es_backtest +from quantkit.var.forecasts import ( + METHODS, + RISKMETRICS_DECAY, + Method, + RiskForecast, + ewma_weights, + rolling_var_es, +) +from quantkit.var.traffic_light import ( + BASEL_DAYS, + BASEL_LEVEL, + TrafficLight, + Zone, + basel_traffic_light, +) + +__all__ = [ + "BASEL_DAYS", + "BASEL_LEVEL", + "METHODS", + "RISKMETRICS_DECAY", + "ESBacktest", + "LikelihoodRatioTest", + "Method", + "RejectionRates", + "RiskForecast", + "TrafficLight", + "TransitionCounts", + "VaRBacktest", + "Zone", + "acerbi_szekely_z2", + "backtest_var", + "basel_traffic_light", + "christoffersen_test", + "conditional_coverage_test", + "es_backtest", + "ewma_weights", + "kupiec_test", + "null_rejection_rates", + "rolling_var_es", + "transition_counts", + "var_exceptions", +] diff --git a/src/quantkit/var/backtests.py b/src/quantkit/var/backtests.py new file mode 100644 index 0000000..b361574 --- /dev/null +++ b/src/quantkit/var/backtests.py @@ -0,0 +1,395 @@ +"""Exception backtests of VaR forecasts: Kupiec, Christoffersen and conditional coverage. + +An *exception* is a day whose loss exceeded that day's VaR forecast. If the +forecasts are correct at confidence level ``c``, exceptions are independent +Bernoulli trials with probability ``p = 1 - c``. Three likelihood-ratio tests +check this (Christoffersen 1998): + +* **Unconditional coverage** (Kupiec 1995, proportion of failures): is the + exception rate consistent with ``p``? With ``x`` exceptions in ``T`` days and + ``pi = x / T``:: + + LR_uc = -2 log[(1 - p)^(T - x) p^x / ((1 - pi)^(T - x) pi^x)] ~ chi2(1) + +* **Independence** (Christoffersen 1998): does an exception today change the + chance of one tomorrow? With ``n_ij`` the number of days in state ``i`` + followed by a day in state ``j`` (1 = exception), the alternative is a + first-order Markov chain with ``pi_01 = n01 / (n00 + n01)`` and + ``pi_11 = n11 / (n10 + n11)``; the null is one common ``pi``:: + + LR_ind = -2 log[L(pi) / L(pi_01, pi_11)] ~ chi2(1) + +* **Conditional coverage** (Christoffersen 1998): the Markov chain against + independent trials with probability exactly ``p``:: + + LR_cc = -2 log[L(p) / L(pi_01, pi_11)] ~ chi2(2) + + ``LR_ind`` and ``LR_cc`` use the ``T - 1`` transitions, so ``LR_cc`` equals + ``LR_ind`` plus the Kupiec statistic of days ``2..T``. + +Likelihood terms use ``0 log 0 = 0``, so every statistic stays finite with no +exceptions, with only exceptions, or with no pair of consecutive exceptions. + +P-values come from the chi-squared limits, which are rough when few exceptions +are expected: at 99% over 250 days a correct model expects only 2.5. +:func:`null_rejection_rates` computes how often each test really rejects a +correct model at a given sample size, by enumerating the exact distribution of +the transition counts. +""" + +from __future__ import annotations + +import math +from dataclasses import dataclass +from statistics import NormalDist + +import numpy as np +import numpy.typing as npt +import pandas as pd +from scipy.stats import binom + +from quantkit._validation import FloatArray, as_1d_array, require_int, require_level +from quantkit.exceptions import InvalidParameterError + +Counts = int | npt.NDArray[np.int64] +"""A count, or an array of counts for many samples at once.""" + +# null_rejection_rates ignores exception counts whose total probability is below this. +_NEGLIGIBLE_MASS = 1e-15 +# Largest dynamic program (days x states) null_rejection_rates will run, about a second. +_MAX_ENUMERATION_WORK = 500_000_000 + + +@dataclass(frozen=True, slots=True) +class LikelihoodRatioTest: + """A likelihood-ratio statistic and its asymptotic p-value. + + Attributes: + statistic: ``-2 log(L_null / L_alternative)``; never negative. + dof: Degrees of freedom of the limiting chi-squared distribution. + p_value: ``P(chi2(dof) >= statistic)``. + """ + + statistic: float + dof: int + p_value: float + + +@dataclass(frozen=True, slots=True) +class TransitionCounts: + """Counts of consecutive day pairs; ``n01`` counts a quiet day followed by an exception. + + Attributes: + n00: No exception, then no exception. + n01: No exception, then an exception. + n10: An exception, then no exception. + n11: An exception, then another exception. + """ + + n00: int + n01: int + n10: int + n11: int + + +@dataclass(frozen=True, slots=True) +class VaRBacktest: + """Exception counts and the three coverage tests for one series of VaR forecasts. + + Attributes: + level: VaR confidence level of the forecasts. + n_observations: Number of forecast days ``T``. + n_exceptions: Number of exceptions ``x``. + transitions: Consecutive-day transition counts. + unconditional: Kupiec's proportion-of-failures test (LR_uc). + independence: Christoffersen's independence test (LR_ind). + conditional: Christoffersen's conditional coverage test (LR_cc). + """ + + level: float + n_observations: int + n_exceptions: int + transitions: TransitionCounts + unconditional: LikelihoodRatioTest + independence: LikelihoodRatioTest + conditional: LikelihoodRatioTest + + @property + def expected_exceptions(self) -> float: + """Exceptions a correct model expects: ``T * (1 - level)``.""" + return self.n_observations * (1.0 - self.level) + + @property + def exception_rate(self) -> float: + """Observed exception rate ``x / T``.""" + return self.n_exceptions / self.n_observations + + +@dataclass(frozen=True, slots=True) +class RejectionRates: + """How often each coverage test rejects a correct VaR model. + + Attributes: + unconditional: Rejection probability of Kupiec's LR_uc. + independence: Rejection probability of Christoffersen's LR_ind. + conditional: Rejection probability of the conditional coverage LR_cc. + """ + + unconditional: float + independence: float + conditional: float + + +def var_exceptions( + returns: pd.Series | npt.ArrayLike, var: pd.Series | npt.ArrayLike +) -> pd.Series: + """Flag the days on which the loss exceeded the VaR forecast: ``return < -var``. + + Args: + returns: Realized returns. + var: VaR forecasts for the same days, as positive loss fractions. Two + Series must share their index; other inputs must have equal lengths. + + Returns: + A boolean Series (with the index of ``returns`` when it is a Series). + """ + r = as_1d_array("returns", returns) + v = as_1d_array("var", var) + if r.shape != v.shape: + raise InvalidParameterError( + f"returns and var must have the same length, got {r.size} and {v.size}" + ) + if ( + isinstance(returns, pd.Series) + and isinstance(var, pd.Series) + and not returns.index.equals(var.index) + ): + raise InvalidParameterError("returns and var must have the same index") + index = returns.index if isinstance(returns, pd.Series) else None + return pd.Series(r < -v, index=index, name="exception") + + +def kupiec_test(exceptions: pd.Series | npt.ArrayLike, *, level: float) -> LikelihoodRatioTest: + """Kupiec's proportion-of-failures test of unconditional coverage (LR_uc). + + The test is two-sided: too few exceptions also count against the model, + because they mean the VaR is set too high. With no exceptions the + statistic is ``-2 T log(level)``, so a 99% model with zero exceptions in 250 + days is rejected at 5%. + + Args: + exceptions: Exception indicators (booleans or 0/1), one per day. + level: VaR confidence level of the forecasts. + """ + hits = _hits(exceptions, min_length=1) + p = 1.0 - require_level("level", level) + n_one = int(hits.sum()) + statistic = float(_uc_statistic(hits.size - n_one, n_one, p)) + return LikelihoodRatioTest(statistic, 1, _chi2_sf(statistic, 1)) + + +def transition_counts(exceptions: pd.Series | npt.ArrayLike) -> TransitionCounts: + """Count the four kinds of consecutive day pairs in an exception series.""" + hits = _hits(exceptions, min_length=2) + before, after = hits[:-1], hits[1:] + n11 = int((before & after).sum()) + n10 = int((before & ~after).sum()) + n01 = int((~before & after).sum()) + return TransitionCounts(n00=before.size - n11 - n10 - n01, n01=n01, n10=n10, n11=n11) + + +def christoffersen_test(exceptions: pd.Series | npt.ArrayLike) -> LikelihoodRatioTest: + """Christoffersen's test that exceptions are independent from day to day (LR_ind). + + With no exceptions, or no quiet days, there is nothing to compare: the + statistic is 0 and the p-value 1. + + Args: + exceptions: Exception indicators (booleans or 0/1), at least two days. + """ + c = transition_counts(exceptions) + statistic = float(_ind_statistic(c.n00, c.n01, c.n10, c.n11)) + return LikelihoodRatioTest(statistic, 1, _chi2_sf(statistic, 1)) + + +def conditional_coverage_test( + exceptions: pd.Series | npt.ArrayLike, *, level: float +) -> LikelihoodRatioTest: + """Christoffersen's joint test of coverage and independence (LR_cc). + + Args: + exceptions: Exception indicators (booleans or 0/1), at least two days. + level: VaR confidence level of the forecasts. + """ + c = transition_counts(exceptions) + p = 1.0 - require_level("level", level) + statistic = float(_cc_statistic(c.n00, c.n01, c.n10, c.n11, p)) + return LikelihoodRatioTest(statistic, 2, _chi2_sf(statistic, 2)) + + +def backtest_var(exceptions: pd.Series | npt.ArrayLike, *, level: float) -> VaRBacktest: + """Count exceptions and run all three coverage tests. + + Args: + exceptions: Exception indicators (booleans or 0/1), at least two days; + for example :attr:`RiskForecast.exceptions` or :func:`var_exceptions`. + level: VaR confidence level of the forecasts. + """ + hits = _hits(exceptions, min_length=2) + level = require_level("level", level) + return VaRBacktest( + level=level, + n_observations=hits.size, + n_exceptions=int(hits.sum()), + transitions=transition_counts(hits), + unconditional=kupiec_test(hits, level=level), + independence=christoffersen_test(hits), + conditional=conditional_coverage_test(hits, level=level), + ) + + +def null_rejection_rates( + n_observations: int, *, level: float, significance: float = 0.05 +) -> RejectionRates: + """Exact probability that each test rejects a correct model, at its asymptotic cutoff. + + Under a correct model the exceptions are independent Bernoulli(p) trials, so + the three statistics have exact distributions that the chi-squared limits + only approximate. A dynamic program over the days tracks the joint + distribution of the first and last day's state, the number of exceptions and + the number of consecutive exception pairs; those four numbers determine the + transition counts and so every statistic. Exception counts with a total + probability below 1e-15 are left out. + + For example, at 99% over 250 days the 5% Kupiec test rejects a correct model + far more often than 5%: no exceptions at all, which happens with probability + 0.99^250 = 8.1%, already gives ``LR_uc = 5.0 > 3.84``. + + Args: + n_observations: Number of backtest days ``T`` (at least 2). + level: VaR confidence level, e.g. 0.99. + significance: Nominal level of the tests, e.g. 0.05. + + Returns: + The probability that each test's asymptotic p-value falls below + ``significance`` when the model is correct. + + Raises: + InvalidParameterError: For invalid inputs, or a sample so long (with so + many expected exceptions) that the enumeration would be too slow. + """ + n = require_int("n_observations", n_observations, minimum=2) + p = 1.0 - require_level("level", level) + significance = require_level("significance", significance) + max_ones = min(n, int(binom.isf(_NEGLIGIBLE_MASS, n, p)) + 1) + if n * (max_ones + 1) ** 2 > _MAX_ENUMERATION_WORK: + raise InvalidParameterError( + f"enumerating {n} days with up to {max_ones} exceptions is too slow; " + "simulate i.i.d. Bernoulli exceptions instead" + ) + + # prob[first, last, ones, pairs]: distribution after the days processed so far. + prob = np.zeros((2, 2, max_ones + 1, max_ones + 1)) + prob[0, 0, 0, 0] = 1.0 - p + prob[1, 1, 1, 0] = p + for _ in range(n - 1): + step = np.zeros_like(prob) + step[:, 0] = (prob[:, 0] + prob[:, 1]) * (1.0 - p) # a quiet day + step[:, 1, 1:, :] = prob[:, 0, :-1, :] * p # an exception after a quiet day + step[:, 1, 1:, 1:] += prob[:, 1, :-1, :-1] * p # an exception after an exception + prob = step + + first, last, ones, n11 = (axis.ravel() for axis in np.indices(prob.shape)) + weight = prob.ravel() + # Exceptions on days 1..T-1 and 2..T give the transition counts. + n10 = ones - last - n11 + n01 = ones - first - n11 + n00 = (n - 1) - n01 - n10 - n11 + possible = (weight > 0) & (np.minimum(np.minimum(n00, n01), n10) >= 0) + first, ones, n00, n01, n10, n11, weight = ( + a[possible] for a in (first, ones, n00, n01, n10, n11, weight) + ) + cutoff_1 = NormalDist().inv_cdf(1.0 - significance / 2.0) ** 2 # chi2(1) quantile + cutoff_2 = -2.0 * math.log(significance) # chi2(2) quantile + return RejectionRates( + unconditional=float(weight[_uc_statistic(n - ones, ones, p) > cutoff_1].sum()), + independence=float(weight[_ind_statistic(n00, n01, n10, n11) > cutoff_1].sum()), + conditional=float(weight[_cc_statistic(n00, n01, n10, n11, p) > cutoff_2].sum()), + ) + + +def _hits(exceptions: pd.Series | npt.ArrayLike, *, min_length: int) -> npt.NDArray[np.bool_]: + """Exception indicators as a boolean array, accepting booleans or 0/1 numbers.""" + array = np.asarray(exceptions) + if array.dtype != np.bool_: + values = as_1d_array("exceptions", array) + if not np.isin(values, (0.0, 1.0)).all(): + raise InvalidParameterError("exceptions must be booleans or 0/1 values") + array = values == 1.0 + if array.ndim != 1: + raise InvalidParameterError(f"exceptions must be one-dimensional, got shape {array.shape}") + if array.size < min_length: + raise InvalidParameterError( + f"exceptions need at least {min_length} observations, got {array.size}" + ) + return array + + +def _xlogy(count: Counts, probability: float | FloatArray) -> FloatArray: + """``count * log(probability)``, with ``0 * log(0) = 0``.""" + count_array = np.asarray(count, dtype=np.float64) + with np.errstate(divide="ignore", invalid="ignore"): + terms = count_array * np.log(probability) + result: FloatArray = np.where(count_array > 0, terms, 0.0) + return result + + +def _loglik(n_zero: Counts, n_one: Counts, p: float) -> FloatArray: + """Log-likelihood of ``n_zero`` zeros and ``n_one`` ones from Bernoulli(p) trials.""" + return _xlogy(n_zero, 1.0 - p) + _xlogy(n_one, p) + + +def _max_loglik(n_zero: Counts, n_one: Counts) -> FloatArray: + """The Bernoulli log-likelihood at its maximum, ``p = n_one / (n_zero + n_one)``. + + With no trials at all the likelihood is an empty product: log-likelihood 0. + """ + total = np.asarray(n_zero + n_one, dtype=np.float64) + trials = np.where(total > 0, total, 1.0) + return _xlogy(n_zero, n_zero / trials) + _xlogy(n_one, n_one / trials) + + +def _likelihood_ratio(null_loglik: FloatArray, alternative_loglik: FloatArray) -> FloatArray: + # The alternative nests the null, so the ratio is >= 0 up to round-off. + ratio: FloatArray = np.maximum(0.0, 2.0 * (alternative_loglik - null_loglik)) + return ratio + + +def _uc_statistic(n_zero: Counts, n_one: Counts, p: float) -> FloatArray: + return _likelihood_ratio(_loglik(n_zero, n_one, p), _max_loglik(n_zero, n_one)) + + +def _markov_loglik(n00: Counts, n01: Counts, n10: Counts, n11: Counts) -> FloatArray: + """Maximized log-likelihood of a first-order Markov chain with these transitions.""" + return _max_loglik(n00, n01) + _max_loglik(n10, n11) + + +def _ind_statistic(n00: Counts, n01: Counts, n10: Counts, n11: Counts) -> FloatArray: + iid = _max_loglik(n00 + n10, n01 + n11) + return _likelihood_ratio(iid, _markov_loglik(n00, n01, n10, n11)) + + +def _cc_statistic(n00: Counts, n01: Counts, n10: Counts, n11: Counts, p: float) -> FloatArray: + null = _loglik(n00 + n10, n01 + n11, p) + return _likelihood_ratio(null, _markov_loglik(n00, n01, n10, n11)) + + +def _chi2_sf(statistic: float, dof: int) -> float: + """Upper tail of the chi-squared distribution with 1 or 2 degrees of freedom. + + chi2(1) is the square of a standard normal, so ``P(X > x) = erfc(sqrt(x / 2))``; + chi2(2) is exponential with mean 2, so ``P(X > x) = exp(-x / 2)``. + """ + if dof == 1: + return math.erfc(math.sqrt(statistic / 2.0)) + return math.exp(-statistic / 2.0) diff --git a/src/quantkit/var/es_backtest.py b/src/quantkit/var/es_backtest.py new file mode 100644 index 0000000..7420301 --- /dev/null +++ b/src/quantkit/var/es_backtest.py @@ -0,0 +1,144 @@ +"""Acerbi and Székely's Z2 backtest of expected shortfall forecasts. + +Expected shortfall (ES) is the mean loss on the days beyond VaR. Acerbi and +Székely (2014) test ES forecasts directly with their second statistic:: + + Z2 = sum_t r[t] I[t] / (T alpha ES[t]) + 1, I[t] = 1 if r[t] < -VaR[t] else 0 + +where ``alpha = 1 - level`` is the tail probability and ``T`` the number of +days. If each day's forecast distribution is right and continuous, then +``E[r I] = -alpha ES`` on every day, so ``E[Z2] = 0``. When ES is +underestimated, the losses beyond VaR are larger or more frequent than +forecast, and ``Z2`` is negative. The test is one-sided: small values reject. + +``Z2`` has no standard distribution. As the paper prescribes, its null +distribution is simulated: draw many return paths from the forecast +distributions themselves (:meth:`~quantkit.var.RiskForecast.simulate`), compute +``Z2`` on each, and locate the observed value among them. The p-value counts the +observed path as one of the draws, ``(1 + #{Z_sim <= Z_obs}) / (N + 1)``, which +makes it a Monte Carlo test with the correct size at any sample length +(Davison & Hinkley 1997, chapter 4). Ties, such as ``Z2 = 1`` on paths without +exceptions, are counted against the model, which keeps the test conservative. +""" + +from __future__ import annotations + +from dataclasses import dataclass + +import numpy as np +import numpy.typing as npt +import pandas as pd + +from quantkit._validation import FloatArray, SeedLike, as_1d_array, require_int, require_level +from quantkit.exceptions import InvalidParameterError +from quantkit.var.forecasts import RiskForecast + +# Simulate at most this many path-days at once to bound memory use (about 16 MB). +_MAX_DRAWS_PER_BATCH = 2_000_000 + + +@dataclass(frozen=True, slots=True) +class ESBacktest: + """The Z2 statistic and its simulated null distribution. + + Attributes: + statistic: Observed ``Z2``; 0 is ideal, negative values mean ES was too low. + p_value: Simulated one-sided p-value, ``P(Z2 <= statistic)`` under the null. + critical_value: The ``significance`` quantile of the simulated ``Z2``; + the test rejects when ``statistic`` is below it. + significance: Significance level of ``critical_value``. + n_simulations: Number of simulated paths. + """ + + statistic: float + p_value: float + critical_value: float + significance: float + n_simulations: int + + +def acerbi_szekely_z2( + returns: pd.Series | npt.ArrayLike, + var: pd.Series | npt.ArrayLike, + es: pd.Series | npt.ArrayLike, + *, + level: float, +) -> float: + """Acerbi and Székely's ``Z2`` statistic for VaR and ES forecasts. + + Args: + returns: Realized returns. + var: VaR forecasts for the same days, as positive loss fractions. + es: ES forecasts for the same days, positive. + level: Confidence level of the forecasts, e.g. 0.975. + + Returns: + ``Z2``: 1 if there are no exceptions, about 0 for correct forecasts on + average, negative when losses beyond VaR exceed the ES forecasts. + """ + r = as_1d_array("returns", returns) + v = as_1d_array("var", var) + e = as_1d_array("es", es) + if not r.shape == v.shape == e.shape: + raise InvalidParameterError( + f"returns, var and es must have the same length, got {r.size}, {v.size} and {e.size}" + ) + _require_positive_es(e) + return float(_z2(r[np.newaxis, :], v, e, 1.0 - require_level("level", level))[0]) + + +def es_backtest( + forecast: RiskForecast, + *, + n_simulations: int = 10_000, + significance: float = 0.05, + seed: SeedLike = None, +) -> ESBacktest: + """Test a forecast's ES with ``Z2`` against its own simulated null distribution. + + Args: + forecast: Rolling forecasts from :func:`~quantkit.var.rolling_var_es`. + n_simulations: Number of null paths to simulate. + significance: Level for the reported critical value, e.g. 0.05. + seed: Seed or :class:`numpy.random.Generator` for reproducibility. + + Returns: + The observed statistic, its simulated p-value and critical value. + """ + n_simulations = require_int("n_simulations", n_simulations, minimum=1) + significance = require_level("significance", significance) + var = forecast.var.to_numpy() + es = forecast.es.to_numpy() + _require_positive_es(es) + alpha = 1.0 - forecast.level + observed = float(_z2(forecast.returns.to_numpy()[np.newaxis, :], var, es, alpha)[0]) + + rng = np.random.default_rng(seed) + batch = max(1, _MAX_DRAWS_PER_BATCH // len(forecast)) + simulated = np.concatenate( + [ + _z2(forecast.simulate(min(batch, n_simulations - start), seed=rng), var, es, alpha) + for start in range(0, n_simulations, batch) + ] + ) + return ESBacktest( + statistic=observed, + p_value=(1 + int(np.count_nonzero(simulated <= observed))) / (n_simulations + 1), + critical_value=float(np.quantile(simulated, significance)), + significance=significance, + n_simulations=n_simulations, + ) + + +def _z2(paths: FloatArray, var: FloatArray, es: FloatArray, alpha: float) -> FloatArray: + """``Z2`` of every row of ``paths`` (one return path per row).""" + tail_terms = np.where(paths < -var, paths / es, 0.0) + z2: FloatArray = tail_terms.sum(axis=1) / (paths.shape[1] * alpha) + 1.0 + return z2 + + +def _require_positive_es(es: FloatArray) -> None: + if not (es > 0.0).all(): + raise InvalidParameterError( + "Z2 divides by the ES forecasts, so every ES forecast must be positive" + ) diff --git a/src/quantkit/var/forecasts.py b/src/quantkit/var/forecasts.py new file mode 100644 index 0000000..8f9f254 --- /dev/null +++ b/src/quantkit/var/forecasts.py @@ -0,0 +1,301 @@ +"""Rolling one-day Value at Risk and expected shortfall forecasts. + +A forecast for day ``t`` is made at the close of day ``t - 1``. It uses the +returns of days ``t - h, ..., t - 1`` and nothing else: never the return +``r[t]`` that it is scored against, and nothing older than its history ``h``. + +Every method describes the next return as a location-scale model:: + + r[t] = mu[t] + sigma[t] * e[t] + +The methods differ in how ``mu[t]`` and ``sigma[t]`` are estimated from the +window ``r[t - window], ..., r[t - 1]`` and in where the innovation ``e[t]`` +comes from: + +========== =========== ==================== =============================== ========== +method mu[t] sigma[t] e[t] history h +========== =========== ==================== =============================== ========== +historical 0 1 resampled window returns window +gaussian window mean window std (ddof=1) standard normal window +ewma 0 EWMA volatility standard normal window +fhs 0 EWMA volatility resampled standardized returns 2 * window +========== =========== ==================== =============================== ========== + +* **EWMA** (RiskMetrics, J.P. Morgan/Reuters 1996) weights the squared return + ``k`` days back by ``decay**k``, normalized over the window so the weights sum + to one: ``sigma[t]**2 = sum_k w_k r[t-1-k]**2``. The mean is taken as zero. +* **Filtered historical simulation** (Barone-Adesi, Giannopoulos & Vosper 1999; + the volatility weighting of Hull & White 1998) divides each window return + ``r[s]`` by the EWMA volatility forecast for its own day ``s``, then rescales + the standardized returns by today's volatility. Each of those volatilities + needs its own window, so FHS needs two windows of history. + +VaR and ES at confidence level ``c`` are positive loss fractions, as in +:mod:`quantkit.portfolio.risk`: + +* standard normal innovations: ``VaR = z_c sigma - mu`` and + ``ES = sigma phi(z_c) / (1 - c) - mu``; +* empirical innovations: the tail is the ``m = ceil(window * (1 - c))`` smallest + innovations, ``VaR = -(mu + sigma e_(m))`` and ``ES = -(mu + sigma mean(tail))``. + For historical simulation these are exactly :func:`~quantkit.portfolio.historical_var` + and :func:`~quantkit.portfolio.historical_cvar` of the window. + +The same representation defines :meth:`RiskForecast.simulate`, which draws +returns from each day's forecast distribution; the expected shortfall backtest +uses it to simulate its null distribution. +""" + +from __future__ import annotations + +from dataclasses import dataclass, replace +from statistics import NormalDist +from typing import Literal + +import numpy as np +import numpy.typing as npt +import pandas as pd +from numpy.lib.stride_tricks import sliding_window_view + +from quantkit._validation import ( + FloatArray, + SeedLike, + as_1d_array, + require_int, + require_level, + tail_count, +) +from quantkit.exceptions import InvalidParameterError +from quantkit.var.backtests import var_exceptions + +Method = Literal["historical", "gaussian", "ewma", "fhs"] +"""A rolling forecasting method; see the module docstring.""" + +METHODS: tuple[Method, ...] = ("historical", "gaussian", "ewma", "fhs") + +RISKMETRICS_DECAY = 0.94 +"""Decay factor that RiskMetrics uses for daily volatility forecasts.""" + +_STANDARD_NORMAL = NormalDist() + + +@dataclass(frozen=True) +class RiskForecast: + """Rolling one-day VaR and ES forecasts and the returns they are scored against. + + Every series is indexed by the day being forecast: ``var.iloc[k]`` was made at + the previous close, before ``returns.iloc[k]`` was known. + + Attributes: + method: The forecasting method. + level: Confidence level, e.g. 0.99. + window: Number of past returns in each estimation window. + decay: EWMA decay factor (used by ``ewma`` and ``fhs`` only). + returns: Realized returns on the forecast days. + var: VaR forecasts, as positive loss fractions. + es: Expected shortfall forecasts; ``es >= var``. + location: ``mu[t]`` of the location-scale model. + scale: ``sigma[t]`` of the location-scale model. + innovations: For ``historical`` and ``fhs``, each day's pool of innovations + (one row per forecast day, ``window`` columns, oldest first). None when + the innovations are standard normal. + """ + + method: Method + level: float + window: int + decay: float + returns: pd.Series + var: pd.Series + es: pd.Series + location: pd.Series + scale: pd.Series + innovations: FloatArray | None + + def __len__(self) -> int: + """Number of forecast days.""" + return len(self.var) + + @property + def exceptions(self) -> pd.Series: + """True on the days whose loss exceeded the VaR forecast (``return < -VaR``).""" + return var_exceptions(self.returns, self.var) + + def last(self, n_days: int) -> RiskForecast: + """The same forecasts restricted to the most recent ``n_days`` days.""" + n_days = require_int("n_days", n_days, minimum=1) + if n_days > len(self): + raise InvalidParameterError( + f"n_days must be at most the number of forecasts ({len(self)}), got {n_days}" + ) + keep = slice(len(self) - n_days, None) + return replace( + self, + returns=self.returns.iloc[keep], + var=self.var.iloc[keep], + es=self.es.iloc[keep], + location=self.location.iloc[keep], + scale=self.scale.iloc[keep], + innovations=None if self.innovations is None else self.innovations[keep], + ) + + def simulate(self, n_paths: int, *, seed: SeedLike = None) -> FloatArray: + """Draw returns from each day's forecast distribution. + + Column ``k`` holds independent draws of ``location[k] + scale[k] * e``, with + ``e`` standard normal or drawn uniformly from that day's innovation pool. + If the forecasts are correct, the realized returns are one such path. + + Args: + n_paths: Number of simulated paths (rows). + seed: Seed or :class:`numpy.random.Generator` for reproducibility. + + Returns: + An array of shape ``(n_paths, len(self))``. + """ + n_paths = require_int("n_paths", n_paths, minimum=1) + rng = np.random.default_rng(seed) + n_days = len(self) + if self.innovations is None: + draws = rng.standard_normal((n_paths, n_days)) + else: + picks = rng.integers(0, self.window, size=(n_paths, n_days)) + draws = self.innovations[np.arange(n_days), picks] + simulated: FloatArray = self.location.to_numpy() + self.scale.to_numpy() * draws + return simulated + + +def rolling_var_es( + returns: pd.Series | npt.ArrayLike, + *, + method: Method = "historical", + level: float = 0.99, + window: int = 250, + decay: float = RISKMETRICS_DECAY, +) -> RiskForecast: + """Forecast one-day VaR and ES for every day that has enough history. + + Args: + returns: Daily returns in time order, e.g. a portfolio's simple returns. + A Series keeps its index (a DatetimeIndex must be sorted and unique); + other sequences get a RangeIndex. + method: ``"historical"``, ``"gaussian"``, ``"ewma"`` or ``"fhs"``. + level: Confidence level strictly between 0 and 1, e.g. 0.99. + window: Number of past returns used by each forecast (at least 2). FHS + uses two windows: see the module docstring. + decay: EWMA decay factor strictly between 0 and 1 (``ewma`` and ``fhs``). + + Returns: + Forecasts for the days ``h, ..., n - 1``, where ``h`` is the method's + history, aligned with the realized returns of those days. + + Raises: + InvalidParameterError: For an unknown method, invalid parameters, too few + returns, or an FHS window whose EWMA volatility is zero. + """ + if method not in METHODS: + raise InvalidParameterError(f"method must be one of {METHODS}, got {method!r}") + level = require_level("level", level) + window = require_int("window", window, minimum=2) + decay = require_level("decay", decay) + series = _return_series(returns) + r = series.to_numpy() + history = 2 * window if method == "fhs" else window + if len(r) <= history: + raise InvalidParameterError( + f"the {method} method with window {window} needs more than {history} returns, " + f"got {len(r)}" + ) + + # Row k holds r[k], ..., r[k + window - 1]: the window for day t = k + window. + windows = sliding_window_view(r, window)[: len(r) - window] + n_days = len(r) - history + location = np.zeros(n_days) + innovations: FloatArray | None = None + if method == "historical": + scale = np.ones(n_days) + innovations = windows + elif method == "gaussian": + location = windows.mean(axis=1) + scale = windows.std(axis=1, ddof=1) + else: + volatility = _ewma_volatility(windows, decay) # forecasts for days window..n-1 + if method == "ewma": + scale = volatility + else: + if (volatility <= 0.0).any(): + day = series.index[window + int(np.argmax(volatility <= 0.0))] + raise InvalidParameterError( + f"the EWMA volatility forecast for {day} is zero, so FHS cannot " + "standardize that return; every window needs a non-zero return" + ) + standardized = r[window:] / volatility + # Row j holds the standardized returns of days window + j, ..., + # 2 * window + j - 1: the pool for day t = 2 * window + j. + innovations = sliding_window_view(standardized, window)[:n_days] + scale = volatility[window:] + + var, es = _var_and_es(location, scale, innovations, level) + index = series.index[history:] + return RiskForecast( + method=method, + level=level, + window=window, + decay=decay, + returns=series.iloc[history:], + var=pd.Series(var, index=index, name="var"), + es=pd.Series(es, index=index, name="es"), + location=pd.Series(location, index=index, name="location"), + scale=pd.Series(scale, index=index, name="scale"), + innovations=innovations, + ) + + +def ewma_weights(window: int, decay: float = RISKMETRICS_DECAY) -> FloatArray: + """EWMA weights for a window of squared returns, oldest first, summing to one. + + The return ``k`` days before the forecast day gets weight proportional to + ``decay**k``. RiskMetrics uses ``(1 - decay) * decay**k`` over an unlimited + history; normalizing over a finite window differs from that by the omitted + weight ``decay**window`` (about 2e-7 for 250 days at 0.94). + """ + window = require_int("window", window, minimum=1) + decay = require_level("decay", decay) + weights: FloatArray = np.power(decay, np.arange(window - 1, -1, -1, dtype=np.float64)) + normalized: FloatArray = weights / weights.sum() + return normalized + + +def _ewma_volatility(windows: FloatArray, decay: float) -> FloatArray: + volatility: FloatArray = np.sqrt(np.square(windows) @ ewma_weights(windows.shape[1], decay)) + return volatility + + +def _var_and_es( + location: FloatArray, scale: FloatArray, innovations: FloatArray | None, level: float +) -> tuple[FloatArray, FloatArray]: + """VaR and ES of ``location + scale * e`` for normal or empirical innovations.""" + if innovations is None: + z = _STANDARD_NORMAL.inv_cdf(level) + var = z * scale - location + es = scale * (_STANDARD_NORMAL.pdf(z) / (1.0 - level)) - location + return var, es + m = tail_count(innovations.shape[1], level) + # After partitioning, column m - 1 holds the m-th smallest innovation and the + # columns before it the m - 1 smaller ones (in no particular order). + tail = np.partition(innovations, m - 1, axis=1)[:, :m] + var = -(location + scale * tail[:, -1]) + es = -(location + scale * tail.mean(axis=1)) + return var, es + + +def _return_series(returns: pd.Series | npt.ArrayLike) -> pd.Series: + """A float copy of ``returns`` as a Series, after checking it is a usable series.""" + values = as_1d_array("returns", returns, min_length=2).copy() + if not isinstance(returns, pd.Series): + return pd.Series(values, name="returns") + index = returns.index + if index.has_duplicates: + raise InvalidParameterError("returns must not have duplicate index labels") + if isinstance(index, pd.DatetimeIndex) and not index.is_monotonic_increasing: + raise InvalidParameterError("returns must be sorted by date") + return pd.Series(values, index=index, name="returns") diff --git a/src/quantkit/var/traffic_light.py b/src/quantkit/var/traffic_light.py new file mode 100644 index 0000000..b3b7aa1 --- /dev/null +++ b/src/quantkit/var/traffic_light.py @@ -0,0 +1,94 @@ +"""The Basel Committee's traffic-light classification of VaR backtests. + +Banks using internal models compare each day's 99% one-day VaR with the next +day's trading outcome over the most recent 250 trading days, and count the +exceptions. The count puts the model in a zone and sets a "plus factor": an +increase in the multiplier, at least 3, that scales the VaR-based market risk +capital charge (Basel Committee on Banking Supervision 1996, Table 2): + +====== ========== =========== ====================== +zone exceptions plus factor cumulative probability +====== ========== =========== ====================== +green 0 to 4 0.00 8.11% to 89.22% +yellow 5 0.40 95.88% +yellow 6 0.50 98.63% +yellow 7 0.65 99.60% +yellow 8 0.75 99.89% +yellow 9 0.85 99.97% +red 10 or more 1.00 99.99% +====== ========== =========== ====================== + +The cumulative probability is the chance that a correct model produces that +many exceptions or fewer: a binomial distribution with 250 trials and +probability 1%. The green zone holds the counts whose cumulative probability is +below 95%, and the red zone starts where it reaches 99.99%. + +Reference: Basel Committee on Banking Supervision (1996), *Supervisory framework +for the use of "backtesting" in conjunction with the internal models approach to +market risk capital requirements*, https://www.bis.org/publ/bcbs22.htm +""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import Literal + +from quantkit._validation import require_int +from quantkit.exceptions import InvalidParameterError + +Zone = Literal["green", "yellow", "red"] +"""A Basel traffic-light zone.""" + +BASEL_DAYS = 250 +"""Number of daily forecasts the traffic light counts exceptions over.""" + +BASEL_LEVEL = 0.99 +"""VaR confidence level the traffic light applies to.""" + +BASE_MULTIPLIER = 3.0 +"""Minimum multiplication factor applied to VaR in the capital charge.""" + +_YELLOW_PLUS_FACTORS = {5: 0.40, 6: 0.50, 7: 0.65, 8: 0.75, 9: 0.85} + + +@dataclass(frozen=True, slots=True) +class TrafficLight: + """The Basel zone and capital plus factor for a count of exceptions. + + Attributes: + zone: ``"green"``, ``"yellow"`` or ``"red"``. + n_exceptions: Exceptions among the last 250 daily 99% VaR forecasts. + plus_factor: Increase in the capital multiplier. + """ + + zone: Zone + n_exceptions: int + plus_factor: float + + @property + def multiplier(self) -> float: + """Capital multiplier: ``3 + plus_factor``.""" + return BASE_MULTIPLIER + self.plus_factor + + +def basel_traffic_light(n_exceptions: int) -> TrafficLight: + """Classify the number of exceptions of 250 daily 99% VaR forecasts. + + Args: + n_exceptions: Exceptions over the most recent 250 trading days. The zones + are defined for that sample size and for 99% VaR only. + + Returns: + The zone and plus factor from Table 2 of the Basel Committee's 1996 + backtesting framework. + """ + count = require_int("n_exceptions", n_exceptions, minimum=0) + if count > BASEL_DAYS: + raise InvalidParameterError( + f"n_exceptions cannot exceed the {BASEL_DAYS} days counted, got {count}" + ) + if count < 5: + return TrafficLight("green", count, 0.0) + if count < 10: + return TrafficLight("yellow", count, _YELLOW_PLUS_FACTORS[count]) + return TrafficLight("red", count, 1.0) diff --git a/tests/test_cli.py b/tests/test_cli.py index 491cd3a..366d0b4 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -84,6 +84,49 @@ def test_frontier_prints_a_table_and_saves_a_plot( assert plot.stat().st_size > 10_000 +def test_var_backtest_compares_every_method_on_common_days( + prices_csv: Path, capsys: pytest.CaptureFixture[str] +) -> None: + args = ["var-backtest", str(prices_csv), "--window", "100", "--es-simulations", "199"] + assert main(args) == 0 + out = capsys.readouterr().out + assert "equal-weight portfolio of 5 assets, rebalanced daily: 399 daily returns" in out + # FHS needs 2 x 100 returns of history, which leaves 199 days for every method. + assert "backtest: 199 days" in out + assert "the period every method covers" in out + header = next(line for line in out.splitlines() if "historical" in line) + assert header.split() == ["historical", "Gaussian", "EWMA", "FHS"] + for row in ("Kupiec LR_uc p-value", "conditional coverage LR_cc p-value", "Z2 simulated"): + assert row in out + assert "n/a (< 250 days)" in out + assert "a correct model fails the 5% tests" in out + + +def test_var_backtest_one_method_one_ticker( + prices_csv: Path, capsys: pytest.CaptureFixture[str] +) -> None: + base = ["var-backtest", str(prices_csv), "--method", "ewma", "--window", "100"] + assert main([*base, "--ticker", "CCC", "--es-simulations", "99"]) == 0 + out = capsys.readouterr().out + assert out.startswith("CCC: 399 daily returns") + assert "backtest: 299 days, " in out + zone = next(line for line in out.splitlines() if line.startswith("Basel zone")) + assert zone.split()[-2] in {"green", "yellow", "red"} + assert main([*base, "--level", "0.975", "--es-simulations", "99"]) == 0 + out = capsys.readouterr().out + assert "97.5% one-day VaR and ES, 100-day window, EWMA decay 0.94" in out + assert "n/a (99% only)" in out + + +def test_var_backtest_reports_bad_inputs( + prices_csv: Path, capsys: pytest.CaptureFixture[str] +) -> None: + assert main(["var-backtest", str(prices_csv), "--ticker", "ZZZ"]) == 1 + assert "ticker 'ZZZ' is not in" in capsys.readouterr().err + assert main(["var-backtest", str(prices_csv), "--method", "fhs"]) == 1 + assert "needs more than 500 returns" in capsys.readouterr().err + + def test_python_dash_m_entry_point( monkeypatch: pytest.MonkeyPatch, capsys: pytest.CaptureFixture[str] ) -> None: diff --git a/tests/test_var_backtests.py b/tests/test_var_backtests.py new file mode 100644 index 0000000..4f3f018 --- /dev/null +++ b/tests/test_var_backtests.py @@ -0,0 +1,304 @@ +"""Exception backtests against hand calculations, SciPy and brute-force enumeration.""" + +import itertools +import math + +import numpy as np +import pandas as pd +import pytest +from hypothesis import given +from hypothesis import strategies as st +from scipy.stats import binom, chi2 + +from quantkit.exceptions import InvalidParameterError +from quantkit.var import ( + BASEL_DAYS, + TransitionCounts, + backtest_var, + basel_traffic_light, + christoffersen_test, + conditional_coverage_test, + kupiec_test, + null_rejection_rates, + transition_counts, + var_exceptions, +) + +hit_sequences = st.lists(st.booleans(), min_size=2, max_size=400) +levels = st.floats(0.8, 0.999) + + +def hits(n_days: int, exception_days: list[int]) -> np.ndarray: + flags = np.zeros(n_days, dtype=bool) + flags[exception_days] = True + return flags + + +# ------------------------------------------------------------- Kupiec (1995) + + +def test_kupiec_by_hand() -> None: + # Six exceptions in 250 days at 99% (2.5 expected), Kupiec's LR_uc: + # 2 [6 ln(6 / 2.5) + 244 ln(244 / 247.5)] = 2 [5.25281 - 3.47514] = 3.5554. + result = kupiec_test(hits(250, [3, 50, 90, 130, 170, 210]), level=0.99) + expected = 2 * (6 * math.log(6 / 2.5) + 244 * math.log(244 / 247.5)) + assert result.statistic == pytest.approx(expected, rel=1e-12) + assert result.statistic == pytest.approx(3.5554, abs=1e-4) + assert result.dof == 1 + assert result.p_value == pytest.approx(float(chi2.sf(expected, 1)), rel=1e-12) + # Seven exceptions cross the 5% critical value 3.84: 2 [7 ln 2.8 + 243 ln(243 / 247.5)]. + seven = kupiec_test(hits(250, [3, 50, 90, 130, 170, 210, 240]), level=0.99) + assert seven.statistic == pytest.approx( + 2 * (7 * math.log(2.8) + 243 * math.log(243 / 247.5)), rel=1e-12 + ) + assert seven.p_value < 0.05 < result.p_value + + +def test_kupiec_with_no_exceptions_or_only_exceptions() -> None: + # pi_hat = 0 or 1 makes the alternative likelihood exactly 1 (0^0 = 1). + none = kupiec_test(np.zeros(250, dtype=bool), level=0.99) + assert none.statistic == pytest.approx(-2 * 250 * math.log(0.99), rel=1e-12) # 5.025 + assert none.p_value == pytest.approx(float(chi2.sf(none.statistic, 1)), rel=1e-12) + assert none.p_value < 0.05 # two-sided: too few exceptions also reject + only = kupiec_test(np.ones(20, dtype=bool), level=0.99) + assert only.statistic == pytest.approx(-2 * 20 * math.log(0.01), rel=1e-12) + assert kupiec_test([True], level=0.5).statistic == pytest.approx(2 * math.log(2)) + + +@given(n_days=st.integers(1, 2000), share=st.floats(0.0, 1.0), level=st.floats(0.5, 0.9999)) +def test_kupiec_matches_binomial_log_likelihoods(n_days: int, share: float, level: float) -> None: + """An independent formula: the binomial coefficient cancels in the ratio.""" + x = round(share * n_days) + result = kupiec_test(hits(n_days, list(range(x))), level=level) + expected = 2 * (binom.logpmf(x, n_days, x / n_days) - binom.logpmf(x, n_days, 1 - level)) + assert result.statistic == pytest.approx(float(expected), rel=1e-9, abs=1e-9) + + +# ------------------------------------------------------- Christoffersen (1998) + +CLUSTERED = [0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 1, 0] + + +def test_transition_counts_by_hand() -> None: + assert transition_counts(CLUSTERED) == TransitionCounts(n00=4, n01=2, n10=2, n11=3) + assert transition_counts([True, False]) == TransitionCounts(n00=0, n01=0, n10=1, n11=0) + + +def test_christoffersen_by_hand() -> None: + # n00 = 4, n01 = 2, n10 = 2, n11 = 3: pi_01 = 2/6, pi_11 = 3/5, pooled pi = 5/11. + markov = 4 * math.log(2 / 3) + 2 * math.log(1 / 3) + 2 * math.log(2 / 5) + 3 * math.log(3 / 5) + iid = 6 * math.log(6 / 11) + 5 * math.log(5 / 11) + independence = christoffersen_test(CLUSTERED) + assert independence.statistic == pytest.approx(2 * (markov - iid), rel=1e-12) # 0.7899 + assert independence.statistic == pytest.approx(0.7899, abs=1e-4) + assert independence.p_value == pytest.approx(float(chi2.sf(2 * (markov - iid), 1))) + # Conditional coverage at 95%: L(p) = 0.95^6 * 0.05^5 over the 11 transitions. + coverage = conditional_coverage_test(CLUSTERED, level=0.95) + expected = 2 * (markov - 6 * math.log(0.95) - 5 * math.log(0.05)) + assert coverage.statistic == pytest.approx(expected, rel=1e-12) # 16.2046 + assert coverage.dof == 2 + assert coverage.p_value == pytest.approx(float(chi2.sf(expected, 2)), rel=1e-12) + + +def direct_markov_loglik(flags: list[bool]) -> float: + """Markov-chain log-likelihood by walking the day pairs (shares no code).""" + pairs = list(itertools.pairwise(flags)) + total = 0.0 + for state in (False, True): + following = [after for before, after in pairs if before == state] + if following: + pi = sum(following) / len(following) + total += sum(math.log(pi if after else 1 - pi) for after in following) + return total + + +@given(flags=hit_sequences) +def test_independence_matches_a_direct_likelihood(flags: list[bool]) -> None: + after = flags[1:] + pi = sum(after) / len(after) + iid = sum(math.log(pi if hit else 1 - pi) for hit in after) + expected = max(0.0, 2 * (direct_markov_loglik(flags) - iid)) + assert christoffersen_test(flags).statistic == pytest.approx(expected, abs=1e-9) + + +@given(flags=hit_sequences, level=levels) +def test_conditional_coverage_is_kupiec_on_days_two_onward_plus_independence( + flags: list[bool], level: float +) -> None: + # Christoffersen (1998): LR_cc = LR_uc + LR_ind when LR_uc conditions on the first day. + combined = conditional_coverage_test(flags, level=level).statistic + parts = kupiec_test(flags[1:], level=level).statistic + christoffersen_test(flags).statistic + assert combined == pytest.approx(parts, rel=1e-9, abs=1e-9) + + +@given(flags=hit_sequences, level=levels) +def test_p_values_match_scipy_chi2(flags: list[bool], level: float) -> None: + for result in ( + kupiec_test(flags, level=level), + christoffersen_test(flags), + conditional_coverage_test(flags, level=level), + ): + assert result.statistic >= 0 + assert result.p_value == pytest.approx( + float(chi2.sf(result.statistic, result.dof)), rel=1e-9, abs=1e-15 + ) + + +def test_edge_cases_have_finite_statistics() -> None: + quiet, busy = np.zeros(250, dtype=bool), np.ones(250, dtype=bool) + for flags, p_tail in ((quiet, 0.99), (busy, 0.01)): + # No transition between states: both likelihoods are equal, LR_ind = 0. + assert christoffersen_test(flags).statistic == 0.0 + assert christoffersen_test(flags).p_value == 1.0 + # LR_cc then reduces to Kupiec on the 249 transitions. + coverage = conditional_coverage_test(flags, level=0.99) + assert coverage.statistic == pytest.approx(-2 * 249 * math.log(p_tail), rel=1e-12) + # Isolated exceptions only (n11 = 0): pi_11 = 0 fits exactly and LR_ind is finite. + isolated = hits(250, [10, 100, 200]) + counts = transition_counts(isolated) + assert counts == TransitionCounts(n00=243, n01=3, n10=3, n11=0) + markov = 243 * math.log(243 / 246) + 3 * math.log(3 / 246) # the n10 term is 3 ln 1 = 0 + iid = 246 * math.log(246 / 249) + 3 * math.log(3 / 249) + assert christoffersen_test(isolated).statistic == pytest.approx(2 * (markov - iid)) # 0.0732 + # The only exception on the last day: no pair starts with an exception at all. + last = hits(10, [9]) + assert transition_counts(last) == TransitionCounts(n00=8, n01=1, n10=0, n11=0) + assert christoffersen_test(last).statistic == 0.0 + assert math.isfinite(conditional_coverage_test(last, level=0.9).statistic) + + +def test_backtest_var_collects_everything() -> None: + flags = pd.Series(hits(250, [5, 6, 100, 180]), index=pd.bdate_range("2024-01-01", periods=250)) + result = backtest_var(flags, level=0.99) + assert result.n_observations == 250 + assert result.n_exceptions == 4 + assert result.expected_exceptions == pytest.approx(2.5) + assert result.exception_rate == pytest.approx(0.016) + assert result.transitions == TransitionCounts(n00=242, n01=3, n10=3, n11=1) + assert result.unconditional == kupiec_test(flags, level=0.99) + assert result.independence == christoffersen_test(flags) + assert result.conditional == conditional_coverage_test(flags, level=0.99) + + +def test_var_exceptions() -> None: + returns = pd.Series([-0.03, -0.02, 0.01], index=pd.bdate_range("2024-01-01", periods=3)) + var = pd.Series(0.02, index=returns.index) + assert list(var_exceptions(returns, var)) == [True, False, False] # -0.02 is not < -0.02 + assert var_exceptions(returns, var).index.equals(returns.index) + assert list(var_exceptions([-0.05, 0.0], [0.02, 0.02])) == [True, False] + with pytest.raises(InvalidParameterError, match="same length"): + var_exceptions([0.01, 0.02], [0.02]) + with pytest.raises(InvalidParameterError, match="same index"): + var_exceptions(returns, var.set_axis(pd.bdate_range("2025-01-01", periods=3))) + + +def test_rejects_invalid_exception_series() -> None: + with pytest.raises(InvalidParameterError, match="0/1"): + kupiec_test([0, 2, 1], level=0.99) + with pytest.raises(InvalidParameterError, match="at least 2"): + christoffersen_test([True]) + with pytest.raises(InvalidParameterError, match="one-dimensional"): + backtest_var(np.zeros((3, 2), dtype=bool), level=0.99) + with pytest.raises(InvalidParameterError, match="level"): + kupiec_test([0, 1], level=1.5) + assert kupiec_test([0, 1, 0], level=0.9) == kupiec_test([False, True, False], level=0.9) + + +# -------------------------------------------------------- exact null sizes + + +@pytest.mark.parametrize("n_days", [2, 3, 7, 12]) +@pytest.mark.parametrize("level", [0.6, 0.9]) +def test_exact_sizes_match_brute_force_enumeration(n_days: int, level: float) -> None: + """Sum the probability of every one of the 2^T exception sequences that is rejected.""" + p = 1 - level + rejected = np.zeros(3) + for flags in itertools.product((False, True), repeat=n_days): + x = sum(flags) + probability = p**x * (1 - p) ** (n_days - x) + tests = ( + kupiec_test(flags, level=level), + christoffersen_test(flags), + conditional_coverage_test(flags, level=level), + ) + rejected += probability * np.array([test.p_value < 0.05 for test in tests]) + rates = null_rejection_rates(n_days, level=level) + np.testing.assert_allclose( + [rates.unconditional, rates.independence, rates.conditional], rejected, atol=1e-12 + ) + + +@pytest.mark.parametrize("n_days", [250, 1000]) +def test_exact_kupiec_size_matches_the_binomial_distribution(n_days: int) -> None: + critical = float(chi2.ppf(0.95, 1)) + counts = np.arange(n_days + 1) + statistics = [ + 2 * (binom.logpmf(x, n_days, x / n_days) - binom.logpmf(x, n_days, 0.01)) for x in counts + ] + exact = float(binom.pmf(counts, n_days, 0.01)[np.array(statistics) > critical].sum()) + assert null_rejection_rates(n_days, level=0.99).unconditional == pytest.approx( + exact, abs=1e-12 + ) + + +def test_nominal_five_percent_tests_are_far_from_five_percent_at_250_days() -> None: + rates = null_rejection_rates(250, level=0.99) + # Zero exceptions alone (probability 0.99^250 = 8.1%) already rejects LR_uc. + assert rates.unconditional == pytest.approx(0.0948, abs=5e-5) + assert rates.unconditional > 0.99**250 + assert rates.independence == pytest.approx(0.0140, abs=5e-5) + assert rates.conditional == pytest.approx(0.0110, abs=5e-5) + + +def test_exact_sizes_refuse_huge_enumerations() -> None: + with pytest.raises(InvalidParameterError, match="too slow"): + null_rejection_rates(20_000, level=0.9) + with pytest.raises(InvalidParameterError, match="n_observations"): + null_rejection_rates(1, level=0.99) + + +# --------------------------------------------- Basel traffic light (BCBS 1996) + +# Table 2 of the Basel Committee's 1996 backtesting framework, 250 observations: +# exceptions -> (zone, plus factor, cumulative probability in percent). +BASEL_TABLE = { + 0: ("green", 0.00, 8.11), + 1: ("green", 0.00, 28.58), + 2: ("green", 0.00, 54.32), + 3: ("green", 0.00, 75.81), + 4: ("green", 0.00, 89.22), + 5: ("yellow", 0.40, 95.88), + 6: ("yellow", 0.50, 98.63), + 7: ("yellow", 0.65, 99.60), + 8: ("yellow", 0.75, 99.89), + 9: ("yellow", 0.85, 99.97), + 10: ("red", 1.00, 99.99), +} + + +@pytest.mark.parametrize(("count", "row"), BASEL_TABLE.items()) +def test_traffic_light_reproduces_the_basel_table( + count: int, row: tuple[str, float, float] +) -> None: + zone, plus_factor, cumulative = row + light = basel_traffic_light(count) + assert (light.zone, light.n_exceptions, light.plus_factor) == (zone, count, plus_factor) + assert light.multiplier == pytest.approx(3.0 + plus_factor) + # The published probabilities are binomial(250, 1%) CDFs rounded to 0.01%. + assert round(100 * float(binom.cdf(count, 250, 0.01)), 2) == cumulative + + +def test_zones_follow_the_cumulative_probability_rule() -> None: + # Green below a 95% cumulative probability, red from 99.99%, yellow between. + for count in range(BASEL_DAYS + 1): + cumulative = float(binom.cdf(count, BASEL_DAYS, 0.01)) + expected = "green" if cumulative < 0.95 else "yellow" if cumulative < 0.9999 else "red" + assert basel_traffic_light(count).zone == expected + + +def test_traffic_light_rejects_impossible_counts() -> None: + for bad in (-1, 251): + with pytest.raises(InvalidParameterError, match="n_exceptions"): + basel_traffic_light(bad) + with pytest.raises(InvalidParameterError, match="integer"): + basel_traffic_light(True) diff --git a/tests/test_var_es_backtest.py b/tests/test_var_es_backtest.py new file mode 100644 index 0000000..d91dce9 --- /dev/null +++ b/tests/test_var_es_backtest.py @@ -0,0 +1,109 @@ +"""Acerbi-Szekely Z2: hand calculations, its zero mean under the null, invariances.""" + +import math +from statistics import NormalDist + +import numpy as np +import pandas as pd +import pytest + +from quantkit.exceptions import InvalidParameterError +from quantkit.var import acerbi_szekely_z2, es_backtest, rolling_var_es + + +@pytest.fixture(scope="module") +def returns() -> pd.Series: + """SYNTHETIC normal daily returns with a volatility that doubles halfway.""" + rng = np.random.default_rng(21) + volatility = np.r_[np.full(300, 0.01), np.full(300, 0.02)] + return pd.Series(volatility * rng.standard_normal(600)) + + +def z2_of_paths(paths: np.ndarray, var: pd.Series, es: pd.Series, level: float) -> np.ndarray: + """Z2 of every row, written out from the paper's formula for many paths at once.""" + exceptions = paths < -var.to_numpy() + tail = np.where(exceptions, paths / es.to_numpy(), 0.0).sum(axis=1) + result: np.ndarray = tail / (paths.shape[1] * (1.0 - level)) + 1.0 + return result + + +def test_z2_by_hand() -> None: + # Exceptions on days 1 and 3 (returns below -VaR = -0.02); alpha = 1 - 0.75. + # Z2 = (-0.03 / 0.03 - 0.05 / 0.04) / (4 * 0.25) + 1 = -2.25 + 1 = -1.25. + z2 = acerbi_szekely_z2( + [-0.03, 0.01, -0.05, 0.02], [0.02, 0.02, 0.02, 0.02], [0.03, 0.03, 0.04, 0.04], level=0.75 + ) + assert z2 == pytest.approx(-1.25, rel=1e-14) + # A loss exactly at VaR is not an exception, and without exceptions Z2 = 1. + assert acerbi_szekely_z2([-0.02, 0.01], [0.02, 0.02], [0.03, 0.03], level=0.9) == 1.0 + + +def test_z2_is_zero_on_average_when_the_forecasts_are_right(returns: pd.Series) -> None: + # If returns follow the forecast distributions, E[r 1{r < -VaR}] = -alpha ES each + # day, so E[Z2] = 0 (Acerbi and Szekely 2014). Check on paths drawn from them. + for method in ("gaussian", "ewma"): + forecast = rolling_var_es(returns, method=method, level=0.975, window=100) + paths = forecast.simulate(20_000, seed=2) + z2 = z2_of_paths(paths, forecast.var, forecast.es, 0.975) + assert z2[0] == pytest.approx( + acerbi_szekely_z2(paths[0], forecast.var, forecast.es, level=0.975), rel=1e-12 + ) + assert abs(z2.mean()) < 4 * z2.std(ddof=1) / math.sqrt(z2.size) + + +def test_null_distribution_does_not_depend_on_the_volatility_path(returns: pd.Series) -> None: + # With zero-mean normal forecasts, r / sigma is standard normal and Z2 depends + # only on it, so two volatility paths share one null distribution. With the + # same seed the simulated statistics are identical up to round-off. + level, n_days = 0.99, 250 + calm = rolling_var_es(returns.iloc[:350], method="ewma", level=level, window=100) + wild = rolling_var_es(3.0 * returns.iloc[-350:], method="ewma", level=level, window=100) + first = es_backtest(calm.last(n_days), n_simulations=4_000, seed=9) + second = es_backtest(wild.last(n_days), n_simulations=4_000, seed=9) + assert first.critical_value == pytest.approx(second.critical_value, abs=1e-12) + # The same quantile from standard normal draws with constant VaR = z and ES = phi(z) / alpha. + z = NormalDist().inv_cdf(level) + draws = np.random.default_rng(9).standard_normal((4_000, n_days)) + direct = [ + acerbi_szekely_z2( + path, np.full(n_days, z), np.full(n_days, NormalDist().pdf(z) / 0.01), level=level + ) + for path in draws + ] + assert first.critical_value == pytest.approx(float(np.quantile(direct, 0.05)), abs=1e-12) + + +def test_es_backtest_reports_a_monte_carlo_p_value(returns: pd.Series) -> None: + forecast = rolling_var_es(returns, method="historical", level=0.975, window=250) + result = es_backtest(forecast, n_simulations=999, seed=4) + assert result.statistic == pytest.approx( + acerbi_szekely_z2(forecast.returns, forecast.var, forecast.es, level=0.975) + ) + assert 1 / 1000 <= result.p_value <= 1.0 + assert result.n_simulations == 999 + assert result.significance == 0.05 + # The volatility doubled on day 300, so a 250-day window badly underestimates ES. + assert result.statistic < result.critical_value + assert result.p_value < 0.05 + assert es_backtest(forecast, n_simulations=999, seed=4) == result + + +def test_es_backtest_batches_long_simulations(returns: pd.Series) -> None: + # 2,000,000 path-days per batch: 30,000 paths of 350 days need several batches, + # and must consume the random stream exactly as one large draw would. + forecast = rolling_var_es(returns, method="ewma", window=250) + batched = es_backtest(forecast, n_simulations=30_000, seed=1) + direct = z2_of_paths(forecast.simulate(30_000, seed=1), forecast.var, forecast.es, 0.99) + assert batched.critical_value == pytest.approx(float(np.quantile(direct, 0.05)), abs=1e-12) + + +def test_rejects_invalid_inputs(returns: pd.Series) -> None: + with pytest.raises(InvalidParameterError, match="same length"): + acerbi_szekely_z2([0.01, 0.02], [0.02], [0.03, 0.03], level=0.975) + with pytest.raises(InvalidParameterError, match="positive"): + acerbi_szekely_z2([0.01, 0.02], [0.02, 0.02], [0.03, 0.0], level=0.975) + forecast = rolling_var_es(returns, method="gaussian", window=100) + with pytest.raises(InvalidParameterError, match="n_simulations"): + es_backtest(forecast, n_simulations=0) + with pytest.raises(InvalidParameterError, match="significance"): + es_backtest(forecast, significance=1.0) diff --git a/tests/test_var_forecasts.py b/tests/test_var_forecasts.py new file mode 100644 index 0000000..11538f4 --- /dev/null +++ b/tests/test_var_forecasts.py @@ -0,0 +1,329 @@ +"""Rolling VaR/ES forecasts: single-window formulas, independent loops, no look-ahead.""" + +import math +from collections.abc import Callable +from statistics import NormalDist + +import numpy as np +import pandas as pd +import pytest +from hypothesis import given +from hypothesis import strategies as st +from scipy.stats import t as student_t + +from quantkit.exceptions import InvalidParameterError +from quantkit.portfolio import historical_cvar, historical_var, parametric_cvar, parametric_var +from quantkit.var import METHODS, Method, RiskForecast, ewma_weights, rolling_var_es + +WINDOW = 60 +Forecaster = Callable[[pd.Series], pd.DataFrame] + + +@pytest.fixture(scope="module") +def returns() -> pd.Series: + """SYNTHETIC daily returns: Student-t(5) shocks with slowly varying volatility.""" + rng = np.random.default_rng(7) + n = 400 + volatility = 0.01 * np.exp(0.5 * np.sin(np.arange(n) / 40.0)) + shocks = rng.standard_t(5, n) * math.sqrt(3.0 / 5.0) + return pd.Series(volatility * shocks, index=pd.bdate_range("2021-01-01", periods=n)) + + +def windows_before(r: np.ndarray, first: int, window: int) -> list[np.ndarray]: + """For each day t >= first, the returns r[t - window], ..., r[t - 1].""" + return [r[t - window : t] for t in range(first, len(r))] + + +@pytest.mark.parametrize("level", [0.95, 0.975, 0.99]) +def test_historical_matches_the_single_window_functions(returns: pd.Series, level: float) -> None: + forecast = rolling_var_es(returns, method="historical", level=level, window=WINDOW) + r = returns.to_numpy() + windows = windows_before(r, WINDOW, WINDOW) + # np.partition (rolling) and np.sort (single window) must pick the same returns. + expected_var = [historical_var(w, level=level) for w in windows] + expected_es = [historical_cvar(w, level=level) for w in windows] + np.testing.assert_array_equal(forecast.var, expected_var) + np.testing.assert_allclose(forecast.es, expected_es, rtol=1e-14) + assert forecast.var.index.equals(returns.index[WINDOW:]) + pd.testing.assert_series_equal(forecast.returns, returns.iloc[WINDOW:], check_names=False) + + +@pytest.mark.parametrize("level", [0.95, 0.99]) +def test_gaussian_matches_the_single_window_functions(returns: pd.Series, level: float) -> None: + forecast = rolling_var_es(returns, method="gaussian", level=level, window=WINDOW) + windows = windows_before(returns.to_numpy(), WINDOW, WINDOW) + np.testing.assert_allclose( + forecast.var, [parametric_var(w, level=level) for w in windows], rtol=1e-12 + ) + np.testing.assert_allclose( + forecast.es, [parametric_cvar(w, level=level) for w in windows], rtol=1e-12 + ) + + +def test_ewma_by_hand() -> None: + # Weights for lags 2, 1, 0 are 0.25, 0.5, 1 before normalizing by their sum 1.75. + returns = [0.01, -0.02, 0.03, 0.0] + forecast = rolling_var_es(returns, method="ewma", level=0.99, window=3, decay=0.5) + sigma = math.sqrt((0.25 * 0.01**2 + 0.5 * 0.02**2 + 1.0 * 0.03**2) / 1.75) + assert forecast.scale.iloc[0] == pytest.approx(sigma, rel=1e-14) + z = NormalDist().inv_cdf(0.99) + assert forecast.var.iloc[0] == pytest.approx(z * sigma, rel=1e-14) + assert forecast.es.iloc[0] == pytest.approx(sigma * NormalDist().pdf(z) / 0.01, rel=1e-14) + np.testing.assert_allclose(ewma_weights(3, 0.5), [0.25 / 1.75, 0.5 / 1.75, 1.0 / 1.75]) + + +def test_ewma_matches_the_riskmetrics_recursion(returns: pd.Series) -> None: + """The RiskMetrics update s2 <- decay * s2 + (1 - decay) * r^2, run over each window.""" + decay = 0.94 + forecast = rolling_var_es(returns, method="ewma", window=WINDOW, decay=decay) + expected = [] + for window in windows_before(returns.to_numpy(), WINDOW, WINDOW): + s2 = 0.0 + for r in window: + s2 = decay * s2 + (1.0 - decay) * r * r + # Started from zero, the recursion's weights sum to 1 - decay**WINDOW. + expected.append(math.sqrt(s2 / (1.0 - decay**WINDOW))) + np.testing.assert_allclose(forecast.scale, expected, rtol=1e-12) + assert (forecast.location == 0.0).all() + + +def test_fhs_matches_a_volatility_weighted_loop(returns: pd.Series) -> None: + """Hull and White's rescaling: r[s] * sigma[t] / sigma[s], then historical VaR.""" + window, level = 40, 0.975 + forecast = rolling_var_es(returns, method="fhs", level=level, window=window) + r = returns.to_numpy() + weights = ewma_weights(window) + + def sigma(day: int) -> float: + return math.sqrt(float(np.square(r[day - window : day]) @ weights)) + + expected_var, expected_es = [], [] + for t in range(2 * window, len(r)): + rescaled = [r[s] * sigma(t) / sigma(s) for s in range(t - window, t)] + expected_var.append(historical_var(rescaled, level=level)) + expected_es.append(historical_cvar(rescaled, level=level)) + np.testing.assert_allclose(forecast.var, expected_var, rtol=1e-12) + np.testing.assert_allclose(forecast.es, expected_es, rtol=1e-12) + assert forecast.var.index[0] == returns.index[2 * window] + + +@given(factor=st.floats(0.05, 20.0), method=st.sampled_from(METHODS)) +def test_forecasts_scale_with_the_returns(factor: float, method: Method) -> None: + returns = np.random.default_rng(3).normal(0.0005, 0.01, 200) + base = rolling_var_es(returns, method=method, window=50) + scaled = rolling_var_es(factor * returns, method=method, window=50) + np.testing.assert_allclose(scaled.var, factor * base.var, rtol=1e-9) + np.testing.assert_allclose(scaled.es, factor * base.es, rtol=1e-9) + + +@pytest.mark.parametrize("method", METHODS) +def test_es_is_at_least_var(returns: pd.Series, method: Method) -> None: + forecast = rolling_var_es(returns, method=method, level=0.975, window=WINDOW) + assert (forecast.es >= forecast.var - 1e-15).all() + assert (forecast.var > 0).all() + + +# --------------------------------------------------------------------- look-ahead + + +def as_forecaster(method: Method, window: int = WINDOW) -> Forecaster: + def forecaster(r: pd.Series) -> pd.DataFrame: + forecast = rolling_var_es(r, method=method, window=window) + return pd.DataFrame({"var": forecast.var, "es": forecast.es}) + + return forecaster + + +def peeking_forecaster(r: pd.Series) -> pd.DataFrame: + """A deliberately broken forecaster whose window ends on the day it forecasts.""" + tomorrow = pd.Series(r.to_numpy()[1:], index=r.index[:-1]) + forecast = rolling_var_es(tomorrow, method="historical", window=WINDOW) + return pd.DataFrame({"var": forecast.var, "es": forecast.es}) + + +def lookahead_days( + forecaster: Forecaster, returns: pd.Series, *, cutoffs: int = 6, seed: int = 0 +) -> list[pd.Timestamp]: + """Cutoff days whose forecast changed when only that day's and later returns changed. + + Like ``check_no_lookahead`` for strategies: the returns from each cutoff day + on are scrambled, and every forecast up to and including the cutoff day + (whose forecast is scored against the scrambled return) must stay identical. + The cutoff day itself becomes a 50% loss, so a forecaster that peeks at it + changes even when its VaR is a single order statistic, such as the worst + return in the window. + """ + rng = np.random.default_rng(seed) + baseline = forecaster(returns) + first = returns.index.get_loc(baseline.index[0]) + assert isinstance(first, int) + violations = [] + for cutoff in sorted(rng.choice(np.arange(first, len(returns)), cutoffs, replace=False)): + scrambled = returns.copy() + later = scrambled.iloc[cutoff:] + scrambled.iloc[cutoff:] = later.to_numpy() * rng.uniform(-3.0, 3.0, len(later)) + scrambled.iloc[cutoff] = -0.5 + rerun = forecaster(scrambled) + through_cutoff = baseline.index <= returns.index[cutoff] + if not np.array_equal( + baseline.to_numpy()[through_cutoff], rerun.to_numpy()[through_cutoff] + ): + violations.append(returns.index[cutoff]) + return violations + + +@pytest.mark.parametrize("method", METHODS) +def test_no_forecast_uses_the_return_it_is_scored_against( + returns: pd.Series, method: Method +) -> None: + assert lookahead_days(as_forecaster(method), returns) == [] + + +def test_lookahead_check_catches_a_peeking_forecaster(returns: pd.Series) -> None: + # Negative control: the check must fail for a forecaster that sees r[t]. + assert len(lookahead_days(peeking_forecaster, returns)) == 6 + + +@pytest.mark.parametrize("method", METHODS) +def test_forecasts_ignore_returns_older_than_their_history( + returns: pd.Series, method: Method +) -> None: + history = 2 * WINDOW if method == "fhs" else WINDOW + baseline = as_forecaster(method)(returns) + for day in (history + 5, len(returns) - 1): + scrambled = returns.copy() + scrambled.iloc[: day - history] = np.random.default_rng(day).normal(0, 0.05, day - history) + rerun = as_forecaster(method)(scrambled) + label = returns.index[day] + np.testing.assert_array_equal(baseline.loc[label], rerun.loc[label]) + + +# ------------------------------------------------- statistical checks on i.i.d. data + + +def exception_rate_and_error( + method: Method, window: int, paths: np.ndarray +) -> tuple[float, float]: + """Exception rate over independent paths, with a standard error from their spread.""" + rates = np.array( + [rolling_var_es(path, method=method, window=window).exceptions.mean() for path in paths] + ) + return float(rates.mean()), float(rates.std(ddof=1) / math.sqrt(len(rates))) + + +def test_historical_exception_probability_is_m_over_window_plus_one() -> None: + # With i.i.d. continuous returns, r[t] is equally likely to take each of the + # window + 1 ranks among r[t - window..t], so it falls below the m-th smallest + # window return with probability m / (window + 1), whatever the distribution. + # At 99% with 250 days, m = 3: the exception probability is 1.195%, not 1%. + paths = np.random.default_rng(11).standard_t(3, (200, 1250)) + rate, error = exception_rate_and_error("historical", 250, paths) + assert abs(rate - 3 / 251) < 4 * error + assert abs(rate - 0.01) > 4 * error + + +def test_gaussian_exception_probability_includes_estimation_error() -> None: + # For i.i.d. normal returns, (r[t] - mean) / (std * sqrt(1 + 1 / w)) has a + # Student-t distribution with w - 1 degrees of freedom (the prediction + # interval), so the exception probability is t_{w-1}(-z / sqrt(1 + 1 / w)). + window, z = 30, NormalDist().inv_cdf(0.99) + exact = float(student_t.cdf(-z / math.sqrt(1 + 1 / window), window - 1)) # 1.479% + paths = np.random.default_rng(12).normal(0.0, 0.01, (400, 1030)) + rate, error = exception_rate_and_error("gaussian", window, paths) + assert abs(rate - exact) < 4 * error + assert abs(rate - 0.01) > 4 * error + + +@pytest.mark.parametrize("method", METHODS) +def test_simulated_returns_follow_each_days_forecast(returns: pd.Series, method: Method) -> None: + level = 0.95 + forecast = rolling_var_es(returns, method=method, level=level, window=WINDOW).last(5) + draws = forecast.simulate(100_000, seed=1) + below = draws <= -forecast.var.to_numpy() + if forecast.innovations is None: + # Normal innovations: P(r <= -VaR) = 1 - level and E[r | r <= -VaR] = -ES. + np.testing.assert_allclose(below.mean(axis=0), 1 - level, atol=0.0035) + tail_mean = (draws * below).sum(axis=0) / below.sum(axis=0) + np.testing.assert_allclose(tail_mean, -forecast.es, rtol=0.02) + else: + # Empirical innovations: every draw is location + scale * a pool value, + # and the m = 3 smallest of the 60 pool values carry probability 3 / 60. + pools = forecast.location.to_numpy() + forecast.scale.to_numpy() * forecast.innovations.T + assert all(np.isin(draws[:, k], pools[:, k]).all() for k in range(len(forecast))) + np.testing.assert_allclose(below.mean(axis=0), 3 / 60, atol=0.0035) + tail_mean = (draws * below).sum(axis=0) / below.sum(axis=0) + np.testing.assert_allclose(tail_mean, -forecast.es, rtol=0.01) + + +def test_simulation_is_reproducible(returns: pd.Series) -> None: + forecast = rolling_var_es(returns, method="fhs", window=WINDOW) + np.testing.assert_array_equal(forecast.simulate(3, seed=5), forecast.simulate(3, seed=5)) + assert forecast.simulate(3, seed=5).shape == (3, len(forecast)) + + +# --------------------------------------------------------------------- the result + + +def test_last_keeps_the_most_recent_days(returns: pd.Series) -> None: + forecast = rolling_var_es(returns, method="historical", window=WINDOW) + tail = forecast.last(10) + assert isinstance(tail, RiskForecast) + assert len(tail) == 10 + pd.testing.assert_series_equal(tail.var, forecast.var.iloc[-10:]) + assert tail.innovations is not None + assert forecast.innovations is not None + np.testing.assert_array_equal(tail.innovations, forecast.innovations[-10:]) + pd.testing.assert_series_equal(tail.exceptions, forecast.exceptions.iloc[-10:]) + with pytest.raises(InvalidParameterError, match="at most"): + forecast.last(len(forecast) + 1) + + +def test_exceptions_use_a_strict_inequality() -> None: + # The fourth return equals the VaR forecast (the worst of the window): not an exception. + forecast = rolling_var_es([-0.02, 0.01, 0.005, -0.02, -0.03], method="historical", window=3) + assert list(forecast.var) == [0.02, 0.02] + assert list(forecast.exceptions) == [False, True] + + +def test_arrays_get_a_range_index() -> None: + forecast = rolling_var_es(np.linspace(-0.02, 0.02, 30), method="gaussian", window=10) + assert list(forecast.var.index) == list(range(10, 30)) + + +# --------------------------------------------------------------------- validation + + +@pytest.mark.parametrize( + ("kwargs", "message"), + [ + ({"method": "garch"}, "method must be one of"), + ({"level": 1.0}, "level"), + ({"window": 1}, "window must be >= 2"), + ({"decay": 1.0}, "decay"), + ({"window": 400}, "needs more than 400 returns"), + ({"method": "fhs", "window": 200}, "needs more than 400 returns"), + ], +) +def test_rejects_invalid_arguments( + returns: pd.Series, kwargs: dict[str, object], message: str +) -> None: + with pytest.raises(InvalidParameterError, match=message): + rolling_var_es(returns, **kwargs) # type: ignore[arg-type] + + +def test_rejects_unusable_return_series(returns: pd.Series) -> None: + with pytest.raises(InvalidParameterError, match="NaN"): + rolling_var_es([0.01, math.nan, 0.02], window=2) + with pytest.raises(InvalidParameterError, match="sorted by date"): + rolling_var_es(returns.iloc[::-1], window=WINDOW) + duplicated = pd.Series(returns.to_numpy()[:100], index=[0, 1] * 50) + with pytest.raises(InvalidParameterError, match="duplicate"): + rolling_var_es(duplicated, window=WINDOW) + + +def test_fhs_rejects_a_window_without_movement() -> None: + flat = np.r_[np.full(30, 0.0), np.random.default_rng(0).normal(0, 0.01, 60)] + with pytest.raises(InvalidParameterError, match=r"volatility forecast .* is zero"): + rolling_var_es(flat, method="fhs", window=20) + with pytest.raises(InvalidParameterError, match="decay"): + ewma_weights(10, 0.0) diff --git a/tests/test_var_size_power.py b/tests/test_var_size_power.py new file mode 100644 index 0000000..6948336 --- /dev/null +++ b/tests/test_var_size_power.py @@ -0,0 +1,113 @@ +"""Size and power of the backtests in simulation, with fixed seeds. + +Size: returns are drawn from the forecast model itself, so every rejection is a +false alarm. Power: Gaussian VaR is applied to fat-tailed Student-t returns. +""" + +import math +from statistics import NormalDist + +import numpy as np +import pytest + +from quantkit.var import backtest_var, es_backtest, null_rejection_rates, rolling_var_es + +WINDOW, DECAY, LEVEL = 250, 0.94, 0.99 + + +def ewma_model_paths(n_paths: int, n_test: int, seed: int) -> tuple[np.ndarray, ...]: + """SYNTHETIC returns drawn from the EWMA model: r[t] = sigma[t] z[t]. + + ``sigma[t]**2`` is the normalized EWMA of the previous WINDOW squared returns, + computed here without quantkit. Returns the paths, sigma and the shocks z for + the last ``n_test`` days. + """ + rng = np.random.default_rng(seed) + weights = DECAY ** np.arange(WINDOW - 1, -1, -1.0) + weights /= weights.sum() + paths = np.empty((n_paths, WINDOW + n_test)) + paths[:, :WINDOW] = 0.01 * rng.standard_normal((n_paths, WINDOW)) + shocks = rng.standard_normal((n_paths, n_test)) + sigma = np.empty((n_paths, n_test)) + for k in range(n_test): + t = WINDOW + k + sigma[:, k] = np.sqrt(np.square(paths[:, t - WINDOW : t]) @ weights) + paths[:, t] = sigma[:, k] * shocks[:, k] + return paths, sigma, shocks + + +def standard_error(rate: float, n: int) -> float: + return math.sqrt(rate * (1 - rate) / n) + + +@pytest.fixture(scope="module") +def model_paths() -> tuple[np.ndarray, ...]: + return ewma_model_paths(n_paths=1_000, n_test=250, seed=2026) + + +def test_ewma_forecasts_recover_the_model_volatility_and_shocks( + model_paths: tuple[np.ndarray, ...], +) -> None: + paths, sigma, shocks = model_paths + z = NormalDist().inv_cdf(LEVEL) + for k in range(20): + forecast = rolling_var_es(paths[k], method="ewma", level=LEVEL, window=WINDOW) + np.testing.assert_allclose(forecast.scale, sigma[k], rtol=1e-12) + # The forecast is the true conditional quantile, so the exceptions are + # exactly the days with a shock below -z: i.i.d. Bernoulli(1%). + np.testing.assert_array_equal(forecast.exceptions, shocks[k] < -z) + + +def test_size_matches_the_exact_null_rejection_rates(model_paths: tuple[np.ndarray, ...]) -> None: + paths, _, _ = model_paths + rejected = np.zeros(4) + for k, path in enumerate(paths): + forecast = rolling_var_es(path, method="ewma", level=LEVEL, window=WINDOW) + result = backtest_var(forecast.exceptions, level=LEVEL) + # 199 simulations make 5% an exact Monte Carlo cutoff: 0.05 * (199 + 1) = 10. + z2 = es_backtest(forecast, n_simulations=199, seed=k) + rejected += [ + result.unconditional.p_value < 0.05, + result.independence.p_value < 0.05, + result.conditional.p_value < 0.05, + z2.p_value <= 0.05, + ] + rates = rejected / len(paths) + exact = null_rejection_rates(250, level=LEVEL) + for measured, expected in zip( + rates, (exact.unconditional, exact.independence, exact.conditional, 0.05), strict=True + ): + assert abs(measured - expected) < 4 * standard_error(expected, len(paths)) + + +def gaussian_rejection_rates(paths: np.ndarray, seed: int) -> np.ndarray: + """Share of paths whose Gaussian VaR/ES forecasts fail LR_uc, LR_cc and Z2 at 5%.""" + rejected = np.zeros(3) + for k, path in enumerate(paths): + forecast = rolling_var_es(path, method="gaussian", level=LEVEL, window=WINDOW) + result = backtest_var(forecast.exceptions, level=LEVEL) + z2 = es_backtest(forecast, n_simulations=199, seed=seed + k) + rejected += [ + result.unconditional.p_value < 0.05, + result.conditional.p_value < 0.05, + z2.p_value <= 0.05, + ] + rates: np.ndarray = rejected / len(paths) + return rates + + +def test_gaussian_var_on_student_t_returns_is_rejected_far_more_often() -> None: + rng = np.random.default_rng(7) + shape = (200, WINDOW + 1_000) + normal = 0.01 * rng.standard_normal(shape) + fat = 0.01 * rng.standard_t(4, shape) * math.sqrt(2 / 4) # unit variance before scaling + size = gaussian_rejection_rates(normal, seed=0) + power = gaussian_rejection_rates(fat, seed=0) + # examples/var_size_power.py (1,000 paths) measures rejection rates of 2.2% (LR_uc), + # 0.9% (LR_cc) and 3.9% (Z2) on normal returns, and 46.9%, 31.3% and 87.3% on + # Student-t returns. Here (200 paths) these bounds leave at least four standard errors. + assert (size < 0.12).all() + assert power[0] > 0.30 + assert power[1] > 0.15 + assert power[2] > 0.65 + assert power[2] > power[0] # Z2 sees the size of tail losses, not only their count