#!/usr/bin/env python3
"""Lower-frequency directional IC on ES OHLCV-1m bars (minutes → days).

After measuring that seconds-horizon ES is structurally unprofitable for a
non-colocated participant (crossing edge ~100x too small vs spread; passive MM
adversely selected — pearl_ofi_edge_uncapturable_by_crossing,
pearl_passive_mm_knife_edge_ofi_conditioning_promising), this tests the remaining
thesis: at minutes-to-days horizons the directional move DWARFS the 1-tick spread,
so crossing cost stops being the killer. The question becomes purely: does a cheap
directional predictor have positive OUT-OF-SAMPLE Information Coefficient?

Tests two regimes on front-month ES 1-min OHLCV (2024-Q1 .. 2026-Q1, ~2y):
  * INTRADAY: momentum/reversion (trailing return over L minutes) vs forward
    return over H minutes; contiguity-guarded (no windows spanning session/roll gaps).
  * DAILY: time-series momentum (Moskowitz-Ooi-Pedersen) — trailing return over
    L days vs forward return over H days.
For each predictor×horizon: Pearson/Spearman IC, sign-accuracy, Newey-West t-stat,
OOS IC (train 2024 / test 2025+), and median |forward move| in ticks (to confirm
moves >> the 1-tick = $12.50 crossing cost).

Usage: python3 scripts/measure_lowfreq_ic.py [dir] [--self-test]
"""
import argparse
import glob
import sys

import numpy as np

TICK = 0.25
PX_SCALE = 1e9


def newey_west_tstat(x, y, lags):
    n = len(x)
    if n < lags + 3:
        return float("nan")
    X = np.column_stack([np.ones(n), x])
    xtx_inv = np.linalg.inv(X.T @ X)
    beta = xtx_inv @ (X.T @ y)
    u = (y - X @ beta)[:, None] * X
    S = u.T @ u
    for l in range(1, lags + 1):
        w = 1.0 - l / (lags + 1.0)
        G = u[l:].T @ u[:-l]
        S += w * (G + G.T)
    cov = xtx_inv @ S @ xtx_inv
    se = np.sqrt(cov[1, 1])
    return float(beta[1] / se) if se > 0 else float("nan")


def load_ohlcv(d):
    import zstandard
    import databento_dbn as dbn
    files = sorted(glob.glob(f"{d}/*.dbn.zst"))
    ts_all, close_all, q_all = [], [], []
    for qi, f in enumerate(files):
        dec = dbn.DBNDecoder(); dctx = zstandard.ZstdDecompressor()
        ts, cl, inst = [], [], []
        with open(f, "rb") as fh:
            r = dctx.stream_reader(fh)
            while True:
                c = r.read(1 << 20)
                if not c:
                    break
                dec.write(c)
                for rec in dec.decode():
                    if "OHLCV" not in type(rec).__name__:
                        continue
                    ts.append(rec.ts_event); cl.append(rec.close); inst.append(rec.instrument_id)
        ts = np.array(ts, np.int64); cl = np.array(cl, np.int64); inst = np.array(inst, np.int32)
        if len(ts) == 0:
            print(f"  {f.split('/')[-1]}: EMPTY (no OHLCV) — skipped")
            continue
        ids, cnts = np.unique(inst, return_counts=True)
        dom = ids[cnts.argmax()]
        m = inst == dom
        order = np.argsort(ts[m], kind="stable")
        ts_all.append(ts[m][order]); close_all.append(cl[m][order] / PX_SCALE)
        q_all.append(np.full(m.sum(), qi))
        print(f"  {f.split('/')[-1]}: {m.sum()} front-month 1m bars (dom={dom})")
    return np.concatenate(ts_all), np.concatenate(close_all), np.concatenate(q_all)


def ic_block(name, ts, close, q, lookbacks, horizons, bar_unit_s, oos_split_ts):
    from scipy.stats import pearsonr, spearmanr
    logc = np.log(close)
    n = len(close)
    print(f"\n=== {name}  ({n} bars) ===")
    print(f"{'L':>4} {'H':>4} {'n':>7} {'IC_pear':>9} {'IC_spear':>9} {'sign%':>7} "
          f"{'NW_t':>7} {'OOS_IC':>8} {'|fwd|ticks':>11}")
    for L in lookbacks:
        for H in horizons:
            # predictor: trailing return over L bars; forward: return over H bars
            # contiguity: bar index window must not span a quarter (roll) boundary
            idx = np.arange(L, n - H)
            mom = logc[idx] - logc[idx - L]
            fwd = logc[idx + H] - logc[idx]
            same_q = (q[idx - L] == q[idx]) & (q[idx + H] == q[idx])
            mom, fwd, ii = mom[same_q], fwd[same_q], idx[same_q]
            if len(mom) < 100 or mom.std() == 0 or fwd.std() == 0:
                continue
            ic = pearsonr(mom, fwd)[0]
            ics = spearmanr(mom, fwd)[0]
            sign = np.mean(np.sign(mom) == np.sign(fwd))
            nwt = newey_west_tstat(mom, fwd, lags=H)
            oos = ts[ii] >= oos_split_ts
            oos_ic = (pearsonr(mom[oos], fwd[oos])[0]
                      if oos.sum() > 50 and mom[oos].std() > 0 else float("nan"))
            fwd_ticks = np.median(np.abs(np.exp(fwd) - 1.0) * close[ii]) / TICK
            print(f"{L:>4} {H:>4} {len(mom):>7} {ic:>9.4f} {ics:>9.4f} {100*sign:>6.2f}% "
                  f"{nwt:>7.2f} {'n/a' if np.isnan(oos_ic) else f'{oos_ic:>8.4f}'} {fwd_ticks:>11.2f}")


def to_daily(ts, close, q):
    day = ts // (86400 * 10**9)
    ud, last = np.unique(day, return_index=False), None
    # last bar per day (ts sorted ascending overall within concat? sort to be safe)
    order = np.argsort(ts, kind="stable")
    ts, close, q, day = ts[order], close[order], q[order], day[order]
    udays, first = np.unique(day, return_index=True)
    last_idx = np.append(first[1:] - 1, len(day) - 1)
    return ts[last_idx], close[last_idx], q[last_idx]


def self_test():
    # varying-momentum series: per-bar return increases linearly, so trailing
    # 1-bar return (= r[t]) and forward 1-bar return (= r[t+1]) are perfectly
    # rank-correlated → IC≈1. (A constant step would be near-constant → undefined.)
    r = 0.0001 * np.arange(1, 52)
    close = np.exp(np.cumsum(np.r_[0.0, r]))
    q = np.zeros(len(close)); ts = np.arange(len(close)) * 60 * 10**9
    logc = np.log(close)
    L = H = 1
    idx = np.arange(L, len(close) - H)
    mom = logc[idx] - logc[idx - L]; fwd = logc[idx + H] - logc[idx]
    from scipy.stats import pearsonr
    ic = pearsonr(mom, fwd)[0]
    assert ic > 0.99, ic
    print(f"lowfreq self-test PASS: trending series momentum IC={ic:.3f}")
    return 0


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("dir", nargs="?", default="test_data/futures-baseline/ES.FUT")
    ap.add_argument("--self-test", action="store_true")
    args = ap.parse_args()
    if args.self_test:
        return self_test()
    ts, close, q = load_ohlcv(args.dir)
    order = np.argsort(ts, kind="stable")
    ts, close, q = ts[order], close[order], q[order]
    span_days = (ts[-1] - ts[0]) / 1e9 / 86400
    # OOS split at 2025-01-01 UTC
    split = 1735689600 * 10**9
    print(f"\nloaded {len(close)} front-month 1m bars, span {span_days:.0f} days, "
          f"{int((ts<split).sum())} train(<2025) / {int((ts>=split).sum())} test(>=2025)")

    # INTRADAY (minutes), contiguity-guarded
    ic_block("INTRADAY momentum (1m bars; L,H in minutes)", ts, close, q,
             lookbacks=[5, 15, 30, 60, 120], horizons=[5, 15, 30, 60], bar_unit_s=60,
             oos_split_ts=split)

    # DAILY time-series momentum
    dts, dclose, dq = to_daily(ts, close, q)
    print(f"\n[daily] {len(dclose)} daily bars")
    ic_block("DAILY time-series momentum (L,H in days)", dts, dclose, dq,
             lookbacks=[1, 2, 5, 10, 20, 60], horizons=[1, 5, 10, 20], bar_unit_s=86400,
             oos_split_ts=split)

    # Daily contrarian strategy backtest (non-overlapping daily P&L = honest test).
    # signal_t = -sign(trailing L-day return); position held, earns next-day return.
    # Daily P&L is non-overlapping by construction, sidestepping the overlap inflation
    # that exaggerates the multi-day IC above.
    print("\n=== DAILY CONTRARIAN BACKTEST (pos = -sign(trailing L-day ret), daily rebalanced) ===")
    dlogc = np.log(dclose)
    dr = np.diff(dlogc, prepend=dlogc[0])               # daily log returns
    same_day_q = np.r_[True, dq[1:] == dq[:-1]]         # not a roll boundary
    split = 1735689600 * 10**9
    print(f"{'L(d)':>5} {'Sharpe_all':>11} {'Sharpe_IS24':>12} {'Sharpe_OOS25':>13} "
          f"{'ann_ret%':>9} {'turnover/yr':>11} {'hitrate':>8}")
    for L in [2, 5, 10, 20]:
        sig = -np.sign(dlogc[L:] - dlogc[:-L])          # contrarian signal at day t (uses past L days)
        # position at day t earns return at day t+1
        pos = sig[:-1]
        ret = dr[L + 1:]                                 # next-day returns aligned to pos
        qok = same_day_q[L + 1:]                         # exclude roll-boundary next-days
        tsd = dts[L + 1:]
        pos, ret, tsd = pos[qok], ret[qok], tsd[qok]
        pnl = pos * ret
        if len(pnl) < 30 or pnl.std() == 0:
            continue
        def sharpe(p):
            return p.mean() / p.std() * np.sqrt(252) if len(p) > 5 and p.std() > 0 else float("nan")
        ism = tsd < split; oosm = tsd >= split
        turn = np.mean(np.abs(np.diff(pos)) > 0) * 252
        hit = np.mean(np.sign(pnl) > 0)
        print(f"{L:>5} {sharpe(pnl):>11.2f} {sharpe(pnl[ism]):>12.2f} {sharpe(pnl[oosm]):>13.2f} "
              f"{pnl.mean()*252*100:>8.2f}% {turn:>11.0f} {100*hit:>6.1f}%")
    print("Cost: 20d-horizon turnover is low; 1-tick round-trip (~0.005% of price) is negligible.")
    print("HONEST READ: trust OOS Sharpe (non-overlapping). >0.5 OOS on thin 2y bull data is")
    print("suggestive, NOT proven — needs full history + purged CV + regime stratification.")

    print("\nNOTE: at these horizons |fwd|ticks >> 1-tick ($12.50) crossing cost, so cost")
    print("is not the binding constraint — the question is purely whether OOS_IC > 0 and stable.")
    return 0


if __name__ == "__main__":
    sys.exit(main())