Probabilistic Electricity Price Forecasting (Part 2)

Series: Probabilistic Electricity Price Forecasting

Part: 2 (Implementation and Results)

This post extends the earlier deterministic analysis using the same data source, Energi Data Service, with a clean date range from 2024-01-01 to 2025-09-30 UTC, excluding subsequent empty and unreported price periods, and incorporating the same post-production exogenous features.

That window covers 639 days, or 15,336 hourly observations.

1. Data and Features

TARGET   = "DK1_EUR/MWh"
DATE_COL = "HourUTC"

FEATURE_COLS = [
    "LocalPowerMWhDK1",
    "LocalPowerSelfConMWhDK1",
    "CentralPowerMWhDK1",
    "CommercialPowerMWhDK1",
    "HydroPowerMWhDK1",
    "OffshoreWindGe100MW_MWhDK1",
    "OffshoreWindLt100MW_MWhDK1",
    "OnshoreWindGe50kW_MWhDK1",
    "OnshoreWindLt50kW_MWhDK1",
    "SolarPowerGe10Lt40kW_MWhDK1",
    "SolarPowerGe40kW_MWhDK1",
    "SolarPowerLt10kW_MWhDK1",
    "SolarPowerSelfConMWhDK1",
    "PowerToHeatMWhDK1",
    "GrossConsumptionMWhDK1",
]

2. Building the QR and QRF Pipeline

The earlier theory post committed to two models: Quantile Regression (QR) and Quantile Regression Forest (QRF). This section walks through the actual implementation, step by step, and evaluates both on the DK1 data above using walk-forward cross-validation rather than a single train/test split. The full pipeline is split across qr_qrf_walkforward_pipeline.py (data, models, evaluation) and forecasting_plots.py (all figures), both included alongside this post.

2.1 Feature Engineering

Each row needs calendar features (hour, day of week, month) and lagged price features. The lags are built relative to the forecast lead time $L$, so a model trained for the 6h-ahead task only ever sees information that would actually be available 6 hours before delivery:

df["hour"] = df[DATE_COL].dt.hour
df["dayofweek"] = df[DATE_COL].dt.dayofweek
df["month"] = df[DATE_COL].dt.month
df["is_weekend"] = (df["dayofweek"] >= 5).astype(int)
df["is_peak"] = df["hour"].between(7, 20).astype(int)

# Lag features relative to lead time L
df[f"price_lag_{L}"] = df[TARGET].shift(L)
df[f"price_lag_{L + 24}"] = df[TARGET].shift(L + 24)
df["price_lag_168"] = df[TARGET].shift(168)
df["price_rolling_24h"] = df[TARGET].shift(L).rolling(24).mean()
df = df.dropna().reset_index(drop=True)

price_lag_168 is last week’s price at the same hour; price_rolling_24h is a 24-hour rolling average ending at the last point actually known before the target.

2.2 Walk-Forward Train / Validate / Test Split

A single chronological split (train on everything before a cutoff, test on what’s after) only tells us how the model does on one slice of time. DK1 prices are seasonal and regime-dependent, so this project instead uses walk-forward, expanding-window cross-validation: four folds, spread across the full 2024 to 2025 timeline, each with its own Train, then Validate, then Test window, in that order:

def generate_walk_forward_folds(min_date, max_date, n_folds, test_len_h, val_len_h, min_train_frac):
    total_hours = (max_date - min_date) / pd.Timedelta(hours=1)
    min_train_len = total_hours * min_train_frac
    remaining = total_hours - min_train_len
    step = remaining / n_folds

    folds = []
    for k in range(n_folds):
        val_start = min_date + pd.Timedelta(hours=min_train_len + k * step)
        val_end = val_start + pd.Timedelta(hours=val_len_h)
        test_start = val_end
        test_end = min(test_start + pd.Timedelta(hours=test_len_h), max_date)
        folds.append(dict(
            fold=k + 1,
            train_start=min_date, train_end=val_start,
            val_start=val_start, val_end=val_end,
            test_start=test_start, test_end=test_end,
        ))
    return folds

Training always starts from the same date and expands; validation and test windows are spaced evenly across the remaining timeline, so each fold lands in a different season:

The four walk-forward folds actually used. Fold 2’s test window falls over Christmas and New Year; the others fall in early autumn, early spring, and midsummer.

The validation window is not just a formality: it’s where QRF’s min_samples_leaf is chosen (between 20 and 40, by whichever gives the lower validation pinball loss) before either model ever sees the test window.

2.3 Fitting Quantile Regression

QR fits one linear model per quantile level, each minimizing the pinball loss directly rather than the squared error. An earlier version of this pipeline used sklearn’s QuantileRegressor, which solves this as a linear program and turned out to scale badly: about 20 seconds per quantile on this project’s larger folds. Switching to statsmodels’ QuantReg, which uses iteratively reweighted least squares instead of a linear program, cut that to roughly 1 to 4 seconds per quantile, for the same result:

import statsmodels.api as sm
from statsmodels.regression.quantile_regression import QuantReg
from sklearn.preprocessing import StandardScaler

QUANTILES = [0.025, 0.05, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, 0.95, 0.975]

def fit_quantile_regression_models(X_train, y_train, quantiles, max_iter=1000, p_tol=1e-6):
    scaler = StandardScaler()
    X_train_s = scaler.fit_transform(X_train)
    X_train_c = sm.add_constant(X_train_s, has_constant="add")

    models = {}
    for a in quantiles:
        models[a] = QuantReg(y_train, X_train_c).fit(q=a, max_iter=max_iter, p_tol=p_tol)
    return {"scaler": scaler, "models": models}


def predict_quantile_regression(fitted, X):
    X_c = sm.add_constant(fitted["scaler"].transform(X), has_constant="add")
    return {a: np.asarray(res.predict(X_c)) for a, res in fitted["models"].items()}

The models are fit once per fold, on the training window only, then reused to predict on both the validation and test windows, rather than refitting from scratch for each.

Because each quantile is fit independently, nothing stops the 0.90 line from occasionally landing below the 0.80 line, a problem called quantile crossing. The standard fix is to sort the predicted quantiles row by row after the fact:

def enforce_monotonicity(preds, quantiles):
    sorted_q = np.sort(quantiles)
    stacked = np.vstack([preds[a] for a in sorted_q])
    before_cross = float(np.mean(np.any(np.diff(stacked, axis=0) < 0, axis=0)))
    stacked_sorted = np.sort(stacked, axis=0)
    fixed = {a: stacked_sorted[k] for k, a in enumerate(sorted_q)}
    return fixed, before_cross

before_cross is kept as a diagnostic: it reports how often crossing actually happened before the fix, which turns out to be a useful signal on its own (see Discussion below).

2.4 Fitting Quantile Regression Forest

QRF follows the Meinshausen algorithm from the theory post exactly: grow a standard Random Forest, but instead of collapsing each leaf to its mean, keep every training sample that lands there. At prediction time, each tree assigns a weight to every training point depending on whether it shares a leaf with the query point; averaging these weights across trees gives the forest-weighted empirical CDF, and inverting it gives the quantile.

The first version of this class predicted one test row at a time, looping over every tree in Python, a real bottleneck once it had to run across four folds and four lead times. It’s now vectorized with sparse matrices: for each tree, leaf membership becomes a one-hot matrix, weights for the whole test set are computed in a single sparse matrix multiplication per tree, and the final quantile lookup runs as one batched cumulative sum over the entire test set instead of a per-row loop:

from scipy import sparse
from sklearn.ensemble import RandomForestRegressor

class QuantileRandomForest:
    def __init__(self, n_estimators=200, min_samples_leaf=30, random_state=42, n_jobs=-1):
        self.rf = RandomForestRegressor(
            n_estimators=n_estimators, min_samples_leaf=min_samples_leaf,
            random_state=random_state, n_jobs=n_jobs,
        )

    def fit(self, X, y):
        self.y_train_ = np.asarray(y)
        self.rf.fit(X, self.y_train_)
        self.train_leaves_ = self.rf.apply(X)  # (n_train, n_trees)
        self.n_trees_ = self.train_leaves_.shape[1]
        self.n_train_ = X.shape[0]
        return self

    def _tree_weight_matrix(self, t, test_leaf_ids):
        train_leaf_ids = self.train_leaves_[:, t]
        leaves, train_inv = np.unique(train_leaf_ids, return_inverse=True)
        train_onehot = sparse.csr_matrix(
            (np.ones(self.n_train_, dtype=np.float32), (np.arange(self.n_train_), train_inv)),
            shape=(self.n_train_, len(leaves)),
        )
        leaf_counts = np.asarray(train_onehot.sum(axis=0)).ravel()
        inv_counts = np.divide(1.0, leaf_counts, out=np.zeros_like(leaf_counts), where=leaf_counts > 0)
        train_onehot_norm = train_onehot.multiply(inv_counts).tocsr()

        leaf_to_col = {leaf: i for i, leaf in enumerate(leaves)}
        test_col = np.fromiter((leaf_to_col.get(l, -1) for l in test_leaf_ids), dtype=int, count=len(test_leaf_ids))
        valid = test_col >= 0
        test_onehot = sparse.csr_matrix(
            (np.ones(int(valid.sum()), dtype=np.float32), (np.where(valid)[0], test_col[valid])),
            shape=(len(test_leaf_ids), len(leaves)),
        )
        return test_onehot @ train_onehot_norm.T  # (n_test, n_train), sparse

    def predict_quantiles(self, X, quantiles):
        test_leaves = self.rf.apply(X)
        W = sparse.csr_matrix((X.shape[0], self.n_train_), dtype=np.float32)
        for t in range(self.n_trees_):
            W = W + self._tree_weight_matrix(t, test_leaves[:, t])
        W = (W / self.n_trees_).toarray().astype(np.float32)

        order = np.argsort(self.y_train_)
        y_sorted, cdf = self.y_train_[order], np.cumsum(W[:, order], axis=1)
        cdf[:, -1] = 1.0  # guard against floating-point residue

        return {a: y_sorted[np.argmax(cdf >= a, axis=1)] for a in quantiles}

Unlike QR, QRF is monotonic by construction: the quantiles are all read off the same weighted CDF, so a higher $\alpha$ can never produce a lower value. min_samples_leaf is selected per fold from {20, 40} using the validation window described in 2.2.

2.5 Evaluation Metrics, Averaged Across Folds

Both models are scored with the same metrics introduced in the theory post: mean pinball loss, the CRPS approximation (twice the average pinball loss over the grid), empirical coverage, and mean interval width for the 80/90/95% intervals. With four folds instead of one split, each metric is now reported as a mean across folds, with a standard deviation showing how much it varies fold to fold:

def pinball_loss(y_true, q_hat, alpha):
    err = y_true - q_hat
    return float(np.mean(np.where(err >= 0, alpha * err, (alpha - 1) * err)))

def crps_from_pinball(y_true, quantile_preds, quantiles):
    losses = [pinball_loss(y_true, quantile_preds[a], a) for a in quantiles]
    return 2.0 * np.mean(losses), np.mean(losses)

def empirical_coverage(y_true, lower, upper):
    return float(np.mean((y_true >= lower) & (y_true <= upper)))

The pipeline runs both models across four lead times (1h, 6h, 12h, 24h) and all four folds, then aggregates.

3. Results

Running qr_qrf_walkforward_pipeline.py on the real DK1 data, across all four walk-forward folds, gives the following comparison. Each number is the mean across the 4 folds, with the standard deviation in parentheses showing how much that metric actually varies from one time period to the next:

The four walk-forward folds used to produce every number above.

QR vs QRF, 1h-ahead, on the most recent fold’s two-week test window, with 80% and 95% prediction intervals.

QR vs QRF, 6h-ahead.

QR vs QRF, 12h-ahead.

QR vs QRF, 24h-ahead.

Reliability diagram (empirical vs nominal coverage) for QR and QRF, averaged across all four folds, for all four lead times.

4. Discussion

QRF beats QR on every accuracy metric, at every lead time, and the gap widens once the model is judged across multiple time periods instead of one. Median accuracy (MAE, RMSE), sharpness (mean pinball loss), and CRPS are all better for the forest than for the linear model, from 1h through 24h. This matches the earlier single-split finding and the theoretical expectation: DK1 prices are nonlinear and heavy-tailed, and a linear model per quantile can only do so much, while a forest handles nonlinear interactions between generation mix, calendar effects, and lagged prices directly.

Calibration is far less stable across time than a single train/test split would suggest. The standard deviations on the coverage numbers are large, often 0.10 to 0.18 on a nominal target of 0.80 to 0.95. Both models do well on some folds and poorly on others. A single split, like the earlier version of this pipeline used, would only ever show one point on that range, and could easily give a falsely confident or falsely pessimistic impression depending on which period it happened to land on.

Two of the four folds are consistently the hardest to calibrate, for both models. Fold 1 (test window 19 September to 3 October 2024) and Fold 3 (30 March to 13 April 2025) calibrate well, close to nominal coverage at every lead time. Fold 2 (24 December 2024 to 7 January 2025, spanning Christmas and New Year) and Fold 4 (4 to 18 July 2025) calibrate noticeably worse, sometimes dropping to 53 to 66% actual coverage against an 80% nominal target for both QR and QRF. A plausible explanation for Fold 2: the calendar features (is_weekend, dayofweek, is_peak) encode a normal working week, and public holidays break that pattern in ways the model hasn’t seen enough of. Walk-forward validation is built to catch exactly this kind of failure; a single split would have hidden it.

Quantile crossing for QR hasn’t gone away either, and it lines up with the worst-calibrated folds. Averaged across folds it ranges from under 1% at 1h ahead to over 20% at 6h ahead in the worst fold, spiking in the same folds where coverage collapses. QRF has no such issue by construction.

All models degrade sharply from 1h to 6h, then roughly plateau. This matches the single-split finding: the lag features carry the most information about the very next hour, and that advantage fades quickly once the horizon exceeds the shortest lag.

QRF remains the stronger model of the two overall, and walk-forward validation shows that conclusion holds across different seasons and market regimes, not just the one three-month window the original single split happened to test on. The bigger takeaway is that both models’ calibration is regime-dependent: any deployment of this pipeline should expect noticeably worse-calibrated intervals during holiday periods until the feature set explicitly accounts for them.

Code

• qr_qrf_walkforward_pipeline.py: data loading, feature engineering, the QR and QRF model classes, and the walk-forward evaluation loop.
• forecasting_plots.py: every figure in this post, including the train/validate/test split diagram and the reliability diagram.
• Results/ folder: per-fold metrics and the cross-fold summary (mean and standard deviation per model per lead time).

Next: Part 3 compares two more approaches, bootstrapped residuals and split conformal prediction, and puts all four methods head to head.