PythagorasGoal/scripts/research/trackI_momentum_reversal.py

"""TRACK I — ALTERNATIVE MOMENTUM FORMULATIONS + LONG-HORIZON REVERSAL (BTC & ETH, >=12h).

Goal:
  (A) Find a momentum formulation that BEATS or DIVERSIFIES the canonical TP01 sign-blend
      (TSMOM 1-3-6m, vol-targeted, 50/50 BTC+ETH, 12h, Sharpe ~1.32).
  (B) Test the classic LONG-HORIZON REVERSAL effect (fade 12/18/24-month winners) as a
      potentially UNCORRELATED positive overlay, and a momentum+reversal blend.

Honest harness (mirrors src/strategies/trend_portfolio.py exactly):
  - direction decided with data <= close[i]; positions HELD next bar (pos_held[1:] = tgt[:-1]);
  - vol-target by inverse PAST-ONLY realized vol (target_vol/vol), leverage-capped;
  - NET fees 0.10% RT (0.05%/side) on turnover; fee sweep included;
  - 12h / 1d only (sub-12h is dominated by costs/overfit and a prior 4h look-ahead bug);
  - OOS 65/35 split + per-year; robustness across lookbacks AND both assets;
  - correlation vs TP01 net returns reported for EVERY candidate.

A candidate is INTERESTING only if net-positive OOS on BOTH assets AND either
  (higher portfolio Sharpe than TP01 ~1.32) OR (|corr to TP01| < ~0.3 and positive).

Run:  uv run python scripts/research/trackI_momentum_reversal.py
"""
from __future__ import annotations

import sys
from pathlib import Path

import numpy as np
import pandas as pd

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

from src.backtest.harness import load
from src.strategies.trend_portfolio import resample_tf, simple_returns, realized_vol

ASSETS = ["BTC", "ETH"]
FEE_SIDE = 0.0005          # 0.05%/side = 0.10% RT
TARGET_VOL = 0.20
LEVERAGE = 2.0
VOL_WIN_DAYS = 30
OOS_FRAC = 0.65
MONTH = 30                 # days per "month" (calendar-consistent across TFs)

# tf -> bars_per_day
TF_BPD = {"12h": 2, "1d": 1}


# ---------------------------------------------------------------------------
# data
# ---------------------------------------------------------------------------
def get_df(asset: str, tf: str) -> pd.DataFrame:
    df = load(asset, "1h")
    rule = {"12h": "12h", "1d": "1D"}[tf]
    return resample_tf(df, rule)


# ---------------------------------------------------------------------------
# vol-target machinery (identical convention to TP01)
# ---------------------------------------------------------------------------
def build_target(direction, vol, long_only):
    d = np.clip(direction, 0, None) if long_only else direction
    scal = np.where((vol > 0) & np.isfinite(vol), TARGET_VOL / vol, 0.0)
    tgt = np.clip(d * scal, -LEVERAGE, LEVERAGE)
    tgt[~np.isfinite(tgt)] = 0.0
    return tgt


def net_from_target(tgt, r, fee_side=FEE_SIDE):
    pos_held = np.zeros(len(tgt))
    pos_held[1:] = tgt[:-1]
    gross = pos_held * r
    turn = np.abs(np.diff(pos_held, prepend=0.0))
    net = gross - fee_side * turn
    net[0] = 0.0
    return np.clip(net, -0.99, None)


# ---------------------------------------------------------------------------
# DIRECTION FORMULATIONS  (each returns array in roughly [-1, 1], causal, decided <= close[i])
# ---------------------------------------------------------------------------
def _log_mom(c, h):
    """log return over h bars; nan before h."""
    m = np.full(len(c), np.nan)
    m[h:] = np.log(c[h:] / c[:-h])
    return m


def dir_signblend(c, bpd, horizons_m=(1, 3, 6)):
    """TP01 baseline: mean of sign(log return) over horizons."""
    n = len(c)
    acc = np.zeros(n); cnt = np.zeros(n)
    for hm in horizons_m:
        h = hm * MONTH * bpd
        s = np.full(n, np.nan)
        s[h:] = np.sign(c[h:] / c[:-h] - 1.0)
        v = np.isfinite(s); acc[v] += s[v]; cnt[v] += 1
    out = np.zeros(n); nz = cnt > 0; out[nz] = acc[nz] / cnt[nz]
    return out


def dir_zscore(c, bpd, horizons_m=(1, 3, 6), std_win_m=12):
    """(i) Continuous momentum: z-scored cumulative log-return, tanh-bounded, multi-horizon avg."""
    n = len(c); w = std_win_m * MONTH * bpd
    acc = np.zeros(n); cnt = np.zeros(n)
    for hm in horizons_m:
        h = hm * MONTH * bpd
        m = _log_mom(c, h)
        s = pd.Series(m)
        sd = s.rolling(w, min_periods=w // 3).std().values
        z = np.where((sd > 0) & np.isfinite(sd), m / sd, np.nan)
        d = np.tanh(z)
        v = np.isfinite(d); acc[v] += d[v]; cnt[v] += 1
    out = np.zeros(n); nz = cnt > 0; out[nz] = acc[nz] / cnt[nz]
    return out


def dir_riskadj(c, bpd, horizons_m=(1, 3, 6)):
    """(ii) Risk-adjusted momentum: h-horizon return / vol-of-that-horizon, tanh, multi-horizon."""
    n = len(c); r = simple_returns(c)
    acc = np.zeros(n); cnt = np.zeros(n)
    for hm in horizons_m:
        h = hm * MONTH * bpd
        ret = np.full(n, np.nan); ret[h:] = c[h:] / c[:-h] - 1.0
        # vol of the h-bar return = per-bar std over last h bars * sqrt(h)
        sd = pd.Series(r).rolling(h, min_periods=h // 2).std().values * np.sqrt(h)
        ra = np.where((sd > 0) & np.isfinite(sd), ret / sd, np.nan)
        d = np.tanh(ra)
        v = np.isfinite(d); acc[v] += d[v]; cnt[v] += 1
    out = np.zeros(n); nz = cnt > 0; out[nz] = acc[nz] / cnt[nz]
    return out


def _ema(c, span):
    return pd.Series(c).ewm(span=span, adjust=False).mean().values


def dir_emacross(c, bpd, pairs_m=((1, 3), (2, 6), (3, 9))):
    """(iii) EMA-cross trend: mean of sign(ema_fast - ema_slow) over calendar-day pairs."""
    n = len(c)
    acc = np.zeros(n); cnt = np.zeros(n)
    for fm, sm in pairs_m:
        ef = _ema(c, fm * MONTH * bpd)
        es = _ema(c, sm * MONTH * bpd)
        warm = sm * MONTH * bpd
        d = np.sign(ef - es)
        d[:warm] = np.nan
        v = np.isfinite(d); acc[v] += d[v]; cnt[v] += 1
    out = np.zeros(n); nz = cnt > 0; out[nz] = acc[nz] / cnt[nz]
    return out


def dir_macd(c, bpd):
    """(iii-b) Classic MACD with calendar spans (fast~1m, slow~2m, signal~0.75m): sign(macd-signal)."""
    n = len(c)
    fast = int(round(1.0 * MONTH * bpd)); slow = int(round(2.0 * MONTH * bpd))
    sig = int(round(0.75 * MONTH * bpd))
    macd = _ema(c, fast) - _ema(c, slow)
    signal = pd.Series(macd).ewm(span=sig, adjust=False).mean().values
    d = np.sign(macd - signal)
    d[:slow] = 0.0
    return d


def dir_donchian(c, bpd, n_m=2):
    """(iv) Donchian breakout (>=12h): +1 if close > prior-N max, -1 if < prior-N min, else hold."""
    n = len(c); N = n_m * MONTH * bpd
    hi = pd.Series(c).rolling(N, min_periods=N).max().shift(1).values
    lo = pd.Series(c).rolling(N, min_periods=N).min().shift(1).values
    d = np.zeros(n); state = 0.0
    for i in range(n):
        if np.isfinite(hi[i]) and c[i] >= hi[i]:
            state = 1.0
        elif np.isfinite(lo[i]) and c[i] <= lo[i]:
            state = -1.0
        d[i] = state
    return d


def dir_accel(c, bpd, horizons_m=(3, 6), lag_m=1):
    """(v) Acceleration: sign of CHANGE in momentum (mom[i] - mom[i-lag]) i.e. 2nd derivative."""
    n = len(c); lag = lag_m * MONTH * bpd
    acc = np.zeros(n); cnt = np.zeros(n)
    for hm in horizons_m:
        h = hm * MONTH * bpd
        m = _log_mom(c, h)
        dm = np.full(n, np.nan)
        dm[lag:] = m[lag:] - m[:-lag]
        d = np.sign(dm)
        v = np.isfinite(d); acc[v] += d[v]; cnt[v] += 1
    out = np.zeros(n); nz = cnt > 0; out[nz] = acc[nz] / cnt[nz]
    return out


def dir_mom12_1(c, bpd, lookbacks_m=(6, 12), skip_m=1):
    """(vi) 12-1 momentum: return from (i-L) to (i-skip), skipping the most-recent `skip` month.
    For index i (>=L): sign( c[i-skip] / c[i-L] - 1 ).  Causal (uses data <= close[i-skip])."""
    n = len(c); skip = skip_m * MONTH * bpd
    acc = np.zeros(n); cnt = np.zeros(n)
    for Lm in lookbacks_m:
        L = Lm * MONTH * bpd
        s = np.full(n, np.nan)
        # i runs L..n-1:  c[i-skip] = c[L-skip : n-skip],  c[i-L] = c[0 : n-L]
        s[L:] = np.sign(c[L - skip:n - skip] / c[:n - L] - 1.0)
        v = np.isfinite(s); acc[v] += s[v]; cnt[v] += 1
    out = np.zeros(n); nz = cnt > 0; out[nz] = acc[nz] / cnt[nz]
    return out


def make_reversal(lookbacks_m):
    """(B) long-horizon reversal: -sign of long-horizon return (short past winners)."""
    def fn(c, bpd):
        n = len(c)
        acc = np.zeros(n); cnt = np.zeros(n)
        for Lm in lookbacks_m:
            L = Lm * MONTH * bpd
            s = np.full(n, np.nan)
            s[L:] = -np.sign(c[L:] / c[:-L] - 1.0)
            v = np.isfinite(s); acc[v] += s[v]; cnt[v] += 1
        out = np.zeros(n); nz = cnt > 0; out[nz] = acc[nz] / cnt[nz]
        return out
    return fn


def make_mom_minus_rev(mom_m, rev_m, rev_w=0.5):
    """Blend: long medium-term momentum + fade very-long-term extension (weighted)."""
    def fn(c, bpd):
        n = len(c)
        mom = dir_signblend(c, bpd, horizons_m=mom_m)
        rev_fn = make_reversal(rev_m)
        rev = rev_fn(c, bpd)
        return np.clip(mom + rev_w * rev, -1.0, 1.0)
    return fn


# ---------------------------------------------------------------------------
# run a formulation -> per-asset net series, combined portfolio series, metrics
# ---------------------------------------------------------------------------
def asset_net_series(asset, tf, dir_fn, long_only, fee_side=FEE_SIDE):
    df = get_df(asset, tf); bpd = TF_BPD[tf]
    c = df["close"].values.astype(float)
    r = simple_returns(c)
    bpy = bpd * 365.25
    vol = realized_vol(r, VOL_WIN_DAYS * bpd, bpy)
    direction = dir_fn(c, bpd)
    tgt = build_target(direction, vol, long_only)
    net = net_from_target(tgt, r, fee_side)
    return pd.Series(net, index=pd.to_datetime(df["datetime"].values))


def portfolio_combo(tf, dir_fn, long_only, fee_side=FEE_SIDE):
    s = {a: asset_net_series(a, tf, dir_fn, long_only, fee_side) for a in ASSETS}
    J = pd.concat(s, axis=1, join="inner").fillna(0.0)
    combo = 0.5 * J[ASSETS[0]].values + 0.5 * J[ASSETS[1]].values
    return pd.Series(combo, index=J.index), s


def sharpe_of(series, bpy):
    r = series.values[np.isfinite(series.values)]
    return float(np.mean(r) / np.std(r) * np.sqrt(bpy)) if len(r) and np.std(r) > 0 else 0.0


def metrics_of(combo: pd.Series, bpy):
    idx = combo.index
    equity = np.cumprod(1.0 + np.clip(combo.values, -0.99, None))
    sharpe = sharpe_of(combo, bpy)
    peak = np.maximum.accumulate(equity)
    dd = float(np.max((peak - equity) / peak))
    years = (idx[-1] - idx[0]).total_seconds() / 86400 / 365.25
    total = equity[-1] / equity[0]
    cagr = total ** (1 / years) - 1 if years > 0 and total > 0 else -1.0
    eq = pd.Series(equity, index=idx)
    yearly = {}
    for y, g in eq.groupby(eq.index.year):
        if len(g) > 1 and g.iloc[0] > 0:
            v = g.values; pk = np.maximum.accumulate(v)
            yearly[int(y)] = (float(g.iloc[-1] / g.iloc[0] - 1), float(np.max((pk - v) / pk)))
    # OOS split
    k = int(len(combo) * OOS_FRAC)
    is_sh = sharpe_of(combo.iloc[:k], bpy)
    oos_sh = sharpe_of(combo.iloc[k:], bpy)
    return dict(sharpe=sharpe, max_dd=dd, cagr=cagr, total=total - 1,
                yearly=yearly, is_sharpe=is_sh, oos_sharpe=oos_sh, equity=eq)


ALL_YEARS = list(range(2018, 2027))


def fmt_yearly(yearly):
    return "".join(("   .  " if y not in yearly else f"{yearly[y][0]*100:>+6.0f}") for y in ALL_YEARS)


# ---------------------------------------------------------------------------
# main
# ---------------------------------------------------------------------------
PART_A = [
    ("baseline signblend 1-3-6m", dir_signblend),
    ("(i)   z-score cum-ret",     dir_zscore),
    ("(ii)  risk-adj momentum",   dir_riskadj),
    ("(iii) EMA-cross trend",     dir_emacross),
    ("(iii-b) MACD",              dir_macd),
    ("(iv)  Donchian breakout",   dir_donchian),
    ("(v)   acceleration",        dir_accel),
    ("(vi)  12-1 skip momentum",  dir_mom12_1),
]


def report_block(title, items, tf, long_only, tp_combo, bpy):
    mode = "LONG-FLAT" if long_only else "LONG-SHORT"
    print(f"\n{'='*112}\n  {title}  |  TF={tf}  mode={mode}\n{'='*112}")
    print(f"  {'formulation':<26s} {'Shrp':>5s} {'IS':>5s} {'OOS':>5s} {'CAGR':>6s} "
          f"{'maxDD':>6s} {'corrTP':>7s} {'aBTC':>5s} {'aETH':>5s}  per-year PnL%")
    print(f"  {'':<26s} {'':>5s} {'':>5s} {'':>5s} {'':>6s} {'':>6s} {'':>7s} {'':>5s} {'':>5s}  "
          + "".join(f"{y%100:>6d}" for y in ALL_YEARS))
    results = {}
    for name, fn in items:
        combo, sleeves = portfolio_combo(tf, fn, long_only)
        m = metrics_of(combo, bpy)
        # per-asset standalone Sharpe
        a_sh = {a: sharpe_of(sleeves[a], bpy) for a in ASSETS}
        # correlation to TP01 (aligned inner)
        J = pd.concat([combo.rename("x"), tp_combo.rename("t")], axis=1, join="inner").dropna()
        corr = float(np.corrcoef(J["x"], J["t"])[0, 1]) if len(J) > 2 else float("nan")
        print(f"  {name:<26s} {m['sharpe']:>5.2f} {m['is_sharpe']:>5.2f} {m['oos_sharpe']:>5.2f} "
              f"{m['cagr']*100:>+5.0f}% {m['max_dd']*100:>5.1f}% {corr:>7.2f} "
              f"{a_sh['BTC']:>5.2f} {a_sh['ETH']:>5.2f}  {fmt_yearly(m['yearly'])}")
        results[name] = dict(metrics=m, corr=corr, combo=combo, a_sh=a_sh)
    return results


def main():
    print("#" * 112)
    print("# TRACK I — alternative momentum formulations + long-horizon reversal (BTC&ETH, >=12h)")
    print("# vol-target 20%, lev cap 2x, fee 0.10% RT, positions +1 bar, 50/50 BTC+ETH. OOS 65/35.")
    print("#" * 112)

    for tf in ("12h", "1d"):
        bpy = TF_BPD[tf] * 365.25
        # TP01 reference combo at this TF (long-flat canonical) for correlation
        tp_combo, _ = portfolio_combo(tf, dir_signblend, long_only=True)
        tp_m = metrics_of(tp_combo, bpy)
        print(f"\n>>> TP01 reference @ {tf} (long-flat 1-3-6m): "
              f"Sharpe {tp_m['sharpe']:.2f}  IS {tp_m['is_sharpe']:.2f}  OOS {tp_m['oos_sharpe']:.2f}  "
              f"CAGR {tp_m['cagr']*100:+.0f}%  maxDD {tp_m['max_dd']*100:.1f}%")

        # PART A — long-flat (fair vs canonical) and long-short
        report_block("PART A — momentum formulations", PART_A, tf, True, tp_combo, bpy)
        if tf == "12h":
            report_block("PART A — momentum formulations (long-short)", PART_A, tf, False, tp_combo, bpy)

    # ----- PART B: reversal + blends, focus 12h -----
    tf = "12h"; bpy = TF_BPD[tf] * 365.25
    tp_combo, _ = portfolio_combo(tf, dir_signblend, long_only=True)

    rev_items = [
        ("reversal 12m",        make_reversal((12,))),
        ("reversal 18m",        make_reversal((18,))),
        ("reversal 24m",        make_reversal((24,))),
        ("reversal 12-18-24m",  make_reversal((12, 18, 24))),
    ]
    print("\n\n" + "#" * 112)
    print("# PART B — LONG-HORIZON REVERSAL (fade past winners). Must be net-positive AND uncorrelated.")
    print("#" * 112)
    revB = report_block("PART B — reversal (long-short)", rev_items, tf, False, tp_combo, bpy)
    # reversal long-flat (long past losers only) for completeness
    report_block("PART B — reversal (long-flat)", rev_items, tf, True, tp_combo, bpy)

    blend_items = [
        ("mom(1-6) - 0.5*rev(12-24)", make_mom_minus_rev((1, 3, 6), (12, 24), 0.5)),
        ("mom(1-6) - 1.0*rev(12-24)", make_mom_minus_rev((1, 3, 6), (12, 24), 1.0)),
        ("mom(1-3) - 0.5*rev(18-24)", make_mom_minus_rev((1, 3), (18, 24), 0.5)),
    ]
    report_block("PART B — momentum + reversal blend", blend_items, tf, True, tp_combo, bpy)

    # ----- COMBINED PORTFOLIO: TP01 + best diversifier -----
    print("\n\n" + "#" * 112)
    print("# COMBINED: TP01 (long-flat) + candidate diversifier, blended on net returns")
    print("#" * 112)
    tp_m = metrics_of(tp_combo, bpy)
    print(f"  TP01 alone: Sharpe {tp_m['sharpe']:.3f}  CAGR {tp_m['cagr']*100:+.0f}%  maxDD {tp_m['max_dd']*100:.1f}%")

    # candidates to try as overlay: the best A formulations + reversal variants
    overlays = {
        "z-score":            (dir_zscore, True),
        "risk-adj":           (dir_riskadj, True),
        "12-1 skip":          (dir_mom12_1, True),
        "reversal 12-18-24 LS": (make_reversal((12, 18, 24)), False),
        "reversal 24m LS":    (make_reversal((24,)), False),
    }
    for name, (fn, lo) in overlays.items():
        cand, _ = portfolio_combo(tf, fn, lo)
        J = pd.concat([tp_combo.rename("t"), cand.rename("c")], axis=1, join="inner").fillna(0.0)
        corr = float(np.corrcoef(J["t"], J["c"])[0, 1])
        for w in (0.5, 0.3, 0.2):
            mix = pd.Series((1 - w) * J["t"].values + w * J["c"].values, index=J.index)
            mm = metrics_of(mix, bpy)
            tag = f"TP01 + {w:.0%} {name}"
            print(f"  {tag:<30s} Sharpe {mm['sharpe']:.3f}  CAGR {mm['cagr']*100:+5.0f}%  "
                  f"maxDD {mm['max_dd']*100:4.1f}%  OOS {mm['oos_sharpe']:.2f}  (corr={corr:+.2f})")

    # ----- FEE SWEEP (robustness): 0.00 .. 0.40% RT -----
    print("\n\n" + "#" * 112)
    print("# FEE SWEEP — portfolio Sharpe @12h across round-trip fees (0.00-0.40% RT)")
    print("#" * 112)
    sweep = [
        ("baseline 1-3-6m (LF)",      dir_signblend, True),
        ("z-score cum-ret (LF)",      dir_zscore, True),
        ("MACD (LF)",                 dir_macd, True),
        ("mom(1-6)-0.5rev(12-24)(LF)", make_mom_minus_rev((1, 3, 6), (12, 24), 0.5), True),
        ("reversal 24m (LS)",         make_reversal((24,)), False),
    ]
    rts = [0.0, 0.0005, 0.0010, 0.0020, 0.0040]
    print(f"  {'formulation':<28s}" + "".join(f"{rt*100:>7.2f}%" for rt in rts) + "  (RT)")
    for name, fn, lo in sweep:
        row = [sharpe_of(portfolio_combo(tf, fn, lo, fee_side=rt / 2)[0], bpy) for rt in rts]
        print(f"  {name:<28s}" + "".join(f"{v:>8.2f}" for v in row))

    print("\nDone. See verdict in the script docstring / diary.")


if __name__ == "__main__":
    main()