quant-trading/research/lab/letf_history.py

"""
Synthetic leveraged ETFs back to 1986, and the trend strategies run on them.

The live ETF panel starts in 2016 and the TrendRider research in 2010, so
neither saw the 2000-2002 Nasdaq crash or 2008 with 3x ETFs. This rebuilds
the TrendRider basket from its underlying indices and replays it.

    PYTHONPATH=. uv run python -m research.lab.fetch_long_history   # once
    PYTHONPATH=. uv run python -m research.lab.letf_history

Synthetic series (daily):
  SPY   S&P 500 total return: SPY (1993-), VFINX before that
  QQQ   Nasdaq-100 total return: QQQ (1999-03-), ^NDX (or FRED NASDAQ100)
        + 0.5%/yr dividends before
  UPRO / TQQQ  3 x underlying - 2 x (T-bill + 0.4%) - 0.95%/yr expense, daily reset
  SHY   2-year Treasury from DGS2 (duration 1.9); TLT from DGS30 (duration 17;
        DGS10 changes fill the 2002-2006 gap in DGS30)
  GLD   gold futures GC=F from 2000-08 (- 0.4%/yr); cash (the SHY series) before
  DBC   DBC from 2006-02; cash before (TrendRider's risk-off pick then has no
        commodity leg)

The synthetic 3x series are checked against the real TQQQ/UPRO. Strategies
are run in two fill modes that bracket the live book:
  close  weights from closes through t-1 traded at the close of t-1 (optimistic)
  lag    the same weights traded one session later (pessimistic; live trades
         the same session on the open's information)
Book: $10k, 5 bp, $2 per order, 10% resize band (as trend_rider_long.py).
"""

from __future__ import annotations

import os
import pickle
import warnings

import numpy as np
import pandas as pd

from main import simulate_live_book, strategy_positions
from strategies.base import Strategy
from trader import RESIZE_BAND, STRATEGY_REGISTRY

warnings.filterwarnings("ignore")
SRC = os.path.join("data", "research", "long_history.pkl")
OUT = os.path.join("data", "research", "letf_history.pkl")
FIN_SPREAD, EXPENSE = 0.004, 0.0095
BLOCKS = [("1990-01-01", "1999-12-31"), ("2000-01-01", "2002-12-31"), ("2003-01-01", "2007-12-31"),
          ("2008-01-01", "2009-12-31"), ("2010-01-01", "2015-12-31"), ("2016-01-01", "2026-09-25"),
          ("1990-01-01", "2026-09-25")]


def _splice(primary: pd.Series, backfill_ret: pd.Series, index: pd.DatetimeIndex) -> pd.Series:
    """Daily returns: the primary series where it exists, backfill before it."""
    rp = primary.reindex(index).pct_change(fill_method=None)
    rb = backfill_ret.reindex(index)
    first = primary.first_valid_index()
    return rb.where(index <= first, rp).fillna(0.0) if first is not None else rb.fillna(0.0)


def build(src: dict) -> tuple[pd.DataFrame, pd.DataFrame]:
    """(synthetic closes, real closes for validation) on the S&P 500 calendar from 1986."""
    c, fred = src["close"], src["fred"].ffill()
    idx = c["^GSPC"].dropna().loc["1986-01-01":].index
    rf = fred["DTB3"].reindex(idx.union(fred.index)).ffill().reindex(idx) / 100
    y2 = fred["DGS2"].reindex(idx.union(fred.index)).ffill().reindex(idx) / 100
    y30 = fred["DGS30"].reindex(idx.union(fred.index)).ffill(limit=5).reindex(idx) / 100
    y10 = fred["DGS10"].reindex(idx.union(fred.index)).ffill().reindex(idx) / 100
    dy30 = y30.diff().where(y30.notna() & y30.shift(1).notna(), y10.diff()).fillna(0.0)
    y30f = y30.fillna(y10)

    spx_bf = c["VFINX"].reindex(idx).pct_change(fill_method=None)
    spx_bf = spx_bf.fillna(c["^GSPC"].reindex(idx).pct_change(fill_method=None) + 0.03 / 252)
    spx = _splice(c["SPY"], spx_bf, idx)
    ndx_px = c["^NDX"] if "^NDX" in c else fred["NASDAQ100"]
    ndx_bf = ndx_px.reindex(idx).pct_change(fill_method=None) + 0.005 / 252
    ndx = _splice(c["QQQ"], ndx_bf, idx)

    def letf(r: pd.Series, lev: float) -> pd.Series:
        return lev * r - (lev - 1) * (rf + FIN_SPREAD) / 252 - EXPENSE / 252

    shy = y2.shift(1) / 252 - 1.9 * y2.diff().fillna(0.0)
    tlt = y30f.shift(1) / 252 - 17.0 * dy30
    gold = c["GC=F"].reindex(idx).ffill().pct_change(fill_method=None) - 0.004 / 252
    gld = gold.where(idx > c["GC=F"].first_valid_index(), shy)
    dbc = c["DBC"].reindex(idx).ffill().pct_change(fill_method=None)
    dbc = dbc.where(idx > c["DBC"].first_valid_index(), shy)
    rets = pd.DataFrame({"SPY": spx, "QQQ": ndx, "UPRO": letf(spx, 3), "TQQQ": letf(ndx, 3),
                         "SSO": letf(spx, 2), "QLD": letf(ndx, 2), "SHY": shy, "TLT": tlt,
                         "GLD": gld, "DBC": dbc}).fillna(0.0)
    syn = 100 * (1 + rets).cumprod()
    real = c.reindex(idx)[["SPY", "QQQ", "UPRO", "TQQQ", "SSO", "QLD", "SHY", "TLT", "GLD", "DBC"]]
    return syn, real


def validate(syn: pd.DataFrame, real: pd.DataFrame) -> None:
    print("synthetic vs real (overlap): daily corr, annual return real / synthetic")
    for t in ["UPRO", "TQQQ", "SSO", "QLD", "SHY", "TLT", "GLD"]:
        r = real[t].pct_change(fill_method=None)
        s = syn[t].pct_change(fill_method=None)
        ok = r.notna() & s.notna()
        r, s = r[ok], s[ok]
        yrs = len(r) / 252
        print(f"  {t:5s} {r.index[0].date()}- corr {r.corr(s):.4f}  real {((1 + r).prod() ** (1 / yrs) - 1) * 100:6.2f}%"
              f"  syn {((1 + s).prod() ** (1 / yrs) - 1) * 100:6.2f}%")


class HoldRotation(Strategy):
    """Gayed & Bilello (2016) style: hold ``on`` while the signal is above its
    ``ma``-day average, else ``off``."""

    def __init__(self, on: str, off: str = "SHY", signal: str = "SPY", ma: int = 200):
        self.on, self.off, self.signal, self.ma = on, off, signal, ma

    def generate_signals(self, data: pd.DataFrame) -> pd.DataFrame:
        s = data[self.signal]
        up = s > s.rolling(self.ma).mean()
        w = pd.DataFrame(0.0, index=data.index, columns=data.columns)
        w[self.on] = up.astype(float)
        w[self.off] = (~up).astype(float)
        w.iloc[:self.ma] = 0.0
        return w.shift(1).fillna(0.0)


class Hold(Strategy):
    def __init__(self, asset: str):
        self.asset = asset

    def generate_signals(self, data: pd.DataFrame) -> pd.DataFrame:
        w = pd.DataFrame(0.0, index=data.index, columns=data.columns)
        w[self.asset] = 1.0
        return w


def strategies() -> dict[str, Strategy]:
    out = {k: STRATEGY_REGISTRY[k]() for k in ("trend_rider_v7_vt24", "trend_rider_v7_vt36")}
    out.update({"ma200 TQQQ/SHY": HoldRotation("TQQQ"), "ma200 UPRO/SHY": HoldRotation("UPRO"),
                "ma200 SSO/SHY": HoldRotation("SSO"), "hold SPY": Hold("SPY"), "hold QQQ": Hold("QQQ"),
                "hold TQQQ": Hold("TQQQ"), "hold UPRO": Hold("UPRO")})
    return out


def run(strategy: Strategy, close: pd.DataFrame, mode: str) -> pd.Series:
    held = strategy_positions(strategy, close)
    if mode == "lag":
        held = held.shift(1).fillna(0.0)
    eq, _ = simulate_live_book(held, close, 10_000, 0.0005, 2.0, 0, 0, 50.0, RESIZE_BAND)
    return eq


def block_stats(eq: pd.Series) -> list[str]:
    r = eq.pct_change().fillna(0.0)
    out = []
    for a, b in BLOCKS:
        x = r.loc[a:b]
        if len(x) < 100:
            out.append("   n/a")
            continue
        e = (1 + x).cumprod()
        cagr = e.iloc[-1] ** (252 / len(x)) - 1
        sh = x.mean() / x.std() * np.sqrt(252) if x.std() > 0 else 0.0
        out.append(f"{cagr*100:6.1f}%/{sh:4.2f}/{((e / e.cummax()) - 1).min()*100:4.0f}%")
    return out


def main() -> None:
    with open(SRC, "rb") as fh:
        src = pickle.load(fh)
    syn, real = build(src)
    validate(syn, real)
    close = syn.loc["1989-01-01":]
    head = " ".join(f"{a[:4]}-{b[2:4]:>2s}".rjust(16) for a, b in BLOCKS)
    print(f"\nCAGR/Sharpe/MaxDD by block (synthetic ETFs; GLD = cash before 2000-09, DBC = cash before 2006-02)")
    print(f"{'':30s}{head}")
    eqs = {}
    for name, strat in strategies().items():
        for mode in ("close", "lag"):
            eq = run(strat, close, mode)
            eqs[f"{name} [{mode}]"] = eq
            print(f"{name + ' [' + mode + ']':30s}" + " ".join(s.rjust(16) for s in block_stats(eq)), flush=True)
    with open(OUT, "wb") as fh:
        pickle.dump({"close": syn, "equity": eqs}, fh)


if __name__ == "__main__":
    main()