Probabilistic Electricity Price Forecasting (Part 3)

Series: Probabilistic Electricity Price Forecasting

Part: 3 (Bootstrapped Residuals and Conformal Prediction)

QR and QRF, covered in Part 2, estimate the conditional distribution directly: QR fits one linear model per quantile, QRF reads every quantile off a single forest’s weighted empirical CDF. This post covers two different, older techniques that instead take a single point-forecast model and turn its residuals into an interval. Both come from skforecast’s user guides (Bootstrapped Residuals, Conformal Prediction), and this write-up adapts the concepts to this project’s own XGBoost pipeline and tests them on real DK1 data rather than skforecast’s bike-sharing/ETT examples.

1. Bootstrapped Residuals

The idea: fit one point-forecast model, collect its residuals $e_t = y_t - \hat y_t$, then simulate alternative future outcomes by resampling from that residual pool and adding the draws back to the point forecast. Repeating this many times (n_boot draws) builds an empirical distribution at each forecast horizon; percentiles of that distribution become the prediction interval.

Repeated resampling of past residuals builds a collection of alternative predictions, whose spread represents forecast uncertainty. Source: skforecast.org, Bootstrapped Residuals.

Taking the $\alpha/2$ and $1-\alpha/2$ percentiles of that spread, at every forecasting horizon, turns the individual bootstrap paths into a single interval band:

Each gray line is one bootstrap path; the red band is the resulting prediction interval once percentiles are taken at each horizon. Source: skforecast.org, Bootstrapped Residuals.

“One of the main advantages of this strategy is that it requires only a single model to estimate any interval. However, performing hundreds or thousands of bootstrapping iterations can be computationally expensive.” (skforecast, Bootstrapped Residuals)

Which residuals: in-sample vs. out-of-sample

This is the part skforecast’s docs spend the most time on, and it turns out to matter enormously for DK1.

• In-sample residuals: computed on the same data the model was trained on. The model has already fit that data, so its errors there systematically understate how wrong it will be on new data. skforecast’s own docs flag this: “this can result in intervals that are too narrow (overly optimistic).”
• Out-of-sample residuals: computed on a held-out calibration/validation split the model never trained on. Wider, but honest.

This project’s forecast.py currently uses in-sample residuals: train_residuals = y_train.values - model.predict(X_train) is fit directly on the training set the model just saw. Testing this on real DK1 data (1h-ahead XGBoost, last walk-forward fold) shows exactly the effect skforecast warns about:

Real DK1 residuals from the same XGBoost model: in-sample (evaluated on training data) vs out-of-sample (evaluated on a held-out validation window).

A 5x gap. The consequence on actual interval coverage, measured on held-out test data:

The in-sample intervals claim 80% coverage and deliver barely half that: they are badly overconfident, and this is what the existing code does today, not a hypothetical failure mode. The fix is to compute residuals on a held-out split rather than the training set; the walk-forward validation split already built in qr_qrf_walkforward_pipeline.py is the natural place to source these residuals from.

Binned residuals (conditioning on heteroscedasticity)

Bootstrapping assumes residuals are exchangeable regardless of the situation, which is rarely true: errors are usually larger when the underlying value is larger or more volatile. skforecast partitions residuals into bins (by default, bins of the predicted value) so that bootstrapping samples from the right slice of the error distribution depending on what’s being predicted.

This project’s BootstrappedPI class already bins, but by hour-of-day, not predicted value, on the reasoning that DK1 volatility is hour-dependent (peak vs. off-peak). Testing both binning strategies against no binning at all, using an 8-week calibration window (large enough for about 56 samples per hour bin) on the same test fold:

Empirical coverage against the nominal target, for three ways of conditioning the residual pool: no binning, hour-of-day, and predicted-value quantile bins.

Hour-binning makes the interval narrower but noticeably worse-calibrated on this fold, 73% actual coverage against an 80% claim, while value-binning performs about the same as the unbinned baseline. A likely explanation: the XGBoost point forecast already includes hour as an input feature, so much of the systematic hour-of-day heteroscedasticity is already absorbed into the point prediction itself, leaving less for hour-binned residuals to explain. Binning by the predicted value instead, which reflects the model’s own sense of how extreme the situation is, captures residual heteroscedasticity that hour alone misses.

For this project, that means switching BootstrappedPI’s conditioning variable from hour-of-day to predicted-value bins, or testing both properly within the walk-forward framework before committing, and always sourcing the residual pool from a held-out split.

2. Conformal Prediction (Split Conformal)

Conformal prediction solves the same problem, turning a point forecast plus residuals into a calibrated interval, without the resampling loop. Skforecast implements Split Conformal Prediction (SCP): reserve a calibration set, compute residuals on it once, then read the $(\alpha/2, 1-\alpha/2)$ empirical quantiles of those residuals directly as a fixed offset added to every future point forecast.

Conformal regression turns point predictions into prediction intervals using a single correction learned once from a calibration set. Source: skforecast.org, Conformal Prediction, citing Christoph Molnar, Introduction To Conformal Prediction With Python.

“Conformal methods can also calibrate prediction intervals generated by other techniques, such as quantile regression or bootstrapped residuals… Skforecast implements Split Conformal Prediction (SCP) due to its simplicity and efficiency.” (skforecast, Conformal Prediction)

No draws, no resampling, just one quantile lookup from the calibration residuals. A first, one-off timing test on a single fold (336 test rows, one 80% interval, fully vectorized) suggested roughly 54x, but that number was optimistic: it timed the narrowest possible case. Running both methods properly, across all 4 walk-forward folds and all 4 lead times, with per-row binned lookups for the full 13-quantile grid (needed for pinball loss and CRPS, not just one interval), gives a more honest and still solid speedup:

The gap narrows from the earlier single-shot estimate because computing 13 quantile levels with per-row bin lookups carries its own Python-level overhead for both methods; conformal’s advantage is the resampling step it skips, not the binning logic, which both methods share. Still a clear, consistent win, just a more modest one than the first quick test suggested.

Results: Walk-Forward Comparison

Both methods were run through the exact same walk-forward folds as QR and QRF (bootstrap_conformal_walkforward_pipeline.py), calibrated on out-of-sample validation residuals, value-binned, wrapping the same XGBoost point-forecast family used in forecast.py:

Bootstrap and Conformal are statistically indistinguishable on point accuracy and sharpness, exactly as expected: both derive from the same calibration-residual quantiles.

Coverage tells a different story. Bootstrap consistently covers more than Conformal, at every lead time and every nominal level, while also being wider:

The single-fold test didn’t reveal this: resampling with replacement from a calibration pool (Bootstrap) introduces extra variance beyond the pool’s raw empirical quantiles, which widens the tails and pushes coverage slightly higher than the direct quantile lookup (Conformal) achieves from the exact same pool. Neither dominates. Conformal is sharper, narrower intervals for near-identical accuracy; Bootstrap is safer, better coverage at the cost of width. Which one to prefer depends on whether the downstream decision cares more about tight intervals or about not being caught outside them.

Bootstrap vs Conformal, 1h ahead, on the most recent fold’s test window.

Empirical vs nominal coverage, averaged across all 4 folds, for both methods at all 4 lead times. Conformal sits consistently below Bootstrap, matching the coverage tables above.

Real timing, averaged across all 16 fold/lead-time combinations.

How do all four methods compare?

Putting Bootstrap and Conformal’s mean pinball loss next to QR and QRF’s from Part 2:

QR, the linear model, is the weakest of the four at every horizon, as it was on its own. At 1h, 6h, and 12h, though, the XGBoost-plus-residuals methods (Bootstrap, Conformal) actually beat QRF on sharpness, sometimes by a wide margin (4.93 vs 5.55 at 6h); only at 24h does QRF pull back ahead. On coverage, the pattern reverses: QRF’s 80% interval coverage (0.726 to 0.862 across horizons) is consistently closer to nominal than Bootstrap’s (0.682 to 0.781) or Conformal’s (0.667 to 0.763). The sharpest method here isn’t the best-calibrated one, so both properties need to be reported together rather than picking a single winner.

The multi-step-ahead caveat

Skforecast is explicit about this: SCP’s coverage guarantee is only rigorous for one-step-ahead predictions. For multi-step horizons, the coverage probability is not theoretically guaranteed, it has to be checked empirically.

“When applied to time series forecasting, their coverage guarantees are only valid for one-step-ahead predictions. For multi-step-ahead predictions, the coverage probability is not guaranteed.” (skforecast, Conformal Prediction)

This matters directly for this project, which forecasts 1h, 6h, 12h, and 24h ahead. There’s no shortcut around it: coverage has to be measured empirically at each horizon, exactly as the walk-forward pipeline already does for QR and QRF. Conformal prediction should be evaluated the same way, per lead time, rather than assumed valid because it “is conformal.”

3. Implementation for This Project

The simplified version below shows the core mechanism. The actual implementation used to produce every result in this post, including both the single-fold binning comparison above and the full walk-forward comparison below, lives in bootstrap_conformal_walkforward_pipeline.py (included alongside this post), alongside a matching BootstrappedPI class built the same way.

import numpy as np

class ConformalPI:
    """
    Split Conformal Prediction (SCP), adapted from skforecast's approach
    (https://skforecast.org/0.15.0/user_guides/probabilistic-forecasting-conformal-prediction).

    Fits once on a held-out calibration set's residuals, no resampling.
    Optionally conditions the residual quantiles on a bin variable
    (predicted value or hour) for adaptive/heteroscedasticity-aware
    intervals, matching skforecast's `use_binned_residuals`.
    """

    def __init__(self, coverage_levels=(0.80, 0.90, 0.95), bin_by=None, n_bins=5, min_bin_size=15):
        self.coverage_levels = coverage_levels
        self.bin_by = bin_by  # None, "value", or "hour"
        self.n_bins = n_bins
        self.min_bin_size = min_bin_size

    def fit(self, cal_residuals, cal_point_preds=None, cal_hours=None):
        self._residuals = np.asarray(cal_residuals)
        self._bin_edges = None

        if self.bin_by == "value":
            self._bin_edges = np.quantile(cal_point_preds, np.linspace(0, 1, self.n_bins + 1))
            bins = np.clip(np.digitize(cal_point_preds, self._bin_edges[1:-1]), 0, self.n_bins - 1)
        elif self.bin_by == "hour":
            bins = np.asarray(cal_hours, dtype=int)
        else:
            bins = np.zeros(len(self._residuals), dtype=int)

        self._offsets = {}
        for b in np.unique(bins):
            pool = self._residuals[bins == b]
            pool = pool if len(pool) >= self.min_bin_size else self._residuals
            for cov in self.coverage_levels:
                alpha = 1 - cov
                self._offsets[(b, cov)] = (
                    np.quantile(pool, alpha / 2),
                    np.quantile(pool, 1 - alpha / 2),
                )
        self._bins_seen = set(np.unique(bins))
        return self

    def predict(self, point_preds, point_or_hour=None):
        point_preds = np.asarray(point_preds)
        if self.bin_by == "value":
            bins = np.clip(np.digitize(point_or_hour, self._bin_edges[1:-1]), 0, self.n_bins - 1)
        elif self.bin_by == "hour":
            bins = np.asarray(point_or_hour, dtype=int)
        else:
            bins = np.zeros(len(point_preds), dtype=int)

        result = {}
        for cov in self.coverage_levels:
            lo, hi = np.empty(len(point_preds)), np.empty(len(point_preds))
            for i, b in enumerate(bins):
                b_use = b if b in self._bins_seen else min(self._bins_seen, key=lambda x: abs(x - b))
                lo_off, hi_off = self._offsets[(b_use, cov)]
                lo[i], hi[i] = point_preds[i] + lo_off, point_preds[i] + hi_off
            result[f"lower_{int(cov * 100)}"] = lo
            result[f"upper_{int(cov * 100)}"] = hi
        return result

Usage, matching this project’s walk-forward folds:

# 1. Fit the point-forecast model on the fold's training window only
model.fit(X_train, y_train)

# 2. Calibrate on the held-out validation window (out-of-sample residuals)
cal_point_preds = model.predict(X_val)
cal_residuals = y_val - cal_point_preds

conformal = ConformalPI(coverage_levels=(0.80, 0.90, 0.95), bin_by="value", n_bins=5)
conformal.fit(cal_residuals, cal_point_preds=cal_point_preds)

# 3. Predict intervals on the test window
point_test = model.predict(X_test)
intervals = conformal.predict(point_test, point_or_hour=point_test)

The same ConformalPI object, with bin_by=None, also reproduces the corrected (out-of-sample) bootstrap-equivalent interval: the two techniques share the same calibration-set requirement, and conformal is simply the fast, closed-form way to get there.

4. Summary: What to Change in This Project

1. Fix BootstrappedPI to use out-of-sample residuals. It currently calibrates on training-set residuals, which this data shows understate the true error spread by roughly 5x, and under-cover its own nominal target by nearly 30 percentage points (52% actual vs 80% claimed) on a single test fold. bootstrap_conformal_walkforward_pipeline.py implements the fix: calibrate on the walk-forward validation split instead.
2. Reconsider the binning variable. Hour-of-day binning, this project’s current choice, was empirically worse-calibrated than no binning at all in the single-fold test; value-binning (skforecast’s default) matched the unbinned baseline and is what the full walk-forward pipeline now uses for both methods.
3. Conformal is meaningfully faster, but the margin is smaller than a quick test suggests. Measured across all 16 fold/lead-time combinations: about 7.4x (5.5 ms vs 40.3 ms per fold), not the ~54x a single fully-vectorized test implied. Still a solid, consistent win for anywhere the existing bootstrap approach’s computational cost has been a bottleneck.
4. Pick Bootstrap or Conformal based on what the decision needs, not blindly. Across all four lead times, Bootstrap consistently achieves better coverage (closer to nominal) at the cost of wider intervals, and Conformal is consistently sharper but less well-calibrated. Neither is a strict upgrade over the other.
5. QRF is not automatically the best model once all four are compared. At 1h, 6h, and 12h ahead, the XGBoost-plus-residuals methods (Bootstrap, Conformal) actually beat QRF on mean pinball loss; only at 24h does QRF pull ahead. QRF does have the most consistently well-calibrated coverage of the four, though, so the right choice again depends on whether sharpness or calibration matters more for the use case.
6. Validate coverage per lead time, not just once. Split conformal’s coverage guarantee formally holds only for one-step-ahead forecasts; at 6h/12h/24h it must be checked empirically, exactly as this project’s walk-forward pipeline now does for all four methods.

References

• Skforecast. Probabilistic Forecasting: Bootstrapped Residuals. skforecast.org/0.15.0/user_guides/probabilistic-forecasting-bootstrapped-residuals
• Skforecast. Probabilistic Forecasting: Conformal Prediction. skforecast.org/0.15.0/user_guides/probabilistic-forecasting-conformal-prediction
• MAPIE. Theoretical Description, Split Conformal Prediction. mapie.readthedocs.io/en/stable/theoretical_description_regression.html
• Molnar, C. Introduction To Conformal Prediction With Python. leanpub.com/conformal-prediction

Code

• bootstrap_conformal_walkforward_pipeline.py: the BootstrappedPI and ConformalPI classes and the walk-forward evaluation loop for both.
• qr_qrf_walkforward_pipeline.py: shared data loading, feature engineering, and fold generation (same file as Part 2, included here so this post runs standalone).
• forecasting_plots.py: every figure in this post.
• Results/ folder: per-fold metrics, the cross-fold summary, and the timing summary.