"""Equity-index Volatility Risk Premium (VRP) — a STRUCTURAL risk premium (not momentum): implied vol
exceeds subsequent realized vol on average because option sellers are paid to bear crash risk. Short vol
harvests IV − RV. The catch (and why the deflated Sharpe matters): VRP returns are negatively skewed and
fat-tailed (vol spikes = big losses), so a raw Sharpe overstates it — the skew/kurtosis-adjusted deflated
Sharpe is the honest measure. Built on XSP (mini-SPX) options, 2013-2026, internal put-call-parity (vrp3).
"""
from __future__ import annotations

import datetime as dt
import glob
import math
import os
from dataclasses import dataclass

import numpy as np

from fxhnt.domain.validation import deflated_sharpe_ratio, excess_kurtosis, sharpe_ratio, skewness

_DAY_NS = 86_400 * 10**9


def _to_epoch_day(yymmdd: str) -> int:
    y, m, d = 2000 + int(yymmdd[:2]), int(yymmdd[2:4]), int(yymmdd[4:6])
    return (dt.date(y, m, d) - dt.date(1970, 1, 1)).days


def _load_xsp(xsp_dir: str):  # type: ignore[no-untyped-def]
    import databento as db
    import pandas as pd

    rows = []
    for p in sorted(glob.glob(os.path.join(xsp_dir, "*.dbn"))):
        try:
            df = db.DBNStore.from_file(p).to_df().reset_index()
        except Exception:
            continue
        if df.empty or "symbol" not in df.columns:
            continue
        df = df[df["close"] > 0]
        rows.append(df[["ts_event", "symbol", "close"]])
    d = pd.concat(rows, ignore_index=True)
    s = d["symbol"].astype(str).str.replace(" ", "", regex=False)   # XSP240119C00450000
    d["right"] = s.str[-9]
    d["strike"] = s.str[-8:].astype(float) / 1000.0
    d["expiry"] = s.str[-15:-9]
    d["day"] = d["ts_event"].astype("int64") // _DAY_NS
    return d


@dataclass
class VrpResult:
    n_periods: int
    mean_vrp: float
    sharpe: float
    skew: float
    excess_kurt: float
    deflated_sharpe: float       # skew/kurtosis-adjusted (the honest measure)
    dsr_pvalue: float
    is_sharpe: float
    oos_sharpe: float
    max_drawdown: float


def evaluate(xsp_dir: str, *, holding_days: int = 21, n_trials: int = 5, oos_fraction: float = 0.40) -> VrpResult:
    import pandas as pd

    d = _load_xsp(xsp_dir)
    d["exp_day"] = d["expiry"].map(_to_epoch_day)
    d["ttm"] = d["exp_day"] - d["day"]
    d = d[(d["ttm"] >= 20) & (d["ttm"] <= 45)]

    forward, iv = {}, {}
    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()]
        near = common[np.argsort(np.abs(common - rough))[:5]]
        f = float(np.median([k + float(calls[k]) - float(puts[k]) for k in near]))
        katm = common[np.abs(common - f).argmin()]
        straddle = float(calls[katm] + puts[katm])
        t_yrs = float(ge["ttm"].iloc[0]) / 365.0
        if f <= 0 or straddle <= 0 or t_yrs <= 0:
            continue
        forward[int(day)] = f
        iv[int(day)] = straddle / (0.8 * f * math.sqrt(t_yrs))   # ATM-straddle implied-vol approximation

    days = np.array(sorted(forward))
    fwd = np.array([forward[x] for x in days])
    iv_arr = np.array([iv[x] for x in days])
    logret = np.zeros(len(days))
    logret[1:] = np.clip(np.log(fwd[1:] / fwd[:-1]), -0.25, 0.25)

    # NON-OVERLAPPING short-vol periods: enter, hold `holding_days`, P&L proxy = IV_entry − RV_realized
    returns = []
    step = holding_days
    for t in range(0, len(days) - step, step):
        rv = float(np.std(logret[t + 1:t + 1 + step]) * math.sqrt(252))
        returns.append(iv_arr[t] - rv)
    r = np.array(returns)
    if len(r) < 8:
        raise ValueError("not enough VRP periods")

    sk, ek = skewness(r), excess_kurtosis(r)
    dsr = deflated_sharpe_ratio(sharpe_ratio(r), n_trials, 0.0, sk, ek, len(r))
    split = int((1.0 - oos_fraction) * len(r))
    eq = np.cumprod(1.0 + r / 100.0)                                # scale vol-points to a rough equity curve
    mdd = float((eq / np.maximum.accumulate(eq) - 1.0).min())
    return VrpResult(
        n_periods=len(r), mean_vrp=float(r.mean()), sharpe=sharpe_ratio(r), skew=sk, excess_kurt=ek,
        deflated_sharpe=1.0 - dsr.pvalue, dsr_pvalue=dsr.pvalue,
        is_sharpe=sharpe_ratio(r[:split]), oos_sharpe=sharpe_ratio(r[split:]), max_drawdown=mdd)
