#!/usr/bin/env python3
"""PEAD gate: does post-earnings-announcement drift survive in small/neglected names, net of cost?

The strongest classic anomaly that persists BECAUSE it's too small/illiquid for funds to arb.
Proxy (no earnings-date data needed): an "earnings-like surprise" = a large abnormal-volume single
-day move (|ret| > K_SIG x trailing-vol AND volume > K_VOL x trailing-avg). PEAD predicts
CONTINUATION: enter at the event-day CLOSE (leak-free, after the full move), hold N days, signed by
the surprise direction. Measure net-of-cost drift, t-stat, long vs short, IS/OOS, by liquidity band.
"""
import math
import os
import sys

import numpy as np

sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
from equity_factor_gate import load, roll  # noqa: E402

K_SIG, K_VOL = 2.0, 2.0          # surprise = |ret|>2sigma AND volume>2x average
LO, HI = 1e6, 30e6               # small-but-tradeable $/day band (where PEAD survives)
HORIZONS = [1, 3, 5, 10, 20]


def tstat(x):
    x = x[np.isfinite(x)]
    return float(x.mean() / (x.std() / math.sqrt(len(x)))) if len(x) > 30 and x.std() > 0 else float("nan")


def main():
    insts, days, close, dvol = load()
    T, N = close.shape
    lc = np.log(close)
    R = np.zeros((T, N)); R[1:] = lc[1:] - lc[:-1]; R = np.where(np.isfinite(R), R, 0.0)
    vol20 = roll(np.std, R, 20)                              # trailing std (excludes t -> causal)
    advol = roll(np.mean, np.nan_to_num(dvol), 20)
    year = (1970 + days / 365.25).astype(int)
    rt_cost = np.clip(40.0 / np.sqrt(np.maximum(advol, 1.0) / 1e6), 5.0, 80.0) / 1e4   # round-trip, illiquidity-scaled

    valid = (np.isfinite(close) & (vol20 > 0) & np.isfinite(advol) & (advol > LO) & (advol < HI))
    surprise = valid & (np.abs(R) > K_SIG * vol20) & (dvol > K_VOL * advol)
    surprise[:21] = False; surprise[T - max(HORIZONS):] = False
    sgn = np.sign(R)
    print(f"PEAD gate (DBEQ, band ${LO/1e6:.0f}-{HI/1e6:.0f}M/day, surprise=|ret|>{K_SIG}sig & vol>{K_VOL}x)")
    print(f"  total surprise events: {int(surprise.sum())}  (avg cost {1e4*rt_cost[surprise].mean():.0f}bp round-trip)")

    oos = np.repeat((year >= 2025)[:, None], N, 1)
    up = surprise & (sgn > 0); dn = surprise & (sgn < 0)
    print(f"\n{'horizon':>8} {'n':>7} {'gross%':>7} {'net%':>7} {'t(net)':>7} {'fracpos':>8} | {'LONG net%':>9} {'SHORT net%':>10} | {'OOS net%':>9}")
    for h in HORIZONS:
        dr = np.full((T, N), np.nan); dr[:T - h] = lc[h:] - lc[:T - h]   # dr[t] = lc[t+h]-lc[t]
        signed = sgn * dr
        net = signed - rt_cost
        g = signed[surprise]; nt = net[surprise]
        ln = net[up]; sh = net[dn]
        oo = net[surprise & oos]
        print(f"{h:>8} {int(np.isfinite(nt).sum()):>7} {100*np.nanmean(g):>+7.2f} {100*np.nanmean(nt):>+7.2f} "
              f"{tstat(nt):>+7.1f} {np.nanmean(nt>0):>8.2f} | {100*np.nanmean(ln):>+9.2f} {100*np.nanmean(sh):>+10.2f} | {100*np.nanmean(oo):>+9.2f}")

    # by surprise strength at the 10d horizon
    h = 10
    dr = np.full((T, N), np.nan); dr[:T - h] = lc[h:] - lc[:T - h]
    net = sgn * dr - rt_cost
    strong = surprise & (np.abs(R) > 3 * vol20)
    print(f"\n  10d net by surprise strength:  2-3sig {100*np.nanmean(net[surprise & ~strong]):+.2f}%  "
          f">3sig {100*np.nanmean(net[strong]):+.2f}%  (PEAD: bigger surprise -> bigger drift)")
    print("\nVERDICT: net drift > 0 with t(net) > 3 (deflated) AND positive OOS AND long-side works = real PEAD edge")
    print("worth a proper book. If net ~0 or t<2 or OOS collapses = eaten by cost / arbed even in small names.")


if __name__ == "__main__":
    main()