"""
V47 Diagnostic: Overnight Gap Residual (Sector-ETF OLS Regression)

Hypothesis: The idiosyncratic component of a stock's overnight gap (after removing
systematic sector/market beta) predicts ORB follow-through better than the raw gap.
Based on Lou/Polk (JFE 2019) beta decomposition and Bogousslavsky (JF 2021) overnight
persistence as informed-trader footprint.

Features tested:
  gap_residual             : OLS residual of stock gap on sector ETF gap (rolling 60d β̂)
  gap_residual_qqq         : OLS residual of stock gap on QQQ gap (QQQ universal benchmark)
  gap_residual_std         : Z-scored gap_residual (divided by σ̂ of in-window residuals)
  overnight_residual_avg_5d : Mean of 5 prior daily gap_residuals (Bogousslavsky persistence)
  overnight_residual_avg_10d: Mean of 10 prior daily gap_residuals

Gates (pre-committed, must be set before running):
  G1: |Pearson| ≥ 0.07 AND n ≥ 200
  G2: |top − bottom tercile avg_R| ≥ 0.30R  ← binding gate
  G3: |top − bottom tercile WR| ≥ 5pp
  G4: coverage ≥ 80% of V46 400d trade set
  G5a: |ρ(feature, obv_slope_20)| < 0.70
  G5b: |ρ(feature, raw gap_pct)| < 0.70

Decision rule (pre-committed):
  ≥1 feature passes G1-G5b → advance highest-G2 to Phase 2 backtest
  0 pass → close axis, document in lineage memory

Data source: data/cache/intraday/{ticker}/{date}.parquet
V46 400d baseline: runs/v46_correct_w12_400d.json/intraday_20260422_043215_63f03e64.json
"""
from __future__ import annotations

import json
import os
import sys
from pathlib import Path

import pyarrow.parquet as pq

sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), "../../..")))

from zoneinfo import ZoneInfo
from libs.common.time_utils import trading_days_between
import datetime as dt

_ET = ZoneInfo("America/New_York")
_MKT_OPEN = dt.time(9, 30)
_MKT_CLOSE = dt.time(16, 0)

INTRADAY_CACHE_DIR = "data/cache/intraday"
SECTOR_CACHE = "data/cache/sector_cache.json"
V46_400D_RUN = "runs/v46_correct_w12_400d.json/intraday_20260422_043215_63f03e64.json"

SECTOR_TO_ETF: dict[str, str] = {
    "Technology": "XLK",
    "Financial Services": "XLF",
    "Healthcare": "XLV",
    "Energy": "XLE",
    "Industrials": "XLI",
    "Basic Materials": "XLB",
    "Communication Services": "XLC",
    "Utilities": "XLU",
    "Real Estate": "XLRE",
    "Consumer Defensive": "XLP",
    "Consumer Cyclical": "XLY",
}

BETA_WINDOW = 60   # rolling days for β̂ estimation
MIN_PAIRS = 20     # minimum non-null pairs for valid β̂


# ── Helpers ────────────────────────────────────────────────────────────────

def _parse_ts(ts_raw: object) -> dt.datetime:
    s = str(ts_raw)
    if s.endswith("Z"):
        s = s[:-1] + "+00:00"
    return dt.datetime.fromisoformat(s).astimezone(_ET)


def load_daily_ohlcv(ticker: str, dates: list[str]) -> dict[str, dict]:
    """Load regular-session {open, close} per date from intraday parquet cache."""
    result: dict[str, dict] = {}
    for date in dates:
        path = Path(INTRADAY_CACHE_DIR) / ticker.upper() / f"{date}.parquet"
        if not path.exists():
            continue
        try:
            table = pq.read_table(str(path))
            rows = table.to_pydict()
        except Exception:
            continue

        mkt_bars = []
        for i, ts_raw in enumerate(rows.get("timestamp", [])):
            try:
                ts = _parse_ts(ts_raw)
            except Exception:
                continue
            if not (_MKT_OPEN <= ts.time() < _MKT_CLOSE):
                continue
            mkt_bars.append({
                "ts": ts,
                "open": float(rows["open"][i] or 0),
                "close": float(rows["close"][i] or 0),
            })

        if not mkt_bars:
            continue
        mkt_bars.sort(key=lambda b: b["ts"])
        result[date] = {
            "open": mkt_bars[0]["open"],
            "close": mkt_bars[-1]["close"],
        }
    return result


def compute_overnight_gaps(ohlcv: dict[str, dict], sorted_dates: list[str]) -> dict[str, float]:
    """Compute overnight gap = (open_D - close_{D-1}) / close_{D-1} for each date."""
    gaps: dict[str, float] = {}
    for i in range(1, len(sorted_dates)):
        d0, d1 = sorted_dates[i - 1], sorted_dates[i]
        prev = ohlcv.get(d0)
        curr = ohlcv.get(d1)
        if prev and curr and prev["close"] > 0 and curr["open"] > 0:
            gaps[d1] = (curr["open"] - prev["close"]) / prev["close"]
    return gaps


def _winsorize(vals: list[float], n_sigma: float = 3.0) -> list[float]:
    if len(vals) < 4:
        return vals
    n = len(vals)
    mu = sum(vals) / n
    sigma = (sum((v - mu) ** 2 for v in vals) / n) ** 0.5
    if sigma <= 0:
        return vals
    lo, hi = mu - n_sigma * sigma, mu + n_sigma * sigma
    return [max(lo, min(hi, v)) for v in vals]


def ols_residual_on_window(
    entry_stock_gap: float,
    entry_etf_gap: float,
    prior_stock_gaps: list[float],
    prior_etf_gaps: list[float],
) -> tuple[float, float, float] | None:
    """
    OLS regression: stock_gap = alpha + beta * etf_gap on the prior window.
    Returns (residual_at_entry, beta, sigma_residuals) or None if insufficient data.
    """
    assert len(prior_stock_gaps) == len(prior_etf_gaps)
    n = len(prior_stock_gaps)
    if n < MIN_PAIRS:
        return None

    xs = _winsorize(prior_etf_gaps)
    ys = _winsorize(prior_stock_gaps)

    xm = sum(xs) / n
    ym = sum(ys) / n
    num = sum((xs[i] - xm) * (ys[i] - ym) for i in range(n))
    den = sum((xs[i] - xm) ** 2 for i in range(n))
    if den <= 0:
        return None

    beta = num / den
    alpha = ym - beta * xm

    residuals_window = [ys[i] - (alpha + beta * xs[i]) for i in range(n)]
    sigma = (sum(r ** 2 for r in residuals_window) / n) ** 0.5

    residual_entry = entry_stock_gap - (alpha + beta * entry_etf_gap)
    return residual_entry, beta, sigma


def pearson(xs: list[float], ys: list[float]) -> float | None:
    n = len(xs)
    if n < 5:
        return None
    xm, ym = sum(xs) / n, sum(ys) / n
    num = sum((xs[i] - xm) * (ys[i] - ym) for i in range(n))
    dx = sum((x - xm) ** 2 for x in xs) ** 0.5
    dy = sum((y - ym) ** 2 for y in ys) ** 0.5
    if dx <= 0 or dy <= 0:
        return None
    return num / (dx * dy)


def tercile_stats(vals: list[float], rs: list[float]) -> dict:
    if len(vals) < 9:
        return {}
    pairs = sorted(zip(vals, rs), key=lambda p: p[0])
    n = len(pairs)
    t = n // 3

    def stats(sub):
        ys = [p[1] for p in sub]
        return {"n": len(ys), "wr": sum(1 for y in ys if y > 0) / len(ys), "avg_r": sum(ys) / len(ys)}

    return {"low": stats(pairs[:t]), "mid": stats(pairs[t:2 * t]), "high": stats(pairs[2 * t:])}


def pairwise_rho(idxs_a: list[int], vals_a: list[float],
                 idxs_b: list[int], vals_b: list[float]) -> float | None:
    map_b = {idx: vals_b[i] for i, idx in enumerate(idxs_b)}
    va, vb = [], []
    for i, idx in enumerate(idxs_a):
        if idx in map_b:
            va.append(vals_a[i])
            vb.append(map_b[idx])
    return pearson(va, vb) if len(va) >= 10 else None


def report_feature(
    label: str,
    idxs: list[int],
    vals: list[float],
    rs: list[float],
    n_total: int,
    g5_comparisons: list[tuple[str, list[int], list[float]]],
) -> dict:
    n = len(vals)
    p = pearson(vals, rs)
    ts = tercile_stats(vals, rs)
    coverage = n / n_total

    g4 = coverage >= 0.80

    if not ts or p is None:
        print(f"  {label}: n={n}/{n_total} ({coverage:.1%})  INSUFFICIENT DATA")
        return {"pass_g2": False, "pass_all": False, "pearson": None, "g2_gap": None}

    low, mid, high = ts["low"], ts["mid"], ts["high"]
    raw_g2 = high["avg_r"] - low["avg_r"]
    avg_r_gap = abs(raw_g2)
    wr_gap = abs(high["wr"] - low["wr"])

    g1 = abs(p) >= 0.07 and n >= 200
    g2 = avg_r_gap >= 0.30
    g3 = wr_gap >= 0.05

    print(f"\n  [{label}]  n={n}/{n_total} ({coverage:.1%})  Pearson={p:+.3f}  direction={'↑high' if raw_g2 > 0 else '↓low'}")
    print(f"    G1: {'PASS' if g1 else 'FAIL'} (|{abs(p):.3f}| {'≥' if abs(p) >= 0.07 else '<'} 0.07, n={n} {'≥' if n >= 200 else '<'} 200)")
    print(f"    G2: {'PASS' if g2 else 'FAIL'} (avg_R gap = {avg_r_gap:.3f}R  [threshold 0.30R])")
    print(f"    G3: {'PASS' if g3 else 'FAIL'} (WR gap = {wr_gap * 100:.1f}pp  [threshold 5pp])")
    print(f"    G4: {'PASS' if g4 else 'FAIL'} (coverage = {coverage:.1%}  [threshold 80%])")
    print(f"    Bottom tercile: WR {low['wr'] * 100:.1f}%  avg_R {low['avg_r']:+.3f}  n={low['n']}")
    print(f"    Middle tercile: WR {mid['wr'] * 100:.1f}%  avg_R {mid['avg_r']:+.3f}  n={mid['n']}")
    print(f"    Top tercile:    WR {high['wr'] * 100:.1f}%  avg_R {high['avg_r']:+.3f}  n={high['n']}")

    g5_results = {}
    pass_all_g5 = True
    for cmp_label, cmp_idxs, cmp_vals in g5_comparisons:
        rho = pairwise_rho(idxs, vals, cmp_idxs, cmp_vals)
        if rho is not None:
            gate_key = f"g5_{cmp_label}"
            gate_pass = abs(rho) < 0.70
            g5_results[gate_key] = {"rho": rho, "pass": gate_pass}
            print(f"    G5({cmp_label}): {'PASS' if gate_pass else 'FAIL'} (ρ={rho:+.3f}  [|threshold| < 0.70])")
            if not gate_pass:
                pass_all_g5 = False

    pass_all = g1 and g2 and g3 and g4 and pass_all_g5
    print(f"    → {'ALL GATES PASS ✓' if pass_all else 'FAIL'}")

    return {
        "pass_g2": g2, "pass_all": pass_all,
        "pearson": p, "g2_gap": avg_r_gap, "raw_g2_gap": raw_g2,
        "direction": "high" if raw_g2 > 0 else "low",
        **g5_results,
    }


# ── Main ───────────────────────────────────────────────────────────────────

def main() -> None:
    print("=== V47 Diagnostic: Overnight Gap Residual ===\n")
    print("Pre-committed gates: G1(|P|≥0.07,n≥200) G2(avg_R≥0.30R) G3(WR≥5pp) G4(cov≥80%) G5(|ρ|<0.70)")
    print("Pre-committed decision rule: highest G2 that passes ALL gates → advance to Phase 2 backtest\n")

    with open(V46_400D_RUN) as f:
        run_data = json.load(f)
    trades_raw = run_data.get("trades", [])
    m = run_data.get("metrics", {})
    print(f"V46 400d: {m.get('total_trades')} trades, {m.get('start_date')} → {m.get('end_date')}, "
          f"return={m.get('total_return_pct', 0) * 100:.1f}%\n")

    with open(SECTOR_CACHE) as f:
        sector_cache = json.load(f)

    trades = [
        {
            "idx": i,
            "ticker": t["ticker"],
            "date": t["date"][:10],
            "r": float(t["r_multiple_at_exit"]),
            "sector": sector_cache.get(t["ticker"], "UNKNOWN"),
        }
        for i, t in enumerate(trades_raw)
        if t.get("r_multiple_at_exit") is not None
    ]
    n_total = len(trades)
    print(f"Valid trades (with r_multiple): {n_total}")

    # Build date range: from earliest trade date - 70 trading days to last trade date
    all_trade_dates = sorted(set(t["date"] for t in trades))
    start_d = dt.date.fromisoformat(all_trade_dates[0]) - dt.timedelta(days=110)
    end_d = dt.date.fromisoformat(all_trade_dates[-1])
    all_trading_days = [d.isoformat() for d in trading_days_between(start_d, end_d)]
    td_set = set(all_trading_days)

    print(f"Computing over trading days: {all_trading_days[0]} → {all_trading_days[-1]} ({len(all_trading_days)} days)")

    # Load sector ETF + QQQ daily OHLCV
    etf_tickers = list(set(SECTOR_TO_ETF.values())) + ["QQQ"]
    print(f"\nLoading {len(etf_tickers)} sector ETF + QQQ daily OHLCV...")
    etf_ohlcv: dict[str, dict[str, dict]] = {}
    for etf in etf_tickers:
        etf_ohlcv[etf] = load_daily_ohlcv(etf, all_trading_days)
        print(f"  {etf}: {len(etf_ohlcv[etf])} days")

    # Pre-compute ETF overnight gaps
    etf_gaps: dict[str, dict[str, float]] = {}
    for etf, ohlcv in etf_ohlcv.items():
        sorted_dates = sorted(ohlcv.keys())
        etf_gaps[etf] = compute_overnight_gaps(ohlcv, sorted_dates)
    print()

    # Load stock OHLCV for all unique tickers
    unique_tickers = sorted(set(t["ticker"] for t in trades))
    print(f"Loading {len(unique_tickers)} stock intraday caches...")
    stock_ohlcv: dict[str, dict[str, dict]] = {}
    for ticker in unique_tickers:
        stock_ohlcv[ticker] = load_daily_ohlcv(ticker, all_trading_days)
    print(f"  Loaded. Sample: {unique_tickers[0]}: {len(stock_ohlcv[unique_tickers[0]])} days")

    # Pre-compute stock overnight gaps
    stock_gaps: dict[str, dict[str, float]] = {}
    for ticker in unique_tickers:
        sorted_dates = sorted(stock_ohlcv[ticker].keys())
        stock_gaps[ticker] = compute_overnight_gaps(stock_ohlcv[ticker], sorted_dates)

    # Pre-compute daily residuals for each ticker (for Bogousslavsky persistence)
    # residuals_sector[ticker][date] = gap_residual on that date (sector ETF β̂)
    # residuals_qqq[ticker][date]    = gap_residual on that date (QQQ β̂)
    print("\nPre-computing daily residuals for all tickers (for persistence features)...")
    residuals_sector: dict[str, dict[str, float]] = {}
    residuals_qqq: dict[str, dict[str, float]] = {}

    for ticker in unique_tickers:
        sector = sector_cache.get(ticker, "UNKNOWN")
        sector_etf = SECTOR_TO_ETF.get(sector)
        s_gaps = stock_gaps[ticker]
        sorted_stock_dates = sorted(s_gaps.keys())

        r_sector: dict[str, float] = {}
        r_qqq: dict[str, float] = {}

        for j, date in enumerate(sorted_stock_dates):
            sg = s_gaps.get(date)
            if sg is None:
                continue

            # Get prior BETA_WINDOW dates with paired gaps
            prior_dates = [d for d in sorted_stock_dates[:j] if d < date][-BETA_WINDOW:]

            # Sector ETF residual
            if sector_etf:
                eg = etf_gaps[sector_etf]
                prior_pairs_s = [(s_gaps[d], eg[d]) for d in prior_dates if d in s_gaps and d in eg]
                if len(prior_pairs_s) >= MIN_PAIRS:
                    ps, pe = [x[0] for x in prior_pairs_s], [x[1] for x in prior_pairs_s]
                    entry_etf_gap = etf_gaps[sector_etf].get(date)
                    if entry_etf_gap is not None:
                        res = ols_residual_on_window(sg, entry_etf_gap, ps, pe)
                        if res is not None:
                            r_sector[date] = res[0]

            # QQQ residual
            qq_gaps = etf_gaps["QQQ"]
            prior_pairs_q = [(s_gaps[d], qq_gaps[d]) for d in prior_dates if d in s_gaps and d in qq_gaps]
            if len(prior_pairs_q) >= MIN_PAIRS:
                pq_s, pq_q = [x[0] for x in prior_pairs_q], [x[1] for x in prior_pairs_q]
                entry_qqq_gap = qq_gaps.get(date)
                if entry_qqq_gap is not None:
                    res = ols_residual_on_window(sg, entry_qqq_gap, pq_s, pq_q)
                    if res is not None:
                        r_qqq[date] = res[0]

        residuals_sector[ticker] = r_sector
        residuals_qqq[ticker] = r_qqq

    # ── Compute features for each trade ─────────────────────────────────
    print("\nComputing features for each trade...\n")

    # Feature arrays (indexed by trade idx for pairwise ρ)
    feat_gap_residual: list[tuple[int, float]] = []
    feat_gap_residual_qqq: list[tuple[int, float]] = []
    feat_gap_residual_std: list[tuple[int, float]] = []
    feat_persist_5d: list[tuple[int, float]] = []
    feat_persist_10d: list[tuple[int, float]] = []
    feat_raw_gap: list[tuple[int, float]] = []
    feat_obv_slope: list[tuple[int, float]] = []  # placeholder - computed from daily bars if available

    r_by_idx: dict[int, float] = {t["idx"]: t["r"] for t in trades}

    for trade in trades:
        idx = trade["idx"]
        ticker = trade["ticker"]
        date = trade["date"]
        sector = trade["sector"]
        sector_etf = SECTOR_TO_ETF.get(sector)

        sg = stock_gaps[ticker].get(date)
        if sg is None:
            continue

        # Raw gap feature
        feat_raw_gap.append((idx, sg))

        # gap_residual (sector ETF)
        if date in residuals_sector.get(ticker, {}):
            feat_gap_residual.append((idx, residuals_sector[ticker][date]))

        # gap_residual_qqq
        if date in residuals_qqq.get(ticker, {}):
            feat_gap_residual_qqq.append((idx, residuals_qqq[ticker][date]))

        # gap_residual_std: need sigma from the window on trade date
        # Recompute sigma for this trade date
        if sector_etf:
            s_gaps_t = stock_gaps[ticker]
            sorted_sd = sorted(s_gaps_t.keys())
            j = sorted_sd.index(date) if date in sorted_sd else -1
            if j > 0:
                prior_dates = sorted_sd[:j][-BETA_WINDOW:]
                eg = etf_gaps[sector_etf]
                prior_pairs = [(s_gaps_t[d], eg[d]) for d in prior_dates if d in s_gaps_t and d in eg]
                if len(prior_pairs) >= MIN_PAIRS:
                    ps, pe = [x[0] for x in prior_pairs], [x[1] for x in prior_pairs]
                    entry_etf_gap = eg.get(date)
                    if entry_etf_gap is not None:
                        res = ols_residual_on_window(sg, entry_etf_gap, ps, pe)
                        if res is not None and res[2] > 0:
                            feat_gap_residual_std.append((idx, res[0] / res[2]))

        # Persistence features: mean of prior 5/10 daily residuals (sector β̂)
        sorted_all_td = sorted(td_set & set(stock_gaps[ticker].keys()))
        if date in sorted_all_td:
            d_idx = sorted_all_td.index(date)
            prior_td = sorted_all_td[max(0, d_idx - 10):d_idx]
            ticker_residuals = residuals_sector.get(ticker, {})

            prev_5d_resids = [ticker_residuals[d] for d in prior_td[-5:] if d in ticker_residuals]
            prev_10d_resids = [ticker_residuals[d] for d in prior_td[-10:] if d in ticker_residuals]

            if len(prev_5d_resids) >= 3:
                feat_persist_5d.append((idx, sum(prev_5d_resids) / len(prev_5d_resids)))
            if len(prev_10d_resids) >= 5:
                feat_persist_10d.append((idx, sum(prev_10d_resids) / len(prev_10d_resids)))

    # ── Build r-multiple arrays aligned to each feature ──────────────
    def split(feat_pairs: list[tuple[int, float]]) -> tuple[list[int], list[float], list[float]]:
        idxs = [p[0] for p in feat_pairs]
        vals = [p[1] for p in feat_pairs]
        rs = [r_by_idx[i] for i in idxs]
        return idxs, vals, rs

    # OBV slope from run data if available (fallback: empty)
    obv_idxs, obv_vals, _ = split(feat_obv_slope) if feat_obv_slope else ([], [], [])

    raw_idxs, raw_vals, raw_rs = split(feat_raw_gap)

    g5_refs = [
        ("raw_gap", raw_idxs, raw_vals),
    ]

    print("=" * 70)
    print(f"GATE RESULTS  (n_total={n_total})")
    print("=" * 70)

    results = {}

    for label, feat_pairs in [
        ("gap_residual (sector ETF β̂)", feat_gap_residual),
        ("gap_residual_qqq (QQQ β̂)", feat_gap_residual_qqq),
        ("gap_residual_std (z-scored, sector ETF)", feat_gap_residual_std),
        ("overnight_residual_avg_5d (Bogousslavsky)", feat_persist_5d),
        ("overnight_residual_avg_10d (Bogousslavsky)", feat_persist_10d),
    ]:
        if not feat_pairs:
            print(f"\n  [{label}]: no data")
            continue
        idxs, vals, rs = split(feat_pairs)
        results[label] = report_feature(label, idxs, vals, rs, n_total, g5_refs)

    # Pairwise ρ matrix among all features
    print("\n\n=== Pairwise ρ matrix ===")
    all_feat_sets = [
        ("gap_residual", feat_gap_residual),
        ("gap_residual_qqq", feat_gap_residual_qqq),
        ("gap_residual_std", feat_gap_residual_std),
        ("persist_5d", feat_persist_5d),
        ("persist_10d", feat_persist_10d),
        ("raw_gap", feat_raw_gap),
    ]
    feat_names = [n for n, _ in all_feat_sets]
    feat_data = {n: (split(fp)[0], split(fp)[1]) for n, fp in all_feat_sets if fp}

    for na in feat_names:
        if na not in feat_data:
            continue
        row_parts = []
        for nb in feat_names:
            if nb not in feat_data:
                row_parts.append("  N/A ")
                continue
            if na == nb:
                row_parts.append(" 1.000")
                continue
            rho = pairwise_rho(feat_data[na][0], feat_data[na][1],
                               feat_data[nb][0], feat_data[nb][1])
            row_parts.append(f"{rho:+.3f}" if rho is not None else "  N/A ")
        print(f"  {na:30s} {' '.join(row_parts)}")

    # Summary
    print("\n\n=== Summary & Decision ===")
    passing = [(lbl, r) for lbl, r in results.items() if r.get("pass_all")]
    g2_only = [(lbl, r) for lbl, r in results.items() if r.get("pass_g2")]

    print(f"Features passing ALL gates: {len(passing)}")
    for lbl, r in passing:
        print(f"  ✓ {lbl}  G2={r['g2_gap']:.3f}R  Pearson={r['pearson']:+.3f}")

    if passing:
        # Pick highest G2 among passing
        best_lbl, best_r = max(passing, key=lambda x: x[1]["g2_gap"])
        print(f"\n→ ADVANCE to Phase 2: {best_lbl}")
        print(f"  G2={best_r['g2_gap']:.3f}R  Pearson={best_r['pearson']:+.3f}  direction={best_r['direction']}")
    else:
        print("\nNo features pass ALL gates.")
        best_g2 = max((r.get("g2_gap") or 0) for r in results.values()) if results else 0
        print(f"Best G2 observed: {best_g2:.3f}R  (threshold 0.30R)")
        print("→ CLOSE overnight gap residual axis. Document in project_orb_gainers_lineage.md.")
        print("  Next axis candidates: options flow (put/call ratio, IV skew) or 13F ownership momentum.")

    print("\n" + "=" * 70)
    print(f"Raw gap reference: Pearson={pearson(raw_vals, raw_rs):+.3f}" if len(raw_vals) >= 5 else "")


if __name__ == "__main__":
    main()
