feat(conformal): CV-calibrated prediction intervals

docs/api/core/conformal.md

@@ -0,0 +1,3 @@
+# `foundationforecast.core.conformal`
+
+::: foundationforecast.core.conformal.ConformalIntervals

docs/conformal-prediction.md

@@ -0,0 +1,89 @@
+# Conformal prediction intervals
+
+FoundationForecast can calibrate **distribution-free prediction intervals** from cross-validation absolute errors on the **point forecast**, instead of using a foundation model’s native quantile head.
+
+## When to use it
+
+- You need **wide levels** (e.g. 95%) from a model whose quantile head stops at fixed knots. Tafsut and Chronos-Bolt predict quantiles 0.1 to 0.9 only, so their native 95% interval is clamped to the 80% one; conformal intervals are not bound by those knots.
+- You want intervals from **empirical forecast errors** on your own history (per series), rather than the model’s built-in uncertainty.
+- You accept the extra compute: about **one point forecast per calibration window per series**, plus the final forecast.
+
+Native quantiles (default when `prediction_intervals` is unset) remain unchanged. **`cross_validation()`** uses conformal intervals too: every window is recalibrated on the data before its cutoff, so backtests show the coverage you would get in production. **`detect_anomalies`** is unchanged and does not read `prediction_intervals`.
+
+## API
+
+Pass [`ConformalIntervals`][foundationforecast.core.conformal.ConformalIntervals] as `prediction_intervals` when creating the model, then pass `level` or `quantiles` to [`forecast()`][foundationforecast.FoundationForecast.forecast]:
+
+```python
+import pandas as pd
+
+from foundationforecast import ConformalIntervals, FoundationForecast
+from foundationforecast.models import Chronos
+
+df = pd.read_csv(
+    "https://timecopilot.s3.amazonaws.com/public/data/air_passengers.csv",
+    parse_dates=["ds"],
+)
+
+model = Chronos(
+    repo_id="amazon/chronos-bolt-tiny",
+    alias="Chronos-Bolt",
+    prediction_intervals=ConformalIntervals(n_windows=3),
+)
+
+ff = FoundationForecast(models=[model])
+fcst = ff.forecast(df=df, h=12, freq="MS", level=[80, 95])
+```
+
+### Parameters
+
+| Field | Meaning |
+| --- | --- |
+| `n_windows` | Rolling backtest windows for residual estimation (minimum **2**). |
+| `h` | Set automatically from `forecast(h=...)` on each call. |
+| `method` | `conformal_distribution` (default) or `conformal_error`. |
+
+### Requirements
+
+- You must pass **`level`** or **`quantiles`** whenever `prediction_intervals` is set.
+- Each series needs at least **`n_windows * h + 1`** observations.
+
+## How it works
+
+Calibration runs separately for every `unique_id`. Let \(T\) be the last observed timestamp of a series, \(h\) the horizon and \(W\) = `n_windows`.
+
+**1. Backtest the point forecast.** For each window \(w = 1, \dots, W\), the model forecasts \(h\) steps from the cutoff \(t_w = T - (W - w + 1)\,h\), using only data up to \(t_w\). The absolute error at horizon step \(j = 1, \dots, h\) is
+
+\[
+e_{w,j} = \left| \, y_{t_w + j} - \hat{y}_{t_w + j \mid t_w} \right|.
+\]
+
+This gives \(W\) calibration errors per horizon step, so later steps, which are usually harder to forecast, get their own (typically larger) errors.
+
+**2. Forecast.** The model produces the point forecast \(\hat{y}_{T+j}\) from the full history, exactly as without conformal intervals.
+
+**3. Build the interval.** For a level \(L\) (in percent), let \(\alpha = 1 - L/100\).
+
+- **`conformal_error`** widens the point forecast symmetrically by the \(L\)% quantile of the calibration errors:
+
+    \[
+    \hat{q}_j = \operatorname{Quantile}_{1-\alpha}\big(e_{1,j}, \dots, e_{W,j}\big),
+    \qquad
+    \big[\, \hat{y}_{T+j} - \hat{q}_j,\; \hat{y}_{T+j} + \hat{q}_j \,\big].
+    \]
+
+- **`conformal_distribution`** (default) forms \(2W\) synthetic outcomes around the point forecast and takes their symmetric quantiles:
+
+    \[
+    S_j = \big\{ \hat{y}_{T+j} - e_{w,j} \big\}_{w=1}^{W} \cup \big\{ \hat{y}_{T+j} + e_{w,j} \big\}_{w=1}^{W},
+    \]
+
+    \[
+    \big[\, \operatorname{Quantile}_{\alpha/2}(S_j),\; \operatorname{Quantile}_{1-\alpha/2}(S_j) \,\big].
+    \]
+
+The point forecast \(\hat{y}_{T+j}\) is never modified; only the interval columns change. When `quantiles` are requested, each quantile \(q \neq 0.5\) is mapped to the level \(L = |100 - 200\,q|\) and read from the matching lower or upper bound.
+
+**Coverage.** If recent forecast errors are representative of the next \(h\) steps, about \(L\)% of actuals fall inside the \(L\)% interval. With few windows the empirical quantiles are coarse (for example, with \(W = 2\) the 95% bound is essentially the largest observed error), so use more windows when you need high levels.
+
+See the [conformal prediction example notebook](examples/conformal-prediction.ipynb) on electricity prices.

docs/examples/index.md

@@ -12,6 +12,7 @@ For model API details, see the [Model Hub](../model-hub.md).
 | [Anomaly Detection](anomaly-detection-forecaster-quickstart.ipynb) | Detect anomalies with Chronos and FlowState | Python 3.11+ for FlowState |
 | [Compare Foundation Models](ts-foundation-models-comparison-quickstart.ipynb) | Benchmark multiple foundation models side by side | Python 3.10+; GPU optional |
 | [Exogenous variables](exogenous-variables.ipynb) | Panel forecast with `X_df`, native exog, and `exog_strategy=False` for mixed models | Python 3.10+; GPU recommended for Chronos-2 |
+| [Conformal prediction](conformal-prediction.ipynb) | CV-calibrated intervals with `prediction_intervals` on electricity prices | Python 3.10+; GPU optional |
 
 ## Foundation Models
 

docs/forecasting-parameters.md

@@ -40,6 +40,10 @@ ff = FoundationForecast(models=[Chronos()])
 fcst_df = ff.forecast(df=df, h=12, freq="MS")
 ```
 
+## Conformal prediction intervals
+
+When a model is constructed with [`ConformalIntervals`][foundationforecast.core.conformal.ConformalIntervals] (`prediction_intervals=...`), `forecast()` calibrates `level` / `quantiles` from cross-validation absolute errors on the point forecast instead of native model quantiles. See [Conformal prediction](conformal-prediction.md).
+
 ## Anomaly detection defaults
 
 When calling `detect_anomalies()` without `h`:

docs/javascripts/mathjax.js

@@ -0,0 +1,19 @@
+window.MathJax = {
+  tex: {
+    inlineMath: [["\\(", "\\)"]],
+    displayMath: [["\\[", "\\]"]],
+    processEscapes: true,
+    processEnvironments: true,
+  },
+  options: {
+    ignoreHtmlClass: ".*|",
+    processHtmlClass: "arithmatex",
+  },
+};
+
+document$.subscribe(() => {
+  MathJax.startup.output.clearCache();
+  MathJax.typesetClear();
+  MathJax.texReset();
+  MathJax.typesetPromise();
+});

foundationforecast/__init__.py

@@ -1,6 +1,7 @@
 import sys
 
 from ._foundation_forecast import FoundationForecast
+from .core.conformal import ConformalIntervals
 from .core.model_weight_cache import set_max_cached_models
 from .models import (
     Chronos,
@@ -15,6 +16,7 @@
 )
 
 __all__ = [
+    "ConformalIntervals",
     "FoundationForecast",
     "set_max_cached_models",
     "Chronos",

foundationforecast/core/conformal.py

@@ -0,0 +1,213 @@
+from __future__ import annotations
+
+from collections.abc import Callable
+from dataclasses import dataclass
+from typing import TYPE_CHECKING
+
+import numpy as np
+import pandas as pd
+import utilsforecast.processing as ufp
+
+if TYPE_CHECKING:
+    from .forecaster import ExogCapableForecaster
+
+
+@dataclass
+class ConformalIntervals:
+    """Cross-validation settings for distribution-free prediction intervals.
+
+    Args:
+        n_windows: Number of rolling backtest windows used to estimate residuals.
+            Must be at least 2.
+        h: Forecast horizon used when estimating residuals. Set automatically
+            from ``forecast(h=...)`` on each call.
+        method: ``conformal_distribution`` builds bounds from mean ± error
+            paths; ``conformal_error`` uses quantiles of absolute errors per
+            horizon step.
+    """
+
+    n_windows: int = 2
+    h: int = 1
+    method: str = "conformal_distribution"
+
+    def __post_init__(self) -> None:
+        if self.n_windows < 2:
+            raise ValueError(
+                "You need at least two windows to compute conformal intervals"
+            )
+        allowed = ("conformal_distribution", "conformal_error")
+        if self.method not in allowed:
+            raise ValueError(f"method must be one of {allowed}")
+
+
+def validate_conformal_samples(n_samples: int, h: int, n_windows: int) -> None:
+    min_samples = 2 * h + 1
+    if n_samples < min_samples:
+        raise ValueError(
+            f"Prediction intervals require at least {min_samples:,} samples; "
+            f"series has {n_samples:,}."
+        )
+    required = n_windows * h + 1
+    if n_samples < required:
+        raise ValueError(
+            f"Prediction intervals require at least {required:,} samples for "
+            f"n_windows={n_windows} and h={h}; series has {n_samples:,}."
+        )
+
+
+def effective_n_windows(n_samples: int, h: int, n_windows: int) -> int:
+    capped = min(n_windows, (n_samples - 1) // h)
+    if capped < 2:
+        raise ValueError(
+            f"Prediction intervals require at least {2 * h + 1:,} samples; "
+            f"series has {n_samples:,}."
+        )
+    return capped
+
+
+def add_conformal_distribution_intervals(
+    mean: np.ndarray,
+    cs: np.ndarray,
+    level: list[int | float],
+) -> dict[str, np.ndarray]:
+    level_sorted = sorted(level)
+    alphas = [100 - lv for lv in level_sorted]
+    cuts = [alpha / 200 for alpha in reversed(alphas)]
+    cuts.extend(1 - alpha / 200 for alpha in alphas)
+    mean_row = mean.reshape(1, -1)
+    scores = np.vstack([mean_row - cs, mean_row + cs])
+    quantiles = np.quantile(scores, cuts, axis=0)
+    quantiles = quantiles.reshape(len(cuts), -1)
+    lo_cols = [f"lo-{lv}" for lv in reversed(level_sorted)]
+    hi_cols = [f"hi-{lv}" for lv in level_sorted]
+    out: dict[str, np.ndarray] = {}
+    for i, col in enumerate(lo_cols + hi_cols):
+        out[col] = quantiles[i]
+    return out
+
+
+def add_conformal_error_intervals(
+    mean: np.ndarray,
+    cs: np.ndarray,
+    level: list[int | float],
+) -> dict[str, np.ndarray]:
+    margin = {lv: np.quantile(cs, lv / 100, axis=0) for lv in level}
+    out: dict[str, np.ndarray] = {}
+    for lv in reversed(sorted(level)):
+        out[f"lo-{lv}"] = mean - margin[lv]
+    for lv in sorted(level):
+        out[f"hi-{lv}"] = mean + margin[lv]
+    return out
+
+
+def _apply_conformal_method(
+    mean: np.ndarray,
+    cs: np.ndarray,
+    level: list[int | float],
+    method: str,
+) -> dict[str, np.ndarray]:
+    if method == "conformal_distribution":
+        return add_conformal_distribution_intervals(mean, cs, level)
+    if method == "conformal_error":
+        return add_conformal_error_intervals(mean, cs, level)
+    raise ValueError(f"Unsupported conformal method: {method}")
+
+
+def add_conformal_columns_to_forecast(
+    fcst_df: pd.DataFrame,
+    alias: str,
+    scores_by_id: dict[object, np.ndarray],
+    level: list[int | float],
+    method: str,
+) -> pd.DataFrame:
+    fcst_df = ufp.copy_if_pandas(fcst_df, deep=False)
+    for uid, cs in scores_by_id.items():
+        mask = fcst_df["unique_id"] == uid
+        if not mask.any():
+            continue
+        sub = fcst_df.loc[mask].sort_values("ds")
+        mean = sub[alias].to_numpy()
+        bands = _apply_conformal_method(mean, cs, level, method)
+        for col_suffix, values in bands.items():
+            side, lv = col_suffix.split("-", 1)
+            out_col = f"{alias}-{side}-{lv}"
+            fcst_df.loc[sub.index, out_col] = values
+    return fcst_df
+
+
+def conformity_scores_abs_errors(
+    y: np.ndarray,
+    h: int,
+    n_windows: int,
+    forecast_values: Callable[[int, int], np.ndarray],
+) -> np.ndarray:
+    n_samples = y.size
+    n_windows_eff = effective_n_windows(n_samples, h, n_windows)
+    test_size = n_windows_eff * h
+    cs = np.empty((n_windows_eff, h), dtype=float)
+    for i_window in range(n_windows_eff):
+        train_end = n_samples - test_size + i_window * h
+        y_test = y[train_end : train_end + h]
+        fcst = forecast_values(train_end, i_window)
+        cs[i_window] = np.abs(fcst - y_test)
+    return cs
+
+
+def compute_panel_conformity_scores(
+    forecaster: ExogCapableForecaster,
+    df: pd.DataFrame,
+    h: int,
+    freq: str,
+    prediction_intervals: ConformalIntervals,
+    *,
+    exog_cols: list[str],
+) -> dict[object, np.ndarray]:
+    """Absolute CV errors of the point forecast, keyed by ``unique_id``.
+
+    Each window is forecast from that series' training slice only, so callers
+    must not forward a precomputed panel or the forecast horizon ``X_df``.
+    """
+    from .exog.covariates import exog_strategy_disabled, sort_exog_panel
+
+    scores: dict[object, np.ndarray] = {}
+    work_df = df
+    if exog_cols and not exog_strategy_disabled(forecaster.exog_strategy):
+        work_df = sort_exog_panel(df)
+
+    for uid, group in work_df.groupby("unique_id", sort=False):
+        series = group.sort_values("ds")
+        y = series["y"].to_numpy(dtype=float)
+        validate_conformal_samples(len(y), h, prediction_intervals.n_windows)
+
+        def forecast_values(
+            train_end: int,
+            _i_window: int,
+            _series: pd.DataFrame = series,
+        ) -> np.ndarray:
+            train = _series.iloc[:train_end]
+            valid = _series.iloc[train_end : train_end + h]
+            horizon_x = valid[["unique_id", "ds", *exog_cols]] if exog_cols else None
+            pred = forecaster._forecast_with_exog(
+                train,
+                h,
+                freq,
+                None,
+                None,
+                None,
+                horizon_x,
+                None,
+                exog_cols or None,
+                lambda: forecaster._forecast_univariate(
+                    train, h, freq, None, None, None
+                ),
+            )
+            pred = pred.sort_values("ds")
+            return pred[forecaster.alias].to_numpy(dtype=float)
+
+        scores[uid] = conformity_scores_abs_errors(
+            y,
+            h,
+            prediction_intervals.n_windows,
+            forecast_values,
+        )
+    return scores