feat(harvest): premium-harvesting tested — simple 60/40 beats the sophisticated system

The strategy class we never tried: harvest structural premia + risk system, not predict.
Diversified futures + risk-parity/vol-target/trend; clean 5-asset-class version (ES/ZN/GC/CL/BTC).
Result: plain 60/40 (Sharpe +0.72, 2010-2026) BEATS risk-parity+vol-target (+0.46) and RP+trend
(+0.33, trend hurts); 5-asset+crypto RP only ties 60/40 and loses to buy-hold equity. Meta-pattern
now complete in BOTH games: simple beats/equals sophisticated in prediction AND harvesting.
Constructive deliverable: a simple premium harvest (60/40 / risk-parity) IS a real deployable
robust strategy (~0.5-0.72 Sharpe, low DD, no prediction, minimal complexity). The engine's
sophistication was never the return-generator -- deployable path is light (harvest + risk overlay),
engine's value is infra/discipline/product. (Also: CAISO intraday gate blocked by OASIS plumbing.)

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

scripts/surfer/energy_intraday_gate.py

@@ -0,0 +1,98 @@
+#!/usr/bin/env python3
+"""Intraday gate: on CAISO real-time prices (less seasonal, forecast-error-driven), does the
+ENGINE finally beat the heuristics — unlike day-ahead, where climatology won?
+
+Tests: (1) is RT more volatile than DA (more battery value)? (2) RT-price forecastability:
+baseline "RT=DA" vs climatology vs ML(walk-fwd, features incl. DA price + DART lags) ->
+capture % of perfect-foresight RT battery value. If ML >> baselines on RT, the engine earns
+its keep in the less-seasonal market. If nothing beats "RT=DA", intraday deviations are noise.
+"""
+import datetime
+import json
+import math
+import os
+import sys
+
+import numpy as np
+
+sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
+from sklearn.ensemble import HistGradientBoostingRegressor  # noqa: E402
+
+ETA = 0.85; HRS = 4; ed = math.sqrt(ETA)
+OUT = "data/surfer/caiso"
+
+
+def load():
+    dam = json.load(open(f"{OUT}/dam.json")); rtm = json.load(open(f"{OUT}/rtm.json"))
+    keys = sorted(set(dam) & set(rtm))                          # "YYYY-MM-DDHHH"
+    # group into days with full 24h
+    byday = {}
+    for k in keys:
+        d, h = k[:10], int(k[-2:])
+        byday.setdefault(d, {})[h] = (dam[k], rtm[k])
+    rows, dord = [], []
+    for d in sorted(byday):
+        if len(byday[d]) == 24:
+            da = [byday[d][h][0] for h in range(1, 25)]
+            rt = [byday[d][h][1] for h in range(1, 25)]
+            rows.append((da, rt)); dord.append(d)
+    DA = np.array([r[0] for r in rows]); RT = np.array([r[1] for r in rows])
+    D = [datetime.date.fromisoformat(x) for x in dord]
+    return DA, RT, D
+
+
+def dispatch(yhat, actual):
+    ch = np.argsort(yhat)[:HRS]; dis = np.argsort(yhat)[-HRS:]
+    return actual[dis].sum() * ed - actual[ch].sum() / ed
+
+
+def main():
+    DA, RT, D = load()
+    ND = len(DA)
+    dow = np.array([d.weekday() for d in D]); month = np.array([d.month for d in D])
+    sprDA = DA.max(1) - DA.min(1); sprRT = RT.max(1) - RT.min(1)
+    foreDA = np.array([np.sort(DA[d])[-HRS:].sum() * ed - np.sort(DA[d])[:HRS].sum() / ed for d in range(ND)])
+    foreRT = np.array([np.sort(RT[d])[-HRS:].sum() * ed - np.sort(RT[d])[:HRS].sum() / ed for d in range(ND)])
+    print(f"CAISO NP15 {ND} days ({D[0]}..{D[-1]})")
+    print(f"daily spread: DA mean ${sprDA.mean():.0f}/MWh  RT mean ${sprRT.mean():.0f}/MWh  (RT/DA = {sprRT.mean()/sprDA.mean():.1f}x)")
+    print(f"perfect-foresight battery: DA ${foreDA.mean()*365/1000:.0f}k/yr/MW  RT ${foreRT.mean()*365/1000:.0f}k/yr/MW")
+
+    # RT forecasters
+    roll7 = np.full_like(RT, np.nan)
+    for d in range(7, ND):
+        roll7[d] = RT[d - 7:d].mean(0)
+    dart = RT - DA
+    def feats(d):
+        X = np.zeros((24, 9))
+        for h in range(24):
+            X[h] = [h, dow[d], month[d], DA[d, h], RT[d - 1, h], RT[d - 7, h],
+                    dart[d - 1, h], DA[d].mean(), roll7[d, h]]
+        return X
+    ml = np.full_like(RT, np.nan); INIT, STEP = 60, 30
+    Xc = {d: feats(d) for d in range(7, ND)}
+    for s in range(INIT, ND, STEP):
+        e = min(s + STEP, ND)
+        Xtr = np.vstack([Xc[d] for d in range(7, s)]); ytr = np.concatenate([RT[d] for d in range(7, s)])
+        gb = HistGradientBoostingRegressor(max_depth=4, max_iter=200, learning_rate=0.05, min_samples_leaf=40).fit(Xtr, ytr)
+        for d in range(s, e):
+            ml[d] = gb.predict(Xc[d])
+    clim = np.full_like(RT, np.nan)
+    for d in range(14, ND):
+        wk = dow[d] >= 5
+        past = [k for k in range(max(0, d - 28), d) if (dow[k] >= 5) == wk]
+        if past:
+            clim[d] = RT[past].mean(0)
+
+    valid = np.arange(INIT, ND)
+    fc = {"baseline RT=DA": DA, "climatology": clim, "ML walk-fwd": ml}
+    print(f"\n===== INTRADAY (RT) CAPTURE GATE — {len(valid)} days, perfect-foresight RT ${foreRT[valid].mean()*365/1000:.0f}k/yr/MW =====")
+    print(f"{'RT forecaster':>16} {'capture%':>8} {'EURk/yr/MW':>11}")
+    for nm, yh in fc.items():
+        vv = np.array([dispatch(yh[d], RT[d]) for d in valid]); cap = vv.sum() / foreRT[valid].sum()
+        print(f"{nm:>16} {100*cap:>7.0f}% {vv.mean()*365/1000:>10.0f}")
+    print("\nVERDICT: ML capture >> 'RT=DA' baseline AND > climatology = the engine FINALLY earns its keep")
+    print("(less-seasonal RT is where forecasting/RL beats heuristics). If ML ~ baselines, intraday is noise too.")
+
+
+if __name__ == "__main__":
+    main()

scripts/surfer/fetch_caiso.py

@@ -0,0 +1,83 @@
+#!/usr/bin/env python3
+"""Fetch CAISO day-ahead (DAM) + real-time (RTM) hourly LMP for one hub (free OASIS, no key).
+
+The RT-vs-DA spread is the forecast-error-driven, less-seasonal intraday opportunity. Monthly
+chunks with backoff (OASIS rate-limits hard). Aggregates RT to hourly. Caches per chunk.
+Saves data/surfer/caiso/{dam,rtm}.json as {hour_key: price}.
+"""
+import io
+import json
+import os
+import time
+import urllib.request
+import zipfile
+
+NODE = "TH_NP15_GEN-APND"
+OUT = "data/surfer/caiso"
+MONTHS = [(2024, m) for m in range(1, 13)]
+
+
+def oasis(qn, market, y, m):
+    s = f"{y}{m:02d}01T08:00-0000"
+    ny, nm = (y + 1, 1) if m == 12 else (y, m + 1)
+    e = f"{ny}{nm:02d}01T08:00-0000"
+    u = (f"http://oasis.caiso.com/oasisapi/SingleZip?queryname={qn}&startdatetime={s}"
+         f"&enddatetime={e}&version=1&market_run_id={market}&node={NODE}&resultformat=6")
+    for a in range(5):
+        try:
+            raw = urllib.request.urlopen(urllib.request.Request(u, headers={"User-Agent": "curl/8"}), timeout=90).read()
+            z = zipfile.ZipFile(io.BytesIO(raw))
+            txt = z.read(z.namelist()[0]).decode()
+            return txt
+        except Exception as ex:
+            time.sleep(12 * (a + 1))
+    return None
+
+
+def parse_hourly(txt):
+    lines = txt.strip().split("\n")
+    hdr = lines[0].split(",")
+    if "LMP_TYPE" not in hdr:
+        return {}                                              # error report / wrong format
+    iT = hdr.index("LMP_TYPE")
+    iV = hdr.index("MW") if "MW" in hdr else (hdr.index("PRC") if "PRC" in hdr else -1)
+    if iV < 0:
+        return {}
+    iD = hdr.index("OPR_DT"); iH = hdr.index("OPR_HR")
+    agg = {}
+    for ln in lines[1:]:
+        c = ln.split(",")
+        if len(c) <= max(iT, iV, iD, iH):
+            continue
+        if c[iT] != "LMP":
+            continue
+        k = f"{c[iD]}H{int(c[iH]):02d}"
+        agg.setdefault(k, []).append(float(c[iV]))
+    return {k: sum(v) / len(v) for k, v in agg.items()}
+
+
+def fetch(market, qn):
+    os.makedirs(OUT, exist_ok=True)
+    out = {}
+    for y, m in MONTHS:
+        cf = f"{OUT}/{market}_{y}{m:02d}.json"
+        if os.path.exists(cf):
+            out.update(json.load(open(cf))); continue
+        txt = oasis(qn, market, y, m)
+        if not txt:
+            print(f"  {market} {y}-{m:02d}: FAILED"); continue
+        h = parse_hourly(txt)
+        json.dump(h, open(cf, "w")); out.update(h)
+        print(f"  {market} {y}-{m:02d}: {len(h)} hrs")
+        time.sleep(6)
+    return out
+
+
+def main():
+    dam = fetch("DAM", "PRC_LMP"); rtm = fetch("RTPD", "PRC_RTPD_LMP")   # 15-min RT pre-dispatch
+    json.dump(dam, open(f"{OUT}/dam.json", "w")); json.dump(rtm, open(f"{OUT}/rtm.json", "w"))
+    print(f"DONE: DAM {len(dam)} hrs, RTM {len(rtm)} hrs")
+
+
+if __name__ == "__main__":
+    main()

scripts/surfer/portfolio_harvest.py

@@ -0,0 +1,90 @@
+#!/usr/bin/env python3
+"""The strategy class we never tested: HARVEST structural premia with the risk system, not predict.
+
+Diversified multi-asset futures (equity/rates/metals/energy/ags/FX). Build a ladder:
+equal-weight -> inverse-vol (risk-parity) -> + vol-target -> + trend-overlay (long-flat by 252d trend,
+the crisis-alpha de-risk). Measure RISK-ADJUSTED (Sharpe, ann-vol, max-DD, Calmar, per-year) vs
+benchmarks (equal-weight, equity-only, 60/40). No prediction/alpha needed — just disciplined
+harvesting + risk control (what the engine is actually good at). Roll-zeroed returns (understates
+carry; fine for comparing the risk-management strategies).
+"""
+import math
+import os
+import sys
+
+import numpy as np
+
+sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
+from signal_sweep import load_panel, build_returns  # noqa: E402
+
+TARGET_VOL = 0.10   # 10% annualized portfolio vol
+
+
+def metrics(r):
+    r = r[np.isfinite(r)]
+    if len(r) < 50 or r.std() == 0:
+        return dict(sr=float("nan"), ann=float("nan"), vol=float("nan"), dd=float("nan"), calmar=float("nan"))
+    ann = r.mean() * 252; vol = r.std() * math.sqrt(252); sr = ann / vol
+    eq = np.cumprod(1 + r); dd = float((eq / np.maximum.accumulate(eq) - 1).min())
+    return dict(sr=sr, ann=ann, vol=vol, dd=dd, calmar=ann / (abs(dd) + 1e-9))
+
+
+def main():
+    roots, days, close, op, inst, weekday = load_panel()
+    lc, R, ov, intr = build_returns(close, op, inst)
+    T, N = close.shape
+    year = (1970 + days / 365.25).astype(int)
+    vol63 = np.full_like(lc, np.nan)
+    for t in range(63, T):
+        vol63[t] = np.nanstd(R[t - 63:t], axis=0)
+    iv = 1.0 / np.where(vol63 > 0, vol63, np.nan)
+    avail = np.isfinite(close) & np.isfinite(vol63) & (vol63 > 0)
+
+    def port(weights):
+        w = np.nan_to_num(weights)
+        w = np.where(avail, w, 0.0)
+        g = np.sum(np.abs(w), axis=1, keepdims=True); g[g == 0] = 1
+        w = w / g                                              # unit gross (fully invested)
+        return np.sum(w[:-1] * R[1:], axis=1)
+
+    def voltarget(r):
+        out = np.zeros_like(r)
+        for t in range(63, len(r)):
+            rv = np.std(r[t - 63:t]) * math.sqrt(252)
+            out[t] = r[t] * min(3.0, TARGET_VOL / (rv + 1e-9))   # scale to target vol, cap 3x leverage
+        return out
+
+    # weight schemes
+    eqw = np.where(avail, 1.0, 0.0)
+    rp = np.where(avail, iv, 0.0)                              # inverse-vol risk parity
+    trend = (np.full_like(lc, np.nan)); trend[252:] = lc[252:] - lc[:-252]
+    rp_trend = np.where(avail & (trend > 0), iv, 0.0)         # trend-overlay: hold only uptrending assets
+
+    books = {
+        "equal-weight": port(eqw),
+        "risk-parity (inv-vol)": port(rp),
+        "risk-parity + vol-target": voltarget(port(rp)),
+        "RP + trend + vol-target": voltarget(port(rp_trend)),
+    }
+    # benchmarks
+    es = roots.index("ES") if "ES" in roots else 0
+    bench = {"equity-only (ES)": R[1:, es]}
+    if "ZN" in roots:
+        zn = roots.index("ZN")
+        w6040 = np.zeros((T, N)); w6040[:, es] = 0.6; w6040[:, zn] = 0.4
+        bench["60/40 (ES/ZN)"] = port(w6040)
+
+    print(f"\n===== PREMIUM HARVEST (diversified futures, {N} assets, {T} days, target {int(TARGET_VOL*100)}% vol) =====")
+    print(f"{'strategy':>26} {'Sharpe':>7} {'ann%':>6} {'vol%':>6} {'maxDD%':>7} {'Calmar':>7}")
+    for nm, r in {**bench, **books}.items():
+        m = metrics(r)
+        print(f"{nm:>26} {m['sr']:>+7.2f} {100*m['ann']:>+6.1f} {100*m['vol']:>5.1f} {100*m['dd']:>+7.1f} {m['calmar']:>7.2f}")
+    print("\nper-year Sharpe (RP + trend + vol-target):")
+    rtv = books["RP + trend + vol-target"]
+    print("  " + " ".join(f"{y}:{metrics(rtv[year[1:]==y])['sr']:+.1f}" for y in range(2010, 2027) if (year[1:] == y).sum() > 100))
+    print("\nVERDICT: if risk-managed harvest beats benchmarks on Sharpe/Calmar (esp lower maxDD) = the engine's")
+    print("real job (risk management of real premia) delivers, no alpha needed. Modest but robust = the honest win.")
+
+
+if __name__ == "__main__":
+    main()

scripts/surfer/portfolio_harvest2.py

@@ -0,0 +1,98 @@
+#!/usr/bin/env python3
+"""Clean premium-harvest: risk-parity across ~5 broad ASSET CLASSES (not 39 instruments).
+
+Equity (ES), bonds (ZN), gold (GC), energy (CL), crypto (BTC) — genuinely uncorrelated premia.
+Risk-parity (inv-vol) + vol-target + optional trend-overlay, vs 60/40 and buy-and-hold equity.
+Runs 4-asset (2010-2026, longer) and 5-asset (+BTC, 2019-2026). The fair test of whether the
+risk system adds value over plain 60/40 when applied to the RIGHT universe.
+Futures = roll-zeroed price returns (understates commodity carry); BTC = clean close-close.
+"""
+import math
+import os
+import sys
+
+import numpy as np
+
+sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
+from signal_sweep import load_panel, build_returns  # noqa: E402
+import pit_sweep  # noqa: E402
+
+TVOL = 0.10
+
+
+def metrics(r):
+    r = np.asarray(r); r = r[np.isfinite(r)]
+    if len(r) < 50 or r.std() == 0:
+        return dict(sr=float("nan"), ann=float("nan"), dd=float("nan"), cal=float("nan"))
+    ann = r.mean() * 252; vol = r.std() * math.sqrt(252)
+    eq = np.cumprod(1 + r); dd = float((eq / np.maximum.accumulate(eq) - 1).min())
+    return dict(sr=ann / vol, ann=ann, dd=dd, cal=ann / (abs(dd) + 1e-9))
+
+
+def build():
+    roots, fdays, close, op, inst = load_panel()[:5]
+    lc, R, ov, intr = build_returns(close, op, inst)
+    fcol = {r: roots.index(r) for r in ["ES", "ZN", "GC", "CL"] if r in roots}
+    fut = {r: {int(fdays[t]): R[t, c] for t in range(len(fdays))} for r, c in fcol.items()}
+    syms, cdays, cc, cqv, cf = pit_sweep.load()
+    j = syms.index("BTCUSDT"); lcb = np.log(cc[:, j])
+    btc = {int(cdays[t]): (lcb[t] - lcb[t - 1]) for t in range(1, len(cdays)) if np.isfinite(lcb[t]) and np.isfinite(lcb[t - 1])}
+    return fut, btc
+
+
+def run(assets, rets, label):
+    days = np.array(sorted(set.intersection(*[set(rets[a]) for a in assets])))
+    M = np.array([[rets[a][int(d)] for a in assets] for d in days])   # [T, K]
+    M = np.where(np.isfinite(M), M, 0.0)
+    T, K = M.shape
+    year = (1970 + days / 365.25).astype(int)
+    vol = np.full((T, K), np.nan)
+    for t in range(63, T):
+        vol[t] = M[t - 63:t].std(0)
+    iv = 1.0 / np.where(vol > 0, vol, np.nan)
+    trend = np.full((T, K), np.nan)
+    L = 126
+    for t in range(L, T):
+        trend[t] = M[t - L:t].sum(0)
+
+    def book(w):
+        w = np.nan_to_num(w); g = np.abs(w).sum(1, keepdims=True); g[g == 0] = 1
+        w = w / g
+        return np.sum(w[:-1] * M[1:], axis=1)
+
+    def vt(r):
+        out = np.zeros_like(r)
+        for t in range(63, len(r)):
+            rv = np.std(r[t - 63:t]) * math.sqrt(252)
+            out[t] = r[t] * min(3.0, TVOL / (rv + 1e-9))
+        return out
+
+    rp = book(iv)
+    rp_t = book(np.where(trend > 0, iv, 0.0))
+    es_i = assets.index("ES") if "ES" in assets else 0
+    bench_eq = M[1:, es_i]
+    w60 = np.zeros((T, K)); w60[:, es_i] = 0.6
+    if "ZN" in assets:
+        w60[:, assets.index("ZN")] = 0.4
+    print(f"\n===== {label}: {assets}  ({T} days) =====")
+    print(f"{'strategy':>24} {'Sharpe':>7} {'ann%':>6} {'maxDD%':>7} {'Calmar':>7}")
+    for nm, r in [("equity-only (ES)", bench_eq), ("60/40", book(w60)),
+                  ("risk-parity", rp), ("RP + vol-target", vt(rp)),
+                  ("RP + trend + vol-target", vt(rp_t))]:
+        m = metrics(r)
+        print(f"{nm:>24} {m['sr']:>+7.2f} {100*m['ann']:>+6.1f} {100*m['dd']:>+7.1f} {m['cal']:>7.2f}")
+    best = vt(rp)
+    print("  per-year Sharpe (RP+vol-target): " + " ".join(
+        f"{y}:{metrics(best[year[1:]==y])['sr']:+.1f}" for y in sorted(set(year)) if (year[1:] == y).sum() > 100))
+
+
+def main():
+    fut, btc = build()
+    run(["ES", "ZN", "GC", "CL"], {**fut}, "4-asset (2010-2026)")
+    run(["ES", "ZN", "GC", "CL", "BTC"], {**fut, "BTC": btc}, "5-asset +crypto (2019-2026)")
+    print("\nVERDICT: if RP/vol-target beats 60/40 on Sharpe AND maxDD = risk system adds real value on the")
+    print("right universe (the engine's true job). If 60/40 still wins, the deployable answer is just simple harvest.")
+
+
+if __name__ == "__main__":
+    main()