From 244ccaaf0b9df8e559bfa3ea7a5f139119171442 Mon Sep 17 00:00:00 2001 From: jgrusewski Date: Sun, 7 Jun 2026 00:42:32 +0200 Subject: [PATCH] =?UTF-8?q?feat(harvest):=20premium-harvesting=20tested=20?= =?UTF-8?q?=E2=80=94=20simple=2060/40=20beats=20the=20sophisticated=20syst?= =?UTF-8?q?em?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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) --- scripts/surfer/energy_intraday_gate.py | 98 ++++++++++++++++++++++++++ scripts/surfer/fetch_caiso.py | 83 ++++++++++++++++++++++ scripts/surfer/portfolio_harvest.py | 90 +++++++++++++++++++++++ scripts/surfer/portfolio_harvest2.py | 98 ++++++++++++++++++++++++++ 4 files changed, 369 insertions(+) create mode 100644 scripts/surfer/energy_intraday_gate.py create mode 100644 scripts/surfer/fetch_caiso.py create mode 100644 scripts/surfer/portfolio_harvest.py create mode 100644 scripts/surfer/portfolio_harvest2.py diff --git a/scripts/surfer/energy_intraday_gate.py b/scripts/surfer/energy_intraday_gate.py new file mode 100644 index 000000000..16297ac11 --- /dev/null +++ b/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() diff --git a/scripts/surfer/fetch_caiso.py b/scripts/surfer/fetch_caiso.py new file mode 100644 index 000000000..2f02b9b7d --- /dev/null +++ b/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() diff --git a/scripts/surfer/portfolio_harvest.py b/scripts/surfer/portfolio_harvest.py new file mode 100644 index 000000000..bfcef4458 --- /dev/null +++ b/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() diff --git a/scripts/surfer/portfolio_harvest2.py b/scripts/surfer/portfolio_harvest2.py new file mode 100644 index 000000000..b5b0de2c5 --- /dev/null +++ b/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()