diff --git a/tests/backtests/test_option_pricer_backtest.py b/tests/backtests/test_option_pricer_backtest.py new file mode 100644 index 0000000..483b84d --- /dev/null +++ b/tests/backtests/test_option_pricer_backtest.py @@ -0,0 +1,54 @@ +"""Unit tests for the IV solver and the option pricer backtest (no network).""" + +import pytest + +from trading_system.backtests.options_simulation import bs_call_implied_vol, bs_call_price +from trading_system.backtests.option_pricer_backtest import ( + occ_symbol, + realized_vol, + run_pricer_backtest, + summarize, + year_fraction, +) + + +@pytest.mark.parametrize("sigma", [0.2, 0.6, 1.2]) +def test_implied_vol_round_trip(sigma): + p = bs_call_price(170, 230, 31 / 365, 0.04, sigma) + assert bs_call_implied_vol(p, 170, 230, 31 / 365, 0.04) == pytest.approx(sigma, abs=1e-4) + + +def test_implied_vol_rejects_out_of_band(): + assert bs_call_implied_vol(0.0, 170, 230, 0.1, 0.04) is None + assert bs_call_implied_vol(200.0, 170, 230, 0.1, 0.04) is None # above spot + assert bs_call_implied_vol(1.0, 170, 230, 0.0, 0.04) is None # expired + + +def test_occ_symbol(): + assert occ_symbol("NBIS", "2026-11-20", 230) == "NBIS261120C00230000" + + +def test_year_fraction(): + assert year_fraction("2026-10-20", "2026-11-20") == pytest.approx(31 / 365) + assert year_fraction("2026-12-01", "2026-11-20") == 0.0 + + +def test_realized_vol_flat_and_short(): + assert realized_vol([100.0] * 25, 20) == 0.0 + assert realized_vol([100.0] * 5, 20) is None + + +def test_constant_iv_market_is_priced_exactly(): + """If the market prices at one fixed sigma, entry_iv and prev_iv are exact.""" + sigma, K, exp = 0.9, 230.0, "2026-11-20" + days = [f"2026-10-{d:02d}" for d in range(1, 11)] + spots = [200 - 3 * i for i in range(len(days))] + under = dict(zip(days, spots)) + opt = {d: bs_call_price(s, K, year_fraction(d, exp), 0.04, sigma) for d, s in under.items()} + + rows = run_pricer_backtest(under, opt, K, exp, vol_window=3) + stats = summarize(rows) + assert rows[0]["prev_iv"] is None # no prior day yet + assert stats["entry_iv"]["mae"] == pytest.approx(0, abs=1e-3) + assert stats["prev_iv"]["mae"] == pytest.approx(0, abs=1e-3) + assert stats["prev_iv"]["n"] == len(days) - 1 diff --git a/tests/backtests/test_option_sensitivity.py b/tests/backtests/test_option_sensitivity.py new file mode 100644 index 0000000..bfd50e1 --- /dev/null +++ b/tests/backtests/test_option_sensitivity.py @@ -0,0 +1,49 @@ +"""Greeks are checked against finite differences of bs_call_price.""" + +import pytest + +from trading_system.backtests.options_simulation import bs_call_price +from trading_system.backtests.option_sensitivity import ( + breakeven_spot, + bs_call_greeks, + scenario_grid, +) + +S, K, T, R, V = 180.0, 230.0, 30 / 365, 0.04, 0.93 + + +def test_greeks_match_finite_differences(): + g = bs_call_greeks(S, K, T, R, V) + h = 0.01 + p = lambda s=S, t=T, r=R, v=V: bs_call_price(s, K, t, r, v) + assert g["price"] == pytest.approx(p()) + assert g["delta"] == pytest.approx((p(s=S + h) - p(s=S - h)) / (2 * h), rel=1e-4) + assert g["gamma"] == pytest.approx((p(s=S + h) - 2 * p() + p(s=S - h)) / h ** 2, rel=1e-2) + assert g["vega"] == pytest.approx((p(v=V + 1e-4) - p(v=V - 1e-4)) / 2e-4 / 100, rel=1e-4) + assert g["theta"] == pytest.approx(p(t=T - 1 / 365) - p(), rel=2e-2) + assert g["rho"] == pytest.approx((p(r=R + 1e-4) - p(r=R - 1e-4)) / 2e-4 / 100, rel=1e-4) + + +def test_expired_greeks(): + assert bs_call_greeks(240, K, 0, R, V)["delta"] == 1.0 + assert bs_call_greeks(200, K, 0, R, V)["price"] == 0.0 + + +def test_scenario_grid_center_and_monotone(): + grid = scenario_grid(S, K, 30, R, V, [-0.1, 0.0, 0.1], [-0.1, 0.0, 0.1]) + assert grid[1][1] == pytest.approx(bs_call_price(S, K, T, R, V)) + assert grid[0][1] < grid[1][1] < grid[2][1] # higher spot -> higher call + assert grid[1][0] < grid[1][1] < grid[1][2] # higher IV -> higher call + later = scenario_grid(S, K, 30, R, V, [0.0], [0.0], days_forward=7) + assert later[0][0] < grid[1][1] # time decay + + +def test_breakeven(): + assert breakeven_spot(230, 35.78) == pytest.approx(265.78) + + +@pytest.mark.parametrize("x", [-3, -2, -1, -0.5, 0, 0.5, 1, 2, 3]) +def test_norm_cdf_matches_erf(x): + import math + from trading_system.backtests.options_simulation import _norm_cdf + assert _norm_cdf(x) == pytest.approx(0.5 * (1 + math.erf(x / math.sqrt(2))), abs=1e-6) diff --git a/trading_system/backtests/option_pricer_backtest.py b/trading_system/backtests/option_pricer_backtest.py new file mode 100644 index 0000000..2946bc0 --- /dev/null +++ b/trading_system/backtests/option_pricer_backtest.py @@ -0,0 +1,170 @@ +""" +Backtest simple Black-Scholes pricing models against a real option's daily closes. + +Models (all price the call at each day's underlying close): + * realized — sigma = trailing N-day realized vol of the underlying + * entry_iv — sigma frozen at the IV implied on the first option bar + * prev_iv — sigma = IV implied by the previous day's option close + (one-day-ahead: uses only information available yesterday) + +Usage: + python -m trading_system.backtests.option_pricer_backtest + python -m trading_system.backtests.option_pricer_backtest \ + --symbol NBIS --strike 230 --expiry 2026-11-20 --start 2026-09-29 +""" + +import argparse +import http.client +import json +import math +import time +import urllib.error +import urllib.request +from datetime import date, datetime, timezone +from typing import Dict, List, Optional + +from trading_system.backtests.options_simulation import ( + bs_call_implied_vol, + bs_call_price, +) + +MODELS = ("realized", "entry_iv", "prev_iv") + + +def year_fraction(day: str, expiry: str) -> float: + """Calendar-time to expiry in years, close-to-close.""" + d0 = date.fromisoformat(day) + d1 = date.fromisoformat(expiry) + return max((d1 - d0).days, 0) / 365.0 + + +def realized_vol(closes: List[float], window: int) -> Optional[float]: + """Annualized close-to-close vol over the last `window` returns.""" + if len(closes) < window + 1: + return None + tail = closes[-(window + 1):] + rets = [math.log(b / a) for a, b in zip(tail, tail[1:])] + mean = sum(rets) / len(rets) + var = sum((x - mean) ** 2 for x in rets) / (len(rets) - 1) + return math.sqrt(var * 252) + + +def run_pricer_backtest( + underlying: Dict[str, float], + option: Dict[str, float], + strike: float, + expiry: str, + r: float = 0.04, + vol_window: int = 20, + start: Optional[str] = None, +) -> List[dict]: + """ + Price the call on every day that has both an underlying and option close. + + `underlying` and `option` map ISO date -> close. Underlying history before + `start` is used only to seed realized vol. + """ + days = sorted(d for d in option if d in underlying and (start is None or d >= start)) + u_days = sorted(underlying) + rows: List[dict] = [] + entry_iv: Optional[float] = None + prev_iv: Optional[float] = None + + for d in days: + S = underlying[d] + actual = option[d] + T = year_fraction(d, expiry) + hist = [underlying[x] for x in u_days if x <= d] + rv = realized_vol(hist, vol_window) + iv_today = bs_call_implied_vol(actual, S, strike, T, r) + if entry_iv is None: + entry_iv = iv_today + + row = {"date": d, "spot": S, "actual": actual, "iv": iv_today, "rv": rv} + for name, sigma in (("realized", rv), ("entry_iv", entry_iv), ("prev_iv", prev_iv)): + row[name] = bs_call_price(S, strike, T, r, sigma) if sigma else None + rows.append(row) + if iv_today is not None: + prev_iv = iv_today + return rows + + +def summarize(rows: List[dict]) -> Dict[str, dict]: + """MAE ($), mean signed error ($), and mean absolute % error per model.""" + out = {} + for m in MODELS: + pairs = [(r[m], r["actual"]) for r in rows if r[m] is not None] + if not pairs: + continue + errs = [p - a for p, a in pairs] + out[m] = { + "n": len(pairs), + "mae": sum(abs(e) for e in errs) / len(errs), + "bias": sum(errs) / len(errs), + "mape": 100 * sum(abs(p - a) / a for p, a in pairs) / len(pairs), + } + return out + + +# ── Data fetch (Yahoo public chart endpoint) ───────────────────────────── + + +def fetch_daily_closes(symbol: str, range_: str = "6mo") -> Dict[str, float]: + url = ( + f"https://query1.finance.yahoo.com/v8/finance/chart/{symbol}" + f"?interval=1d&range={range_}" + ) + req = urllib.request.Request(url, headers={"User-Agent": "Mozilla/5.0"}) + for attempt in range(5): + try: + with urllib.request.urlopen(req, timeout=20) as resp: + res = json.load(resp)["chart"]["result"][0] + break + except (http.client.IncompleteRead, urllib.error.URLError): + if attempt == 4: + raise + time.sleep(1 + attempt) + out = {} + for ts, c in zip(res.get("timestamp", []), res["indicators"]["quote"][0]["close"]): + if c is not None: + out[datetime.fromtimestamp(ts, tz=timezone.utc).date().isoformat()] = c + return out + + +def occ_symbol(symbol: str, expiry: str, strike: float, right: str = "C") -> str: + return f"{symbol}{date.fromisoformat(expiry):%y%m%d}{right}{int(round(strike * 1000)):08d}" + + +def main(argv=None): + ap = argparse.ArgumentParser(description=__doc__.splitlines()[1]) + ap.add_argument("--symbol", default="NBIS") + ap.add_argument("--strike", type=float, default=230.0) + ap.add_argument("--expiry", default="2026-11-20") + ap.add_argument("--start", default="2026-09-29") + ap.add_argument("--rate", type=float, default=0.04) + ap.add_argument("--vol-window", type=int, default=20) + args = ap.parse_args(argv) + + contract = occ_symbol(args.symbol, args.expiry, args.strike) + rows = run_pricer_backtest( + fetch_daily_closes(args.symbol), + fetch_daily_closes(contract), + args.strike, args.expiry, args.rate, args.vol_window, args.start, + ) + + fmt = lambda v, w=7: f"{v:{w}.2f}" if v is not None else " " * (w - 1) + "-" + print(f"{contract} r={args.rate} rv_window={args.vol_window}") + print(f"{'date':10} {'spot':>7} {'actual':>7} {'iv%':>6} {'rv%':>6} " + f"{'realized':>8} {'entry_iv':>8} {'prev_iv':>8}") + for r in rows: + iv = r["iv"] * 100 if r["iv"] else None + rv = r["rv"] * 100 if r["rv"] else None + print(f"{r['date']:10} {fmt(r['spot'])} {fmt(r['actual'])} {fmt(iv, 6)} {fmt(rv, 6)} " + f"{fmt(r['realized'], 8)} {fmt(r['entry_iv'], 8)} {fmt(r['prev_iv'], 8)}") + print() + for m, s in summarize(rows).items(): + print(f"{m:9} n={s['n']:3} MAE ${s['mae']:.2f} bias ${s['bias']:+.2f} MAPE {s['mape']:.1f}%") + + +if __name__ == "__main__": + main() diff --git a/trading_system/backtests/option_sensitivity.py b/trading_system/backtests/option_sensitivity.py new file mode 100644 index 0000000..7437d8c --- /dev/null +++ b/trading_system/backtests/option_sensitivity.py @@ -0,0 +1,119 @@ +""" +Sensitivity model for a European call: Greeks plus scenario grids. + +Greeks (per-unit conventions, matching how brokers display them): + * delta — $ change per $1 move in the underlying + * gamma — change in delta per $1 move in the underlying + * vega — $ change per 1 vol point (0.01 sigma) + * theta — $ change per calendar day of time decay + * rho — $ change per 1% (0.01) move in rates + +Scenario grid: reprices the call across spot moves x IV shifts x days forward, +so you can read off e.g. "NBIS +5%, IV -10 pts, 3 days later". + +Usage: + python -m trading_system.backtests.option_sensitivity \ + --spot 180 --price 3.80 --ref-spot 171.45 --days 30 +""" + +import argparse +import math +from typing import Dict, List, Optional, Sequence + +from trading_system.backtests.options_simulation import ( + _norm_cdf, + bs_call_implied_vol, + bs_call_price, +) + + +def _norm_pdf(x: float) -> float: + return math.exp(-0.5 * x * x) / math.sqrt(2 * math.pi) + + +def bs_call_greeks(S: float, K: float, T: float, r: float, sigma: float) -> Dict[str, float]: + """Analytic Black-Scholes call Greeks. At/after expiry returns intrinsic-only values.""" + if T <= 0 or sigma <= 0: + itm = 1.0 if S > K else 0.0 + return {"price": max(S - K, 0.0), "delta": itm, "gamma": 0.0, + "vega": 0.0, "theta": 0.0, "rho": 0.0} + sq = math.sqrt(T) + d1 = (math.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * sq) + d2 = d1 - sigma * sq + disc = math.exp(-r * T) + theta_yr = -S * _norm_pdf(d1) * sigma / (2 * sq) - r * K * disc * _norm_cdf(d2) + return { + "price": S * _norm_cdf(d1) - K * disc * _norm_cdf(d2), + "delta": _norm_cdf(d1), + "gamma": _norm_pdf(d1) / (S * sigma * sq), + "vega": S * _norm_pdf(d1) * sq / 100.0, + "theta": theta_yr / 365.0, + "rho": K * T * disc * _norm_cdf(d2) / 100.0, + } + + +def scenario_grid( + S: float, K: float, days: int, r: float, sigma: float, + spot_moves: Sequence[float], iv_shifts: Sequence[float], days_forward: int = 0, +) -> List[List[float]]: + """ + Call prices on a grid: rows = spot moves (fractions, e.g. 0.05), + cols = IV shifts (absolute vol, e.g. -0.10), `days_forward` days later. + """ + T = max(days - days_forward, 0) / 365.0 + return [ + [bs_call_price(S * (1 + m), K, T, r, max(sigma + s, 0.01)) for s in iv_shifts] + for m in spot_moves + ] + + +def breakeven_spot(K: float, premium: float) -> float: + """Underlying price at expiry where the long call breaks even.""" + return K + premium + + +def main(argv=None): + ap = argparse.ArgumentParser(description="Call option Greeks + scenario grid") + ap.add_argument("--spot", type=float, default=180.0, help="underlying to center the grid on") + ap.add_argument("--strike", type=float, default=230.0) + ap.add_argument("--days", type=int, default=30, help="calendar days to expiry") + ap.add_argument("--rate", type=float, default=0.04) + ap.add_argument("--iv", type=float, help="vol (e.g. 0.93); else implied from --price/--ref-spot") + ap.add_argument("--price", type=float, default=3.80, help="last option price, to imply IV") + ap.add_argument("--ref-spot", type=float, default=171.45, help="underlying when --price printed") + ap.add_argument("--ref-days", type=int, default=31, help="days to expiry when --price printed") + ap.add_argument("--cost", type=float, default=35.78, help="entry premium per share") + ap.add_argument("--contracts", type=int, default=1) + ap.add_argument("--days-forward", type=int, nargs="+", default=[0, 7, 14]) + args = ap.parse_args(argv) + + sigma: Optional[float] = args.iv or bs_call_implied_vol( + args.price, args.ref_spot, args.strike, args.ref_days / 365.0, args.rate) + if sigma is None: + raise SystemExit("could not imply vol from --price; pass --iv") + mult = 100 * args.contracts + T = args.days / 365.0 + g = bs_call_greeks(args.spot, args.strike, T, args.rate, sigma) + + print(f"K={args.strike:g} spot={args.spot:g} {args.days}d IV={sigma*100:.1f}% r={args.rate}") + print(f"price ${g['price']:.2f} (position ${g['price']*mult:,.0f}, " + f"P&L ${(g['price']-args.cost)*mult:+,.0f} vs ${args.cost:.2f} cost)") + print(f"delta {g['delta']:.3f} gamma {g['gamma']:.4f} vega ${g['vega']:.3f}/pt " + f"theta ${g['theta']:.3f}/day rho ${g['rho']:.3f}/1%") + print(f"position: ${g['delta']*mult:,.0f} per $1 NBIS, ${g['vega']*mult:,.0f} per vol pt, " + f"${g['theta']*mult:,.0f}/day decay") + print(f"breakeven at expiry ${breakeven_spot(args.strike, args.cost):.2f} " + f"({(breakeven_spot(args.strike, args.cost)/args.spot-1)*100:+.1f}% from spot)") + + moves = [-0.15, -0.10, -0.05, 0.0, 0.05, 0.10, 0.20, 0.30] + shifts = [-0.30, -0.20, -0.10, 0.0, 0.10] + for fwd in args.days_forward: + grid = scenario_grid(args.spot, args.strike, args.days, args.rate, sigma, moves, shifts, fwd) + print(f"\n+{fwd}d rows=spot move, cols=IV shift (pts)") + print(f"{'':>14}" + "".join(f"{s*100:+7.0f}" for s in shifts)) + for m, row in zip(moves, grid): + print(f"{m*100:+4.0f}% {args.spot*(1+m):7.2f} " + "".join(f"{p:7.2f}" for p in row)) + + +if __name__ == "__main__": + main() diff --git a/trading_system/backtests/options_simulation.py b/trading_system/backtests/options_simulation.py index 2fa97bd..281500a 100644 --- a/trading_system/backtests/options_simulation.py +++ b/trading_system/backtests/options_simulation.py @@ -26,10 +26,11 @@ def _norm_cdf(x: float) -> float: a5 = 1.061405429 p = 0.3275911 + # 7.1.26 approximates erf(z); Phi(x) = 0.5 * (1 + erf(x / sqrt(2))). sign = 1.0 if x >= 0 else -1.0 - x = abs(x) - t = 1.0 / (1.0 + p * x) - y = 1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * math.exp(-x * x / 2.0) + z = abs(x) / math.sqrt(2.0) + t = 1.0 / (1.0 + p * z) + y = 1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * math.exp(-z * z) return 0.5 * (1.0 + sign * y) @@ -51,6 +52,34 @@ def bs_call_price(S: float, K: float, T: float, r: float, sigma: float) -> float return S * _norm_cdf(d1) - K * math.exp(-r * T) * _norm_cdf(d2) +def bs_call_implied_vol( + price: float, S: float, K: float, T: float, r: float, + lo: float = 0.01, hi: float = 5.0, tol: float = 1e-6, +) -> Optional[float]: + """ + Implied volatility of a European call by bisection. + + Returns None when the price is outside the no-arbitrage band + [intrinsic, S] or T <= 0 — there is no sigma that reproduces it. + """ + if T <= 0 or price <= 0: + return None + intrinsic = max(S - K * math.exp(-r * T), 0.0) + if price <= intrinsic or price >= S: + return None + if bs_call_price(S, K, T, r, hi) < price: + return None + for _ in range(200): + mid = 0.5 * (lo + hi) + if bs_call_price(S, K, T, r, mid) < price: + lo = mid + else: + hi = mid + if hi - lo < tol: + break + return 0.5 * (lo + hi) + + # ── Data structures ──────────────────────────────────────────────────────