#!/usr/bin/env python
"""Final robustness: outlier sensitivity, side symmetry by era, overnight extension."""
import os
import numpy as np
import pandas as pd

HERE = os.path.dirname(os.path.abspath(__file__))
sess = pd.read_parquet(os.path.join(HERE, "sessions.parquet"))
sess = sess.dropna(subset=["p1500", "p1530", "p1600", "prev_p1600"]).copy()
sess["year"] = sess["date"].dt.year
sess = sess.set_index("date").sort_index()

r_sig = np.log(sess["p1530"] / sess["p1500"])
label = np.log(sess["p1600"] / sess["p1530"])
pts = sess["p1600"] - sess["p1530"]
r_day = np.log(sess["p1530"] / sess["prev_p1600"])
# overnight: this close -> next session's 09:30 open
next_open = sess["p0930"].shift(-1)
r_on = np.log(next_open / sess["p1600"])
pts_on = next_open - sess["p1600"]
COST = 1.0

mask = np.abs(r_sig) > 0.003
side = np.sign(r_sig)

def line(pl_r, pl_p, name):
    n = len(pl_r)
    if n < 20:
        return f"{name:52s} n={n} (too few)"
    m = pl_r.mean(); t = m / (pl_r.std(ddof=1) / np.sqrt(n))
    med = np.median(pl_p)
    return (f"{name:52s} n={n:4d} {m*1e4:+6.2f}bps t={t:+5.2f} hit={(pl_r>0).mean():.3f} "
            f"net={pl_p.mean()-COST:+6.2f}pt med={med:+5.1f}pt")

d_r, d_p = (side * label)[mask], (side * pts)[mask]
print("=== outlier sensitivity (|r|>0.3%) ===")
print(line(d_r, d_p, "raw"))
w = d_r.clip(d_r.quantile(0.01), d_r.quantile(0.99))
wp = d_p.clip(d_p.quantile(0.01), d_p.quantile(0.99))
print(line(w, wp, "winsorized 1/99"))
trim_mask = (d_r > d_r.quantile(0.02)) & (d_r < d_r.quantile(0.98))
print(line(d_r[trim_mask], d_p[trim_mask], "trimmed 2/98"))

print("\n=== side symmetry by era (|r|>0.3%) ===")
for lo, hi in ((2010, 2014), (2015, 2019), (2020, 2022), (2023, 2026)):
    era = mask & (sess.year >= lo) & (sess.year <= hi)
    for s, nm in ((1, "long"), (-1, "short")):
        m2 = era & (side == s)
        print(line((side * label)[m2], (side * pts)[m2], f"{lo}-{hi} {nm}s"))

print("\n=== overnight extension: hold signal direction 16:00 -> next 09:30 ===")
m = mask & r_on.notna()
print(line((side * r_on)[m], (side * pts_on)[m], "momentum days: hold overnight"))
m2 = m & (np.sign(r_day) == side)
print(line((side * r_on)[m2], (side * pts_on)[m2], "momentum+agree days: hold overnight"))
# and unconditional overnight (sanity: known long drift?)
print(line(r_on.dropna(), pts_on.dropna(), "all days: long overnight (sanity)"))

print("\n=== final candidate rule, per-era summary ===")
# RULE: |r_1500_1530|>0.3% AND sign agrees with day move AND day_range>1.5%... but 2025 broke it.
# Alternative final: same but LONG side only in 2023+? -> report both eras for honesty
day_range_pct = (sess["day_hi_1530"] - sess["day_lo_1530"]) / sess["p1530"]
ruleB = mask & (np.sign(r_day) == side) & (day_range_pct > 0.015)
print(line((side * label)[ruleB], (side * pts)[ruleB], "rule B all years"))
for lo, hi in ((2010, 2019), (2020, 2026)):
    m3 = ruleB & (sess.year >= lo) & (sess.year <= hi)
    print(line((side * label)[m3], (side * pts)[m3], f"rule B {lo}-{hi}"))

print("\n=== trade count reality check ===")
print(f"base signal days/yr: {mask.sum() / sess.year.nunique():.1f}")
print(f"rule B days/yr: {ruleB.sum() / sess.year.nunique():.1f}")
print(f"2026 YTD base: n={int((mask & (sess.year==2026)).sum())}, ruleB n={int((ruleB & (sess.year==2026)).sum())}")