Single horizon predictive modeling

Environment setup

We need to install some extra dependencies for this notebook if needed (when running jupyterlite).

%pip install -q skrub altair holidays plotly nbformat polars

Note: you may need to restart the kernel to use updated packages.

import warnings
from pathlib import Path

import altair
import cloudpickle
import numpy as np
import polars as pl
import pyarrow  # noqa: F401
import skrub
import tzdata  # noqa: F401
from plotly.io import read_json, write_json  # noqa: F401
from sklearn.ensemble import HistGradientBoostingRegressor

from feature_engineering_lib import feature_engineering_outputs
from tutorial_helpers import (
    collect_cv_predictions,
    plot_binned_residuals,
    plot_lorenz_curve,
    plot_reliability_diagram,
    plot_residuals_vs_predicted,
)

# Ignore warnings from pkg_resources triggered by Python 3.13's multiprocessing.
warnings.filterwarnings("ignore", category=UserWarning, module="pkg_resources")

In this notebook we focus on predicting the electricity demand at a single, fixed horizon (12 hours). In the next one we will alter our pipeline to predict multiple horizons.

Cross-validation splitter

In the previous notebook we built predictive features. Here we will complete the pipeline and we will evaluate it with cross-validation. Therefore, we need to define a cross-validation splitter that will divide our data into train and test sets.

We define a custom time-based cross-validation splitter. We do not use scikit-learn’s TimeSeriesSplit because it is based on positional indices, while here we can have an irregular grid after dropping rows with missing ground truth. It is also easier to inspect and debug splits based on actual dates and a datetime column than on row positions.

In this implementation, the test windows are 3-month blocks. We start from the end of the historical data and go back in regular 3-month steps, so that the last validation block never has an arbitrary length. We stop when the training set would be shorter than an initial training period of two years. Each fold keeps a 7-day gap between the end of training and the start of the test window. This is an example of a Backtesting with intermittent refit including gap setup that is common in forecasting problems, where split boundaries must mimic operational constraints in deployment.

import datetime

from dateutil.relativedelta import relativedelta


def _split_indices(X, test_start_date, test_end_date, gap_days=7):
    train = (
        X.with_row_index()
        .filter(pl.col("prediction_time") < test_start_date - datetime.timedelta(days=gap_days))[
            "index"
        ]
        .to_numpy()
    )
    test = (
        X.with_row_index()
        .filter(
            (pl.col("prediction_time") >= test_start_date)
            & (pl.col("prediction_time") < test_end_date)
        )["index"]
        .to_numpy()
    )
    return train, test


class TimeSeriesSplitter:
    train_test_gap_days = 7
    test_blocks = 3

    def split(self, X, y=None, groups=None, blocks=None):
        if blocks is None:
            blocks = self.test_blocks
        min_train_days = 365 * 2  # Initial train period: 2 years
        min_date = X["prediction_time"].min()
        max_date = X["prediction_time"].max()

        first_allowed = (
            min_date
            + relativedelta(days=min_train_days)
            + datetime.timedelta(days=self.train_test_gap_days)
        )

        # Start from the end of the data and go back by `blocks` months at a time.
        current_test_end = max_date.replace(
            day=1, hour=0, minute=0, second=0, microsecond=0
        ) + relativedelta(months=1)
        folds = []
        while True:
            test_start = current_test_end - relativedelta(months=blocks)
            if test_start < first_allowed:
                break
            folds.append((test_start, current_test_end))
            current_test_end = test_start

        for test_start, test_end in reversed(folds):
            train, test = _split_indices(X, test_start, test_end, gap_days=self.train_test_gap_days)
            if len(train) and len(test):
                yield train, test

    def get_n_splits(self, X, y=None, groups=None):
        return len(list(self.split(X, y)))

TIME_HORIZON = 12  # Focus on 12-hour horizon
features, y = feature_engineering_outputs(horizons=TIME_HORIZON, cv_splitter=TimeSeriesSplitter())

Adding a supervised regressor

To complete our pipeline, we add a HistGradientBoostingRegressor that will use the features defined in the previous notebook to predict the electricity demand.

def get_hgb_params(max_iter=100):
    return dict(
        random_state=0,
        loss=skrub.choose_from(["squared_error", "poisson", "gamma"], name="loss"),
        early_stopping=True,
        n_iter_no_change=50,
        max_depth=None,
        learning_rate=skrub.choose_float(0.01, 0.5, default=0.1, log=True, name="learning_rate"),
        max_leaf_nodes=skrub.choose_int(2, 128, default=31, log=True, name="max_leaf_nodes"),
        min_samples_leaf=skrub.choose_int(2, 128, default=20, log=True, name="min_samples_leaf"),
    )


def get_regressor(max_iter=100):
    # max_iter = 100 for fast notebook but should set to 1000 or more for best results
    return HistGradientBoostingRegressor(**get_hgb_params(max_iter=max_iter))


pred = features.skb.apply(get_regressor(), y=y).skb.with_scoring(
    ["neg_mean_absolute_percentage_error", "r2"]
)
pred

<Scoring <Apply HistGradientBoostingRegressor> (2 scorers)>

    This DataOp will be scored with:
      - 'neg_mean_absolute_percentage_error'
      - 'r2'
    Use .skb.cross_validate(…) or .skb.make_learner(…).score(…) to compute scores.

Show/Hide graph

Result (preview of the learner's output, not the scores):

Please enable javascript

The skrub table reports need javascript to display correctly. If you are displaying a report in a Jupyter notebook and you see this message, you may need to re-execute the cell or to trust the notebook (button on the top right or "File > Trust notebook").

Assessing the model performance via cross-validation

pred.skb.cross_validate()

	fit_time	score_time	test_neg_mean_absolute_percentage_error	test_r2
0	0.599063	0.186611	-0.028182	0.927435
1	0.660780	0.192731	-0.021063	0.970250
2	0.686824	0.188721	-0.023043	0.949349
3	0.696022	0.187387	-0.026650	0.949349
4	0.728424	0.189730	-0.016591	0.970999
5	0.746868	0.189618	-0.019243	0.972456
6	0.776612	0.192372	-0.023370	0.930070
7	0.793216	0.187886	-0.022071	0.971033

For further inspection of predictions, we will collect the cross-validated prediction into a dataframe.

def get_cv_results(pred, return_train_score=False):
    """Collect cross-validated predictions and scores.

    The predictions can be a numpy array (one column per target in ``y_test``)
    or a polars dataframe (e.g. several quantiles per horizon, indicated by the
    column names). In both cases they are returned as ``pred_<column>`` columns
    next to ``X_test`` and ``y_test``.
    """
    predictions = []
    scores = []
    for i, split in enumerate(pred.skb.iter_cv_splits()):
        learner = pred.skb.make_learner().fit(split["train"])

        split_scores, split_predictions = learner.score(split["test"], return_predictions=True)
        if return_train_score:
            split_scores.update({f"train_{k}": v for k, v in learner.score(split["train"]).items()})
        scores.append(split_scores | {"split": i})
        y_test = pl.DataFrame(split["y_test"])
        pred_values = split_predictions["predict"]
        if not isinstance(pred_values, pl.DataFrame):
            pred_values = np.asarray(pred_values)
            if pred_values.ndim == 1:
                pred_values = pred_values[:, None]
            pred_values = pl.DataFrame(pred_values, schema=y_test.columns)
        predictions.append(
            pl.concat(
                [
                    split["X_test"],
                    y_test,
                    pred_values.rename("pred_{}".format),
                ],
                how="horizontal",
            ).with_columns(split=pl.lit(i))
        )
        print(f"split {i}:", split["X_test"]["prediction_time"].min().isoformat())
        print(split_scores)

    return pl.concat(predictions, how="vertical"), pl.DataFrame(scores)


cv_predictions, cv_scores = get_cv_results(pred)

split 0: 2023-06-01T00:00:00+00:00
{'neg_mean_absolute_percentage_error': -0.02818191714011722, 'r2': 0.92743480539659}

split 1: 2023-09-01T00:00:00+00:00
{'neg_mean_absolute_percentage_error': -0.021062619514067016, 'r2': 0.970249951860939}

split 2: 2023-12-01T00:00:00+00:00
{'neg_mean_absolute_percentage_error': -0.02304311321256355, 'r2': 0.9493488394475555}

split 3: 2024-03-01T00:00:00+00:00
{'neg_mean_absolute_percentage_error': -0.026650471567178303, 'r2': 0.949349241518662}

split 4: 2024-06-01T00:00:00+00:00
{'neg_mean_absolute_percentage_error': -0.016590871382778576, 'r2': 0.9709993319737483}

split 5: 2024-09-01T00:00:00+00:00
{'neg_mean_absolute_percentage_error': -0.019242750207894135, 'r2': 0.9724555377162379}

split 6: 2024-12-01T00:00:00+00:00
{'neg_mean_absolute_percentage_error': -0.023370467183172653, 'r2': 0.9300704289002145}

split 7: 2025-03-01T00:00:00+00:00
{'neg_mean_absolute_percentage_error': -0.022071030361231234, 'r2': 0.9710326349135896}

As a sanity check, we will take a look at the predictions on the first fold and plot the observed values and the prediction values from the model. We limit the visualization to the last 7 days of the fold.

altair.Chart(cv_predictions.tail(100)).transform_fold(
    ["12h", "pred_12h"],
).mark_line(
    tooltip=True
).encode(x="prediction_time:T", y="value:Q", color="key:N").interactive()

Now, let’s check the performance of our models.

The first curve is called the Lorenz curve. It shows on the x-axis the fraction of observations sorted by predicted values and on the y-axis the cumulative observed load proportion.

plot_lorenz_curve(cv_predictions, TIME_HORIZON).interactive()

The diagonal on the plot corresponds to a model predicting a constant value that is therefore not an informative model. The oracle model corresponds to the “perfect” model that would provide the an output identical to the observed values. Thus, the ranking of such hypothetical model is the best possible ranking. However, you should note that the oracle model is not the line passing through the right-hand corner of the plot. Instead, this curvature is defined by the distribution of the observations. Indeed, more the observations are composed of small values and a couple of large values, the more the oracle model is closer to the right-hand corner of the plot.

A true model is navigating between the diagonal and the oracle model. The area between the diagonal and the Lorenz curve of a model is called the Gini index.

For our usecase, we observe that each oracle model is not far from the diagonal. It means that the observed values do not contain a couple of large values with high variability. Therefore, it informs us that the complexity of our problem at hand is not too high. Looking at the Lorenz curve of each model, we observe that it is quite close to the oracle model. Therefore, the gradient boosting regressor is discriminative for our task.

Then, we have a look at the reliability diagram. This diagram shows on the x-axis the mean predicted load and on the y-axis the mean observed load.

plot_reliability_diagram(cv_predictions, TIME_HORIZON).interactive().properties(
    title="Reliability diagram from cross-validation predictions"
)

The diagonal on the reliability diagram corresponds to the best possible model: for a level of predicted load that fall into a bin, then the mean observed load is also in the same bin. If the line is above the diagonal, it means that our model is predicted a value too low in comparison to the observed values. If the line is below the diagonal, it means that our model is predicted a value too high in comparison to the observed values.

For our cross-validated model, we observe that each reliability curve is close to the diagonal. We only observe a mis-calibration for the extremum values.

plot_residuals_vs_predicted(cv_predictions, TIME_HORIZON).interactive().properties(
    title="Residuals vs Predicted Values from cross-validation predictions"
)

plot_binned_residuals(cv_predictions, TIME_HORIZON, by="hour").interactive().properties(
    title="Residuals by hour of the day from cross-validation predictions"
)

plot_binned_residuals(cv_predictions, TIME_HORIZON, by="month").interactive().properties(
    title="Residuals by hour of the day from cross-validation predictions"
)

Exercise: non-linear feature engineering coupled with linear predictive model

Now, it is your turn to make a predictive model. Towards this end, we request you to preprocess the input features with non-linear feature engineering:

• the first step is to impute the missing values using a SimpleImputer. Make sure to include the indicator of missing values in the feature set (i.e. look at the add_indicator parameter);

• use a SplineTransformer to create non-linear features. Use the default parameters but make sure to set sparse_output=True since it subsequent processing will be faster and more memory efficient with such data structure;

• use a VarianceThreshold to remove features with potential constant features;

• use a SelectKBest to select the most informative features. Set k to be chosen from a log-uniform distribution between 100 and 1,000 (i.e. use skrub.choose_int);

• use a Nystroem to approximate an RBF kernel. Set n_components to be chosen from a log-uniform distribution between 10 and 200 (i.e. use skrub.choose_int).

• finally, use a Ridge as the final predictive model. Set alpha to be chosen from a log-uniform distribution between 1e-6 and 1e3 (i.e. use skrub.choose_float).

Use a scikit-learn Pipeline using make_pipeline to chain the steps together.

Chaining several .skb.apply(...) calls can be useful to inspect intermediate outputs in reports, but for this exercise we keep a single pipeline. With the current behavior of SplineTransformer(sparse_output=True), step-by-step .skb.apply(...) would require no_wrap=True to avoid wrapping sparse output into a dataframe.

Once the predictive model is defined, apply it on X and pass y as the target.

# Here we provide all the imports for creating the predictive model.
from functools import partial

from sklearn.feature_selection import SelectKBest, VarianceThreshold, f_regression
from sklearn.impute import SimpleImputer
from sklearn.kernel_approximation import Nystroem
from sklearn.linear_model import Ridge
from sklearn.metrics import get_scorer, make_scorer, mean_absolute_percentage_error
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import SplineTransformer

# Write your code here.
#
#
#
#
#
#
#
#
#
#
#

def make_ridge_predictions(features, y):
    return features.skb.apply(
        make_pipeline(
            SimpleImputer(add_indicator=True),
            SplineTransformer(sparse_output=True),
            VarianceThreshold(threshold=1e-6),
            SelectKBest(
                score_func=partial(f_regression, force_finite=True),
                k=skrub.choose_int(100, 400, log=True, name="n_selected_splines"),
            ),
            Nystroem(
                n_components=skrub.choose_int(10, 200, log=True, name="n_components", default=150)
            ),
            Ridge(alpha=skrub.choose_float(1e-6, 1e3, log=True, name="alpha", default=1e-2)),
        ),
        y=y,
    ).skb.with_scoring(["neg_mean_absolute_percentage_error", "r2"])


predictions_ridge = make_ridge_predictions(features, y)
predictions_ridge

<Scoring <Apply Pipeline> (2 scorers)>

    This DataOp will be scored with:
      - 'neg_mean_absolute_percentage_error'
      - 'r2'
    Use .skb.cross_validate(…) or .skb.make_learner(…).score(…) to compute scores.

Show/Hide graph

Result (preview of the learner's output, not the scores):

Please enable javascript

The skrub table reports need javascript to display correctly. If you are displaying a report in a Jupyter notebook and you see this message, you may need to re-execute the cell or to trust the notebook (button on the top right or "File > Trust notebook").

Now that you defined the predictive model, let’s evaluate the performance of the model using cross-validation. Use the time-based cross-validation splitter defined earlier in TimeSeriesSplitter. Make sure to compute the R2 score and the mean absolute percentage error. Return the training scores as well as the fitted pipeline such that we can make additional analysis.

# Write your code here.
#
#
#
#
#
#
#
#
#
#
#

cv_predictions_ridge, cv_scores_ridge = get_cv_results(predictions_ridge, return_train_score=True)

split 0: 2023-06-01T00:00:00+00:00
{'neg_mean_absolute_percentage_error': -0.03480123919990411, 'r2': 0.8815176903197758, 'train_neg_mean_absolute_percentage_error': -0.029757130541505247, 'train_r2': 0.9642633304808311}

split 1: 2023-09-01T00:00:00+00:00
{'neg_mean_absolute_percentage_error': -0.03236115011847307, 'r2': 0.9308733249225162, 'train_neg_mean_absolute_percentage_error': -0.029723022120445634, 'train_r2': 0.966202468017059}

split 2: 2023-12-01T00:00:00+00:00
{'neg_mean_absolute_percentage_error': -0.033093632109240094, 'r2': 0.8998589720086219, 'train_neg_mean_absolute_percentage_error': -0.028436629483585517, 'train_r2': 0.9677636462093347}

split 3: 2024-03-01T00:00:00+00:00
{'neg_mean_absolute_percentage_error': -0.03636323692486444, 'r2': 0.9024017722335035, 'train_neg_mean_absolute_percentage_error': -0.030307498697119205, 'train_r2': 0.9630811428996459}

split 4: 2024-06-01T00:00:00+00:00
{'neg_mean_absolute_percentage_error': -0.0258822568910276, 'r2': 0.9296507612016761, 'train_neg_mean_absolute_percentage_error': -0.029480311817503217, 'train_r2': 0.964913473023066}

split 5: 2024-09-01T00:00:00+00:00
{'neg_mean_absolute_percentage_error': -0.028220508746968062, 'r2': 0.9398203750821706, 'train_neg_mean_absolute_percentage_error': -0.028993711573717766, 'train_r2': 0.9652314571448289}

split 6: 2024-12-01T00:00:00+00:00
{'neg_mean_absolute_percentage_error': -0.028205507424844938, 'r2': 0.8910346747240206, 'train_neg_mean_absolute_percentage_error': -0.028738571051164216, 'train_r2': 0.9650984467678659}

split 7: 2025-03-01T00:00:00+00:00
{'neg_mean_absolute_percentage_error': -0.03256502938753288, 'r2': 0.9332220704850676, 'train_neg_mean_absolute_percentage_error': -0.02890865692975494, 'train_r2': 0.9657723082509972}

Do a sanity check by plotting the observed values and predictions for the first fold as we did earlier.

Then, make an analysis of the cross-validated metrics. Does this model perform better or worse than the previous model? Is it underfitting or overfitting?

# Write your code here.
#
#
#
#
#
#
#
#
#
#
#

cv_scores_ridge

shape: (8, 5)

neg_mean_absolute_percentage_error	r2	train_neg_mean_absolute_percentage_error	train_r2	split
f64	f64	f64	f64	i64
-0.034801	0.881518	-0.029757	0.964263	0
-0.032361	0.930873	-0.029723	0.966202	1
-0.033094	0.899859	-0.028437	0.967764	2
-0.036363	0.902402	-0.030307	0.963081	3
-0.025882	0.929651	-0.02948	0.964913	4
-0.028221	0.93982	-0.028994	0.965231	5
-0.028206	0.891035	-0.028739	0.965098	6
-0.032565	0.933222	-0.028909	0.965772	7

altair.Chart(cv_predictions_ridge.tail(24 * 7)).transform_fold(
    ["12h", "pred_12h"],
).mark_line(
    tooltip=True
).encode(x="prediction_time:T", y="value:Q", color="key:N").interactive()

Compute the Lorenz curve and the reliability diagram for this pipeline.

# Write your code here.
#
#
#
#
#
#
#
#
#
#
#

plot_lorenz_curve(cv_predictions_ridge, TIME_HORIZON).interactive()

plot_reliability_diagram(cv_predictions_ridge, TIME_HORIZON).interactive().properties(
    title="Reliability diagram from cross-validation predictions"
)

Hyperparameter search

To keep the notebook fast we reload an existing search that has 100 iterations (specified by storage and study_name), and resume it to add 2 more iterations.

randomized_search_ridge = predictions_ridge.skb.make_randomized_search(
    backend="optuna",
    storage="sqlite:///single_horizon_prediction_ridge.optuna",
    study_name="randomized_search",
    refit="r2",
    n_iter=2,
    fitted=True,
    verbose=1,
    n_jobs=-1,
)

Running optuna search for study randomized_search in storage sqlite:///single_horizon_prediction_ridge.optuna

[I 2026-11-09 15:50:32,178] Using an existing study with name 'randomized_search' instead of creating a new one.

[I 2026-11-09 15:50:34,483] Using an existing study with name 'randomized_search' instead of creating a new one.
[I 2026-11-09 15:50:34,501] Using an existing study with name 'randomized_search' instead of creating a new one.

[Parallel(n_jobs=1)]: Done   8 out of   8 | elapsed:   12.4s finished
[I 2026-11-09 15:50:47,448] Trial 128 finished with value: 0.9262830894636531 and parameters: {'0:temperature_only': '0:True', '1:cities': "0:'all'", '2:n_selected_splines': 328, '3:n_components': 132, '4:alpha': 4.454586997246425e-05}. Best is trial 87 with value: 0.935208416218609.

[Parallel(n_jobs=1)]: Done   8 out of   8 | elapsed:   13.0s finished
[I 2026-11-09 15:50:48,083] Trial 129 finished with value: 0.9365480129633736 and parameters: {'0:temperature_only': '0:True', '1:cities': "0:'all'", '2:n_selected_splines': 326, '3:n_components': 153, '4:alpha': 0.05368520820586315}. Best is trial 129 with value: 0.9365480129633736.
[I 2026-11-09 15:50:48,110] Using an existing study with name 'randomized_search' instead of creating a new one.

randomized_search_ridge.plot_results().update_layout(margin=dict(l=200))

We observe that the default values of the hyper-parameters are in the optimal region explored by the randomized search. This is a good sign that the model is well-specified and that the hyper-parameters are not too sensitive to small changes of those values.

We could further assess the stability of those optimal hyper-parameters by running a nested cross-validation, where we would perform a randomized search for each fold of the outer cross-validation loop as below but this is computationally expensive.

# nested_cv_results_ridge = skrub.cross_validate(
#      environment=predictions_ridge.skb.get_data(),
#      learner=randomized_search_ridge,
#      cv=TimeSeriesSplitter(),
#      scoring={
#          "r2": get_scorer("r2"),
#          "mape": make_scorer(mean_absolute_percentage_error),
#      },
#      n_jobs=-1,
#      return_learner=True,
#  ).round(3)

# nested_cv_results_ridge.round(3)