#!/usr/bin/env python3
"""
Research backtests: **new** option strategies from academic asset-pricing ideas,
using **raw** Theta 15:45 SPY chains + Yahoo SPY/VIX/VVIX — **not** VRP trade logs.

Execution: conservative bid/ask (short: bid in / ask out; long: ask in / bid out).
Sequential positions (no overlap). Sharpe on daily returns of an equity curve that
steps only on realized exit PnL (zeros on other days).

Run::

    cd /Users/robzingale/trading_bot && .venv/bin/python RenTech/strategy_stack/research_literature_theta_strategies.py
"""
from __future__ import annotations

import argparse
import json
import math
import sys
import time
from collections import defaultdict
from pathlib import Path
from typing import Any, Callable, Sequence

import numpy as np
import pandas as pd

_REPO_ROOT = Path(__file__).resolve().parents[2]
if str(_REPO_ROOT) not in sys.path:
    sys.path.insert(0, str(_REPO_ROOT))

from RenTech.core.options_backtest import find_contract_in_chain
from RenTech.core.options_data_loader import OptionChain, OptionContract
from RenTech.core.theta_chunks_loader import ThetaChunksLoader, theta_chunks_date_bounds
from RenTech.strategy_stack.vrp_backtester import (
    load_spy_vix_from_yfinance,
    normalize_spy_df,
    trading_days_intersecting_spy,
)

MULT = 100.0
_DEFAULT_THETA = _REPO_ROOT / "RenTech" / "data" / "theta_chunks"


def _norm(d: pd.Timestamp) -> pd.Timestamp:
    return pd.Timestamp(d).normalize()


def _short_open_px(c: OptionContract) -> float:
    return float(c.bid)


def _short_close_px(c: OptionContract) -> float:
    return float(c.ask)


def _long_open_px(c: OptionContract) -> float:
    return float(c.ask)


def _long_close_px(c: OptionContract) -> float:
    return float(c.bid)


def _settle_short(c: OptionContract, spy_close: float) -> float:
    k, ot = float(c.strike), str(c.option_type).upper()
    if ot == "C":
        return max(spy_close - k, 0.0)
    return max(k - spy_close, 0.0)


def _settle_long(c: OptionContract, spy_close: float) -> float:
    if str(c.option_type).upper() == "C":
        return max(spy_close - float(c.strike), 0.0)
    return max(float(c.strike) - spy_close, 0.0)


def sharpe_daily_returns(daily_ret: pd.Series) -> float:
    r = daily_ret.dropna().astype(float)
    if len(r) < 50:
        return float("nan")
    sd = float(r.std(ddof=1))
    if sd < 1e-12:
        return float("nan")
    return float(r.mean() / sd) * math.sqrt(252.0)


def daily_pnl_series(
    exit_dates: Sequence[pd.Timestamp],
    pnls: Sequence[float],
    all_days: Sequence[pd.Timestamp],
) -> pd.Series:
    """Per-session realized PnL (zeros on non-exit days), aligned to ``all_days``."""
    idx = pd.DatetimeIndex([_norm(d) for d in all_days])
    s = pd.Series(0.0, index=idx)
    for d, p in zip(exit_dates, pnls, strict=False):
        s.loc[_norm(d)] += float(p)
    return s


def equity_curve_from_realized(
    exit_dates: Sequence[pd.Timestamp],
    pnls: Sequence[float],
    all_days: Sequence[pd.Timestamp],
    capital: float,
) -> tuple[pd.Series, float]:
    idx = pd.DatetimeIndex([_norm(d) for d in all_days])
    s = daily_pnl_series(exit_dates, pnls, all_days)
    eq = float(capital) + s.cumsum()
    eq = pd.Series(eq.values, index=idx)
    r = eq.pct_change().replace([np.inf, -np.inf], np.nan).fillna(0.0)
    return eq, sharpe_daily_returns(r)


def find_atm_straddle(
    chain: OptionChain, spot: float, target_dte: int
) -> tuple[OptionContract, OptionContract] | None:
    """Same expiry (|DTE−target| minimal), strike closest to spot with both C and P."""
    if not chain.contracts:
        return None
    as_of = _norm(chain.as_of)
    strikes_c: dict[pd.Timestamp, set[float]] = defaultdict(set)
    strikes_p: dict[pd.Timestamp, set[float]] = defaultdict(set)
    for c in chain.contracts:
        exp = _norm(c.expiration)
        dte = int((exp - as_of).days)
        if dte < 5:
            continue
        if c.option_type == "C":
            strikes_c[exp].add(float(c.strike))
        elif c.option_type == "P":
            strikes_p[exp].add(float(c.strike))
    best: tuple[int, float, pd.Timestamp, float] | None = None
    for exp in strikes_c.keys() & strikes_p.keys():
        dte = int((exp - as_of).days)
        d_err = abs(dte - int(target_dte))
        for strike in strikes_c[exp] & strikes_p[exp]:
            sk_err = abs(float(strike) - float(spot))
            tup = (d_err, sk_err, exp, float(strike))
            if best is None or tup < best:
                best = tup
    if best is None:
        return None
    _, _, exp, strike = best
    call = find_contract_in_chain(chain, exp, strike, "C")
    put = find_contract_in_chain(chain, exp, strike, "P")
    if call is None or put is None:
        return None
    return call, put


def atm_iv_straddle(chain: OptionChain, spot: float, target_dte: int) -> float | None:
    pair = find_atm_straddle(chain, spot, target_dte)
    if pair is None:
        return None
    c, p = pair
    if not (math.isfinite(c.iv) and math.isfinite(p.iv) and c.iv > 0 and p.iv > 0):
        return None
    return 0.5 * (float(c.iv) + float(p.iv))


def find_target_leg_safe(
    chain: OptionChain, target_dte: int, target_delta: float, option_type: str
) -> OptionContract | None:
    try:
        return chain.find_target_leg(
            target_dte=target_dte, target_delta=float(target_delta), option_type=option_type
        )
    except ValueError:
        return None


def short_straddle_pnl(
    entry_chain: OptionChain,
    exit_chain: OptionChain,
    entry_spy: float,
    exit_spy: float,
    target_dte: int,
) -> float | None:
    pair_e = find_atm_straddle(entry_chain, entry_spy, target_dte)
    if pair_e is None:
        return None
    ce, pe = pair_e
    exp = _norm(ce.expiration)
    strike = float(ce.strike)
    credit = (_short_open_px(ce) + _short_open_px(pe)) * MULT
    c_x = find_contract_in_chain(exit_chain, exp, strike, "C")
    p_x = find_contract_in_chain(exit_chain, exp, strike, "P")
    exit_d = _norm(exit_chain.as_of)
    if c_x is not None and p_x is not None:
        debit = (_short_close_px(c_x) + _short_close_px(p_x)) * MULT
        return float(credit - debit)
    if exit_d >= exp:
        settle = (_settle_short(ce, exit_spy) + _settle_short(pe, exit_spy)) * MULT
        return float(credit - settle)
    return None


def long_straddle_pnl(
    entry_chain: OptionChain,
    exit_chain: OptionChain,
    entry_spy: float,
    exit_spy: float,
    target_dte: int,
) -> float | None:
    pair_e = find_atm_straddle(entry_chain, entry_spy, target_dte)
    if pair_e is None:
        return None
    ce, pe = pair_e
    exp = _norm(ce.expiration)
    strike = float(ce.strike)
    cost = (_long_open_px(ce) + _long_open_px(pe)) * MULT
    c_x = find_contract_in_chain(exit_chain, exp, strike, "C")
    p_x = find_contract_in_chain(exit_chain, exp, strike, "P")
    exit_d = _norm(exit_chain.as_of)
    if c_x is not None and p_x is not None:
        val = (_long_close_px(c_x) + _long_close_px(p_x)) * MULT
        return float(val - cost)
    if exit_d >= exp:
        val = (_settle_long(ce, exit_spy) + _settle_long(pe, exit_spy)) * MULT
        return float(val - cost)
    return None


def short_put_pnl(
    entry_chain: OptionChain,
    exit_chain: OptionChain,
    exit_spy: float,
    target_dte: int,
    target_delta: float,
) -> float | None:
    p0 = find_target_leg_safe(entry_chain, target_dte, target_delta, "P")
    if p0 is None:
        return None
    exp = _norm(p0.expiration)
    strike = float(p0.strike)
    credit = _short_open_px(p0) * MULT
    p_x = find_contract_in_chain(exit_chain, exp, strike, "P")
    exit_d = _norm(exit_chain.as_of)
    if p_x is not None:
        return float(credit - _short_close_px(p_x) * MULT)
    if exit_d >= exp:
        return float(credit - _settle_short(p0, exit_spy) * MULT)
    return None


def short_strangle_pnl(
    entry_chain: OptionChain,
    exit_chain: OptionChain,
    exit_spy: float,
    target_dte: int,
    put_delta: float,
    call_delta: float,
) -> float | None:
    p0 = find_target_leg_safe(entry_chain, target_dte, put_delta, "P")
    c0 = find_target_leg_safe(entry_chain, target_dte, call_delta, "C")
    if p0 is None or c0 is None:
        return None
    if _norm(p0.expiration) != _norm(c0.expiration):
        return None
    exp = _norm(p0.expiration)
    credit = (_short_open_px(p0) + _short_open_px(c0)) * MULT
    p_x = find_contract_in_chain(exit_chain, exp, float(p0.strike), "P")
    c_x = find_contract_in_chain(exit_chain, exp, float(c0.strike), "C")
    exit_d = _norm(exit_chain.as_of)
    if p_x is not None and c_x is not None:
        return float(credit - (_short_close_px(p_x) + _short_close_px(c_x)) * MULT)
    if exit_d >= exp:
        return float(credit - (_settle_short(p0, exit_spy) + _settle_short(c0, exit_spy)) * MULT)
    return None


def short_risk_reversal_pnl(
    entry_chain: OptionChain,
    exit_chain: OptionChain,
    exit_spy: float,
    target_dte: int,
    put_delta: float,
    call_delta: float,
) -> float | None:
    p0 = find_target_leg_safe(entry_chain, target_dte, put_delta, "P")
    c0 = find_target_leg_safe(entry_chain, target_dte, call_delta, "C")
    if p0 is None or c0 is None or _norm(c0.expiration) != _norm(p0.expiration):
        return None
    exp = _norm(p0.expiration)
    premium = (_short_open_px(p0) - _long_open_px(c0)) * MULT
    p_x = find_contract_in_chain(exit_chain, exp, float(p0.strike), "P")
    c_x = find_contract_in_chain(exit_chain, exp, float(c0.strike), "C")
    exit_d = _norm(exit_chain.as_of)
    if p_x is not None and c_x is not None:
        return float(premium - (_short_close_px(p_x) - _long_close_px(c_x)) * MULT)
    if exit_d >= exp:
        return float(premium - (_settle_short(p0, exit_spy) - _settle_long(c0, exit_spy)) * MULT)
    return None


def short_vertical_put_pnl(
    entry_chain: OptionChain,
    exit_chain: OptionChain,
    exit_spy: float,
    target_dte: int,
    short_delta: float,
    wing_width: float,
) -> float | None:
    ps = find_target_leg_safe(entry_chain, target_dte, short_delta, "P")
    if ps is None:
        return None
    exp = _norm(ps.expiration)
    long_strike = float(ps.strike) - float(wing_width)
    pl = find_contract_in_chain(entry_chain, exp, long_strike, "P")
    if pl is None:
        return None
    credit = (_short_open_px(ps) - _long_open_px(pl)) * MULT
    p_s_x = find_contract_in_chain(exit_chain, exp, float(ps.strike), "P")
    p_l_x = find_contract_in_chain(exit_chain, exp, long_strike, "P")
    exit_d = _norm(exit_chain.as_of)
    if p_s_x is not None and p_l_x is not None:
        return float(credit - (_short_close_px(p_s_x) - _long_close_px(p_l_x)) * MULT)
    if exit_d >= exp:
        return float(credit - (_settle_short(ps, exit_spy) - _settle_long(pl, exit_spy)) * MULT)
    return None


def short_vertical_call_pnl(
    entry_chain: OptionChain,
    exit_chain: OptionChain,
    exit_spy: float,
    target_dte: int,
    short_call_delta: float,
    wing_width: float,
) -> float | None:
    """Bear call credit spread: short lower-strike call, long higher-strike call (same expiry)."""
    cs = find_target_leg_safe(entry_chain, target_dte, short_call_delta, "C")
    if cs is None:
        return None
    exp = _norm(cs.expiration)
    long_k = float(cs.strike) + float(wing_width)
    cl = find_contract_in_chain(entry_chain, exp, long_k, "C")
    if cl is None:
        return None
    credit = (_short_open_px(cs) - _long_open_px(cl)) * MULT
    c_s_x = find_contract_in_chain(exit_chain, exp, float(cs.strike), "C")
    c_l_x = find_contract_in_chain(exit_chain, exp, long_k, "C")
    exit_d = _norm(exit_chain.as_of)
    if c_s_x is not None and c_l_x is not None:
        return float(credit - (_short_close_px(c_s_x) - _long_close_px(c_l_x)) * MULT)
    if exit_d >= exp:
        return float(credit - (_settle_short(cs, exit_spy) - _settle_long(cl, exit_spy)) * MULT)
    return None


SignalFn = Callable[[int, pd.Series, OptionChain, float], bool]
TradeFn = Callable[[OptionChain, OptionChain, float, float, tuple], float | None]


def run_signal_backtest(
    days: list[pd.Timestamp],
    get_chain: Callable[[pd.Timestamp], OptionChain],
    panel: pd.DataFrame,
    signal: SignalFn,
    hold: int,
    trade_fn: TradeFn,
    trade_params: tuple,
) -> tuple[list[pd.Timestamp], list[float], int]:
    exits: list[pd.Timestamp] = []
    pnls: list[float] = []
    pending: tuple[int, float, tuple] | None = None  # entry_idx, entry_spy, params

    for i, d in enumerate(days):
        row = panel.reindex([_norm(d)]).iloc[0]
        spy_c = float(row["close"])
        ch = get_chain(d)
        if pending is not None:
            ent_i, spy_e, tp = pending
            if i >= ent_i + hold:
                ch_e = get_chain(days[ent_i])
                spy_x = float(panel.reindex([_norm(d)]).iloc[0]["close"])
                pnl = trade_fn(ch_e, ch, spy_e, spy_x, tp)
                pending = None
                if pnl is not None and math.isfinite(pnl):
                    exits.append(_norm(d))
                    pnls.append(float(pnl))
        if pending is not None:
            continue
        if not ch.contracts:
            continue
        if signal(i, row, ch, spy_c):
            pending = (i, spy_c, trade_params)

    return exits, pnls, len(pnls)


def run_signal_backtest_stack_while_signal(
    days: list[pd.Timestamp],
    get_chain: Callable[[pd.Timestamp], OptionChain],
    panel: pd.DataFrame,
    signal: SignalFn,
    hold: int,
    trade_fn: TradeFn,
    trade_params: tuple,
) -> tuple[list[pd.Timestamp], list[float], int, int]:
    """
    Open **one new position every session** the signal is true (chains non-empty), while
    allowing unlimited overlap. Each leg uses the same ``trade_fn`` / ``trade_params`` and
    exits ``hold`` **session indices** after its entry index (same convention as
    :func:`run_signal_backtest`).

    A new leg is only opened if ``entry_i + hold < len(days)`` so its exit session exists
    in ``days`` (mirrors the original engine's inability to close beyond the sample end).

    Returns ``(exit_dates, pnls, n_trades, max_concurrent_open)``.
    """
    exits: list[pd.Timestamp] = []
    pnls: list[float] = []
    # (entry_session_index, entry_spy_close, trade_params) — same trade_params for all legs here
    open_legs: list[tuple[int, float, tuple]] = []
    max_open = 0
    n_days = len(days)

    for i, d in enumerate(days):
        row = panel.reindex([_norm(d)]).iloc[0]
        spy_c = float(row["close"])
        ch = get_chain(d)

        max_open = max(max_open, len(open_legs))

        still: list[tuple[int, float, tuple]] = []
        for ent_i, spy_e, tp in open_legs:
            if i >= ent_i + hold:
                ch_e = get_chain(days[ent_i])
                spy_x = float(panel.reindex([_norm(d)]).iloc[0]["close"])
                pnl = trade_fn(ch_e, ch, spy_e, spy_x, tp)
                if pnl is not None and math.isfinite(pnl):
                    exits.append(_norm(d))
                    pnls.append(float(pnl))
            else:
                still.append((ent_i, spy_e, tp))
        open_legs = still

        if ch.contracts and signal(i, row, ch, spy_c):
            if i + hold < n_days:
                open_legs.append((i, spy_c, trade_params))
        max_open = max(max_open, len(open_legs))

    return exits, pnls, len(pnls), int(max_open)


def _exp_str(exp: pd.Timestamp) -> str:
    return _norm(exp).strftime("%Y-%m-%d")


def _short_leg_row(
    c0: OptionContract,
    c_x: OptionContract | None,
    exit_spy: float,
    exit_d: pd.Timestamp,
    exp: pd.Timestamp,
) -> dict[str, Any]:
    exp_n = _norm(exp)
    settled = c_x is None and exit_d >= exp_n
    return {
        "position": "short",
        "right": str(c0.option_type).upper()[:1],
        "strike": float(c0.strike),
        "expiry": _exp_str(exp),
        "entry_bid": float(c0.bid),
        "entry_ask": float(c0.ask),
        "exit_bid": float(c_x.bid) if c_x is not None else None,
        "exit_ask": float(c_x.ask) if c_x is not None else None,
        "exit_settled": settled,
        "exit_intrinsic_per_sh": float(_settle_short(c0, exit_spy)) if settled else None,
        "entry_iv": float(c0.iv),
        "delta_entry": float(c0.delta),
    }


def _long_leg_row(
    c0: OptionContract,
    c_x: OptionContract | None,
    exit_spy: float,
    exit_d: pd.Timestamp,
    exp: pd.Timestamp,
) -> dict[str, Any]:
    exp_n = _norm(exp)
    settled = c_x is None and exit_d >= exp_n
    return {
        "position": "long",
        "right": str(c0.option_type).upper()[:1],
        "strike": float(c0.strike),
        "expiry": _exp_str(exp),
        "entry_bid": float(c0.bid),
        "entry_ask": float(c0.ask),
        "exit_bid": float(c_x.bid) if c_x is not None else None,
        "exit_ask": float(c_x.ask) if c_x is not None else None,
        "exit_settled": settled,
        "exit_intrinsic_per_sh": float(_settle_long(c0, exit_spy)) if settled else None,
        "entry_iv": float(c0.iv),
        "delta_entry": float(c0.delta),
    }


def collect_trade_legs(
    trade_kind: str,
    trade_params: tuple,
    ch_e: OptionChain,
    ch_x: OptionChain,
    entry_spy: float,
    exit_spy: float,
) -> list[dict[str, Any]]:
    """Point-in-time leg snapshot at entry + exit (or settlement), mirroring *_pnl helpers."""
    tp = trade_params
    exit_d = _norm(ch_x.as_of)
    legs: list[dict[str, Any]] = []

    if trade_kind == "ss":
        dte = int(tp[0])
        pair_e = find_atm_straddle(ch_e, entry_spy, dte)
        if pair_e is None:
            return []
        ce, pe = pair_e
        exp = _norm(ce.expiration)
        strike = float(ce.strike)
        c_x = find_contract_in_chain(ch_x, exp, strike, "C")
        p_x = find_contract_in_chain(ch_x, exp, strike, "P")
        legs.append(_short_leg_row(ce, c_x, exit_spy, exit_d, exp))
        legs.append(_short_leg_row(pe, p_x, exit_spy, exit_d, exp))
        return legs

    if trade_kind == "sl":
        dte = int(tp[0])
        pair_e = find_atm_straddle(ch_e, entry_spy, dte)
        if pair_e is None:
            return []
        ce, pe = pair_e
        exp = _norm(ce.expiration)
        strike = float(ce.strike)
        c_x = find_contract_in_chain(ch_x, exp, strike, "C")
        p_x = find_contract_in_chain(ch_x, exp, strike, "P")
        legs.append(_long_leg_row(ce, c_x, exit_spy, exit_d, exp))
        legs.append(_long_leg_row(pe, p_x, exit_spy, exit_d, exp))
        return legs

    if trade_kind == "sg":
        dte, pdel, cdel = int(tp[0]), float(tp[1]), float(tp[2])
        p0 = find_target_leg_safe(ch_e, dte, pdel, "P")
        c0 = find_target_leg_safe(ch_e, dte, cdel, "C")
        if p0 is None or c0 is None or _norm(p0.expiration) != _norm(c0.expiration):
            return []
        exp = _norm(p0.expiration)
        p_x = find_contract_in_chain(ch_x, exp, float(p0.strike), "P")
        c_x = find_contract_in_chain(ch_x, exp, float(c0.strike), "C")
        legs.append(_short_leg_row(p0, p_x, exit_spy, exit_d, exp))
        legs.append(_short_leg_row(c0, c_x, exit_spy, exit_d, exp))
        return legs

    if trade_kind == "put":
        dte, delt = int(tp[0]), float(tp[1])
        p0 = find_target_leg_safe(ch_e, dte, delt, "P")
        if p0 is None:
            return []
        exp = _norm(p0.expiration)
        p_x = find_contract_in_chain(ch_x, exp, float(p0.strike), "P")
        legs.append(_short_leg_row(p0, p_x, exit_spy, exit_d, exp))
        return legs

    if trade_kind == "rr":
        dte, pdel, cdel = int(tp[0]), float(tp[1]), float(tp[2])
        p0 = find_target_leg_safe(ch_e, dte, pdel, "P")
        c0 = find_target_leg_safe(ch_e, dte, cdel, "C")
        if p0 is None or c0 is None or _norm(c0.expiration) != _norm(p0.expiration):
            return []
        exp = _norm(p0.expiration)
        p_x = find_contract_in_chain(ch_x, exp, float(p0.strike), "P")
        c_x = find_contract_in_chain(ch_x, exp, float(c0.strike), "C")
        legs.append(_short_leg_row(p0, p_x, exit_spy, exit_d, exp))
        legs.append(_long_leg_row(c0, c_x, exit_spy, exit_d, exp))
        return legs

    if trade_kind == "vert":
        dte, short_delta, wing = int(tp[0]), float(tp[1]), float(tp[2])
        ps = find_target_leg_safe(ch_e, dte, short_delta, "P")
        if ps is None:
            return []
        exp = _norm(ps.expiration)
        long_strike = float(ps.strike) - float(wing)
        pl = find_contract_in_chain(ch_e, exp, long_strike, "P")
        if pl is None:
            return []
        p_s_x = find_contract_in_chain(ch_x, exp, float(ps.strike), "P")
        p_l_x = find_contract_in_chain(ch_x, exp, long_strike, "P")
        legs.append(_short_leg_row(ps, p_s_x, exit_spy, exit_d, exp))
        legs.append(_long_leg_row(pl, p_l_x, exit_spy, exit_d, exp))
        return legs

    if trade_kind == "vtc":
        dte, short_delt, wing = int(tp[0]), float(tp[1]), float(tp[2])
        cs = find_target_leg_safe(ch_e, dte, short_delt, "C")
        if cs is None:
            return []
        exp = _norm(cs.expiration)
        long_k = float(cs.strike) + float(wing)
        cl = find_contract_in_chain(ch_e, exp, long_k, "C")
        if cl is None:
            return []
        c_s_x = find_contract_in_chain(ch_x, exp, float(cs.strike), "C")
        c_l_x = find_contract_in_chain(ch_x, exp, long_k, "C")
        legs.append(_short_leg_row(cs, c_s_x, exit_spy, exit_d, exp))
        legs.append(_long_leg_row(cl, c_l_x, exit_spy, exit_d, exp))
        return legs

    return []


def run_signal_backtest_trades(
    days: list[pd.Timestamp],
    get_chain: Callable[[pd.Timestamp], OptionChain],
    panel: pd.DataFrame,
    signal: SignalFn,
    hold: int,
    trade_fn: TradeFn,
    trade_params: tuple,
    trade_kind: str,
    *,
    sid: str,
    family: str,
    description: str,
) -> list[dict[str, Any]]:
    """Same sequencing as :func:`run_signal_backtest`, plus per-trade leg detail rows."""
    rows: list[dict[str, Any]] = []
    pending: tuple[int, float, tuple] | None = None

    for i, d in enumerate(days):
        row = panel.reindex([_norm(d)]).iloc[0]
        spy_c = float(row["close"])
        ch = get_chain(d)
        if pending is not None:
            ent_i, spy_e, tp = pending
            if i >= ent_i + hold:
                ch_e = get_chain(days[ent_i])
                spy_x = float(panel.reindex([_norm(d)]).iloc[0]["close"])
                pnl = trade_fn(ch_e, ch, spy_e, spy_x, tp)
                pending = None
                if pnl is not None and math.isfinite(pnl):
                    ent_day = _norm(days[ent_i])
                    ex_day = _norm(d)
                    erow = panel.reindex([ent_day]).iloc[0]
                    xrow = row
                    vix_e = float(erow["vix_close"]) if pd.notna(erow.get("vix_close")) else float("nan")
                    vix_x = float(xrow["vix_close"]) if pd.notna(xrow.get("vix_close")) else float("nan")
                    legs = collect_trade_legs(trade_kind, tp, ch_e, ch, spy_e, spy_x)
                    rows.append(
                        {
                            "sid": sid,
                            "family": family,
                            "description": description,
                            "trade_kind": trade_kind,
                            "entry_date": ent_day.strftime("%Y-%m-%d"),
                            "exit_date": ex_day.strftime("%Y-%m-%d"),
                            "hold_sessions": int(hold),
                            "calendar_days_in_trade": int((ex_day - ent_day).days),
                            "entry_spy": float(spy_e),
                            "exit_spy": float(spy_x),
                            "vix_entry": vix_e,
                            "vix_exit": vix_x,
                            "pnl_usd": float(pnl),
                            "legs_json": json.dumps(legs, separators=(",", ":")),
                            "n_legs": len(legs),
                        }
                    )
        if pending is not None:
            continue
        if not ch.contracts:
            continue
        if signal(i, row, ch, spy_c):
            pending = (i, spy_c, trade_params)

    return rows


def build_panel(spy_df: pd.DataFrame, days: list[pd.Timestamp]) -> pd.DataFrame:
    o = spy_df.copy()
    o["ret_1"] = o["close"].pct_change()
    o["rv21"] = o["ret_1"].rolling(21).std() * math.sqrt(252.0)
    o["rv5"] = o["ret_1"].rolling(5).std() * math.sqrt(252.0)
    o["sma_50"] = o["close"].rolling(50, min_periods=1).mean()
    o["vix_ma20"] = o["vix_close"].rolling(20).mean()
    if "vvix_close" in o.columns:
        o["vvix_med20"] = o["vvix_close"].rolling(20, min_periods=5).median()
    else:
        o["vvix_med20"] = np.nan
    idx = [_norm(d) for d in days]
    return o.reindex(idx).ffill()


def _wrap_straddle_short(tp: tuple) -> TradeFn:
    dte = int(tp[0])

    def fn(
        ch0: OptionChain, ch1: OptionChain, se: float, sx: float, _tp: tuple
    ) -> float | None:
        return short_straddle_pnl(ch0, ch1, se, sx, dte)

    return fn


def _wrap_straddle_long(tp: tuple) -> TradeFn:
    dte = int(tp[0])

    def fn(ch0: OptionChain, ch1: OptionChain, se: float, sx: float, _tp: tuple) -> float | None:
        return long_straddle_pnl(ch0, ch1, se, sx, dte)

    return fn


def _wrap_strangle(tp: tuple) -> TradeFn:
    dte, pdel, cdel = int(tp[0]), float(tp[1]), float(tp[2])

    def fn(ch0: OptionChain, ch1: OptionChain, _se: float, sx: float, _tp: tuple) -> float | None:
        return short_strangle_pnl(ch0, ch1, sx, dte, pdel, cdel)

    return fn


def _wrap_rr(tp: tuple) -> TradeFn:
    dte, pdel, cdel = int(tp[0]), float(tp[1]), float(tp[2])

    def fn(ch0: OptionChain, ch1: OptionChain, _se: float, sx: float, _tp: tuple) -> float | None:
        return short_risk_reversal_pnl(ch0, ch1, sx, dte, pdel, cdel)

    return fn


def _wrap_put(tp: tuple) -> TradeFn:
    dte, delt = int(tp[0]), float(tp[1])

    def fn(ch0: OptionChain, ch1: OptionChain, _se: float, sx: float, _tp: tuple) -> float | None:
        return short_put_pnl(ch0, ch1, sx, dte, delt)

    return fn


def _wrap_vert(tp: tuple) -> TradeFn:
    dte, delt, w = int(tp[0]), float(tp[1]), float(tp[2])

    def fn(ch0: OptionChain, ch1: OptionChain, _se: float, sx: float, _tp: tuple) -> float | None:
        return short_vertical_put_pnl(ch0, ch1, sx, dte, delt, w)

    return fn


def _wrap_vert_call(tp: tuple) -> TradeFn:
    dte, delt, w = int(tp[0]), float(tp[1]), float(tp[2])

    def fn(ch0: OptionChain, ch1: OptionChain, _se: float, sx: float, _tp: tuple) -> float | None:
        return short_vertical_call_pnl(ch0, ch1, sx, dte, delt, w)

    return fn


def prepare_theta_research_context(
    *,
    theta_dir: Path,
    capital: float,
    start: str = "",
    end: str = "",
    max_days: int = 0,
) -> tuple[
    list[pd.Timestamp],
    pd.DataFrame,
    Callable[[pd.Timestamp], OptionChain],
    dict[tuple[int, int], float | None],
    dict[tuple[int, int, float, float], float | None],
    list[int],
    pd.DataFrame,
]:
    """
    Shared precompute for literature backtests: trading days, panel, chain getter,
    ATM IV cache, skew IV-diff cache, contract counts per session, raw spy_wide (VVIX check).
    """
    theta_dir = theta_dir.expanduser().resolve()
    d0, d1 = theta_chunks_date_bounds(theta_dir)
    yf_start = (d0 - pd.Timedelta(days=400)).strftime("%Y-%m-%d")
    yf_end = (d1 + pd.Timedelta(days=14)).strftime("%Y-%m-%d")
    spy_wide = normalize_spy_df(load_spy_vix_from_yfinance(yf_start, yf_end))
    ld = ThetaChunksLoader(theta_dir, spy_wide)
    days = trading_days_intersecting_spy(ld, spy_wide.index, d0, d1)
    if int(max_days) > 0:
        days = days[: int(max_days)]
    if str(start).strip():
        t0 = pd.Timestamp(str(start).strip())
        days = [d for d in days if _norm(d) >= t0]
    if str(end).strip():
        t1 = pd.Timestamp(str(end).strip())
        days = [d for d in days if _norm(d) <= t1]

    panel = build_panel(spy_wide, days)
    chain_cache: dict[pd.Timestamp, OptionChain] = {}
    _CHAIN_KW = dict(strike_pct_lo=0.76, strike_pct_hi=1.12, min_dte=5, max_dte=60)

    def get_chain(d: pd.Timestamp) -> OptionChain:
        t = _norm(d)
        if t not in chain_cache:
            chain_cache[t] = ld.get_chain_for_date(t, **_CHAIN_KW)
        return chain_cache[t]

    n = len(days)
    print(f"  precompute {n} sessions (chains + IV + skew) …", flush=True)
    dte_universe = [7, 10, 14, 18, 20, 21, 25, 30, 35, 40, 45, 55]
    iv_atm: dict[tuple[int, int], float | None] = {}
    n_contracts: list[int] = [0] * n
    for i in range(n):
        if i % 50 == 0:
            print(f"  precompute chains/iv {i + 1}/{n} …", flush=True)
        ch = get_chain(days[i])
        n_contracts[i] = len(ch.contracts)
        sp = float(panel.iloc[i]["close"])
        for dt in dte_universe:
            iv_atm[(i, dt)] = atm_iv_straddle(ch, sp, dt)

    skew_keys = [
        (30, -0.12, 0.08),
        (30, -0.15, 0.10),
        (30, -0.18, 0.12),
        (30, -0.20, 0.12),
        (40, -0.12, 0.08),
        (40, -0.15, 0.10),
        (40, -0.18, 0.12),
        (40, -0.20, 0.12),
        (45, -0.12, 0.08),
        (45, -0.15, 0.10),
        (45, -0.18, 0.12),
        (45, -0.20, 0.12),
    ]
    skew_put_minus_call_iv: dict[tuple[int, int, float, float], float | None] = {}
    for i in range(n):
        if i % 50 == 0:
            print(f"  precompute skew {i + 1}/{n} …", flush=True)
        ch = get_chain(days[i])
        for dte, pdel, cdel in skew_keys:
            pleg = find_target_leg_safe(ch, dte, pdel, "P")
            cleg = find_target_leg_safe(ch, dte, cdel, "C")
            key = (i, dte, pdel, cdel)
            if (
                pleg is not None
                and cleg is not None
                and math.isfinite(float(pleg.iv))
                and math.isfinite(float(cleg.iv))
            ):
                skew_put_minus_call_iv[key] = float(pleg.iv) - float(cleg.iv)
            else:
                skew_put_minus_call_iv[key] = None

    _ = float(capital)  # reserved for API symmetry with callers that log capital
    return days, panel, get_chain, iv_atm, skew_put_minus_call_iv, n_contracts, spy_wide


def main() -> None:
    ap = argparse.ArgumentParser()
    ap.add_argument("--theta-dir", type=Path, default=_DEFAULT_THETA)
    ap.add_argument("--capital", type=float, default=1_000_000.0)
    ap.add_argument("--min-trades", type=int, default=35)
    ap.add_argument("--sharpe-min", type=float, default=1.5)
    ap.add_argument("--max-days", type=int, default=0)
    ap.add_argument("--start", type=str, default="", help="YYYY-MM-DD inclusive lower bound on sessions")
    ap.add_argument("--end", type=str, default="", help="YYYY-MM-DD inclusive upper bound on sessions")
    ap.add_argument("--list-top", type=int, default=25)
    args = ap.parse_args()

    theta_dir = args.theta_dir.expanduser().resolve()
    t_build0 = time.perf_counter()
    print("Building narrow chains + IV cache (one pass) …", flush=True)
    days, panel, get_chain, iv_atm, skew_put_minus_call_iv, n_contracts, spy_wide = (
        prepare_theta_research_context(
            theta_dir=theta_dir,
            capital=float(args.capital),
            start=str(args.start),
            end=str(args.end),
            max_days=int(args.max_days),
        )
    )
    n = len(days)
    print(
        f"Chain+IV+skew precompute done in {(time.perf_counter() - t_build0)/60:.1f} min ({n} sessions).",
        flush=True,
    )

    scored: list[tuple[float, int, str, str]] = []

    def try_add(name: str, lit: str, sig: SignalFn, hold: int, tfn: TradeFn, tp: tuple) -> None:
        ex, pn, n = run_signal_backtest(days, get_chain, panel, sig, hold, tfn, tp)
        _, sh = equity_curve_from_realized(ex, pn, days, float(args.capital))
        if not math.isfinite(sh):
            return
        scored.append((sh, n, name, lit))

    # ---- 1) Variance risk premium: short ATM straddle when IV > RV + k·RV ----
    for thr in (0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 2.0):
        for dte in (20, 25, 30, 35, 40, 45):
            for h in (2, 3, 5, 7, 10):

                def make_vrp_sig(threshold: float, target_dte: int):
                    def sig(i: int, row: pd.Series, ch: OptionChain, _s: float) -> bool:
                        iv = iv_atm.get((i, target_dte))
                        rv = float(row["rv21"]) if pd.notna(row.get("rv21")) else float("nan")
                        if iv is None or not math.isfinite(rv) or rv <= 0:
                            return False
                        return iv > rv + threshold * rv

                    return sig

                try_add(
                    f"VRP short straddle IV>RV+{thr}·RV DTE≈{dte} H={h}d",
                    "Bollerslev–Tauchen–Zhou; Carr–Wu: short variance when implied exceeds realized.",
                    make_vrp_sig(thr, dte),
                    h,
                    _wrap_straddle_short((dte,)),
                    (dte,),
                )

    # ---- 2) Calm regime short strangle (low VIX) ----
    for vx in (11.0, 12.0, 13.0, 14.0, 15.0, 16.0):
        for dte in (30, 40, 45, 55):
            for h in (3, 5, 7, 10, 12):
                for pdel, cdel in ((-0.18, 0.12), (-0.20, 0.15), (-0.22, 0.15), (-0.25, 0.18), (-0.30, 0.20)):

                    def make_lowv(vlim: float):
                        def sig(i: int, row: pd.Series, ch: OptionChain, _s: float) -> bool:
                            return float(row["vix_close"]) < float(vlim) and n_contracts[i] > 80

                        return sig

                    try_add(
                        f"Low-VIX strangle VIX<{vx} ΔP={pdel} ΔC={cdel} DTE≈{dte} H={h}d",
                        "Index short-vol benchmarks (BXM/PUT analogues): harvest in low fear.",
                        make_lowv(float(vx)),
                        h,
                        _wrap_strangle((dte, pdel, cdel)),
                        (dte, pdel, cdel),
                    )

    # ---- 3) VIX spike → long straddle (vol mean reversion) ----
    for chg in (0.08, 0.10, 0.12, 0.15, 0.18, 0.22):
        for dte in (7, 10, 14, 18, 21):
            for h in (2, 3, 4, 5, 7):

                def make_spike(pct: float):
                    def sig(i: int, row: pd.Series, ch: OptionChain, _s: float) -> bool:
                        if i < 5:
                            return False
                        v0 = float(panel.iloc[i - 5]["vix_close"])
                        v1 = float(row["vix_close"])
                        return v1 > v0 * (1.0 + float(pct))

                    return sig

                try_add(
                    f"Post-shock long straddle VIX>{int(chg*100)}%/5d DTE≈{dte} H={h}d",
                    "Equity index vol spikes mean-revert; lottery demand lifts OTM vol (Bondarenko, others).",
                    make_spike(float(chg)),
                    h,
                    _wrap_straddle_long((dte,)),
                    (dte,),
                )

    # ---- 4) VIX decline → short straddle (Stein / Poteshman) ----
    for dv in (-0.75, -1.0, -1.25, -1.5, -2.0, -2.5):
        for dte in (25, 30, 35, 40):
            for h in (3, 5, 7, 10):

                def make_fall(dvx: float):
                    def sig(i: int, row: pd.Series, ch: OptionChain, _s: float) -> bool:
                        if i < 5:
                            return False
                        return float(row["vix_close"]) - float(panel.iloc[i - 5]["vix_close"]) < float(
                            dvx
                        )

                    return sig

                try_add(
                    f"VIX 5d drop<{dv} short straddle DTE≈{dte} H={h}d",
                    "Implied vol mean-reversion after declines.",
                    make_fall(float(dv)),
                    h,
                    _wrap_straddle_short((dte,)),
                    (dte,),
                )

    # ---- 5) Bull trend + short OTM put (PUTWrite-style) ----
    for dte in (35, 45, 55):
        for delt in (-0.18, -0.22, -0.25, -0.30):
            for h in (5, 8, 10, 15, 21):

                def sig_bull(i: int, row: pd.Series, ch: OptionChain, spy: float) -> bool:
                    sma = float(row["sma_200"]) if pd.notna(row.get("sma_200")) else float("nan")
                    return math.isfinite(sma) and spy > sma * 1.0 and float(row["vix_close"]) < 22

                try_add(
                    f"PUTWrite-like SPY>SMA200 δ={delt} DTE≈{dte} H={h}d",
                    "Hill–Jain–Kan (PUT index); Bollen–Whaley crash insurance pricing.",
                    sig_bull,
                    h,
                    _wrap_put((dte, delt)),
                    (dte, delt),
                )

    # ---- 6) Short skew / risk reversal when IV(RM) steep vs calls ----
    for dte in (30, 40, 45):
        for pdel, cdel in ((-0.12, 0.08), (-0.15, 0.10), (-0.18, 0.12), (-0.20, 0.12)):
            for h in (5, 8, 10, 15):
                thr_iv = 0.04

                def make_skew(target_dte: int, put_d: float, call_d: float, gap_iv: float):
                    def sig(i: int, row: pd.Series, ch: OptionChain, _s: float) -> bool:
                        diff = skew_put_minus_call_iv.get((i, target_dte, put_d, call_d))
                        return diff is not None and diff > float(gap_iv)

                    return sig

                try_add(
                    f"Short RR IVput−IVcall>{thr_iv} ΔP={pdel} ΔC={cdel} DTE≈{dte} H={h}d",
                    "Index skew premium (Bakshi–Kapadia–Madan; Carr–Wu skew).",
                    make_skew(int(dte), float(pdel), float(cdel), float(thr_iv)),
                    h,
                    _wrap_rr((dte, pdel, cdel)),
                    (dte, pdel, cdel),
                )

    # ---- 7) VVIX filter + short straddle (tail risk off) ----
    if "vvix_close" in spy_wide.columns:
        for vvr in (1.0, 1.05, 1.1, 1.15):
            for dte in (30, 40):
                for h in (5, 8):

                    def make_vv(vvix_vix_cap: float, target_dte: int):
                        def sig(i: int, row: pd.Series, ch: OptionChain, _s: float) -> bool:
                            vv = row.get("vvix_close")
                            vx = row.get("vix_close")
                            if pd.isna(vv) or pd.isna(vx) or float(vx) <= 0:
                                return False
                            if float(vv) / float(vx) > float(vvix_vix_cap):
                                return False
                            iv = iv_atm.get((i, target_dte))
                            rv = float(row["rv21"]) if pd.notna(row.get("rv21")) else float("nan")
                            return iv is not None and math.isfinite(rv) and rv > 0 and iv > rv * 1.05

                        return sig

                    try_add(
                        f"VVIX/VIX<{vvr} & mild VRP short straddle DTE≈{dte} H={h}d",
                        "Stylized tail-risk scaling: avoid shorting vol when vol-of-vol is rich.",
                        make_vv(float(vvr), int(dte)),
                        h,
                        _wrap_straddle_short((dte,)),
                        (dte,),
                    )

    # ---- 8) Realized vol compression → short straddle ----
    for gap in (0.02, 0.04, 0.06, 0.08):
        for dte in (30, 40):
            for h in (5, 8, 10):

                def make_rvcomp(rv_gap: float, target_dte: int):
                    def sig(i: int, row: pd.Series, ch: OptionChain, _s: float) -> bool:
                        r5 = float(row["rv5"]) if pd.notna(row.get("rv5")) else float("nan")
                        r21 = float(row["rv21"]) if pd.notna(row.get("rv21")) else float("nan")
                        if not (math.isfinite(r5) and math.isfinite(r21)) or r21 <= 0:
                            return False
                        if r21 - r5 <= float(rv_gap):
                            return False
                        iv = iv_atm.get((i, target_dte))
                        return iv is not None and iv > r21

                    return sig

                try_add(
                    f"RV21−RV5>{gap} short straddle DTE≈{dte} H={h}d",
                    "Vol clustering: after spot calms, short-dated RV falls vs longer window.",
                    make_rvcomp(float(gap), int(dte)),
                    h,
                    _wrap_straddle_short((dte,)),
                    (dte,),
                )

    # ---- 9) Credit put spread: mild VIX, defined risk ----
    for dte in (25, 35, 45):
        for delt in (-0.22, -0.28, -0.32):
            for wing in (3.0, 5.0, 7.0, 10.0):
                for h in (5, 10, 15):

                    def sig_midvx(i: int, row: pd.Series, ch: OptionChain, _s: float) -> bool:
                        v = float(row["vix_close"])
                        return 14.0 < v < 22.0

                    try_add(
                        f"Mild-VIX put spread δ={delt} wing={wing} DTE≈{dte} H={h}d",
                        "Collateralized short put risk in normal regimes (CBOE PPUT / defined-risk).",
                        sig_midvx,
                        h,
                        _wrap_vert((dte, delt, wing)),
                        (dte, delt, wing),
                    )

    # ---- 10) High VIX + still upward VRP → short straddle (contrarian) ----
    for vx0 in (22.0, 25.0, 28.0):
        for thr in (0.1, 0.2, 0.3):
            for dte in (35, 45):
                for h in (5, 8):

                    def make_hivx(vix_floor: float, rv_mult: float, target_dte: int):
                        def sig(i: int, row: pd.Series, ch: OptionChain, _s: float) -> bool:
                            if float(row["vix_close"]) < float(vix_floor):
                                return False
                            iv = iv_atm.get((i, target_dte))
                            rv = float(row["rv21"]) if pd.notna(row.get("rv21")) else float("nan")
                            if iv is None or not math.isfinite(rv) or rv <= 0:
                                return False
                            return iv > rv + float(rv_mult) * rv

                        return sig

                    try_add(
                        f"High-VIX>{vx0} VRP short straddle IV>RV+{thr}·RV DTE≈{dte} H={h}d",
                        "Even in stress, documented positive variance risk premium on indices.",
                        make_hivx(float(vx0), float(thr), int(dte)),
                        h,
                        _wrap_straddle_short((dte,)),
                        (dte,),
                    )

    ok = [t for t in scored if t[1] >= int(args.min_trades) and t[0] >= float(args.sharpe_min)]
    ok.sort(key=lambda x: -x[0])

    def _family(name: str) -> str:
        for prefix in (
            "VRP short",
            "Low-VIX",
            "Post-shock",
            "VIX 5d",
            "PUTWrite-like",
            "Short RR",
            "VVIX/VIX",
            "RV21−RV5",
            "Mild-VIX",
            "High-VIX",
        ):
            if name.startswith(prefix):
                return prefix
        return name.split()[0] if name else "?"

    best_by_family: dict[str, tuple[float, int, str, str]] = {}
    for row in ok:
        fam = _family(row[2])
        if fam not in best_by_family or row[0] > best_by_family[fam][0]:
            best_by_family[fam] = row

    print(f"Sessions: {len(days)}  ({days[0].date()} → {days[-1].date()})")
    print(f"Candidates evaluated: {len(scored)}")
    print(f"Pass Sharpe≥{args.sharpe_min} & trades≥{args.min_trades}: {len(ok)}\n")
    for sh, n, name, lit in ok[: int(args.list_top)]:
        print(f"Sharpe={sh:.3f}  trades={n}  {name}")
        print(f"  └─ {lit}\n")
    if len(ok) < 10:
        print("--- Top 15 by Sharpe (relax min) ---")
        scored.sort(key=lambda x: -x[0])
        for sh, n, name, lit in scored[:15]:
            print(f"Sharpe={sh:.3f}  trades={n}  {name}")
    else:
        fam_sorted = sorted(best_by_family.values(), key=lambda x: -x[0])
        print("--- Best qualifying spec per strategy family ---")
        for sh, n, name, lit in fam_sorted:
            print(f"{name}  |  Sharpe={sh:.3f}  n={n}")


if __name__ == "__main__":
    main()