[codex] Add calibrated forecast intervals

CHANGELOG.md

@@ -2,6 +2,13 @@
 
 All notable changes to this project are documented in this file.
 
+## Unreleased
+
+### Added
+
+- Add horizon-specific prediction intervals calibrated on leakage-free,
+  recursive rolling-origin forecast errors.
+
 ## [0.1.1] - 2026-07-28
 
 ### Fixed

README.md

@@ -189,10 +189,19 @@ can also be configured from the command line:
 aep-forecast `
   --horizon 48 `
   --estimators 400 `
+  --interval-coverage 0.9 `
+  --calibration-days 30 `
   --output reports\forecast_next48h.csv `
   --figure reports\figures\forecast_next48h.png
 ```
 
+Forecast CSVs include the point estimate plus
+`forecast_xgb_lower_MW` and `forecast_xgb_upper_MW`. The interval width is
+calibrated separately for every forecast hour using recursive rolling-origin
+errors from the trailing calibration window. The calibration model sees only
+earlier observations; after calibration, the exported production model is
+refitted on all available feature rows.
+
 ## Demo (Streamlit)
 
 ```powershell
@@ -201,8 +210,9 @@ streamlit run streamlit_app.py
 
 The app loads the bundled `assets/forecast_next24h.csv` by default, or lets you
 upload your own forecast CSV. Uploads must contain hourly `Datetime` values and
-a numeric `forecast_xgb_MW` column; `baseline_blend_MW` is optional. Invalid
-files produce a clear error in the app instead of a chart or metric failure.
+a numeric `forecast_xgb_MW` column; `baseline_blend_MW` and the paired
+prediction-interval columns are optional. Invalid files produce a clear error
+in the app instead of a chart or metric failure.
 
 ## Data validation
 
@@ -239,7 +249,9 @@ Python 3.10 and 3.12 for every pull request.
 
 ## Limitations
 
-- Single region (AEP) and a **point** forecast — no uncertainty intervals.
+- Single region (AEP); interval coverage is empirical and can drift when the
+  load process changes. Time-series dependence means the calibrated intervals
+  are not a formal distribution-free coverage guarantee.
 - The 24-hour forecast is **recursive**, so errors can compound over the horizon.
 - Features are calendar + lags only — no weather or holiday signals.
 - Reported metrics come from a single 30-day out-of-time test window.

src/aep_load_forecasting/demo_data.py

@@ -13,7 +13,14 @@
 DATETIME_COLUMN = "Datetime"
 FORECAST_COLUMN = "forecast_xgb_MW"
 BASELINE_COLUMN = "baseline_blend_MW"
-PLOT_COLUMNS = (FORECAST_COLUMN, BASELINE_COLUMN)
+LOWER_COLUMN = "forecast_xgb_lower_MW"
+UPPER_COLUMN = "forecast_xgb_upper_MW"
+PLOT_COLUMNS = (
+    FORECAST_COLUMN,
+    LOWER_COLUMN,
+    UPPER_COLUMN,
+    BASELINE_COLUMN,
+)
 HOURLY_STEP = pd.Timedelta(hours=1)
 
 
@@ -60,6 +67,12 @@ def load_forecast_csv(source: CsvSource) -> pd.DataFrame:
         raise ForecastDataError(
             f"Forecast CSV is missing required column {FORECAST_COLUMN!r}."
         )
+    interval_columns = {LOWER_COLUMN, UPPER_COLUMN}
+    present_interval_columns = interval_columns.intersection(frame.columns)
+    if present_interval_columns and present_interval_columns != interval_columns:
+        raise ForecastDataError(
+            "Forecast CSV must contain both prediction-interval columns."
+        )
 
     timestamps = pd.to_datetime(frame[timestamp_column], errors="coerce")
     invalid_timestamps = int(timestamps.isna().sum())
@@ -90,6 +103,16 @@ def load_forecast_csv(source: CsvSource) -> pd.DataFrame:
             )
         frame[column] = values.astype(float)
 
+    if present_interval_columns:
+        outside_interval = (
+            (frame[LOWER_COLUMN] > frame[FORECAST_COLUMN])
+            | (frame[FORECAST_COLUMN] > frame[UPPER_COLUMN])
+        )
+        if outside_interval.any():
+            raise ForecastDataError(
+                "Prediction intervals must contain the point forecast."
+            )
+
     frame = frame.drop(columns=timestamp_column)
     frame.index = pd.DatetimeIndex(timestamps, name=DATETIME_COLUMN)
     frame = frame.sort_index()

src/aep_load_forecasting/demo_pipeline.py

@@ -22,10 +22,15 @@
 from aep_load_forecasting.evaluation import HOURS_PER_DAY
 from aep_load_forecasting.forecast_24h import (
     DEFAULT_ESTIMATORS,
+    DEFAULT_INTERVAL_COVERAGE,
     save_forecast_plot,
     train_final_model,
 )
-from aep_load_forecasting.forecasting import recursive_forecast
+from aep_load_forecasting.forecasting import (
+    add_prediction_intervals,
+    calibrate_recursive_intervals,
+    recursive_forecast,
+)
 from aep_load_forecasting.make_features import (
     load_hourly_series,
     make_feature_table,
@@ -181,6 +186,7 @@ def _validate_run_settings(
     plot_days: int,
     horizon: int,
     n_estimators: int,
+    interval_coverage: float,
 ) -> None:
     settings = {
         "days": days,
@@ -194,6 +200,13 @@ def _validate_run_settings(
             raise ValueError(f"{name} must be a positive integer.")
     if plot_days > evaluation_days:
         raise ValueError("plot_days must not exceed evaluation_days.")
+    if evaluation_days * HOURS_PER_DAY < horizon:
+        raise ValueError(
+            "The evaluation window must contain at least one complete "
+            "forecast horizon for interval calibration."
+        )
+    if isinstance(interval_coverage, bool) or not 0.0 < interval_coverage < 1.0:
+        raise ValueError("interval_coverage must be between zero and one.")
 
     required_hours = (
         LAG_HISTORY_HOURS
@@ -218,6 +231,7 @@ def run_demo_pipeline(
     plot_days: int = DEFAULT_PLOT_DAYS,
     horizon: int = DEFAULT_HORIZON,
     n_estimators: int = DEFAULT_ESTIMATORS,
+    interval_coverage: float = DEFAULT_INTERVAL_COVERAGE,
 ) -> DemoArtifacts:
     """Generate data, evaluate models, and export a future forecast."""
 
@@ -227,6 +241,7 @@ def run_demo_pipeline(
         plot_days=plot_days,
         horizon=horizon,
         n_estimators=n_estimators,
+        interval_coverage=interval_coverage,
     )
     artifacts = demo_artifacts(output_dir)
     evaluation_hours = evaluation_days * HOURS_PER_DAY
@@ -263,13 +278,35 @@ def run_demo_pipeline(
     )
 
     history = load_hourly_series(artifacts.source)
+    calibration_features = features.iloc[-evaluation_hours:]
+    calibration_model, feature_columns = train_final_model(
+        features.iloc[:-evaluation_hours],
+        n_estimators=n_estimators,
+    )
+    interval_half_widths = calibrate_recursive_intervals(
+        calibration_model,
+        history,
+        calibration_features.index,
+        feature_columns,
+        horizon=horizon,
+        coverage=interval_coverage,
+    )
+
     model, feature_columns = train_final_model(
         features,
         n_estimators=n_estimators,
     )
     artifacts.model.parent.mkdir(parents=True, exist_ok=True)
     joblib.dump(
-        {"model": model, "features": list(feature_columns)},
+        {
+            "model": model,
+            "features": list(feature_columns),
+            "prediction_interval": {
+                "coverage": interval_coverage,
+                "calibration_days": evaluation_days,
+                "half_widths_MW": interval_half_widths.tolist(),
+            },
+        },
         artifacts.model,
     )
 
@@ -279,6 +316,7 @@ def run_demo_pipeline(
         feature_columns,
         horizon=horizon,
     )
+    forecast = add_prediction_intervals(forecast, interval_half_widths)
     artifacts.forecast.parent.mkdir(parents=True, exist_ok=True)
     forecast.to_csv(artifacts.forecast)
     save_forecast_plot(
@@ -297,6 +335,7 @@ def run_demo_pipeline(
             "plot_days": plot_days,
             "horizon": horizon,
             "n_estimators": n_estimators,
+            "interval_coverage": interval_coverage,
         },
     )
     return artifacts
@@ -324,6 +363,11 @@ def parse_args(argv: list[str] | None = None) -> argparse.Namespace:
     parser.add_argument("--plot-days", type=int, default=DEFAULT_PLOT_DAYS)
     parser.add_argument("--horizon", type=int, default=DEFAULT_HORIZON)
     parser.add_argument("--estimators", type=int, default=DEFAULT_ESTIMATORS)
+    parser.add_argument(
+        "--interval-coverage",
+        type=float,
+        default=DEFAULT_INTERVAL_COVERAGE,
+    )
     return parser.parse_args(argv)
 
 
@@ -338,6 +382,7 @@ def main(argv: list[str] | None = None) -> int:
         plot_days=args.plot_days,
         horizon=args.horizon,
         n_estimators=args.estimators,
+        interval_coverage=args.interval_coverage,
     )
 
     print("Demo pipeline completed:")

src/aep_load_forecasting/forecast_24h.py

@@ -11,8 +11,11 @@
 from xgboost import XGBRegressor
 
 from aep_load_forecasting.cli import add_version_argument
+from aep_load_forecasting.evaluation import HOURS_PER_DAY
 from aep_load_forecasting.forecasting import (
     FORECAST_FEATURES,
+    add_prediction_intervals,
+    calibrate_recursive_intervals,
     recursive_forecast,
     validate_feature_columns,
 )
@@ -26,6 +29,8 @@
 DEFAULT_FIGURE = Path("reports") / "figures" / "forecast_next24h.png"
 DEFAULT_MODEL = Path("models") / "aep_xgb.joblib"
 DEFAULT_ESTIMATORS = 800
+DEFAULT_INTERVAL_COVERAGE = 0.9
+DEFAULT_CALIBRATION_DAYS = 30
 
 
 def load_feature_table(csv_path: str | Path) -> pd.DataFrame:
@@ -124,6 +129,17 @@ def save_forecast_plot(
         forecast["forecast_xgb_MW"],
         label="XGBoost forecast",
     )
+    if {
+        "forecast_xgb_lower_MW",
+        "forecast_xgb_upper_MW",
+    }.issubset(forecast.columns):
+        axis.fill_between(
+            forecast.index,
+            forecast["forecast_xgb_lower_MW"],
+            forecast["forecast_xgb_upper_MW"],
+            alpha=0.2,
+            label="Calibrated interval",
+        )
     axis.set_title(f"AEP Load Forecast: next {len(forecast)} hours")
     axis.set_xlabel("Time")
     axis.set_ylabel("MW")
@@ -160,6 +176,24 @@ def parse_args(argv: list[str] | None = None) -> argparse.Namespace:
         default=DEFAULT_ESTIMATORS,
         help=f"Number of boosting rounds (default: {DEFAULT_ESTIMATORS})",
     )
+    parser.add_argument(
+        "--interval-coverage",
+        type=float,
+        default=DEFAULT_INTERVAL_COVERAGE,
+        help=(
+            "Target coverage for calibrated prediction intervals "
+            f"(default: {DEFAULT_INTERVAL_COVERAGE})"
+        ),
+    )
+    parser.add_argument(
+        "--calibration-days",
+        type=int,
+        default=DEFAULT_CALIBRATION_DAYS,
+        help=(
+            "Trailing days reserved for interval calibration "
+            f"(default: {DEFAULT_CALIBRATION_DAYS})"
+        ),
+    )
     parser.add_argument(
         "--show",
         action="store_true",
@@ -172,14 +206,45 @@ def main(argv: list[str] | None = None) -> int:
     args = parse_args(argv)
     features = load_feature_table(args.features)
     history = load_hourly_series(args.input)
+    if args.calibration_days <= 0:
+        raise ValueError("calibration-days must be greater than zero.")
+    calibration_hours = args.calibration_days * HOURS_PER_DAY
+    if len(features) <= calibration_hours:
+        raise ValueError(
+            "Feature data needs at least one training row in addition to "
+            f"{calibration_hours} calibration rows."
+        )
+
+    calibration_features = features.iloc[-calibration_hours:]
+    calibration_model, feature_columns = train_final_model(
+        features.iloc[:-calibration_hours],
+        n_estimators=args.estimators,
+    )
+    interval_half_widths = calibrate_recursive_intervals(
+        calibration_model,
+        history,
+        calibration_features.index,
+        feature_columns,
+        horizon=args.horizon,
+        coverage=args.interval_coverage,
+    )
+
     model, feature_columns = train_final_model(
         features,
         n_estimators=args.estimators,
     )
 
     args.model.parent.mkdir(parents=True, exist_ok=True)
     joblib.dump(
-        {"model": model, "features": list(feature_columns)},
+        {
+            "model": model,
+            "features": list(feature_columns),
+            "prediction_interval": {
+                "coverage": args.interval_coverage,
+                "calibration_days": args.calibration_days,
+                "half_widths_MW": interval_half_widths.tolist(),
+            },
+        },
         args.model,
     )
     print(f"Model saved: {args.model}")
@@ -190,6 +255,10 @@ def main(argv: list[str] | None = None) -> int:
         feature_columns,
         horizon=args.horizon,
     )
+    forecast = add_prediction_intervals(
+        forecast,
+        interval_half_widths,
+    )
     args.output.parent.mkdir(parents=True, exist_ok=True)
     forecast.to_csv(args.output)
     print(f"Forecast saved: {args.output}")