PythagorasGoal/scripts/research/r0725_statarb_demean_skeptic.py

"""SCETTICO su STATARB-DEMEAN (V3) — prima di crederci.

IL CANDIDATO. `r0725_statarb_basket_gate.py` ha prodotto la prima cosa che somiglia a uno sleeve
nuovo da molte ondate: posizioni del meccanismo congelato (W=45, sgn=+1) sulle 50 coppie alt/BTC,
DEMEANATE cross-sezionalmente ogni giorno. Sharpe 1.79, maxDD -2.5%, ampiezza effettiva 37.4,
marginale ADDS, deflated-Sharpe 0.985 (PASS), corr ~0 a tutti e 5 gli sleeve attivi.

Numeri cosi' puliti sono, nella storia di questo progetto, il momento in cui si e' quasi sempre
scoperto un artefatto (feed testnet, look-ahead ffill, CC01 Sharpe 11, celle 0-perdite). Quindi
prima di scrivere "trovata", si prova a UCCIDERLA.

⚠ OSSERVAZIONE STRUTTURALE CHE MOTIVA IL TEST PRINCIPALE. Demeanare annulla algebricamente la
gamba BTC:  sum_i (p_i - p̄)(r_i - r_BTC) = sum_i (p_i - p̄) r_i , perche' sum_i (p_i - p̄) = 0.
Quindi V3 NON e' piu' un paniere di coppie: e' una strategia CROSS-SECTIONAL sui 50 alt con pesi
(p_i - p̄). Va quindi testata come tale, e in particolare contro il sospetto che il merito sia
della COSTRUZIONE (demean + vol-target + 50 asset) e non del SEGNALE.

I TEST (un candidato deve sopravvivere a tutti):
  T1 NULL DI PERMUTAZIONE CROSS-SEZIONALE — ogni giorno si permutano le etichette-asset del vettore
     di posizione. Distribuzione di posizioni, demeaning, vol-target, costi: tutto identico; sparisce
     SOLO l'abbinamento segnale<->asset. Se il candidato non batte questo null, l'edge e' costruzione.
  T2 NULL CASUALE — posizioni casuali tanh(N(0,1)) con la stessa scalatura e lo stesso demeaning.
  T3 NULL DI PERMUTAZIONE TEMPORALE a blocchi (20g) — rompe l'allineamento nel tempo.
  T4 LAG — ritardare la posizione di 1 e 2 giorni. Un segnale lento degrada dolcemente; un
     look-ahead crolla a zero immediatamente.
  T5 SOTTOPERIODI anno per anno — 2.6 anni sono un regime solo (bear alt): dove vive il rendimento?
  T6 COSTI — la fee attuale addebita 2 gambe per coppia, ma dopo il demeaning la gamba BTC si
     annulla: e' quindi SOVRASTIMATA. Si verifica comunque a 0.05 / 0.10 / 0.20% per gamba.
  T7 RIDONDANZA con un momentum cross-sectional semplice sugli stessi 50 (e' XS01 travestito?).
  T8 STRUTTURA — esposizione lorda/netta, bilanciamento long/short, turnover, gambe che si muovono.

    uv run python scripts/research/r0725_statarb_demean_skeptic.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 FEE_LEG, SEED, _dd, _sh
from r0725_statarb_basket_gate import pair_frames

ANN = np.sqrt(365.0)
N_NULL = 300
BLOCK = 20


def ret_from_pos(P: pd.DataFrame, S: pd.DataFrame, fee_leg: float = FEE_LEG,
                 lag: int = 0, demean: bool = True) -> pd.Series:
    """Rendimento EW del paniere da una matrice di posizioni [data x asset]. `lag` giorni in piu'."""
    Q = P.sub(P.mean(axis=1), axis=0) if demean else P
    q = Q.to_numpy(float)
    s = S.to_numpy(float)
    held = np.zeros_like(q)
    k = 1 + lag
    if k < len(q):
        held[k:] = q[:-k]
    turn = np.abs(np.diff(held, axis=0, prepend=np.zeros((1, q.shape[1]))))
    net = held * s - 2.0 * fee_leg * turn
    return pd.Series(np.nanmean(net, axis=1), index=P.index).dropna()


def main() -> None:
    print("=" * 100)
    print("  SCETTICO — STATARB-DEMEAN (V3): si riesce a ucciderla?")
    print("=" * 100)

    P, S = pair_frames()
    P = P.fillna(0.0)
    S = S.reindex(P.index)[P.columns]
    base = ret_from_pos(P, S)
    sh0 = _sh(base.values)
    print(f"\n  candidato: {P.shape[1]} asset, {len(P)} barre "
          f"({P.index[0].date()} -> {P.index[-1].date()})")
    print(f"  Sharpe {sh0:.2f}   maxDD {_dd(base.values)*100:.1f}%   "
          f"vol ann {base.std()*ANN*100:.2f}%   ret ann {base.mean()*365*100:.2f}%")

    rng = np.random.default_rng(SEED)
    q = P.to_numpy(float)

    # ⚠ I NULL SI CONFRONTANO A FEE ZERO. Permutare le posizioni ne fa esplodere il turnover
    # (un asset riceve ogni giorno una posizione scorrelata dalla precedente), quindi un null a
    # fee piene sarebbe battuto per motivi di COSTO e non di segnale: darebbe un p-value
    # artificialmente trionfale. A fee zero il confronto isola l'informazione.
    sh0_gross = _sh(ret_from_pos(P, S, fee_leg=0.0).values)
    print(f"  Sharpe LORDA (fee 0) {sh0_gross:.2f}  <- e' questa che i null devono battere")

    # ---------------- T1 permutazione cross-sezionale
    print("\n" + "-" * 100)
    print("  T1  NULL DI PERMUTAZIONE CROSS-SEZIONALE — il test decisivo")
    print("      (stessa distribuzione di posizioni, stesso demeaning, stesso vol-target e costi;")
    print("       si rompe SOLO l'abbinamento segnale<->asset)")
    print("-" * 100)
    nulls = []
    for _ in range(N_NULL):
        qq = np.take_along_axis(q, rng.permuted(np.tile(np.arange(q.shape[1]), (q.shape[0], 1)),
                                                axis=1), axis=1)
        nulls.append(_sh(ret_from_pos(pd.DataFrame(qq, index=P.index, columns=P.columns), S, fee_leg=0.0).values))
    nulls = np.array(nulls)
    p1 = float((nulls >= sh0_gross).mean())
    print(f"  null: media {nulls.mean():>6.2f}  mediana {np.median(nulls):>6.2f}  "
          f"p95 {np.percentile(nulls,95):>6.2f}  max {nulls.max():>6.2f}")
    print(f"  candidato LORDO {sh0_gross:.2f}  ->  p = {p1:.4f}   "
          f"{'SOPRAVVIVE' if p1 < 0.05 else '*** UCCISO: l edge e la COSTRUZIONE, non il segnale ***'}")

    # ---------------- T2 posizioni casuali
    print("\n" + "-" * 100)
    print("  T2  NULL CASUALE — posizioni tanh(N(0,1)) con stessa scala, demeaning e costi")
    print("-" * 100)
    scale = np.abs(q).mean()
    r2 = []
    for _ in range(N_NULL):
        qq = np.tanh(rng.normal(0, 1, size=q.shape)) * scale / max(np.abs(np.tanh(1.0)), 1e-9)
        qq[q == 0.0] = 0.0                      # stessa maschera di "non ancora attivo"
        r2.append(_sh(ret_from_pos(pd.DataFrame(qq, index=P.index, columns=P.columns), S, fee_leg=0.0).values))
    r2 = np.array(r2)
    p2 = float((r2 >= sh0_gross).mean())
    print(f"  null: media {r2.mean():>6.2f}  p95 {np.percentile(r2,95):>6.2f}  max {r2.max():>6.2f}"
          f"   ->  p = {p2:.4f}  {'SOPRAVVIVE' if p2 < 0.05 else 'UCCISO'}")

    # ---------------- T3 permutazione temporale a blocchi
    print("\n" + "-" * 100)
    print("  T3  NULL DI PERMUTAZIONE TEMPORALE (blocchi 20g)")
    print("-" * 100)
    nb = int(np.ceil(len(q) / BLOCK))
    r3 = []
    for _ in range(N_NULL):
        order = rng.permutation(nb)
        qq = np.concatenate([q[i * BLOCK:(i + 1) * BLOCK] for i in order])[:len(q)]
        r3.append(_sh(ret_from_pos(pd.DataFrame(qq, index=P.index, columns=P.columns), S, fee_leg=0.0).values))
    r3 = np.array(r3)
    p3 = float((r3 >= sh0_gross).mean())
    print(f"  null: media {r3.mean():>6.2f}  p95 {np.percentile(r3,95):>6.2f}   ->  p = {p3:.4f}"
          f"  {'SOPRAVVIVE' if p3 < 0.05 else 'UCCISO'}")

    # ---------------- T4 lag
    print("\n" + "-" * 100)
    print("  T4  LAG — un segnale lento degrada dolcemente, un look-ahead crolla")
    print("-" * 100)
    for lag in (0, 1, 2, 3):
        r = ret_from_pos(P, S, lag=lag)
        print(f"    lag +{lag}g:  Sharpe {_sh(r.values):>5.2f}   maxDD {_dd(r.values)*100:>5.1f}%")

    # ---------------- T5 sottoperiodi
    print("\n" + "-" * 100)
    print("  T5  SOTTOPERIODI — 2.6 anni sono UN regime; dove vive il rendimento?")
    print("-" * 100)
    for y in sorted(set(base.index.year)):
        w = base[base.index.year == y]
        if len(w) < 40:
            continue
        print(f"    {y}:  n={len(w):>4}  Sharpe {_sh(w.values):>5.2f}  "
              f"ret {(np.prod(1+w.values)-1)*100:>6.2f}%  maxDD {_dd(w.values)*100:>5.1f}%")
    half = len(base) // 2
    print(f"    prima meta' Sharpe {_sh(base.values[:half]):>5.2f}   "
          f"seconda meta' {_sh(base.values[half:]):>5.2f}")

    # ---------------- T6 costi
    print("\n" + "-" * 100)
    print("  T6  COSTI — nota: dopo il demeaning la gamba BTC si annulla, quindi addebitare")
    print("      2 gambe per coppia SOVRASTIMA i costi. Si verifica comunque al rialzo.")
    print("-" * 100)
    for f in (0.0, 0.0005, 0.0010, 0.0020, 0.0040):
        r = ret_from_pos(P, S, fee_leg=f)
        print(f"    fee {f*100:>5.2f}%/gamba:  Sharpe {_sh(r.values):>5.2f}   "
              f"ret ann {r.mean()*365*100:>6.2f}%")

    # ---------------- T7 ridondanza con XS momentum semplice
    print("\n" + "-" * 100)
    print("  T7  RIDONDANZA — e' un momentum cross-sectional travestito?")
    print("-" * 100)
    from r0725_statarb_multi import load_hl, universe
    px = {}
    for s_ in universe():
        try:
            px[s_] = load_hl(s_)
        except FileNotFoundError:
            pass
    PX = pd.concat(px, axis=1, sort=True).reindex(P.index).ffill()
    for L in (30, 45, 90):
        mom = np.log(PX / PX.shift(L))
        z = mom.sub(mom.mean(axis=1), axis=0).div(mom.std(axis=1).replace(0, np.nan), axis=0)
        w = np.tanh(z.fillna(0.0))
        w = w.sub(w.mean(axis=1), axis=0)
        rx = ret_from_pos(w[P.columns].fillna(0.0), S, demean=False)
        j = pd.concat({"a": base, "b": rx}, axis=1, sort=True).dropna()
        print(f"    XS-mom L={L:>3}: Sharpe {_sh(rx.values):>5.2f}   "
              f"corr al candidato {j['a'].corr(j['b']):>6.3f}")

    # ---------------- T8 struttura
    print("\n" + "-" * 100)
    print("  T8  STRUTTURA — esposizione, bilanciamento, turnover, gambe operative")
    print("-" * 100)
    Q = P.sub(P.mean(axis=1), axis=0)
    gross = Q.abs().sum(axis=1)
    net = Q.sum(axis=1)
    nlong = (Q > 0).sum(axis=1)
    turn = Q.diff().abs().sum(axis=1)
    print(f"    esposizione LORDA media {gross.mean():>6.2f}  (mediana {gross.median():.2f})")
    print(f"    esposizione NETTA media {net.mean():>8.2e}  (deve essere ~0 per costruzione)")
    print(f"    gambe long/giorno  media {nlong.mean():>5.1f} su {P.shape[1]}  "
          f"(bilanciamento atteso ~50%)")
    print(f"    turnover lordo/giorno {turn.mean():>6.2f}  = {turn.mean()/max(gross.mean(),1e-9)*100:>5.1f}% del lordo")
    n_ass = max(P.shape[1], 1)
    # il rendimento e' la MEDIA sugli asset -> il peso di portafoglio della gamba i e' held_i/N.
    # Quindi il lordo del paniere e' gross/N, e il nozionale per gamba e' cap*|held_i|/N.
    gross_frac = gross.mean() / n_ass
    print(f"    lordo del paniere {gross_frac*100:>5.1f}% del capitale  (leva {gross_frac:.2f}x)")
    for cap in (600, 2000, 5000, 20000):
        tk = cap * gross_frac / n_ass
        print(f"    cap ${cap:>6}: ticket medio per gamba ${tk:>7.2f}  "
              f"{'-> sotto il min-order $5: NON eseguibile' if tk < 5 else '-> sopra il min-order'}")

    print("\n" + "=" * 100)
    verdict = (p1 < 0.05 and p2 < 0.05 and p3 < 0.05)
    print(f"  ESITO SCETTICO: {'i tre null NON la uccidono' if verdict else 'UCCISA da almeno un null'}")
    print("=" * 100)


if __name__ == "__main__":
    main()