fithia2/apps/intraday_bt/scripts/diag_orb_quality_features.py

"""
V23 Quality Feature Diagnostic: Hurst / OU-θ / OBV-Slope

Hypothesis: Daily-prior return-series statistics (Hurst exponent, OU mean-reversion
speed, OBV accumulation slope) carry incremental per-trade edge when overlaid on V23's
candidate set — similar to how entropy_20d works.

Method:
  1. Run V23 simulation over 200d window (baseline, full regime filter active)
  2. For each completed trade, compute three new features on prev_bars at trade date:
       - hurst_60    : Hurst exponent (R/S, 60-day lookback) — H>0.5 = trending
       - ou_theta_60 : OU mean-reversion speed (AR(1) β→-ln(β), 60d) — high = fast reversion
       - obv_slope_20: OBV accumulation slope (normalized, 20d) — positive = accumulation
  3. Report per-feature:
       - Pearson correlation with r_multiple_at_exit
       - Tercile-split WR, avg_pnl, sample counts

Gates (for promotion to Phase 2):
  G1: |Pearson| ≥ 0.07 with r_multiple_at_exit
  G2: top-tercile avg_r − bottom-tercile avg_r ≥ 0.3R
  G3: top-tercile WR ≥ bottom-tercile WR + 5pp
"""
from __future__ import annotations

import datetime as dt
import json
import math
import os
import sys
from concurrent.futures import ThreadPoolExecutor
from pathlib import Path

import pyarrow.parquet as pq
import yaml
from zoneinfo import ZoneInfo

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

from libs.common.time_utils import trading_days_between
from libs.intraday.domain import ORBStrategyParams
from libs.intraday.features import enrich_daily_bars
from libs.intraday.orb_simulator import ORBSimulationState, run_orb_simulation_with_state
from libs.intraday.screener import orb_pre_screen_candidates

# ── Config ─────────────────────────────────────────────────────────────────
V23_CONFIG = "configs/intraday/strategies/orb_gainers_v23.yaml"
UNIVERSE_FILE = "configs/symbols_midlarge_snapshot_exact.yaml"
INTRADAY_CACHE_DIR = "data/cache/intraday"
LOOKBACK_DAYS = 200
FEATURE_LOOKBACK_TRADING = 70  # need 60 trading days; add buffer → 70

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


# ── Feature Computations (ported from libs/features/market_features.py) ────

def compute_hurst_approx(prev_bars: list[dict], lookback: int = 60) -> float | None:
    """Hurst exponent via R/S rescaled-range. H>0.5 trending, H<0.5 mean-reverting."""
    if len(prev_bars) < lookback + 2:
        return None
    sorted_bars = sorted(prev_bars, key=lambda b: b["date"])
    recent = sorted_bars[-(lookback + 1):]
    returns: list[float] = []
    for i in range(1, len(recent)):
        pc = recent[i - 1].get("close")
        cc = recent[i].get("close")
        if not pc or not cc or pc <= 0:
            continue
        returns.append((cc - pc) / pc)
    if len(returns) < 30:
        return None

    def rs_stat(series: list[float]) -> float:
        n = len(series)
        mean = sum(series) / n
        devs = [x - mean for x in series]
        cumdev: list[float] = []
        s = 0.0
        for d in devs:
            s += d
            cumdev.append(s)
        r = max(cumdev) - min(cumdev)
        std = (sum(d ** 2 for d in devs) / n) ** 0.5
        return r / std if std > 0 else 0.0

    window_sizes = [w for w in [8, 12, 16, 24, 32] if w <= len(returns) // 2]
    if len(window_sizes) < 2:
        return None
    log_n: list[float] = []
    log_rs: list[float] = []
    for w in window_sizes:
        rs_vals = [
            rs_stat(returns[start:start + w])
            for start in range(0, len(returns) - w + 1, w)
            if len(returns[start:start + w]) == w
        ]
        if rs_vals:
            avg_rs = sum(rs_vals) / len(rs_vals)
            if avg_rs > 0:
                log_n.append(math.log(w))
                log_rs.append(math.log(avg_rs))
    if len(log_n) < 2:
        return None
    n = len(log_n)
    xm = sum(log_n) / n
    ym = sum(log_rs) / n
    num = sum((log_n[i] - xm) * (log_rs[i] - ym) for i in range(n))
    den = sum((log_n[i] - xm) ** 2 for i in range(n))
    if den <= 0:
        return 0.5
    return num / den


def compute_ou_theta_approx(prev_bars: list[dict], lookback: int = 60) -> float | None:
    """OU mean-reversion speed θ via AR(1). θ = -ln(β). High = fast reversion."""
    if len(prev_bars) < lookback + 2:
        return None
    sorted_bars = sorted(prev_bars, key=lambda b: b["date"])
    recent = sorted_bars[-(lookback + 1):]
    log_prices: list[float] = []
    for bar in recent:
        c = bar.get("close")
        if not c or c <= 0:
            return None
        log_prices.append(math.log(c))
    if len(log_prices) < lookback:
        return None
    y = log_prices[1:]
    x = log_prices[:-1]
    n = len(y)
    xm = sum(x) / n
    ym = sum(y) / n
    num = sum((x[i] - xm) * (y[i] - ym) for i in range(n))
    den = sum((x[i] - xm) ** 2 for i in range(n))
    if den <= 0:
        return None
    beta = num / den
    if beta <= 0 or beta >= 1.0:
        return 0.0
    return -math.log(beta)


def compute_obv_slope_approx(prev_bars: list[dict], lookback: int = 20) -> float | None:
    """OBV accumulation slope, normalized by avg_volume. Positive = accumulation."""
    if len(prev_bars) < lookback + 2:
        return None
    sorted_bars = sorted(prev_bars, key=lambda b: b["date"])
    recent = sorted_bars[-(lookback + 1):]
    obv_series = [0.0]
    total_vol = 0.0
    for i in range(1, len(recent)):
        pc = recent[i - 1].get("close")
        cc = recent[i].get("close")
        vol = float(recent[i].get("volume") or 0)
        total_vol += vol
        if not pc or pc <= 0:
            continue
        if cc > pc:
            obv_series.append(obv_series[-1] + vol)
        elif cc < pc:
            obv_series.append(obv_series[-1] - vol)
        else:
            obv_series.append(obv_series[-1])
    avg_vol = total_vol / lookback if lookback > 0 else 1.0
    if avg_vol <= 0:
        return None
    n = len(obv_series)
    xm = (n - 1) / 2.0
    ym = sum(obv_series) / n
    num = sum((i - xm) * (obv_series[i] - ym) for i in range(n))
    den = sum((i - xm) ** 2 for i in range(n))
    if den <= 0:
        return 0.0
    return (num / den) / avg_vol


# ── Daily Bar Builder ────────────────────────────────────────────────────────

def _build_daily_bar_from_intraday(path: Path, date: str) -> dict | None:
    try:
        table = pq.read_table(str(path))
        rows = table.to_pydict()
    except Exception:
        return None
    opens, highs, lows, closes, vols = [], [], [], [], []
    for i, ts_raw in enumerate(rows.get("timestamp", [])):
        try:
            if ts_raw.endswith("Z"):
                ts_raw = ts_raw[:-1] + "+00:00"
            ts = dt.datetime.fromisoformat(ts_raw).astimezone(_ET)
        except Exception:
            continue
        if _MKT_OPEN <= ts.time() < _MKT_CLOSE:
            opens.append(float(rows["open"][i] or 0))
            highs.append(float(rows["high"][i] or 0))
            lows.append(float(rows["low"][i] or 0))
            closes.append(float(rows["close"][i] or 0))
            vols.append(float(rows["volume"][i] or 0))
    if not opens:
        return None
    return {
        "date": date, "open": opens[0], "high": max(highs),
        "low": min(lows), "close": closes[-1], "volume": sum(vols),
    }


def build_daily_bars(tickers: list[str], dates: list[str], workers: int = 8) -> dict[str, list[dict]]:
    root = Path(INTRADAY_CACHE_DIR)

    def _load(ticker: str) -> tuple[str, list[dict]]:
        d_path = root / ticker
        if not d_path.is_dir():
            return ticker, []
        bars: list[dict] = []
        for date in dates:
            p = d_path / f"{date}.parquet"
            if not p.exists():
                continue
            bar = _build_daily_bar_from_intraday(p, date)
            if bar and bar["close"] > 0:
                bars.append(bar)
        return ticker, bars

    result: dict[str, list[dict]] = {}
    with ThreadPoolExecutor(max_workers=workers) as ex:
        for ticker, bars in ex.map(_load, tickers):
            if bars:
                result[ticker] = bars
    return result


def load_intraday_bulk(candidates: dict[str, list[str]]) -> dict[str, dict[str, list[dict]]]:
    import pandas as pd
    result: dict[str, dict[str, list[dict]]] = {}
    for date, tickers in candidates.items():
        day_bars: dict[str, list[dict]] = {}
        for ticker in tickers:
            p = Path(INTRADAY_CACHE_DIR) / ticker / f"{date}.parquet"
            if not p.exists():
                continue
            try:
                df = pd.read_parquet(str(p))
                if not df.empty and len(df) >= 5:
                    day_bars[ticker] = df.to_dict("records")
            except Exception:
                pass
        if day_bars:
            result[date] = day_bars
    return result


# ── Stats Helpers ─────────────────────────────────────────────────────────────

def pearson(xs: list[float], ys: list[float]) -> float | None:
    if len(xs) != len(ys) or len(xs) < 2:
        return None
    n = len(xs)
    xm = sum(xs) / n
    ym = 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], outcomes_r: list[float]) -> dict:
    if len(vals) < 6:
        return {}
    pairs = sorted(zip(vals, outcomes_r), key=lambda p: p[0])
    n = len(pairs)
    t = n // 3
    def stats(pairs_sub):
        ys = [p[1] for p in pairs_sub]
        wins = [y for y in ys if y > 0]
        wr = len(wins) / len(ys) if ys else 0.0
        avg = sum(ys) / len(ys) if ys else 0.0
        return {"n": len(ys), "wr": wr, "avg_r": avg}
    return {
        "low": stats(pairs[:t]),
        "mid": stats(pairs[t:2*t]),
        "high": stats(pairs[2*t:]),
    }


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

def main() -> None:
    print("=== V23 Quality Feature Diagnostic (Hurst / OU-θ / OBV-Slope) ===\n")

    # 1. Load V23 config
    with open(V23_CONFIG) as f:
        raw = yaml.safe_load(f)
    params = ORBStrategyParams(**raw["orb_strategy"])
    print(f"V23 regime threshold: {params.market_regime_spy_threshold}")

    # 2. Determine 200d trading window (ending today)
    today = dt.date(2026, 4, 21)
    all_td = trading_days_between(today - dt.timedelta(days=400), today)
    trading_days = [d.isoformat() for d in all_td[-LOOKBACK_DAYS:]]
    print(f"Window: {trading_days[0]} → {trading_days[-1]} ({len(trading_days)} trading days)")

    # 3. Build date range: need FEATURE_LOOKBACK_TRADING + LOOKBACK_DAYS trading days
    extended_td = [d.isoformat() for d in all_td[-(LOOKBACK_DAYS + FEATURE_LOOKBACK_TRADING + 10):]]
    first_cal = dt.date.fromisoformat(extended_td[0])
    last_cal = dt.date.fromisoformat(trading_days[-1])

    # Build full calendar range (including weekends)
    needed_dates: list[str] = []
    d = first_cal
    while d <= last_cal:
        needed_dates.append(d.isoformat())
        d += dt.timedelta(days=1)
    print(f"Daily bar date range: {needed_dates[0]} → {needed_dates[-1]} ({len(needed_dates)} calendar days)")

    # 4. Load universe
    with open(UNIVERSE_FILE) as f:
        udata = yaml.safe_load(f)
    universe = udata.get("symbols", udata) if isinstance(udata, dict) else udata
    if "QQQ" not in universe:
        universe = list(universe) + ["QQQ"]
    print(f"\nBuilding daily bars from intraday cache for {len(universe)} tickers...")

    daily_bars = build_daily_bars(universe, needed_dates)
    print(f"Built daily bars for {len(daily_bars)} tickers")

    # 5. Compute enrichment (only over trading_days window)
    print("Computing enrichment...")
    enrichment = enrich_daily_bars(daily_bars, trading_days)
    print(f"Enrichment for {len(enrichment)} tickers")

    # 6. Pre-screen candidates + load intraday bars
    candidates = orb_pre_screen_candidates(
        daily_bars, trading_days, enrichment,
        min_price=params.min_price,
        min_atr=params.min_atr_14,
        min_avg_dollar_vol=params.min_avg_dollar_volume,
        max_per_day=None,
    )
    total_pairs = sum(len(v) for v in candidates.values())
    print(f"Pre-screened: {total_pairs} ticker-days ({total_pairs/len(trading_days):.0f} avg/day)")
    print(f"Loading intraday bars...")
    all_intraday = load_intraday_bulk(candidates)
    intraday_pairs = sum(len(v) for v in all_intraday.values())
    print(f"Loaded: {intraday_pairs} ticker-days")

    # 7. Run V23 simulation (full regime filter active — baseline)
    print("\nRunning V23 simulation...")
    state = ORBSimulationState(equity=params.initial_capital)
    day_results, _ = run_orb_simulation_with_state(
        all_intraday, trading_days, params, enrichment, state=state,
    )
    all_trades = [t for dr in day_results for t in dr.trades]
    trade_days = sum(1 for dr in day_results if dr.trades)
    print(f"Total trades: {len(all_trades)} across {trade_days} trade-days")

    # Filter trades with r_multiple_at_exit
    trades_with_r = [t for t in all_trades if getattr(t, "r_multiple_at_exit", None) is not None]
    print(f"Trades with r_multiple_at_exit: {len(trades_with_r)}/{len(all_trades)}")

    if len(trades_with_r) < 40:
        print("ABORT: fewer than 40 trades with r_multiple — insufficient sample")
        return

    # 8. Compute features for each trade
    print("\nComputing features for each trade...")
    annotated: list[dict] = []
    missing = {"hurst": 0, "ou_theta": 0, "obv_slope": 0}

    # Pre-sort daily_bars by ticker
    sorted_daily: dict[str, list[dict]] = {
        t: sorted(bars, key=lambda b: b["date"])
        for t, bars in daily_bars.items()
    }

    for trade in trades_with_r:
        ticker = trade.ticker
        date = trade.date
        r = trade.r_multiple_at_exit

        bars_for_ticker = sorted_daily.get(ticker, [])
        # Find bars strictly before trade date
        prev_bars = [b for b in bars_for_ticker if b["date"][:10] < date]

        h = compute_hurst_approx(prev_bars, lookback=60)
        ou = compute_ou_theta_approx(prev_bars, lookback=60)
        obv = compute_obv_slope_approx(prev_bars, lookback=20)

        if h is None:
            missing["hurst"] += 1
        if ou is None:
            missing["ou_theta"] += 1
        if obv is None:
            missing["obv_slope"] += 1

        annotated.append({
            "ticker": ticker, "date": date, "r": float(r),
            "hurst": h, "ou_theta": ou, "obv_slope": obv,
            "win": r > 0,
        })

    total = len(annotated)
    print(f"Annotated: {total} trades")
    print(f"Missing values: hurst={missing['hurst']}, ou_theta={missing['ou_theta']}, obv_slope={missing['obv_slope']}")

    # 9. Report per-feature
    feature_defs = [
        ("hurst_60", "hurst", "H>0.5 trending → breakout follow-through", False),
        ("ou_theta_60", "ou_theta", "Low θ = slow reversion = PEAD/breakout friendly", True),
        ("obv_slope_20", "obv_slope", "Positive = accumulation = smart-money pre-positioning", False),
    ]

    print("\n" + "=" * 90)
    print("FEATURE ANALYSIS — V23 200d trade set")
    print("=" * 90)

    results = {}
    for feat_name, feat_key, description, invert in feature_defs:
        valid = [(t[feat_key], t["r"]) for t in annotated if t[feat_key] is not None]
        if len(valid) < 30:
            print(f"\n{feat_name}: SKIP — only {len(valid)} valid trades (need ≥30)")
            results[feat_name] = None
            continue

        vals = [v[0] for v in valid]
        rs = [v[1] for v in valid]
        wins = [v for v in valid if v[1] > 0]
        overall_wr = len(wins) / len(valid)

        rho = pearson(vals, rs)
        tstat = tercile_stats(vals, rs)

        effective_high = "high" if not invert else "low"
        effective_low = "low" if not invert else "high"

        print(f"\n{'─'*60}")
        print(f"FEATURE: {feat_name}")
        print(f"  Description: {description}")
        print(f"  n={len(valid)}, overall WR={overall_wr*100:.1f}%")
        print(f"  Pearson(feature, r_multiple) = {rho:.4f}" if rho is not None else "  Pearson = n/a")

        if tstat:
            h = tstat["high"]
            m = tstat["mid"]
            l = tstat["low"]
            print(f"  Tercile breakdown (low→high feature value):")
            print(f"    Bottom: n={l['n']}, WR={l['wr']*100:.1f}%, avg_R={l['avg_r']:+.3f}")
            print(f"    Middle: n={m['n']}, WR={m['wr']*100:.1f}%, avg_R={m['avg_r']:+.3f}")
            print(f"    Top:    n={h['n']}, WR={h['wr']*100:.1f}%, avg_R={h['avg_r']:+.3f}")

            # Gate evaluation
            best = tstat[effective_high]
            worst = tstat[effective_low]
            rho_abs = abs(rho) if rho is not None else 0.0
            g1 = rho_abs >= 0.07
            g2 = best["avg_r"] - worst["avg_r"] >= 0.30
            g3 = best["wr"] >= worst["wr"] + 0.05

            print(f"  Gates ('{effective_high}' vs '{effective_low}' tercile for {feat_name}):")
            print(f"    G1 |Pearson| ≥ 0.07: {rho_abs:.4f}  → {'PASS ✓' if g1 else 'FAIL ✗'}")
            print(f"    G2 avg_R gap ≥ 0.30R: {best['avg_r'] - worst['avg_r']:+.3f}  → {'PASS ✓' if g2 else 'FAIL ✗'}")
            print(f"    G3 WR gap ≥ 5pp:     {(best['wr'] - worst['wr'])*100:+.1f}pp  → {'PASS ✓' if g3 else 'FAIL ✗'}")
            overall_pass = g1 and g2 and g3
            print(f"  VERDICT: {'ALL GATES PASS → PROCEED TO PHASE 2' if overall_pass else f'FAIL ({sum([not g1, not g2, not g3])} gates failed)'}")
            results[feat_name] = {"pass": overall_pass, "pearson": rho, "stats": tstat}

    # 10. Pairwise feature correlation
    h_vals = [t["hurst"] for t in annotated if t["hurst"] is not None and t["ou_theta"] is not None]
    ou_vals = [t["ou_theta"] for t in annotated if t["hurst"] is not None and t["ou_theta"] is not None]
    obv_h = [t["obv_slope"] for t in annotated if t["hurst"] is not None and t["obv_slope"] is not None]
    h_for_obv = [t["hurst"] for t in annotated if t["hurst"] is not None and t["obv_slope"] is not None]
    obv_ou = [t["obv_slope"] for t in annotated if t["ou_theta"] is not None and t["obv_slope"] is not None]
    ou_for_obv = [t["ou_theta"] for t in annotated if t["ou_theta"] is not None and t["obv_slope"] is not None]

    print(f"\n{'─'*60}")
    print("PAIRWISE FEATURE CORRELATIONS:")
    rho_h_ou = pearson(h_vals, ou_vals)
    rho_h_obv = pearson(h_for_obv, obv_h)
    rho_ou_obv = pearson(ou_for_obv, obv_ou)
    print(f"  Hurst vs OU-θ:     {rho_h_ou:.4f}" if rho_h_ou is not None else "  Hurst vs OU-θ:  n/a")
    print(f"  Hurst vs OBV:      {rho_h_obv:.4f}" if rho_h_obv is not None else "  Hurst vs OBV:   n/a")
    print(f"  OU-θ  vs OBV:      {rho_ou_obv:.4f}" if rho_ou_obv is not None else "  OU-θ  vs OBV:   n/a")
    print("  (|ρ| > 0.7 between any pair → use only stronger one in Phase 2)")

    # 11. Summary
    passing = [name for name, r in results.items() if r is not None and r["pass"]]
    print(f"\n{'='*90}")
    print(f"SUMMARY")
    print(f"{'='*90}")
    if passing:
        print(f"Features passing all gates: {', '.join(passing)}")
        print("VERDICT: PROCEED TO PHASE 2")
        print(f"  → Add {passing} as scoring weights in ORBStrategyParams")
        print("  → Config: orb_gainers_v24_quality_overlay.yaml (parent V23)")
        print("  → Weight allocation: 0.15 total split proportionally by |Pearson|")
    else:
        print("No features passed all gates.")
        print("VERDICT: ABORT — Hurst/OU-θ/OBV statistically neutral on V23 200d trade set")
        print("  → V23 remains terminal. Entropy-sibling axis exhausted.")
        print("  → Next Ralph iteration: FINRA short-volume axis or tape-ignition microstructure")


if __name__ == "__main__":
    main()