PythagorasGoal/scripts/research/options_vrp_lab.py

"""OPTIONS VRP LAB — sleeve income: vendita put settimanali (CSP) che incassa il VRP (IV>RV).

Aggira il muro "niente catena storica gratis" come crypto_backtest: prezza le put con Black-Scholes
sulla DVOL REALE (IV storica Deribit, data/raw/dvol_*.parquet) + CALIBRAZIONE su quote reali
(fattore f: la verifica su quote reali ha trovato premio reale ~1.29x il modellato a IV-ATM per via
dello skew, al netto dello spread). Payoff sul path REALIZZATO dei prezzi certificati. Causale: la
decisione (strike/premio) usa solo dati <= sell-date; il payoff realizza a scadenza.

Onesto: e' SHORT-VOL, il rischio vero e' la CODA (crash). Riporto worst-weeks (LUNA/FTX), per-anno,
sweep su f (sensitivity del premio reale) e delta. NON e' un deploy: e' la prima validazione del lead.

    uv run python scripts/research/options_vrp_lab.py
"""
from __future__ import annotations
import sys
from pathlib import Path
PROJECT_ROOT = Path(__file__).resolve().parents[2]
sys.path.insert(0, str(PROJECT_ROOT))
import numpy as np, pandas as pd
from scipy.stats import norm
from scripts.analysis.research_lab import load_tf

HOLDOUT = pd.Timestamp("2025-01-01", tz="UTC")
WK_PER_YEAR = 365.25 / 7.0


def bs_put(S, K, T, sig):
    if T <= 0 or sig <= 0:
        return max(K - S, 0.0)
    d1 = (np.log(S / K) + 0.5 * sig ** 2 * T) / (sig * np.sqrt(T))
    d2 = d1 - sig * np.sqrt(T)
    return K * norm.cdf(-d2) - S * norm.cdf(-d1)        # r=0


def strike_from_delta(S, T, sig, target_delta=-0.28):
    # delta_put = -N(-d1) = target -> d1 = -N^{-1}(-target)
    d1 = -norm.ppf(-target_delta)
    return S * np.exp(0.5 * sig ** 2 * T - d1 * sig * np.sqrt(T))


def load_series(asset):
    px = load_tf(asset, "1d")
    s = pd.Series(px["close"].values.astype(float), index=pd.to_datetime(px["timestamp"], unit="ms", utc=True))
    dv = pd.read_parquet(PROJECT_ROOT / "data" / "raw" / f"dvol_{asset.lower()}.parquet")
    d = pd.Series(dv["close"].values.astype(float), index=pd.to_datetime(dv["timestamp"], unit="ms", utc=True))
    J = pd.concat({"px": s, "dvol": d}, axis=1, join="inner").sort_index().dropna()
    return J


def put_sell_weekly(asset, delta=-0.28, f=1.0, tenor_d=7):
    """Vendita CSP settimanale. Ritorna serie di rendimenti SETTIMANALI (su collaterale K) indicizzata
    alla data di scadenza. Causale: strike/premio da DVOL e prezzo a sell-date; payoff a scadenza."""
    J = load_series(asset)
    px = J["px"].values; dv = J["dvol"].values / 100.0; idx = J.index
    n = len(px); T = tenor_d / 365.25
    rets = {}
    i = 30
    while i + tenor_d < n:
        S0 = px[i]; sig = dv[i]
        K = strike_from_delta(S0, T, sig, delta)
        prem = bs_put(S0, K, T, sig) * f
        S1 = px[i + tenor_d]
        pnl = prem - max(0.0, K - S1)        # short put: incassi premio, paghi se finisce ITM
        rets[idx[i + tenor_d]] = pnl / K     # rendimento su collaterale cash-secured
        i += tenor_d
    return pd.Series(rets)


def m_weekly(r):
    r = r.dropna()
    if len(r) < 3 or r.std() == 0:
        return dict(sh=0, cagr=0, dd=0, n=len(r))
    eq = np.cumprod(1 + r.values); pk = np.maximum.accumulate(eq)
    yrs = len(r) / WK_PER_YEAR
    return dict(sh=float(r.mean() / r.std() * np.sqrt(WK_PER_YEAR)),
                cagr=float(eq[-1] ** (1 / yrs) - 1) if yrs > 0 and eq[-1] > 0 else 0,
                dd=float(np.max((pk - eq) / pk)), n=len(r))


def per_year(r):
    out = {}
    for y, g in r.groupby(r.index.year):
        eq = np.cumprod(1 + g.values)
        out[int(y)] = float(eq[-1] - 1)
    return out


def main():
    print("=" * 96)
    print("  OPTIONS VRP LAB — vendita put settimanali (CSP), premio BS su DVOL reale + calibrazione f")
    print("=" * 96)

    # contesto VRP: IV (DVOL) vs RV realizzata
    for a in ("BTC", "ETH"):
        J = load_series(a)
        rv = J["px"].pct_change().rolling(30).std() * np.sqrt(365.25) * 100
        vrp = (J["dvol"] - rv).dropna()
        print(f"  {a}: DVOL media {J['dvol'].mean():.0f}% | RV30 media {rv.mean():.0f}% | VRP media {vrp.mean():+.1f} pt, >0 nel {100*(vrp>0).mean():.0f}% del tempo")

    print("\n  (1) SWEEP CALIBRAZIONE f (delta -0.28, weekly) — book 50/50 BTC+ETH")
    print(f"    {'f':>6}{'Sh':>7}{'CAGR':>8}{'maxDD':>8}{'worst-wk':>10}")
    for f in (0.70, 0.85, 1.0, 1.15, 1.29):
        rB = put_sell_weekly("BTC", f=f); rE = put_sell_weekly("ETH", f=f)
        book = pd.concat({"B": rB, "E": rE}, axis=1, join="inner").mean(axis=1)
        mm = m_weekly(book); worst = book.min()
        tag = "  <- reale(calm)" if f == 1.29 else ("  <- conservativo" if f == 1.0 else "")
        print(f"    {f:>6.2f}{mm['sh']:>7.2f}{mm['cagr']*100:>+7.0f}%{mm['dd']*100:>7.1f}%{worst*100:>+9.1f}%{tag}")

    print("\n  (2) SWEEP DELTA (f=1.0 conservativo) — book 50/50")
    print(f"    {'delta':>7}{'Sh':>7}{'CAGR':>8}{'maxDD':>8}")
    for dl in (-0.15, -0.28, -0.40):
        rB = put_sell_weekly("BTC", delta=dl); rE = put_sell_weekly("ETH", delta=dl)
        book = pd.concat({"B": rB, "E": rE}, axis=1, join="inner").mean(axis=1)
        mm = m_weekly(book)
        print(f"    {dl:>7.2f}{mm['sh']:>7.2f}{mm['cagr']*100:>+7.0f}%{mm['dd']*100:>7.1f}%")

    # config centrale: delta -0.28, f=1.0 (conservativo) e f=1.29 (reale misurato)
    print("\n  (3) PER ANNO + WORST WEEKS (delta -0.28, book 50/50) — il rischio e' la CODA")
    for f in (1.0, 1.29):
        rB = put_sell_weekly("BTC", f=f); rE = put_sell_weekly("ETH", f=f)
        book = pd.concat({"B": rB, "E": rE}, axis=1, join="inner").mean(axis=1)
        py = per_year(book)
        worst = book.nsmallest(5)
        print(f"\n    f={f}: per-anno " + "  ".join(f"{y}:{v*100:+.0f}%" for y, v in py.items()))
        print(f"           worst weeks: " + "  ".join(f"{d.date()}:{v*100:.0f}%" for d, v in worst.items()))
        full = m_weekly(book); ho = m_weekly(book[book.index >= HOLDOUT])
        print(f"           FULL Sh {full['sh']:.2f} CAGR {full['cagr']*100:+.0f}% DD {full['dd']*100:.0f}% | HOLD-OUT Sh {ho['sh']:.2f}")

    # correlazione e contributo vs TP01 (resampling settimanale)
    print("\n  (4) CORRELAZIONE + CONTRIBUTO vs TP01 (settimanale; f=1.0 conservativo)")
    from src.portfolio.sleeves import tp01_sleeve
    tp = tp01_sleeve().daily()
    tp_wk = (1 + tp).resample("7D").prod() - 1
    rB = put_sell_weekly("BTC"); rE = put_sell_weekly("ETH")
    opt = pd.concat({"B": rB, "E": rE}, axis=1, join="inner").mean(axis=1)
    opt_wk = opt.copy(); opt_wk.index = opt_wk.index.to_period("W").to_timestamp()
    tp_wk2 = tp_wk.copy(); tp_wk2.index = tp_wk2.index.to_period("W").to_timestamp()
    Jc = pd.concat({"tp": tp_wk2, "opt": opt_wk}, axis=1, join="inner").dropna()
    corr = float(Jc["tp"].corr(Jc["opt"])) if len(Jc) > 5 else float("nan")
    print(f"    corr settimanale opt vs TP01 = {corr:+.2f}  (atteso ~0.2)")
    for w in (0.3, 0.5):
        comb = (1 - w) * Jc["tp"] + w * Jc["opt"]
        mt = m_weekly(Jc["tp"]); mc = m_weekly(comb)
        print(f"    TP01 {1-w:.0%} + OPT {w:.0%}: Sh {mc['sh']:.2f} (TP01-solo {mt['sh']:.2f}) DD {mc['dd']*100:.0f}%")
    print("\n  NB onesto: short-vol -> guarda i worst-weeks e gli anni di crash. Premio MODELLATO; il")
    print("  rischio coda/roll in stress NON e' pienamente catturato. Lead, non deploy.")


if __name__ == "__main__":
    main()