244ccaaf0b
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>
99 lines
4.2 KiB
Python
99 lines
4.2 KiB
Python
#!/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()
|