PythagorasGoal/scripts/research/r0725_statarb_eq.py

"""STATARB-EQ — il meccanismo di STATARB-RESID su coppie ETF con 20-30 ANNI di storia.

PERCHE'. STATARB-RESID (relative-momentum del residuo, W=45, sgn=+1) e' il miglior lead del
progetto e ha un gate di deploy al 2026-09-27, ma la sua debolezza non e' il segno dei numeri:
e' la STATISTICA. Su crypto ha ~2,6 anni di storia; il test multi-coppia su Hyperliquid
(r0725_statarb_multi.py) ha mostrato che le 50 coppie alt/BTC valgono solo ~4,5 scommesse
indipendenti (corr media 0,204), che il paniere fa Sharpe 0,82 ma con IC95% [-0,12, 1,72] —
cioe' NON distinguibile da zero — e che il t apparente 5,05 e' gonfiato dalla dipendenza.

Le coppie ETF risolvono esattamente questo: 20-30 anni di storia certificata (IB ADJUSTED_LAST,
gia' su disco), classi d'attivo diverse fra loro (quindi ampiezza effettiva vera), e 3-4 regimi
di mercato completi (dot-com, GFC, ZIRP, 2022) invece di un solo mini-ciclo alt.

DOMANDA: il "relative-momentum del residuo" e' un fenomeno reale e persistente, o e' un artefatto
della finestra crypto 2024-2026?

DISCIPLINA
  * Meccanismo CONGELATO, importato da r0725_statarb_multi.py: W=45, sgn=+1, beta OLS rolling
    causale, z-score su W, tanh, vol-target 20%, cap 2x. ZERO rifit, zero griglia.
  * Coppie definite A PRIORI dentro la stessa classe d'attivo (relazione economica, non
    data-mining di cointegrazione: cercare le coppie piu' cointegrate SU QUESTI STESSI DATI
    sarebbe selezione). 12 coppie, elencate sotto con la loro ragione economica.
  * Costi: 2 bps/lato per gamba (ETF liquidi) PIU' un costo di prestito sullo short, testato a
    0 / 30 / 100 bps annui — sui long/short azionari il borrow non e' un dettaglio.
  * Null: (a) statica CAUSALE (segno = media espandente del segnale), (b) permutazione a blocchi.
  * Stabilita' per DECADE: e' il test che i 2,6 anni crypto non possono fare.

    uv run python scripts/research/r0725_statarb_eq.py
"""
from __future__ import annotations

import sys
from pathlib import Path

import numpy as np
import pandas as pd

ROOT = Path(__file__).resolve().parents[2]
sys.path.insert(0, str(ROOT))
sys.path.insert(0, str(ROOT / "scripts" / "research"))

from r0725_statarb_multi import (BLOCK, CAP, N_PERM, SEED, TARGET_VOL, VOL_WIN, W, _dd,
                                 _zscore, SGN)
from eqlib import load_eq

ANN_EQ = np.sqrt(252.0)
FEE_LEG = 0.0002                       # 2 bps/lato per gamba (ETF liquidi su IB)
BORROW_GRID = (0.0, 0.0030, 0.0100)    # costo prestito annuo sulla gamba corta

# --- 12 coppie A PRIORI: stessa classe d'attivo, relazione economica dichiarata ------------------
PAIRS = [
    ("DIA", "SPY", "US large cap: due proxy dello stesso indice"),
    ("IWM", "SPY", "small vs large cap US (size)"),
    ("QQQ", "SPY", "tech vs mercato ampio US"),
    ("EEM", "EFA", "emergenti vs sviluppati ex-US"),
    ("EWJ", "EFA", "Giappone vs sviluppati ex-US"),
    ("FXI", "EEM", "Cina vs emergenti"),
    ("EFA", "SPY", "ex-US vs US (home bias)"),
    ("SLV", "GLD", "argento vs oro (metalli preziosi)"),
    ("IEF", "TLT", "duration intermedia vs lunga (curva)"),
    ("HYG", "LQD", "credito HY vs IG (rischio di credito)"),
    ("LQD", "AGG", "credito IG vs aggregato (spread)"),
    ("USO", "DBC", "petrolio vs paniere commodity"),
]


def signal_eq(base_px: pd.Series, tgt_px: pd.Series):
    """Meccanismo CONGELATO, identico a r0725_statarb_multi.signal ma su calendario di borsa."""
    x, y = np.log(base_px.values), np.log(tgt_px.values)
    sx, sy = pd.Series(x), pd.Series(y)
    mx = sx.rolling(W, min_periods=W).mean()
    my = sy.rolling(W, min_periods=W).mean()
    cov = (sx * sy).rolling(W, min_periods=W).mean() - mx * my
    var = (sx * sx).rolling(W, min_periods=W).mean() - mx * mx
    beta = cov / var.replace(0, np.nan)
    resid = (sy - (my - beta * mx) - beta * sx).values
    z = _zscore(np.nan_to_num(resid), W)
    d = SGN * np.tanh(np.nan_to_num(z))

    r_b = np.concatenate([[0.0], np.diff(base_px.values) / base_px.values[:-1]])
    r_t = np.concatenate([[0.0], np.diff(tgt_px.values) / tgt_px.values[:-1]])
    spread = r_t - r_b
    vol = pd.Series(spread).rolling(VOL_WIN, min_periods=VOL_WIN).std().values * ANN_EQ
    scal = np.where((vol > 0) & np.isfinite(vol), TARGET_VOL / vol, 0.0)
    pos = np.clip(np.nan_to_num(d) * np.nan_to_num(scal), -CAP, CAP)
    pos[~np.isfinite(pos)] = 0.0
    return pos, spread, np.nan_to_num(scal)


def pnl_eq(pos: np.ndarray, spread: np.ndarray, borrow: float = 0.0) -> np.ndarray:
    held = np.concatenate([[0.0], pos[:-1]])
    turn = np.abs(np.diff(held, prepend=0.0))
    # il prestito si paga sulla gamba corta: |posizione| di nozionale a prestito, ogni giorno
    return held * spread - 2.0 * FEE_LEG * turn - borrow / 252.0 * np.abs(held)


def _sh(r: np.ndarray) -> float:
    r = r[np.isfinite(r)]
    return float(r.mean() / r.std() * ANN_EQ) if len(r) > 60 and r.std() > 0 else 0.0


def main() -> None:
    print("=" * 108)
    print("  STATARB-EQ — meccanismo congelato (W=45, sgn=+1) su coppie ETF, 20-30 anni di storia")
    print("=" * 108)

    rng = np.random.default_rng(SEED + 3)
    rows = []
    for tgt, base, why in PAIRS:
        try:
            b_all, t_all = load_eq(base)["close"].astype(float), load_eq(tgt)["close"].astype(float)
        except FileNotFoundError:
            print(f"  [{tgt}/{base}] parquet mancante — salto")
            continue
        ix = b_all.index.intersection(t_all.index)
        if len(ix) < 1000:
            continue
        b, t = b_all[ix], t_all[ix]
        pos, spread, scal = signal_eq(b, t)
        r = pnl_eq(pos, spread)

        live = pos[pos != 0.0]
        mono = float(max((live > 0).mean(), (live < 0).mean())) if len(live) else 0.0
        # null statica CAUSALE
        cum, cnt = np.cumsum(pos), np.arange(1, len(pos) + 1)
        run_mean = np.concatenate([[0.0], (cum / cnt)[:-1]])
        sc = np.where(run_mean >= 0, 1.0, -1.0)
        sc[:VOL_WIN + W] = 0.0
        r_stat = pnl_eq(np.clip(sc * scal, -CAP, CAP), spread)
        # null permutazione a blocchi
        nb = int(np.ceil(len(pos) / BLOCK))
        perm = []
        for _ in range(N_PERM):
            order = rng.permutation(nb)
            pp = np.concatenate([pos[k * BLOCK:(k + 1) * BLOCK] for k in order])[:len(pos)]
            perm.append(_sh(pnl_eq(pp, spread)))
        perm = np.array(perm)
        gross = np.concatenate([[0.0], pos[:-1]]) * spread          # senza fee ne' borrow
        turn = np.abs(np.diff(np.concatenate([[0.0], pos[:-1]]), prepend=0.0))
        rows.append(dict(pair=f"{tgt}/{base}", why=why, n=len(ix), idx=ix, ret=r,
                         sh_gross=_sh(gross), turn=float(turn.mean()),
                         sh=_sh(r), dd=_dd(r), mono=mono, sh_stat=_sh(r_stat),
                         up=_sh(r) - _sh(r_stat), pperm=float((perm >= _sh(r)).mean()),
                         yrs=(ix[-1] - ix[0]).days / 365.25))

    if not rows:
        print("  nessuna coppia valutabile")
        return

    print("\n" + "-" * 108)
    print("  (1) COPPIE A PRIORI — meccanismo congelato, netto fee (borrow 0 in questa tabella)")
    print("-" * 108)
    print(f"  {'coppia':<11}{'anni':>6}{'Sh LORDA':>10}{'Sh netta':>10}{'turn/g':>8}{'maxDD':>9}"
          f"{'stat.caus':>11}{'uplift':>8}{'p perm':>9}  ragione economica")
    for d in sorted(rows, key=lambda x: -x["sh_gross"]):
        print(f"  {d['pair']:<11}{d['yrs']:>6.1f}{d['sh_gross']:>10.2f}{d['sh']:>10.2f}"
              f"{d['turn']:>8.2f}{d['dd']*100:>8.1f}%"
              f"{d['sh_stat']:>11.2f}{d['up']:>8.2f}{d['pperm']:>9.3f}  {d['why']}")

    sh = np.array([d["sh"] for d in rows])
    shg = np.array([d["sh_gross"] for d in rows])
    up = np.array([d["up"] for d in rows])
    pp = np.array([d["pperm"] for d in rows])
    n = len(rows)
    print(f"\n  coppie {n}   Sharpe LORDA media {shg.mean():>5.2f}  mediana {np.median(shg):>5.2f}   "
          f"frazione > 0: {(shg>0).mean()*100:.0f}%")
    print(f"              Sharpe NETTA media {sh.mean():>5.2f}  mediana {np.median(sh):>5.2f}   "
          f"frazione > 0: {(sh>0).mean()*100:.0f}%")
    print(f"  turnover medio {np.mean([d['turn'] for d in rows]):.2f}/giorno -> drag di fee "
          f"{np.mean(shg-sh):.2f} Sharpe. LETTURA: se la LORDA e' ~0 l'edge non esiste ed e' inutile")
    print("  cercare il segno opposto (anche lo specchio sarebbe ~0 lordo e negativo netto).")
    print(f"  uplift vs statica causale: media {up.mean():>5.2f}  frazione > 0: {(up>0).mean()*100:.0f}%")
    print(f"  p permutazione < 0.05: {(pp<0.05).mean()*100:.0f}% delle coppie (atteso 5%)")

    # ---------------- paniere + IC
    print("\n" + "-" * 108)
    print("  (2) PANIERE EW DELLE COPPIE — con l'ampiezza e la storia che al crypto mancano")
    print("-" * 108)
    M = pd.concat({d["pair"]: pd.Series(d["ret"], index=d["idx"]) for d in rows},
                  axis=1, sort=True).sort_index()
    bask = M.mean(axis=1, skipna=True).dropna()
    yrs = len(bask) / 252.0
    C = M.corr().values
    off = C[~np.eye(len(C), dtype=bool)]
    rbar = float(np.nanmean(off))
    n_eff = len(C) / (1.0 + (len(C) - 1) * rbar) if rbar > -1 / (len(C) - 1) else float(len(C))
    print(f"  paniere: n={len(bask)} ({yrs:.1f} anni)  Sharpe {_sh(bask.values):>5.2f}  "
          f"maxDD {_dd(bask.values)*100:>5.1f}%")
    print(f"  corr media fra coppie {rbar:>5.3f}  ->  ampiezza EFFETTIVA ~{n_eff:.1f} "
          f"(vs ~4.5 delle 50 coppie crypto)")
    rng2 = np.random.default_rng(SEED + 11)
    v = bask.values
    nb2 = int(np.ceil(len(v) / BLOCK))
    boot = np.array([_sh(np.concatenate([v[i * BLOCK:(i + 1) * BLOCK]
                                         for i in rng2.integers(0, nb2, size=nb2)])[:len(v)])
                     for _ in range(2000)])
    print(f"  block-bootstrap IC95% [{np.percentile(boot,2.5):>5.2f}, {np.percentile(boot,97.5):>5.2f}]"
          f"   P(Sh>0) {(boot>0).mean()*100:.0f}%   t ~ {_sh(v)*np.sqrt(yrs):.2f}")

    # ---------------- stabilita' per decade
    print("\n" + "-" * 108)
    print("  (3) STABILITA' PER DECADE — il test che i 2,6 anni crypto NON possono fare")
    print("-" * 108)
    decs = [("1998-2004", 1998, 2004), ("2005-2011", 2005, 2011),
            ("2012-2018", 2012, 2018), ("2019-2026", 2019, 2026)]
    print(f"  {'periodo':<12}{'Sh paniere':>12}{'maxDD':>9}{'coppie Sh>0':>14}")
    for nm, y0, y1 in decs:
        w = bask[(bask.index.year >= y0) & (bask.index.year <= y1)]
        if len(w) < 200:
            print(f"  {nm:<12}{'(storia corta)':>12}")
            continue
        cnt = sum(1 for d in rows
                  if len(pd.Series(d['ret'], index=d['idx'])[
                      (d['idx'].year >= y0) & (d['idx'].year <= y1)]) > 200
                  and _sh(pd.Series(d['ret'], index=d['idx'])[
                      (d['idx'].year >= y0) & (d['idx'].year <= y1)].values) > 0)
        tot = sum(1 for d in rows if ((d['idx'].year >= y0) & (d['idx'].year <= y1)).sum() > 200)
        print(f"  {nm:<12}{_sh(w.values):>12.2f}{_dd(w.values)*100:>8.1f}%{f'{cnt}/{tot}':>14}")

    # ---------------- costo di prestito
    print("\n" + "-" * 108)
    print("  (4) SENSIBILITA' AL COSTO DI PRESTITO (gamba corta) — su un long/short non e' un dettaglio")
    print("-" * 108)
    print(f"  {'borrow ann':>12}{'Sh paniere':>13}{'coppie Sh>0':>14}")
    for bw in BORROW_GRID:
        rr = {}
        for tgt, base, _ in PAIRS:
            d = next((x for x in rows if x["pair"] == f"{tgt}/{base}"), None)
            if d is None:
                continue
            b_all, t_all = load_eq(base)["close"].astype(float), load_eq(tgt)["close"].astype(float)
            ix = b_all.index.intersection(t_all.index)
            pos, spread, _ = signal_eq(b_all[ix], t_all[ix])
            rr[d["pair"]] = pd.Series(pnl_eq(pos, spread, borrow=bw), index=ix)
        bb = pd.concat(rr, axis=1, sort=True).sort_index().mean(axis=1, skipna=True).dropna()
        pos_cnt = sum(1 for s in rr.values() if _sh(s.values) > 0)
        print(f"  {bw*100:>11.2f}%{_sh(bb.values):>13.2f}{f'{pos_cnt}/{len(rr)}':>14}")

    print("\n" + "=" * 108)


if __name__ == "__main__":
    main()