#!/usr/bin/env python3 """Harder gate: does the engine beat the climatology heuristic when given the REAL price drivers (weather -> renewables)? And specifically on the high-volatility days where the value concentrates? Adds free Open-Meteo weather (wind@100m, solar radiation, temp) for 3 German points to the day-ahead forecaster. Compares: climatology (the 83% champ) vs ML-price-only (76%) vs ML-with-weather. Capture % overall AND on the top-20% highest-spread days (where wind-drought spikes live and fundamentals should matter). Leak-free walk-forward. Weather@day-d is the forecast you'd have at decision time (day-ahead weather forecasts are ~90%+ accurate). """ import datetime import json import math import os import sys import time import urllib.request import numpy as np sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) from energy_battery_gate import fetch # cached prices # noqa: E402 from sklearn.ensemble import HistGradientBoostingRegressor # noqa: E402 ETA = 0.85; HRS = 4; ed = math.sqrt(ETA) WCACHE = "data/surfer/energy/weather" PTS = {"north": (53.55, 9.99), "central": (50.1, 8.68), "south": (48.14, 11.58)} def wget(u): for a in range(6): try: return json.load(urllib.request.urlopen(urllib.request.Request(u, headers={"User-Agent": "curl/8"}), timeout=90)) except urllib.error.HTTPError as e: if e.code == 429: time.sleep(10 * (a + 1)); continue raise raise RuntimeError("rate-limited") def weather(): os.makedirs(WCACHE, exist_ok=True) acc = {} for nm, (la, lo) in PTS.items(): f = f"{WCACHE}/{nm}.json" if os.path.exists(f): d = json.load(open(f)) else: u = (f"https://archive-api.open-meteo.com/v1/archive?latitude={la}&longitude={lo}" f"&start_date=2019-01-01&end_date=2025-09-29" f"&hourly=wind_speed_100m,shortwave_radiation,temperature_2m&timezone=UTC") d = wget(u)["hourly"]; json.dump(d, open(f, "w")); time.sleep(3) acc[nm] = d return acc 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(): ts, px = fetch() ok = np.isfinite(px); ts, px = ts[ok], px[ok] pday = np.array([datetime.datetime.utcfromtimestamp(int(t)).date().toordinal() for t in ts]) phour = np.array([datetime.datetime.utcfromtimestamp(int(t)).hour for t in ts]) # weather aligned to (ordinal_day, hour), averaged over 3 points W = weather() wind = {}; sol = {}; tmp = {} for nm, d in W.items(): for i, tstr in enumerate(d["time"]): dt = datetime.datetime.strptime(tstr, "%Y-%m-%dT%H:%M") k = (dt.date().toordinal(), dt.hour) wind.setdefault(k, []).append(d["wind_speed_100m"][i] or 0.0) sol.setdefault(k, []).append(d["shortwave_radiation"][i] or 0.0) tmp.setdefault(k, []).append(d["temperature_2m"][i] or 0.0) rows, dord = [], [] for dd in np.unique(pday): h = px[pday == dd] if len(h) == 24 and all((dd, hr) in wind for hr in range(24)): rows.append(h); dord.append(dd) P = np.array(rows); D = np.array(dord); ND = len(P) Wd = np.array([[np.mean(wind[(d, h)]) for h in range(24)] for d in D]) Sd = np.array([[np.mean(sol[(d, h)]) for h in range(24)] for d in D]) Td = np.array([[np.mean(tmp[(d, h)]) for h in range(24)] for d in D]) dow = np.array([datetime.date.fromordinal(int(x)).weekday() for x in D]) month = np.array([datetime.date.fromordinal(int(x)).month for x in D]) yr = np.array([datetime.date.fromordinal(int(x)).year for x in D]) fore = np.array([np.sort(P[d])[-HRS:].sum() * ed - np.sort(P[d])[:HRS].sum() / ed for d in range(ND)]) spread = P.max(1) - P.min(1) print(f"DE prices+weather aligned: {ND} days") roll7 = np.full_like(P, np.nan) for d in range(7, ND): roll7[d] = P[d - 7:d].mean(0) def feats(d, with_w): y1, y7 = P[d - 1], P[d - 7] X = np.zeros((24, 12 if with_w else 8)) for h in range(24): f = [h, dow[d], month[d], y1[h], y7[h], y1.mean(), y7.mean(), roll7[d, h]] if with_w: f += [Wd[d, h], Sd[d, h], Td[d, h], Wd[d].mean()] X[h] = f return X INIT, STEP = 365, 90 yhat = {"price_only": np.full_like(P, np.nan), "with_weather": np.full_like(P, np.nan)} for with_w, key in [(False, "price_only"), (True, "with_weather")]: Xc = {d: feats(d, with_w) 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([P[d] for d in range(7, s)]) gb = HistGradientBoostingRegressor(max_depth=4, max_iter=200, learning_rate=0.05, min_samples_leaf=50).fit(Xtr, ytr) for d in range(s, e): yhat[key][d] = gb.predict(Xc[d]) # climatology champ clim = np.full_like(P, 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] = P[past].mean(0) yhat["climatology"] = clim valid = np.arange(INIT, ND) hi = valid[spread[valid] >= np.quantile(spread[valid], 0.80)] # top-20% volatile days print(f"\n===== HARDER CAPTURE GATE (fundamentals + volatile-day focus), {len(valid)} days =====") print(f"perfect-foresight e{fore.mean()*365/1000:.1f}k/yr/MW | top-20%-vol days hold {100*fore[hi].sum()/fore[valid].sum():.0f}% of value") print(f"{'forecaster':>14} {'capture%':>8} {'EURk/yr':>8} {'capture% on HI-VOL days':>24}") for nm in ["climatology", "price_only", "with_weather"]: yh = yhat[nm] vv = np.array([dispatch(yh[d], P[d]) for d in valid]); cap = vv.sum() / fore[valid].sum() vh = np.array([dispatch(yh[d], P[d]) for d in hi]); caph = vh.sum() / fore[hi].sum() print(f"{nm:>14} {100*cap:>7.0f}% {vv.mean()*365/1000:>7.1f} {100*caph:>23.0f}%") print("\nVERDICT: with_weather > climatology (esp on HI-VOL days) = engine+fundamentals earns its keep.") print("If climatology still wins, even fundamentals don't beat the heuristic for day-ahead -> engine's home is elsewhere (intraday/multi-market).") if __name__ == "__main__": main()