"""Fase 2 raffinata: extra-rendimento vs baseline + significativita'.

Per ogni (coppia, direzione) calcola:
  - rendimento medio grezzo a 1/3/5/10/20 g
  - EXTRA-rendimento = evento - baseline nello stesso regime (toglie il bias di periodo)
  - effetto in sigma del movimento tipico a h giorni
  - quota di eventi che battono la baseline
  - p-value (block bootstrap) e q-value (Benjamini-Hochberg)

Uso:
  python excess_analysis.py [SYMBOL] [YEAR_MIN]
Output:
  ../results/summary_<SYM>_<YEAR>.csv  + stampa dei top a 5/10/20 g
"""
import os, sys, numpy as np, pandas as pd
from _common import load_pair, add_regime, regime_edges, block_bootstrap_p, benjamini_hochberg, HZ, HERE

SYM = sys.argv[1] if len(sys.argv) > 1 else "EURUSD"
YEAR = int(sys.argv[2]) if len(sys.argv) > 2 else 1999
RESULTS = os.path.normpath(os.path.join(HERE, "..", "results"))

cr, ba = load_pair(SYM, YEAR)
print(f"[{SYM} {YEAR}+] incroci={len(cr)}  baseline={len(ba)}")

edges = regime_edges(ba)
cr = add_regime(cr, edges)
ba = add_regime(ba, edges)

base_mean = {h: ba.groupby("regime")[f"cret_{h}"].mean() for h in HZ}
base_glob = {h: ba[f"cret_{h}"].mean() for h in HZ}
base_sd   = {h: ba[f"cret_{h}"].std()  for h in HZ}

rows = []
for (pair, d), g in cr.groupby(["pair", "dir"]):
    n = len(g)
    if n < 30:
        continue
    rec = {"pair": pair, "dir": int(d), "n": n}
    for h in HZ:
        ev = g[f"cret_{h}"].values
        bexp = g["regime"].map(base_mean[h]).fillna(base_glob[h]).values
        exc = ev - bexp
        rec[f"raw_{h}"] = np.nanmean(ev)
        rec[f"exc_{h}"] = np.nanmean(exc)
        rec[f"eff_{h}"] = np.nanmean(exc) / base_sd[h]
        rec[f"pos_{h}"] = (exc > 0).mean()
        rec[f"p_{h}"]   = block_bootstrap_p(exc)
    rows.append(rec)

res = pd.DataFrame(rows)
for h in HZ:
    res[f"q_{h}"] = benjamini_hochberg(res[f"p_{h}"].values)

os.makedirs(RESULTS, exist_ok=True)
outp = os.path.join(RESULTS, f"summary_{SYM}_{YEAR}.csv")
res.to_csv(outp, index=False)

print(f"vol giornaliera prezzo (sd cret_1) = {ba['cret_1'].std():.3f}%")
print(f"baseline drift 5/10/20g = {base_glob[5]:.3f}% / {base_glob[10]:.3f}% / {base_glob[20]:.3f}%")
for h in [5, 10, 20]:
    sig = res[res[f"q_{h}"] < 0.10].copy()
    sig["ae"] = sig[f"eff_{h}"].abs()
    sig = sig.sort_values("ae", ascending=False).head(12)
    print(f"\n=== TOP {h}g (q<0.10) ===")
    for _, r in sig.iterrows():
        print(f"  {r['pair']:13s} dir={int(r['dir']):+d} n={int(r['n']):5d} "
              f"exc={r[f'exc_{h}']:+.3f}% eff={r[f'eff_{h}']:+.2f}sd "
              f"pos={r[f'pos_{h}']*100:.0f}% q={r[f'q_{h}']:.3f}")
print(f"\nsalvato: {outp}")