diff --git a/CHANGELOG.md b/CHANGELOG.md index aaa5512..f7ee5dd 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/README.md b/README.md index 4909a68..cc95a71 100644 --- a/README.md +++ b/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. diff --git a/src/aep_load_forecasting/demo_data.py b/src/aep_load_forecasting/demo_data.py index 8c7e81f..58aa547 100644 --- a/src/aep_load_forecasting/demo_data.py +++ b/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() diff --git a/src/aep_load_forecasting/demo_pipeline.py b/src/aep_load_forecasting/demo_pipeline.py index 256a8f6..80d185b 100644 --- a/src/aep_load_forecasting/demo_pipeline.py +++ b/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:") diff --git a/src/aep_load_forecasting/forecast_24h.py b/src/aep_load_forecasting/forecast_24h.py index 546c92a..8f2dbad 100644 --- a/src/aep_load_forecasting/forecast_24h.py +++ b/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,6 +206,29 @@ 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, @@ -179,7 +236,15 @@ def main(argv: list[str] | None = None) -> int: 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}") diff --git a/src/aep_load_forecasting/forecasting.py b/src/aep_load_forecasting/forecasting.py index 22d499c..20924ba 100644 --- a/src/aep_load_forecasting/forecasting.py +++ b/src/aep_load_forecasting/forecasting.py @@ -3,6 +3,7 @@ from __future__ import annotations from collections.abc import Sequence +from math import ceil from typing import Protocol import numpy as np @@ -173,3 +174,102 @@ def recursive_forecast( }, index=future_index, ) + + +def calibrate_recursive_intervals( + model: Regressor, + history: pd.Series, + calibration_index: pd.DatetimeIndex, + feature_columns: Sequence[str] = FORECAST_FEATURES, + *, + horizon: int = 24, + coverage: float = 0.9, +) -> pd.Series: + """Calibrate lead-specific interval widths on recursive forecasts. + + The model must be trained only on observations before ``calibration_index``. + Each non-overlapping calibration block is forecast recursively from the + history available at its origin. Absolute errors are converted to the + finite-sample split-conformal order statistic independently for every lead. + """ + + if isinstance(horizon, bool) or not isinstance(horizon, int) or horizon < 1: + raise ValueError("Forecast horizon must be a positive integer.") + if isinstance(coverage, bool) or not 0.0 < coverage < 1.0: + raise ValueError("Interval coverage must be between zero and one.") + if not isinstance(calibration_index, pd.DatetimeIndex): + raise TypeError("Calibration timestamps must use a DatetimeIndex.") + if calibration_index.has_duplicates: + raise ValueError("Calibration timestamps contain duplicates.") + if not calibration_index.is_monotonic_increasing: + raise ValueError("Calibration timestamps must be sorted.") + if len(calibration_index) < horizon: + raise ValueError( + "Calibration data must contain at least one complete forecast " + f"horizon ({horizon} rows)." + ) + + steps = calibration_index.to_series().diff().dropna() + if not steps.eq(HOURLY_STEP).all(): + raise ValueError("Calibration timestamps must be consecutive hours.") + + columns = validate_feature_columns(feature_columns) + validated_history = _validated_history(history) + missing_timestamps = calibration_index.difference(validated_history.index) + if not missing_timestamps.empty: + raise ValueError("Calibration timestamps are missing from history.") + + complete_blocks = len(calibration_index) // horizon + errors_by_lead: list[list[float]] = [[] for _ in range(horizon)] + for block_number in range(complete_blocks): + start = block_number * horizon + block = calibration_index[start : start + horizon] + past = validated_history.loc[: block[0] - HOURLY_STEP] + forecast = recursive_forecast( + model, + past, + columns, + horizon=horizon, + ) + if not forecast.index.equals(block): + raise ValueError( + "Calibration timestamps must begin immediately after the " + "available history." + ) + + actual = validated_history.loc[block].to_numpy(dtype=float) + predicted = forecast["forecast_xgb_MW"].to_numpy(dtype=float) + for lead, error in enumerate(np.abs(actual - predicted)): + errors_by_lead[lead].append(float(error)) + + half_widths = [] + for errors in errors_by_lead: + rank = min(len(errors), ceil((len(errors) + 1) * coverage)) + half_widths.append(float(np.partition(errors, rank - 1)[rank - 1])) + + return pd.Series( + half_widths, + index=pd.RangeIndex(1, horizon + 1, name="horizon_hour"), + name="interval_half_width_MW", + ) + + +def add_prediction_intervals( + forecast: pd.DataFrame, + half_widths: Sequence[float], +) -> pd.DataFrame: + """Add physically bounded symmetric intervals to a point forecast.""" + + if "forecast_xgb_MW" not in forecast.columns: + raise ValueError("Forecast is missing the forecast_xgb_MW column.") + widths = np.asarray(half_widths, dtype=float).reshape(-1) + if widths.size != len(forecast): + raise ValueError("One interval width is required per forecast row.") + if not np.isfinite(widths).all() or (widths < 0).any(): + raise ValueError("Interval widths must be finite and non-negative.") + + result = forecast.copy() + point = result["forecast_xgb_MW"].to_numpy(dtype=float) + result["forecast_xgb_lower_MW"] = np.maximum(0.0, point - widths) + result["forecast_xgb_upper_MW"] = point + widths + return result diff --git a/tests/test_demo_data.py b/tests/test_demo_data.py index 70196b5..b2a8e74 100644 --- a/tests/test_demo_data.py +++ b/tests/test_demo_data.py @@ -8,6 +8,8 @@ from aep_load_forecasting.demo_data import ( BASELINE_COLUMN, FORECAST_COLUMN, + LOWER_COLUMN, + UPPER_COLUMN, ForecastDataError, forecast_plot_columns, load_forecast_csv, @@ -56,6 +58,24 @@ def test_load_forecast_csv_sorts_named_timestamps() -> None: assert forecast_plot_columns(forecast) == [FORECAST_COLUMN] +def test_load_forecast_csv_accepts_prediction_intervals() -> None: + source = csv_source( + """ +Datetime,forecast_xgb_MW,forecast_xgb_lower_MW,forecast_xgb_upper_MW +2025-01-01 00:00:00,120,115,126 +2025-01-01 01:00:00,121,116,127 +""" + ) + + forecast = load_forecast_csv(source) + + assert forecast_plot_columns(forecast) == [ + FORECAST_COLUMN, + LOWER_COLUMN, + UPPER_COLUMN, + ] + + @pytest.mark.parametrize( ("contents", "message"), [ @@ -103,6 +123,20 @@ def test_load_forecast_csv_sorts_named_timestamps() -> None: """, "'Datetime' column", ), + ( + """ +Datetime,forecast_xgb_MW,forecast_xgb_lower_MW +2025-01-01 00:00:00,120,115 +""", + "both prediction-interval", + ), + ( + """ +Datetime,forecast_xgb_MW,forecast_xgb_lower_MW,forecast_xgb_upper_MW +2025-01-01 00:00:00,120,125,130 +""", + "must contain the point forecast", + ), ], ) def test_load_forecast_csv_rejects_invalid_data( diff --git a/tests/test_demo_pipeline.py b/tests/test_demo_pipeline.py index c855a40..dc038f0 100644 --- a/tests/test_demo_pipeline.py +++ b/tests/test_demo_pipeline.py @@ -57,12 +57,25 @@ def test_command_writes_complete_reproducible_demo(tmp_path) -> None: model_artifact = joblib.load(expected[6]) assert tuple(model_artifact["features"]) == FORECAST_FEATURES + assert model_artifact["prediction_interval"]["coverage"] == 0.9 + assert model_artifact["prediction_interval"]["calibration_days"] == 2 + assert len(model_artifact["prediction_interval"]["half_widths_MW"]) == 3 forecast = pd.read_csv(expected[7]) assert len(forecast) == 3 assert { "forecast_xgb_MW", "baseline_blend_MW", + "forecast_xgb_lower_MW", + "forecast_xgb_upper_MW", }.issubset(forecast.columns) + assert ( + forecast["forecast_xgb_lower_MW"] + <= forecast["forecast_xgb_MW"] + ).all() + assert ( + forecast["forecast_xgb_MW"] + <= forecast["forecast_xgb_upper_MW"] + ).all() manifest = json.loads(expected[9].read_text(encoding="utf-8")) assert manifest["schema_version"] == 1 @@ -74,6 +87,7 @@ def test_command_writes_complete_reproducible_demo(tmp_path) -> None: "plot_days": 1, "horizon": 3, "n_estimators": 5, + "interval_coverage": 0.9, } assert manifest["runtime"]["python"] assert manifest["runtime"]["packages"]["xgboost"] @@ -104,3 +118,21 @@ def test_pipeline_rejects_too_little_history_before_writing(tmp_path) -> None: ) assert not output_dir.exists() + + +def test_pipeline_rejects_calibration_window_shorter_than_horizon( + tmp_path, +) -> None: + output_dir = tmp_path / "demo" + + with pytest.raises(ValueError, match="complete forecast horizon"): + run_demo_pipeline( + output_dir, + days=13, + evaluation_days=1, + plot_days=1, + horizon=25, + n_estimators=5, + ) + + assert not output_dir.exists() diff --git a/tests/test_forecast_24h.py b/tests/test_forecast_24h.py index b13196a..12ea4cf 100644 --- a/tests/test_forecast_24h.py +++ b/tests/test_forecast_24h.py @@ -14,9 +14,10 @@ def test_command_forwards_estimator_count_and_writes_outputs( ) -> None: features = pd.DataFrame( { - "y": [100.0], - **{name: [1.0] for name in FORECAST_FEATURES}, - } + "y": [100.0] * 25, + **{name: [1.0] * 25 for name in FORECAST_FEATURES}, + }, + index=pd.date_range("2025-12-01", periods=25, freq="h"), ) history = pd.Series( [100.0], @@ -29,8 +30,9 @@ def test_command_forwards_estimator_count_and_writes_outputs( }, index=pd.date_range("2026-01-01 01:00", periods=1, freq="h"), ) - model = object() - captured: dict[str, int] = {} + calibration_model = object() + final_model = object() + captured: dict[str, object] = {"training_rows": []} monkeypatch.setattr(command, "load_feature_table", lambda _: features) monkeypatch.setattr(command, "load_hourly_series", lambda _: history) @@ -40,21 +42,44 @@ def fake_train_final_model( *, n_estimators: int, ) -> tuple[object, tuple[str, ...]]: - assert received_features is features captured["n_estimators"] = n_estimators + training_rows = captured["training_rows"] + assert isinstance(training_rows, list) + training_rows.append(len(received_features)) + model = calibration_model if len(training_rows) == 1 else final_model return model, FORECAST_FEATURES monkeypatch.setattr(command, "train_final_model", fake_train_final_model) + monkeypatch.setattr( + command, + "calibrate_recursive_intervals", + lambda received_model, + received_history, + calibration_index, + columns, + *, + horizon, + coverage: pd.Series([2.0]), + ) monkeypatch.setattr( command, "recursive_forecast", - lambda received_model, received_history, columns, *, horizon: forecast, + lambda received_model, + received_history, + columns, + *, + horizon: forecast, ) def fake_dump(payload: object, path: str | Path) -> None: assert payload == { - "model": model, + "model": final_model, "features": list(FORECAST_FEATURES), + "prediction_interval": { + "coverage": 0.9, + "calibration_days": 1, + "half_widths_MW": [2.0], + }, } Path(path).write_bytes(b"model") @@ -68,7 +93,9 @@ def fake_save_plot( show: bool = False, ) -> Path: assert received_history is history - assert received_forecast is forecast + assert received_forecast["forecast_xgb_MW"].tolist() == [101.0] + assert received_forecast["forecast_xgb_lower_MW"].tolist() == [99.0] + assert received_forecast["forecast_xgb_upper_MW"].tolist() == [103.0] assert show is False path = Path(output_path) path.parent.mkdir(parents=True, exist_ok=True) @@ -96,11 +123,17 @@ def fake_save_plot( "1", "--estimators", "17", + "--calibration-days", + "1", ] ) assert result == 0 assert captured["n_estimators"] == 17 + assert captured["training_rows"] == [1, 25] assert model_path.is_file() assert forecast_path.is_file() assert figure_path.is_file() + saved_forecast = pd.read_csv(forecast_path) + assert saved_forecast["forecast_xgb_lower_MW"].tolist() == [99.0] + assert saved_forecast["forecast_xgb_upper_MW"].tolist() == [103.0] diff --git a/tests/test_forecasting.py b/tests/test_forecasting.py index 90cf1bc..221d4e9 100644 --- a/tests/test_forecasting.py +++ b/tests/test_forecasting.py @@ -6,6 +6,8 @@ from aep_load_forecasting.forecasting import ( FORECAST_FEATURES, + add_prediction_intervals, + calibrate_recursive_intervals, make_forecast_row, recursive_forecast, ) @@ -29,6 +31,11 @@ def predict(self, features: pd.DataFrame) -> np.ndarray: return features["lag_1"].to_numpy() + 1.0 +class PersistenceModel: + def predict(self, features: pd.DataFrame) -> np.ndarray: + return features["lag_1"].to_numpy() + + def test_make_forecast_row_uses_only_past_values() -> None: history = hourly_history() timestamp = history.index[-1] + pd.Timedelta(hours=1) @@ -101,3 +108,58 @@ def predict(self, features: pd.DataFrame) -> np.ndarray: def test_recursive_forecast_rejects_non_finite_prediction() -> None: with pytest.raises(ValueError, match="non-finite"): recursive_forecast(NonFiniteModel(), hourly_history()) + + +def test_recursive_interval_calibration_is_lead_specific() -> None: + history = hourly_history() + calibration_index = history.index[-6:] + + widths = calibrate_recursive_intervals( + PersistenceModel(), + history, + calibration_index, + horizon=3, + coverage=0.5, + ) + + assert widths.index.tolist() == [1, 2, 3] + assert widths.tolist() == [1.0, 2.0, 3.0] + + +@pytest.mark.parametrize("coverage", [0.0, 1.0, -0.1, 1.1, True]) +def test_recursive_interval_calibration_rejects_invalid_coverage( + coverage: object, +) -> None: + history = hourly_history() + + with pytest.raises(ValueError, match="between zero and one"): + calibrate_recursive_intervals( + PersistenceModel(), + history, + history.index[-3:], + horizon=3, + coverage=coverage, # type: ignore[arg-type] + ) + + +def test_add_prediction_intervals_clips_negative_load_bound() -> None: + forecast = pd.DataFrame( + {"forecast_xgb_MW": [2.0, 5.0]}, + index=pd.date_range("2026-01-01", periods=2, freq="h"), + ) + + result = add_prediction_intervals(forecast, [3.0, 1.0]) + + assert result["forecast_xgb_lower_MW"].tolist() == [0.0, 4.0] + assert result["forecast_xgb_upper_MW"].tolist() == [5.0, 6.0] + assert "forecast_xgb_lower_MW" not in forecast.columns + + +@pytest.mark.parametrize("widths", [[1.0], [1.0, -1.0], [1.0, np.nan]]) +def test_add_prediction_intervals_rejects_invalid_widths( + widths: list[float], +) -> None: + forecast = pd.DataFrame({"forecast_xgb_MW": [2.0, 5.0]}) + + with pytest.raises(ValueError, match="interval width|finite"): + add_prediction_intervals(forecast, widths)