#!/usr/bin/env python3
"""Equity-index VRP, robust internal-parity version (self-consistent, no external underlying).

Per day: implied forward F = MEDIAN over the 5 near-ATM strikes of (K + C - P) [put-call parity,
median-smoothed to kill per-strike noise]. Use F consistently for ATM, implied vol, AND the
realized-vol series. Internally consistent -> avoids both the single-strike-parity noise and the
SPY!=XSP divergence that produced spurious negative VRP. XSP 2013-2026.
"""
import math
import os
import sys

import numpy as np

sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
from equity_vrp import load_xsp, to_epoch_day  # noqa: E402
from signal_sweep import sharpe_t  # noqa: E402
import torch  # noqa: E402

DEV = "cuda" if torch.cuda.is_available() else "cpu"


def main():
    d = load_xsp()
    d["exp_day"] = d["expiry"].map(to_epoch_day)
    d["ttm"] = d["exp_day"] - d["day"]
    d = d[(d["ttm"] >= 20) & (d["ttm"] <= 45)]
    F, IVm = {}, {}
    for day, g in d.groupby("day"):
        exp = g.iloc[(g["ttm"] - 30).abs().argsort()].iloc[0]["exp_day"]
        ge = g[g["exp_day"] == exp]
        calls = ge[ge["right"] == "C"].groupby("strike")["close"].mean()
        puts = ge[ge["right"] == "P"].groupby("strike")["close"].mean()
        common = np.array(sorted(set(calls.index) & set(puts.index)))
        if len(common) < 5:
            continue
        rough = common[np.abs(calls[common].values - puts[common].values).argmin()]   # min|C-P| ~ ATM
        near = common[np.argsort(np.abs(common - rough))[:5]]                          # 5 nearest strikes
        f = float(np.median([k + float(calls[k]) - float(puts[k]) for k in near]))     # parity forward (median)
        katm = common[np.abs(common - f).argmin()]
        straddle = float(calls[katm] + puts[katm])
        T = float(ge["ttm"].iloc[0]) / 365.0
        if f <= 0 or straddle <= 0 or T <= 0:
            continue
        F[int(day)] = f; IVm[int(day)] = straddle / (0.8 * f * math.sqrt(T))

    days = np.array(sorted(F))
    S = np.array([F[x] for x in days]); IV = np.array([IVm[x] for x in days])
    r = np.zeros(len(days)); r[1:] = np.log(S[1:] / S[:-1])
    r = np.clip(r, -0.25, 0.25)                                    # guard residual parity glitches
    H = 21
    rv = np.full(len(days), np.nan)
    for t in range(len(days) - H):
        rv[t] = np.std(r[t + 1:t + 1 + H]) * math.sqrt(252)
    vrp = IV - rv; fin = np.isfinite(vrp)
    iv_lag = np.concatenate([[np.nan], IV[:-1]])
    svol = (iv_lag ** 2) / 252.0 - r ** 2
    year = np.array([1970 + x / 365.25 for x in days]).astype(int)
    T_ = lambda x: torch.tensor(np.asarray(x)[np.isfinite(np.asarray(x))], device=DEV, dtype=torch.float64)

    print(f"\n===== EQUITY-INDEX VRP (robust internal parity) — {len(days)} days ({days.min()}..{days.max()}) =====")
    print(f"implied-forward level {S[0]:.0f}->{S[-1]:.0f}  (should track SPX/10 ~ 400->650)")
    print(f"mean implied {np.nanmean(IV):.1%}  mean fwd-realized {np.nanmean(rv):.1%}  "
          f"mean VRP {np.nanmean(vrp[fin]):+.1%}  frac>0 {np.mean(vrp[fin]>0):.2f}")
    sk = float(((svol[np.isfinite(svol)] - np.nanmean(svol)) ** 3).mean() / np.nanstd(svol) ** 3)
    print(f"short-vol daily P&L: Sharpe {sharpe_t(T_(svol)):+.2f}  skew {sk:+.2f}  worst {np.nanmin(svol)/np.nanstd(svol):+.1f}σ")
    print("per-year: " + " ".join(f"{y}:{sharpe_t(T_(svol[year==y])):+.1f}" for y in range(2013, 2027) if (year == y).sum() > 60))
    rising = np.zeros(len(days))
    for t in range(6, len(days)):
        rising[t] = 0.0 if IV[t - 1] > IV[t - 6] else 1.0
    tm = rising * svol
    print(f"tail-managed: Sharpe {sharpe_t(T_(tm)):+.2f}  worst {np.nanmin(tm)/np.nanstd(tm):+.1f}σ")
    print("per-year TM: " + " ".join(f"{y}:{sharpe_t(T_(tm[year==y])):+.1f}" for y in range(2013, 2027) if (year == y).sum() > 60))
    print("\nVERDICT: implied-forward tracking ~SPX/10 + mean RV ~13-16% + VRP>0 => measurement clean & premium real.")


if __name__ == "__main__":
    main()