PythagorasGoal/scripts/research/trackF_seasonality.py

"""TRACK F — CALENDAR SEASONALITY on BTC & ETH (hour-of-day, day-of-week, interactions).

Honest test of whether there is a SYSTEMATIC, TRADEABLE calendar edge on the certified
Deribit-mainnet BTC/ETH feeds. Seasonality is the easiest place on earth to overfit
(24 hours x 7 weekdays = 168 buckets => you WILL find "significant" cells by chance), so
every claim here is held to the project's anti-look-ahead, OOS, per-year, both-assets bar.

METHODOLOGY (no shortcuts):
  - ret[i] = close[i]/close[i-1]-1 is known at close[i]. A position decided at close[i]
    earns ret[i+1]. We NEVER include the bar being traded (or any future bar) in the
    statistic that decides the trade.
  - DESCRIPTIVE tables (per-hour / per-weekday mean returns) are split IS(65%)/OOS(35%).
    They are diagnostics, not trades.
  - TRADEABLE rule = ADAPTIVE EXPANDING sign: at close[i] we look up the calendar bucket
    of bar i+1 (the clock is known with zero look-ahead) and take the SIGN of that bucket's
    mean return computed ONLY on bars <= i (expanding, warmup-gated). Long-flat or
    long-short. Fees charged only on |Δposition| (turnover-aware). This lets the data pick
    each bucket's sign LIVE — the honest analogue of "trade the seasonal".
  - Also an in-sample-optimised discrete rule (enter at hour H, hold W bars, best dir) is
    shown ONLY to demonstrate the overfit gap IS->OOS.
  - NET fees fee_side baseline 0.0005 (=0.10% RT); swept 0.0005/0.00075/0.001.
  - A survivor must be net-positive OOS AND across years AND on BOTH BTC & ETH.

Run:  uv run python scripts/research/trackF_seasonality.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  # noqa: E402

ASSETS = ["BTC", "ETH"]
TF = "1h"
FEE_SIDE = 0.0005          # 0.05%/side = 0.10% round-trip
BARS_PER_DAY = 24
BPY = BARS_PER_DAY * 365.25


# ---------------------------------------------------------------------------
# helpers
# ---------------------------------------------------------------------------
def prep(asset: str, tf: str = TF):
    df = load(asset, tf)
    c = df["close"].values.astype(float)
    ret = np.empty(len(c))
    ret[0] = 0.0
    ret[1:] = c[1:] / c[:-1] - 1.0
    dt = pd.to_datetime(df["datetime"])
    return dict(
        df=df, ret=ret,
        hour=dt.dt.hour.values.astype(int),
        dow=dt.dt.dayofweek.values.astype(int),   # 0=Mon..6=Sun
        ts=dt,
    )


def metrics_from_pnl(pnl: np.ndarray, ts: pd.Series):
    """pnl[i] = realized per-bar net return of the strategy (already fee-adjusted)."""
    eq = np.cumprod(1.0 + np.clip(pnl, -0.99, None))
    r = pnl[np.isfinite(pnl)]
    sharpe = float(np.mean(r) / np.std(r) * np.sqrt(BPY)) if np.std(r) > 0 else 0.0
    peak = np.maximum.accumulate(eq)
    maxdd = float(np.max((peak - eq) / peak)) if len(eq) else 0.0
    span_days = (ts.iloc[-1] - ts.iloc[0]).total_seconds() / 86400
    years = span_days / 365.25 if span_days > 0 else 1.0
    total = eq[-1] / eq[0] if len(eq) else 1.0
    cagr = total ** (1 / years) - 1 if years > 0 and total > 0 else -1.0
    daily_2k = (2000 * total - 2000) / span_days if span_days > 0 else 0.0
    return dict(sharpe=sharpe, maxdd=maxdd, cagr=cagr, total=total - 1.0,
                daily_2k=daily_2k, eq=eq)


def per_year_pnl(pnl: np.ndarray, ts: pd.Series):
    s = pd.Series(pnl, index=ts.values)
    out = {}
    for y, g in s.groupby(s.index.year):
        eq = np.cumprod(1.0 + np.clip(g.values, -0.99, None))
        out[int(y)] = float(eq[-1] - 1.0)
    return out


# ---------------------------------------------------------------------------
# 1. DESCRIPTIVE seasonality tables (diagnostics, IS vs OOS)
# ---------------------------------------------------------------------------
def descriptive(data, frac=0.65):
    n = len(data["ret"])
    cut = int(n * frac)
    ret, hour, dow = data["ret"], data["hour"], data["dow"]
    rows_h, rows_d = {}, {}
    for h in range(24):
        m_is = ret[:cut][hour[:cut] == h]
        m_oos = ret[cut:][hour[cut:] == h]
        rows_h[h] = (m_is.mean() * 1e4, m_oos.mean() * 1e4,
                     np.sign(m_is.mean()) == np.sign(m_oos.mean()))
    for d in range(7):
        m_is = ret[:cut][dow[:cut] == d]
        m_oos = ret[cut:][dow[cut:] == d]
        rows_d[d] = (m_is.mean() * 1e4, m_oos.mean() * 1e4,
                     np.sign(m_is.mean()) == np.sign(m_oos.mean()))
    return rows_h, rows_d


# ---------------------------------------------------------------------------
# 2. ADAPTIVE EXPANDING-sign seasonal strategy (the honest tradeable test)
# ---------------------------------------------------------------------------
def adaptive_seasonal(data, bucket="hour", mode="longshort",
                      warmup=200, fee_side=FEE_SIDE):
    """Position at close[i] = sign of the EXPANDING past mean return of bar (i+1)'s
    calendar bucket, using only bars <= i. earns ret[i+1]. Fee on |Δposition|."""
    ret = data["ret"]
    key = data[bucket]
    n = len(ret)
    nbuck = int(key.max()) + 1
    sums = np.zeros(nbuck)
    counts = np.zeros(nbuck)
    pos = np.zeros(n)
    for i in range(1, n - 1):
        b = key[i]
        sums[b] += ret[i]
        counts[b] += 1
        nb = key[i + 1]
        if counts[nb] >= warmup:
            m = sums[nb] / counts[nb]
            if m > 0:
                pos[i] = 1.0
            else:
                pos[i] = -1.0 if mode == "longshort" else 0.0
    # pnl[i] earned over bar i+1
    pnl = np.zeros(n)
    prev = 0.0
    for i in range(1, n - 1):
        turn = abs(pos[i] - prev)
        pnl[i] = pos[i] * ret[i + 1] - fee_side * turn
        prev = pos[i]
    return pnl, pos


def adaptive_hourxdow(data, mode="longshort", warmup=120, fee_side=FEE_SIDE):
    ret, hour, dow = data["ret"], data["hour"], data["dow"]
    key = hour * 7 + dow                      # 168 buckets
    n = len(ret)
    sums = np.zeros(168)
    counts = np.zeros(168)
    pos = np.zeros(n)
    for i in range(1, n - 1):
        b = key[i]
        sums[b] += ret[i]
        counts[b] += 1
        nb = key[i + 1]
        if counts[nb] >= warmup:
            m = sums[nb] / counts[nb]
            if m > 0:
                pos[i] = 1.0
            else:
                pos[i] = -1.0 if mode == "longshort" else 0.0
    pnl = np.zeros(n)
    prev = 0.0
    for i in range(1, n - 1):
        turn = abs(pos[i] - prev)
        pnl[i] = pos[i] * ret[i + 1] - fee_side * turn
        prev = pos[i]
    return pnl, pos


# ---------------------------------------------------------------------------
# 3. In-sample-optimised DISCRETE rule (to expose the overfit gap)
# ---------------------------------------------------------------------------
def discrete_hour_rule_scan(data, frac=0.65, fee_side=FEE_SIDE):
    """Scan IS for best (entry_hour, hold_window, direction) by IS Sharpe; report OOS.

    A trade: enter at close of bar whose hour==H (decided with data<=close[i]), hold W
    bars, exit at close. One trade per day. Fee charged round-trip on each trade.
    """
    ret, hour, ts = data["ret"], data["hour"], data["ts"]
    n = len(ret)
    cut = int(n * frac)

    def rule_pnl(H, W, direction, lo, hi):
        pnl = np.zeros(n)
        i = lo
        last_exit = lo - 1
        while i < hi:
            if hour[i] == H and i > last_exit:
                # cumulative return over the next W bars: prod(1+ret[i+1..i+W]) - 1
                end = min(i + W, n - 1)
                gross = np.prod(1.0 + ret[i + 1:end + 1]) - 1.0
                pnl[i] = direction * gross - 2 * fee_side
                last_exit = end
                i = end
            else:
                i += 1
        return pnl

    best = None
    n_tested = 0
    for H in range(24):
        for W in (1, 2, 3, 4, 6, 8, 12, 24):
            for direction in (+1, -1):
                n_tested += 1
                pnl_is = rule_pnl(H, W, direction, 1, cut)
                r = pnl_is[pnl_is != 0.0]
                if len(r) < 50:
                    continue
                sh = np.mean(r) / np.std(r) * np.sqrt(BPY) if np.std(r) > 0 else 0.0
                if best is None or sh > best[0]:
                    best = (sh, H, W, direction)
    sh, H, W, direction = best
    pnl_oos = rule_pnl(H, W, direction, cut, n)
    r_oos = pnl_oos[pnl_oos != 0.0]
    sh_oos = (np.mean(r_oos) / np.std(r_oos) * np.sqrt(BPY)) if (len(r_oos) and np.std(r_oos) > 0) else 0.0
    return dict(n_tested=n_tested, H=H, W=W, dir=direction, sh_is=sh,
                sh_oos=sh_oos, n_is=int((rule_pnl(H, W, direction, 1, cut) != 0).sum()),
                n_oos=len(r_oos), oos_mean_bp=r_oos.mean() * 1e4 if len(r_oos) else 0.0)


# ---------------------------------------------------------------------------
# reporting
# ---------------------------------------------------------------------------
def split_metrics(pnl, ts, frac=0.65):
    n = len(pnl)
    cut = int(n * frac)
    m_is = metrics_from_pnl(pnl[:cut], ts.iloc[:cut])
    m_oos = metrics_from_pnl(pnl[cut:], ts.iloc[cut:])
    m_all = metrics_from_pnl(pnl, ts)
    return m_is, m_oos, m_all


def turnover_per_year(pos, ts):
    s = pd.Series(np.abs(np.diff(pos, prepend=0.0)), index=ts.values)
    return s.groupby(s.index.year).sum().to_dict()


def main():
    print("=" * 100)
    print("# TRACK F — CALENDAR SEASONALITY (hour-of-day / day-of-week / hour×weekday)")
    print("# certified Deribit-mainnet BTC & ETH, 1h UTC. fee_side=0.0005 (0.10% RT).")
    print("# No look-ahead: bucket stats use only bars <= i; position earns ret[i+1].")
    print("=" * 100)

    data = {a: prep(a) for a in ASSETS}

    # --- DESCRIPTIVE ---------------------------------------------------------
    print("\n" + "#" * 100)
    print("# 1. DESCRIPTIVE per-bucket mean returns (basis points/bar). IS=first 65%, OOS=last 35%.")
    print("#    'sign?' = IS and OOS agree on sign. Diagnostics only (NOT trades, no fees).")
    print("#" * 100)
    for a in ASSETS:
        rows_h, rows_d = descriptive(data[a])
        print(f"\n  ── {a}  HOUR-OF-DAY (UTC)  mean bp/hr ─────────────────────────────")
        print("     hr :  IS_bp   OOS_bp  sign?")
        agree_h = 0
        for h in range(24):
            iv, ov, ag = rows_h[h]
            agree_h += int(ag)
            flag = "  <-- US open" if h in (13, 14) else ("  <-- US close" if h in (20, 21) else "")
            print(f"     {h:>2d} : {iv:>+6.2f}  {ov:>+6.2f}   {'Y' if ag else '.'}{flag}")
        print(f"     hour sign-agreement IS/OOS: {agree_h}/24")
        print(f"\n  ── {a}  DAY-OF-WEEK  mean bp/bar (0=Mon..6=Sun) ──────────────────")
        names = ["Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun"]
        agree_d = 0
        for d in range(7):
            iv, ov, ag = rows_d[d]
            agree_d += int(ag)
            print(f"     {names[d]} : {iv:>+6.3f}  {ov:>+6.3f}   {'Y' if ag else '.'}")
        print(f"     weekday sign-agreement IS/OOS: {agree_d}/7")

    # --- ADAPTIVE EXPANDING-SIGN (the honest tradeable test) ----------------
    print("\n" + "#" * 100)
    print("# 2. ADAPTIVE EXPANDING-SIGN seasonal strategies (HONEST tradeable test).")
    print("#    sign of bucket's PAST-ONLY mean decides position; fee on turnover.")
    print("#" * 100)
    configs = [
        ("HOUR  long-short", "hour", "longshort", 200),
        ("HOUR  long-flat ", "hour", "longflat", 200),
        ("DOW   long-short", "dow", "longshort", 60),
        ("DOW   long-flat ", "dow", "longflat", 60),
    ]
    for label, bucket, mode, warmup in configs:
        print(f"\n  ── {label} ────────────────────────────────────────────────────")
        for a in ASSETS:
            pnl, pos = adaptive_seasonal(data[a], bucket=bucket, mode=mode, warmup=warmup)
            ts = data[a]["ts"]
            m_is, m_oos, m_all = split_metrics(pnl, ts)
            py = per_year_pnl(pnl, ts)
            yrs = "".join(f"{py.get(y, float('nan'))*100:>+6.0f}" for y in range(2019, 2027))
            print(f"     {a}: ALL Sh={m_all['sharpe']:>+5.2f} CAGR={m_all['cagr']*100:>+6.1f}% "
                  f"DD={m_all['maxdd']*100:>4.1f}% €/d={m_all['daily_2k']:>+5.2f} | "
                  f"IS Sh={m_is['sharpe']:>+5.2f}  OOS Sh={m_oos['sharpe']:>+5.2f}")
            print(f"        per-year %: {yrs}   (2019..2026)")

    # buy-and-hold benchmark — the key control: does any 'seasonal' beat just being long?
    print(f"\n  ── BUY-AND-HOLD benchmark (the control for long-bias) ──")
    for a in ASSETS:
        ret = data[a]["ret"].copy()
        ret[0] = 0.0
        m = metrics_from_pnl(ret, data[a]["ts"])
        print(f"     {a}: Sh={m['sharpe']:>+5.2f} CAGR={m['cagr']*100:>+6.1f}% DD={m['maxdd']*100:>4.1f}% "
              f"  <- compare to DOW long-flat above (it's nearly identical = no edge, just long)")

    # hour x weekday interaction (168 buckets — extreme overfit risk)
    print(f"\n  ── HOUR×WEEKDAY long-short (168 buckets, warmup 120) — overfit canary ──")
    for a in ASSETS:
        pnl, pos = adaptive_hourxdow(data[a], mode="longshort", warmup=120)
        ts = data[a]["ts"]
        m_is, m_oos, m_all = split_metrics(pnl, ts)
        print(f"     {a}: ALL Sh={m_all['sharpe']:>+5.2f} CAGR={m_all['cagr']*100:>+6.1f}% "
              f"DD={m_all['maxdd']*100:>4.1f}% | IS Sh={m_is['sharpe']:>+5.2f}  OOS Sh={m_oos['sharpe']:>+5.2f}")

    # --- FEE SWEEP on the best adaptive config -------------------------------
    print("\n" + "#" * 100)
    print("# 3. FEE SWEEP — HOUR long-short adaptive (turnover-aware). Are survivors fee-robust?")
    print("#" * 100)
    for fee in (0.0, 0.0005, 0.00075, 0.001):
        line = f"  fee_side={fee:.5f} (RT {fee*2*100:.2f}%): "
        for a in ASSETS:
            pnl, _ = adaptive_seasonal(data[a], bucket="hour", mode="longshort",
                                       warmup=200, fee_side=fee)
            m = metrics_from_pnl(pnl, data[a]["ts"])
            line += f"{a} Sh={m['sharpe']:>+5.2f} CAGR={m['cagr']*100:>+6.1f}%   "
        print(line)

    # --- TURNOVER (fees are first-order for hour strategies) -----------------
    print("\n" + "#" * 100)
    print("# 4. TURNOVER (HOUR long-short adaptive): position flips/year (each flip costs ~fee).")
    print("#" * 100)
    for a in ASSETS:
        _, pos = adaptive_seasonal(data[a], bucket="hour", mode="longshort", warmup=200)
        tpy = turnover_per_year(pos, data[a]["ts"])
        s = "  ".join(f"{y}:{int(v)}" for y, v in sorted(tpy.items()))
        print(f"     {a} turnover units/yr: {s}")

    # --- IN-SAMPLE-OPTIMISED DISCRETE RULE (overfit demonstration) ----------
    print("\n" + "#" * 100)
    print("# 5. IN-SAMPLE-OPTIMISED discrete rule (enter hour H, hold W, best dir).")
    print("#    Picked by IS Sharpe, reported OOS. Demonstrates the multiple-testing trap.")
    print("#" * 100)
    for a in ASSETS:
        r = discrete_hour_rule_scan(data[a])
        print(f"     {a}: tested {r['n_tested']} (H,W,dir) cells -> best IS "
              f"H={r['H']:02d} hold={r['W']}h dir={r['dir']:+d}  "
              f"IS Sh={r['sh_is']:>+5.2f} (n={r['n_is']})  ->  OOS Sh={r['sh_oos']:>+5.2f} "
              f"(n={r['n_oos']}, mean {r['oos_mean_bp']:>+.1f} bp/trade)")

    # --- VERDICT -------------------------------------------------------------
    print("\n" + "#" * 100)
    print("# MULTIPLE-TESTING CAVEAT")
    print("#" * 100)
    print("""
  Buckets examined: 24 hours + 7 weekdays + 168 hour×weekday = 199 calendar cells PER ASSET,
  each tested IS and OOS, plus discrete grid = 24×8×2 = 384 (H,W,dir) cells per asset.
  With that many cells, spurious 'significant' buckets are GUARANTEED. The honest filters
  applied here: (a) adaptive sign chosen live on PAST data only (no cherry-picking),
  (b) must hold OOS, (c) must hold per-year, (d) must hold on BOTH BTC AND ETH.
  Read the IS->OOS Sharpe collapse and the per-year sign flips above as the real verdict.
""")


if __name__ == "__main__":
    main()