"""Forecast skills v1: load / PV / price with quantile intervals (docs/05 §2.1). Baseline method, deliberately simple and fully deterministic: * point forecast (P50): exponentially-weighted mean of recent history days of the same day-type (weekday vs weekend), per 15-min interval; * interval (P10/P90): empirical quantiles of *backtest residuals* — the same method is re-run on each history day using only the days before it, and the pooled (actual − predicted) residuals give the band. This is split-conformal calibration in its simplest form, so coverage on held-out days tracks the nominal 80% instead of being an optimistic guess (docs/12 L2: 区间校准单独考核). Better models (statsforecast, gradient boosting) plug in behind the same contract; the eval harness (packages/evals) is what decides whether they win. Exogenous inputs are accepted by the contract but unused by v1 — recorded here so nobody mistakes v1 for weather-aware. """ from __future__ import annotations from datetime import date, datetime, timezone from typing import Callable import numpy as np from vpp_contracts.forecast_bundle import ForecastBundle from vpp_contracts.forecast_request import ForecastRequest from . import SKILL_VERSIONS from .numeric import INTERVALS, SCALE_MW, SCALE_PRICE, curve, to_array MODEL_NAME = {"LOAD": "load-forecast", "PV": "pv-forecast", "PRICE": "price-forecast"} UNIT_FOR_KIND = {"LOAD": "mw", "PV": "mw", "PRICE": "yuan_per_mwh"} LOOKBACK_DAYS = 14 HALF_LIFE_DAYS = 3.0 MIN_BACKTEST_HISTORY = 3 NOMINAL_LOWER, NOMINAL_UPPER = 0.10, 0.90 def _is_weekend(d: str) -> bool: return date.fromisoformat(d).weekday() >= 5 def _point_forecast(days: list[tuple[str, np.ndarray]], target: str) -> np.ndarray: """EWM over the most recent LOOKBACK_DAYS same-day-type days (fallback: all days).""" same = [(d, v) for d, v in days if _is_weekend(d) == _is_weekend(target)] pool = (same if len(same) >= 2 else days)[-LOOKBACK_DAYS:] ages = np.array([(date.fromisoformat(target) - date.fromisoformat(d)).days for d, _ in pool], float) w = np.power(0.5, ages / HALF_LIFE_DAYS) stack = np.stack([v for _, v in pool]) return (w[:, None] * stack).sum(axis=0) / w.sum() def _residual_quantiles(days: list[tuple[str, np.ndarray]], relative: bool) -> tuple[float, float]: """Pooled backtest residuals. Relative (actual/pred − 1) for quantities whose spread scales with level (load, PV — a night-time zero must not shrink the noon band); additive for prices, which can sit near or below zero.""" pairs: list[tuple[np.ndarray, np.ndarray]] = [] for i in range(MIN_BACKTEST_HISTORY, len(days)): d, actual = days[i] pairs.append((actual, _point_forecast(days[:i], d))) if not pairs: # Too little history to backtest: fall back to day-to-day differences. pairs = [(days[i][1], days[i - 1][1]) for i in range(1, len(days))] actual = np.concatenate([a for a, _ in pairs]) pred = np.concatenate([p for _, p in pairs]) if relative: floor = 0.05 * max(float(pred.max()), 1e-9) mask = pred > floor pooled = actual[mask] / pred[mask] - 1.0 if mask.any() else np.zeros(1) else: pooled = actual - pred lo, hi = np.quantile(pooled, [NOMINAL_LOWER, NOMINAL_UPPER]) return float(lo), float(hi) def forecast( req: ForecastRequest, clock: Callable[[], datetime] = lambda: datetime.now(timezone.utc), ) -> ForecastBundle: kind = str(req.kind.value) if req.unit.value != UNIT_FOR_KIND[kind]: raise ValueError(f"{kind} forecast unit must be {UNIT_FOR_KIND[kind]}, got {req.unit.value}") days = [(c.date, to_array(v.root for v in c.values)) for c in req.history] dates = [d for d, _ in days] if dates != sorted(dates) or len(set(dates)) != len(dates): raise ValueError("history must be ascending by date with no duplicates") if dates[-1] >= req.market_date: raise ValueError("history must be strictly before market_date") p50 = _point_forecast(days, req.market_date) relative = kind in ("LOAD", "PV") lo, hi = _residual_quantiles(days, relative) if relative: p10, p90 = p50 * (1.0 + lo), p50 * (1.0 + hi) else: p10, p90 = p50 + lo, p50 + hi if kind in ("LOAD", "PV"): p10, p50, p90 = (np.clip(x, 0.0, None) for x in (p10, p50, p90)) if kind == "PV": # Intervals that were zero on every history day (night) stay exactly zero. night = np.all(np.stack([v for _, v in days]) == 0.0, axis=0) for x in (p10, p50, p90): x[night] = 0.0 scale = SCALE_PRICE if kind == "PRICE" else SCALE_MW assert p50.shape == (INTERVALS,) bundle = { "id": f"fc-{kind.lower()}-{req.market_date}", "kind": kind, "market_date": req.market_date, "unit": req.unit.value, "quantiles": { "p10": curve(p10, req.market_date, scale), "p50": curve(p50, req.market_date, scale), "p90": curve(p90, req.market_date, scale), }, "model": {"name": MODEL_NAME[kind], "version": SKILL_VERSIONS[MODEL_NAME[kind]]}, "features_snapshot_ref": req.features_snapshot_ref, "generated_at": clock().isoformat().replace("+00:00", "Z"), } return ForecastBundle.model_validate(bundle)