Fix _norm_cdf scaling; add option pricer backtest + sensitivity model

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

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)

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

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