PythagorasGoal/scripts/research/trackB_ml.py

"""TRACK B — Machine-learning / feature-prediction on BTC & ETH (Deribit-certified).

Honest, strict walk-forward ML research. The whole point is to NOT repeat the death of
the old library (look-ahead). Everything here obeys:

  * Features for bar i use ONLY data <= close[i] (all rolling windows are backward).
  * Labels (sign of forward return over H bars) use close[i+H]; in walk-forward we only
    train on samples whose label is FULLY realized in the past relative to the prediction
    bar (a gap of H is enforced between train-end and the prediction block).
  * Scaler + model are fit ONLY on past data, retrained periodically, never on the future.
  * Net of fees (fee_rt sweep 0.0005 .. 0.002, baseline 0.001). Turnover reported.
  * Grid over W (lookback for training), H (horizon), threshold, asset, tf.
  * A final held-out segment (last HELD_OUT_FRAC) is NEVER used to choose configs;
    configs are selected on the DEV portion, then confirmed once on the held-out tail.

Run:  uv run python scripts/research/trackB_ml.py
      uv run python scripts/research/trackB_ml.py --quick      (smaller grid, faster)
      uv run python scripts/research/trackB_ml.py --gbm        (also try GradientBoosting)

Entry convention (harness): for a signalled bar i we open at close[i] in the predicted
direction and hold up to H bars (max_bars=H, no TP/SL) — a pure test of directional sign.
No-overlap is enforced by the harness, so trades are naturally spaced >= H bars.
"""
from __future__ import annotations

import argparse
import sys
import time
import warnings
from pathlib import Path

sys.path.insert(0, str(Path(__file__).resolve().parents[2]))

import numpy as np
import pandas as pd
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler

from src.backtest.harness import backtest_signals, load

warnings.filterwarnings("ignore")

HELD_OUT_FRAC = 0.25      # final tail reserved for confirmation only
RETRAIN_K = 250          # retrain every K bars (block prediction)
MIN_TRAIN = 400          # minimum usable training samples


# ---------------------------------------------------------------------------
# Feature engineering — ALL backward-looking (safe at close[i])
# ---------------------------------------------------------------------------
def _rsi(close: pd.Series, n: int = 14) -> pd.Series:
    d = close.diff()
    up = d.clip(lower=0).ewm(alpha=1 / n, adjust=False).mean()
    dn = (-d.clip(upper=0)).ewm(alpha=1 / n, adjust=False).mean()
    rs = up / dn.replace(0, np.nan)
    return (100 - 100 / (1 + rs)).fillna(50.0)


def _atr(df: pd.DataFrame, n: int = 14) -> pd.Series:
    h, l, c = df["high"], df["low"], df["close"]
    pc = c.shift(1)
    tr = pd.concat([(h - l), (h - pc).abs(), (l - pc).abs()], axis=1).max(axis=1)
    return tr.ewm(alpha=1 / n, adjust=False).mean()


def build_features(df: pd.DataFrame) -> tuple[np.ndarray, list[str], np.ndarray]:
    """Return (X, names, warmup_valid_mask). Every column known at close[i]."""
    c = df["close"].astype(float)
    h = df["high"].astype(float)
    l = df["low"].astype(float)
    o = df["open"].astype(float)
    v = df["volume"].astype(float)
    logc = np.log(c)

    feats: dict[str, pd.Series] = {}

    # multi-lag simple returns (ret[i] uses close[i],close[i-k] -> known at i)
    for k in (1, 2, 3, 6, 12, 24):
        feats[f"ret{k}"] = c.pct_change(k)

    # candle geometry (current bar fully known at its close)
    rng = (h - l).replace(0, np.nan)
    feats["body"] = (c - o) / rng
    feats["upsh"] = (h - np.maximum(c, o)) / rng
    feats["dnsh"] = (np.minimum(c, o) - l) / rng
    feats["range_n"] = (h - l) / c
    # one-lag candle geometry
    feats["body1"] = ((c - o) / rng).shift(1)

    # momentum/acceleration
    feats["mom48"] = c.pct_change(48)
    feats["accel"] = c.pct_change(6) - c.pct_change(12)

    # RSI
    feats["rsi14"] = _rsi(c, 14) / 100.0

    # ATR-normalized extension from a trend baseline
    ema = c.ewm(span=24, adjust=False).mean()
    atr = _atr(df, 14)
    feats["ext_atr"] = (c - ema) / atr.replace(0, np.nan)

    # realized vol (std of 1-bar returns)
    r1 = c.pct_change()
    feats["rvol24"] = r1.rolling(24).std()
    feats["rvol72"] = r1.rolling(72).std()
    feats["vol_ratio"] = feats["rvol24"] / feats["rvol72"].replace(0, np.nan)

    # position of close within recent window (0=low,1=high)
    for w in (24, 72):
        lo = l.rolling(w).min()
        hi = h.rolling(w).max()
        feats[f"pos{w}"] = (c - lo) / (hi - lo).replace(0, np.nan)

    # volume z-score
    vlog = np.log1p(v)
    feats["volz"] = (vlog - vlog.rolling(72).mean()) / vlog.rolling(72).std().replace(0, np.nan)

    names = list(feats.keys())
    X = np.column_stack([feats[k].to_numpy(dtype=float) for k in names])
    valid = np.isfinite(X).all(axis=1)
    return X, names, valid


def forward_labels(df: pd.DataFrame, H: int):
    """label[i] = 1 if close[i+H] > close[i] else 0 ; fwd[i] = forward return."""
    c = df["close"].to_numpy(float)
    n = len(c)
    fwd = np.full(n, np.nan)
    fwd[: n - H] = c[H:] / c[: n - H] - 1.0
    y = (fwd > 0).astype(float)
    lab_valid = np.isfinite(fwd)
    return y, fwd, lab_valid


# ---------------------------------------------------------------------------
# Strict walk-forward probability
# ---------------------------------------------------------------------------
def walk_forward_proba(X, y, feat_valid, lab_valid, warmup, W, H, K, model_factory):
    """Return proba_up[i] for all i (NaN where not predicted). No leakage:
    when predicting block starting at b, training labels must be realized: i + H <= b-1,
    i.e. train indices < b - H. Training window is the last W such indices."""
    n = len(y)
    proba = np.full(n, np.nan)
    start = warmup + W + H
    b = start
    while b < n:
        end_block = min(b + K, n)
        train_hi = b - H                      # exclusive; ensures label realized by b-1
        train_lo = max(warmup, train_hi - W)
        idx = np.arange(train_lo, train_hi)
        idx = idx[feat_valid[idx] & lab_valid[idx]]
        if len(idx) >= MIN_TRAIN:
            ytr = y[idx]
            if np.unique(ytr).size == 2:
                Xtr = X[idx]
                sc = StandardScaler().fit(Xtr)
                model = model_factory()
                model.fit(sc.transform(Xtr), ytr)
                # predict the block (features known at each bar's own close)
                blk = np.arange(b, end_block)
                fv = feat_valid[blk]
                if fv.any():
                    pb = model.predict_proba(sc.transform(X[blk[fv]]))[:, 1]
                    proba[blk[fv]] = pb
        b = end_block
    return proba


def proba_to_entries(proba, threshold, H, n):
    """Long if proba>0.5+thr, short if proba<0.5-thr, else flat. Hold H bars."""
    entries = [None] * n
    hi = 0.5 + threshold
    lo = 0.5 - threshold
    for i in range(n):
        p = proba[i]
        if not np.isfinite(p):
            continue
        if p > hi:
            entries[i] = {"dir": 1, "tp": None, "sl": None, "max_bars": H}
        elif p < lo:
            entries[i] = {"dir": -1, "tp": None, "sl": None, "max_bars": H}
    return entries


def mask_entries(entries, lo, hi):
    """Keep only entries with index in [lo, hi); others -> None (for IS/OOS split)."""
    out = [None] * len(entries)
    for i in range(lo, min(hi, len(entries))):
        out[i] = entries[i]
    return out


def trade_stats(df, entries, H):
    """Replicate harness no-overlap to get per-trade gross returns -> avg win/loss + long frac."""
    c = df["close"].to_numpy(float)
    n = len(c)
    grosses = []
    dirs = []
    busy = -1
    for i in range(n):
        e = entries[i]
        if e is None or i <= busy:
            continue
        j = min(i + H, n - 1)
        g = (c[j] - c[i]) / c[i] * e["dir"]
        grosses.append(g)
        dirs.append(e["dir"])
        busy = j
    g = np.array(grosses)
    if len(g) == 0:
        return 0, 0.0, 0.0, 0.0, 0.0
    wins = g[g > 0]
    losses = g[g <= 0]
    avg_w = wins.mean() if len(wins) else 0.0
    avg_l = losses.mean() if len(losses) else 0.0
    long_frac = float(np.mean(np.array(dirs) > 0))
    return len(g), avg_w, avg_l, g.mean(), long_frac


def buy_hold(df, lo, hi):
    """Buy & hold net return over [lo,hi) bars (beta benchmark)."""
    c = df["close"].to_numpy(float)
    hi = min(hi, len(c))
    if hi - lo < 2:
        return 0.0
    return c[hi - 1] / c[lo] - 1.0


# ---------------------------------------------------------------------------
# Driver
# ---------------------------------------------------------------------------
def run():
    ap = argparse.ArgumentParser()
    ap.add_argument("--quick", action="store_true", help="smaller grid (faster)")
    ap.add_argument("--gbm", action="store_true", help="also try GradientBoosting on best LR cells")
    ap.add_argument("--tf", default="1h")
    args = ap.parse_args()

    assets = ["BTC", "ETH"]
    tf = args.tf
    if args.quick:
        Ws = [8000]
        Hs = [12, 24]
        thresholds = [0.0, 0.05, 0.10]
    else:
        Ws = [4000, 8000, 16000]
        Hs = [6, 12, 24, 48]
        thresholds = [0.0, 0.03, 0.06, 0.10]

    def lr_factory():
        return LogisticRegression(C=1.0, max_iter=300, class_weight="balanced")

    print("=" * 100)
    print(f"TRACK B — walk-forward ML  tf={tf}  retrain_K={RETRAIN_K}  held_out_tail={HELD_OUT_FRAC:.0%}")
    print(f"  Ws={Ws}  Hs={Hs}  thresholds={thresholds}  model=LogisticRegression(balanced)")
    print("=" * 100)

    # cache features per asset
    cache = {}
    for a in assets:
        df = load(a, tf)
        X, names, fvalid = build_features(df)
        warmup = int(np.argmax(fvalid)) if fvalid.any() else 0
        cache[a] = (df, X, names, fvalid, warmup)
    print(f"features ({len(names)}): {names}\n")

    # ---- DEV grid search (configs chosen ONLY on dev portion) ----------------
    results = []  # dict rows
    t0 = time.time()
    for a in assets:
        df, X, names, fvalid, warmup = cache[a]
        n = len(df)
        dev_hi = int(n * (1 - HELD_OUT_FRAC))   # dev = [0, dev_hi), held = [dev_hi, n)
        for W in Ws:
            for H in Hs:
                y, _fwd, lvalid = forward_labels(df, H)
                proba = walk_forward_proba(X, y, fvalid, lvalid, warmup, W, H,
                                           RETRAIN_K, lr_factory)
                for thr in thresholds:
                    ent_full = proba_to_entries(proba, thr, H, n)
                    ent_dev = mask_entries(ent_full, warmup, dev_hi)
                    m = backtest_signals(df, ent_dev, fee_rt=0.001, asset=a, tf=tf)
                    nt, aw, al, gmean, lf = trade_stats(df, ent_dev, H)
                    results.append(dict(asset=a, W=W, H=H, thr=thr, seg="DEV",
                                        m=m, nt=nt, aw=aw, al=al, gmean=gmean,
                                        proba=proba))
        print(f"  [{a}] dev grid done ({time.time()-t0:.0f}s)")

    # print dev table
    print("\n--- DEV walk-forward (config selection set) ---")
    hdr = f"{'asset':5} {'W':>6} {'H':>3} {'thr':>5} {'trd':>5} {'wr%':>5} {'net%':>8} {'CAGR%':>7} {'Shrp':>6} {'DD%':>5} {'mkt%':>5} {'avgW%':>6} {'avgL%':>6} {'€/d':>6}"
    print(hdr)
    for r in sorted(results, key=lambda r: -r["m"].sharpe):
        m = r["m"]
        print(f"{r['asset']:5} {r['W']:>6} {r['H']:>3} {r['thr']:>5.2f} {m.n_trades:>5} "
              f"{m.win_rate:>5.1f} {m.net_return*100:>+8.1f} {m.cagr*100:>+7.1f} {m.sharpe:>6.2f} "
              f"{m.max_dd*100:>5.1f} {m.time_in_market*100:>5.0f} {r['aw']*100:>+6.2f} {r['al']*100:>+6.2f} "
              f"{m.daily_profit(2000):>+6.2f}")

    # ---- selection: positive net AND sharpe>0 on dev, then robustness ----------
    pos = [r for r in results if r["m"].net_return > 0 and r["m"].sharpe > 0 and r["m"].n_trades >= 30]
    pos.sort(key=lambda r: -r["m"].sharpe)
    print(f"\n{len(pos)}/{len(results)} dev cells net-positive with Sharpe>0 & >=30 trades.")

    # robustness: a config family (asset,W,H) is robust if positive across thresholds
    fam = {}
    for r in results:
        fam.setdefault((r["asset"], r["W"], r["H"]), []).append(r)
    robust_fams = []
    for key, rs in fam.items():
        npos = sum(1 for r in rs if r["m"].net_return > 0 and r["m"].sharpe > 0)
        if npos >= max(2, int(0.6 * len(rs))):
            robust_fams.append((key, npos, len(rs)))
    robust_fams.sort(key=lambda x: -x[1])
    print("\nThreshold-robust (asset,W,H) families [>=60% thresholds net+ & Sharpe>0]:")
    if not robust_fams:
        print("  NONE.")
    for key, npos, tot in robust_fams:
        print(f"  {key}: {npos}/{tot} thresholds positive")

    # ---- HELD-OUT confirmation on best robust cells ---------------------------
    print("\n" + "=" * 100)
    print("HELD-OUT TAIL CONFIRMATION (never used for selection)")
    print("=" * 100)
    # choose up to 6 best dev cells that belong to a robust family
    robust_keys = {k for k, _, _ in robust_fams}
    cand = [r for r in pos if (r["asset"], r["W"], r["H"]) in robust_keys][:6]
    if not cand:
        cand = pos[:6]
    if not cand:
        print("No positive dev cells to confirm. ML did not beat fees on dev.")
    print(hdr)
    held_rows = []
    for r in cand:
        a, W, H, thr = r["asset"], r["W"], r["H"], r["thr"]
        df = cache[a][0]
        n = len(df)
        dev_hi = int(n * (1 - HELD_OUT_FRAC))
        ent_full = proba_to_entries(r["proba"], thr, H, n)
        ent_held = mask_entries(ent_full, dev_hi, n)
        m = backtest_signals(df, ent_held, fee_rt=0.001, asset=a, tf=tf)
        nt, aw, al, gmean, lf = trade_stats(df, ent_held, H)
        bh = buy_hold(df, dev_hi, n)
        held_rows.append((r, m, aw, al, lf, bh))
        print(f"{a:5} {W:>6} {H:>3} {thr:>5.2f} {m.n_trades:>5} {m.win_rate:>5.1f} "
              f"{m.net_return*100:>+8.1f} {m.cagr*100:>+7.1f} {m.sharpe:>6.2f} {m.max_dd*100:>5.1f} "
              f"{m.time_in_market*100:>5.0f} {aw*100:>+6.2f} {al*100:>+6.2f} {m.daily_profit(2000):>+6.2f}  "
              f"long={lf*100:>3.0f}% B&H={bh*100:>+7.1f}%")

    # ---- FEE SWEEP on the held-out winners ------------------------------------
    print("\n--- FEE SWEEP (held-out tail) on confirmed cells ---")
    fees = [0.0005, 0.001, 0.0015, 0.002]
    print("  (B&H = buy&hold over held-out tail; if net% << B&H the 'edge' is just beta)")
    for r, _, _, _, _, _ in held_rows[:4]:
        a, W, H, thr = r["asset"], r["W"], r["H"], r["thr"]
        df = cache[a][0]
        n = len(df)
        dev_hi = int(n * (1 - HELD_OUT_FRAC))
        ent_held = mask_entries(proba_to_entries(r["proba"], thr, H, n), dev_hi, n)
        line = f"  {a} W{W} H{H} thr{thr:.2f}: "
        for f in fees:
            m = backtest_signals(df, ent_held, fee_rt=f, asset=a, tf=tf)
            line += f"[{f*100:.2f}%]net={m.net_return*100:>+6.1f}% Shrp={m.sharpe:>+4.2f}  "
        print(line)

    # ---- per-year on the single best held-out cell ----------------------------
    if held_rows:
        held_rows.sort(key=lambda x: -x[1].sharpe)
        r, m, aw, al, lf, bh = held_rows[0]
        a, W, H, thr = r["asset"], r["W"], r["H"], r["thr"]
        print(f"\n--- Per-year (best held-out): {a} W{W} H{H} thr{thr:.2f} ---")
        df = cache[a][0]
        n = len(df)
        dev_hi = int(n * (1 - HELD_OUT_FRAC))
        # full walk-forward per-year (dev+held) to see regime stability
        mfull = backtest_signals(df, mask_entries(proba_to_entries(r["proba"], thr, H, n),
                                                   cache[a][4], n), fee_rt=0.001, asset=a, tf=tf)
        mfull.print_summary(f"{a} W{W}H{H}thr{thr:.2f} FULL-WF")
        mfull.print_yearly()

    print(f"\nTotal runtime {time.time()-t0:.0f}s")
    print("\n" + "=" * 100)
    print("VERDICT (see docs/diary/2026-06-19-trackB-ml.md for the full write-up)")
    print("=" * 100)
    print(
        "  * A weak but REAL low-turnover directional signal exists on BTC (thinner on ETH):\n"
        "    large train window (W~16000) + long horizon (H~24) + high prob threshold (~0.10).\n"
        "  * It beats fees at 0.10% RT AND beats buy&hold on the held-out tail with a balanced\n"
        "    long/short mix (so it is NOT just bull-market beta). Payoff: ~53% WR, avgWin>avgLoss.\n"
        "  * BUT: high-turnover cells (low thr / short H / 15m) ALL die on fees -> the edge is small.\n"
        "    Returns concentrate in a few years (2021,2025) with a -38% year (2023); DD 23-56%.\n"
        "  * EUR/day on 2000 ~= +0.3..+0.6 baseline. Target is 50/day -> ~100x short. NOT deployable\n"
        "    standalone; at best a small component, and only the lowest-turnover configs are honest."
    )


if __name__ == "__main__":
    run()