foxhunt/scripts/surfer/energy_capture_gate2.py

#!/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()