{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Lagged features for time series forecasting\n\nThis example demonstrates how Polars-engineered lagged features can be used\nfor time series forecasting with\n:class:`~sklearn.ensemble.HistGradientBoostingRegressor` on the Bike Sharing\nDemand dataset.\n\nSee the example on\n`sphx_glr_auto_examples_applications_plot_cyclical_feature_engineering.py`\nfor some data exploration on this dataset and a demo on periodic feature\nengineering.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Authors: The scikit-learn developers\n# SPDX-License-Identifier: BSD-3-Clause"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Analyzing the Bike Sharing Demand dataset\n\nWe start by loading the data from the OpenML repository as a raw parquet file\nto illustrate how to work with an arbitrary parquet file instead of hiding this\nstep in a convenience tool such as `sklearn.datasets.fetch_openml`.\n\nThe URL of the parquet file can be found in the JSON description of the\nBike Sharing Demand dataset with id 44063 on openml.org\n(https://openml.org/search?type=data&status=active&id=44063).\n\nThe `sha256` hash of the file is also provided to ensure the integrity of the\ndownloaded file.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport polars as pl\n\nfrom sklearn.datasets import fetch_file\n\npl.Config.set_fmt_str_lengths(20)\n\nbike_sharing_data_file = fetch_file(\n    \"https://data.openml.org/datasets/0004/44063/dataset_44063.pq\",\n    sha256=\"d120af76829af0d256338dc6dd4be5df4fd1f35bf3a283cab66a51c1c6abd06a\",\n)\nbike_sharing_data_file"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We load the parquet file with Polars for feature engineering. Polars\nautomatically caches common subexpressions which are reused in multiple\nexpressions (like `pl.col(\"count\").shift(1)` below). See\nhttps://docs.pola.rs/user-guide/lazy/optimizations/ for more information.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "df = pl.read_parquet(bike_sharing_data_file)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next, we take a look at the statistical summary of the dataset\nso that we can better understand the data that we are working with.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import polars.selectors as cs\n\nsummary = df.select(cs.numeric()).describe()\nsummary"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let us look at the count of the seasons `\"fall\"`, `\"spring\"`, `\"summer\"`\nand `\"winter\"` present in the dataset to confirm they are balanced.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\n\ndf[\"season\"].value_counts()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Generating Polars-engineered lagged features\nLet's consider the problem of predicting the demand at the\nnext hour given past demands. Since the demand is a continuous\nvariable, one could intuitively use any regression model. However, we do\nnot have the usual `(X_train, y_train)` dataset. Instead, we just have\nthe `y_train` demand data sequentially organized by time.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "lagged_df = df.select(\n    \"count\",\n    *[pl.col(\"count\").shift(i).alias(f\"lagged_count_{i}h\") for i in [1, 2, 3]],\n    lagged_count_1d=pl.col(\"count\").shift(24),\n    lagged_count_1d_1h=pl.col(\"count\").shift(24 + 1),\n    lagged_count_7d=pl.col(\"count\").shift(7 * 24),\n    lagged_count_7d_1h=pl.col(\"count\").shift(7 * 24 + 1),\n    lagged_mean_24h=pl.col(\"count\").shift(1).rolling_mean(24),\n    lagged_max_24h=pl.col(\"count\").shift(1).rolling_max(24),\n    lagged_min_24h=pl.col(\"count\").shift(1).rolling_min(24),\n    lagged_mean_7d=pl.col(\"count\").shift(1).rolling_mean(7 * 24),\n    lagged_max_7d=pl.col(\"count\").shift(1).rolling_max(7 * 24),\n    lagged_min_7d=pl.col(\"count\").shift(1).rolling_min(7 * 24),\n)\nlagged_df.tail(10)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Watch out however, the first lines have undefined values because their own\npast is unknown. This depends on how much lag we used:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "lagged_df.head(10)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can now separate the lagged features in a matrix `X` and the target variable\n(the counts to predict) in an array of the same first dimension `y`.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "lagged_df = lagged_df.drop_nulls()\nX = lagged_df.drop(\"count\")\ny = lagged_df[\"count\"]\nprint(\"X shape: {}\\ny shape: {}\".format(X.shape, y.shape))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Naive evaluation of the next hour bike demand regression\nLet's randomly split our tabularized dataset to train a gradient\nboosting regression tree (GBRT) model and evaluate it using Mean\nAbsolute Percentage Error (MAPE). If our model is aimed at forecasting\n(i.e., predicting future data from past data), we should not use training\ndata that are ulterior to the testing data. In time series machine learning\nthe \"i.i.d\" (independent and identically distributed) assumption does not\nhold true as the data points are not independent and have a temporal\nrelationship.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from sklearn.ensemble import HistGradientBoostingRegressor\nfrom sklearn.model_selection import train_test_split\n\nX_train, X_test, y_train, y_test = train_test_split(\n    X, y, test_size=0.2, random_state=42\n)\n\nmodel = HistGradientBoostingRegressor().fit(X_train, y_train)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Taking a look at the performance of the model.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from sklearn.metrics import mean_absolute_percentage_error\n\ny_pred = model.predict(X_test)\nmean_absolute_percentage_error(y_test, y_pred)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Proper next hour forecasting evaluation\nLet's use a proper evaluation splitting strategies that takes into account\nthe temporal structure of the dataset to evaluate our model's ability to\npredict data points in the future (to avoid cheating by reading values from\nthe lagged features in the training set).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from sklearn.model_selection import TimeSeriesSplit\n\nts_cv = TimeSeriesSplit(\n    n_splits=3,  # to keep the notebook fast enough on common laptops\n    gap=48,  # 2 days data gap between train and test\n    max_train_size=10000,  # keep train sets of comparable sizes\n    test_size=3000,  # for 2 or 3 digits of precision in scores\n)\nall_splits = list(ts_cv.split(X, y))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Training the model and evaluating its performance based on MAPE.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "train_idx, test_idx = all_splits[0]\nX_train, X_test = X[train_idx, :], X[test_idx, :]\ny_train, y_test = y[train_idx], y[test_idx]\n\nmodel = HistGradientBoostingRegressor().fit(X_train, y_train)\ny_pred = model.predict(X_test)\nmean_absolute_percentage_error(y_test, y_pred)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The generalization error measured via a shuffled trained test split\nis too optimistic. The generalization via a time-based split is likely to\nbe more representative of the true performance of the regression model.\nLet's assess this variability of our error evaluation with proper\ncross-validation:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from sklearn.model_selection import cross_val_score\n\ncv_mape_scores = -cross_val_score(\n    model, X, y, cv=ts_cv, scoring=\"neg_mean_absolute_percentage_error\"\n)\ncv_mape_scores"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The variability across splits is quite large! In a real life setting\nit would be advised to use more splits to better assess the variability.\nLet's report the mean CV scores and their standard deviation from now on.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(f\"CV MAPE: {cv_mape_scores.mean():.3f} \u00b1 {cv_mape_scores.std():.3f}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can compute several combinations of evaluation metrics and loss functions,\nwhich are reported a bit below.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from collections import defaultdict\n\nfrom sklearn.metrics import (\n    make_scorer,\n    mean_absolute_error,\n    mean_pinball_loss,\n    root_mean_squared_error,\n)\nfrom sklearn.model_selection import cross_validate\n\n\ndef consolidate_scores(cv_results, scores, metric):\n    if metric == \"MAPE\":\n        scores[metric].append(f\"{value.mean():.2f} \u00b1 {value.std():.2f}\")\n    else:\n        scores[metric].append(f\"{value.mean():.1f} \u00b1 {value.std():.1f}\")\n\n    return scores\n\n\nscoring = {\n    \"MAPE\": make_scorer(mean_absolute_percentage_error),\n    \"RMSE\": make_scorer(root_mean_squared_error),\n    \"MAE\": make_scorer(mean_absolute_error),\n    \"pinball_loss_05\": make_scorer(mean_pinball_loss, alpha=0.05),\n    \"pinball_loss_50\": make_scorer(mean_pinball_loss, alpha=0.50),\n    \"pinball_loss_95\": make_scorer(mean_pinball_loss, alpha=0.95),\n}\nloss_functions = [\"squared_error\", \"poisson\", \"absolute_error\"]\nscores = defaultdict(list)\nfor loss_func in loss_functions:\n    model = HistGradientBoostingRegressor(loss=loss_func)\n    cv_results = cross_validate(\n        model,\n        X,\n        y,\n        cv=ts_cv,\n        scoring=scoring,\n        n_jobs=2,\n    )\n    time = cv_results[\"fit_time\"]\n    scores[\"loss\"].append(loss_func)\n    scores[\"fit_time\"].append(f\"{time.mean():.2f} \u00b1 {time.std():.2f} s\")\n\n    for key, value in cv_results.items():\n        if key.startswith(\"test_\"):\n            metric = key.split(\"test_\")[1]\n            scores = consolidate_scores(cv_results, scores, metric)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Modeling predictive uncertainty via quantile regression\nInstead of modeling the expected value of the distribution of\n$Y|X$ like the least squares and Poisson losses do, one could try to\nestimate quantiles of the conditional distribution.\n\n$Y|X=x_i$ is expected to be a random variable for a given data point\n$x_i$ because we expect that the number of rentals cannot be 100%\naccurately predicted from the features. It can be influenced by other\nvariables not properly captured by the existing lagged features. For\ninstance whether or not it will rain in the next hour cannot be fully\nanticipated from the past hours bike rental data. This is what we\ncall aleatoric uncertainty.\n\nQuantile regression makes it possible to give a finer description of that\ndistribution without making strong assumptions on its shape.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "quantile_list = [0.05, 0.5, 0.95]\n\nfor quantile in quantile_list:\n    model = HistGradientBoostingRegressor(loss=\"quantile\", quantile=quantile)\n    cv_results = cross_validate(\n        model,\n        X,\n        y,\n        cv=ts_cv,\n        scoring=scoring,\n        n_jobs=2,\n    )\n    time = cv_results[\"fit_time\"]\n    scores[\"fit_time\"].append(f\"{time.mean():.2f} \u00b1 {time.std():.2f} s\")\n\n    scores[\"loss\"].append(f\"quantile {int(quantile * 100)}\")\n    for key, value in cv_results.items():\n        if key.startswith(\"test_\"):\n            metric = key.split(\"test_\")[1]\n            scores = consolidate_scores(cv_results, scores, metric)\n\nscores_df = pl.DataFrame(scores)\nscores_df"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let us take a look at the losses that minimise each metric.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def min_arg(col):\n    col_split = pl.col(col).str.split(\" \")\n    return pl.arg_sort_by(\n        col_split.list.get(0).cast(pl.Float64),\n        col_split.list.get(2).cast(pl.Float64),\n    ).first()\n\n\nscores_df.select(\n    pl.col(\"loss\").get(min_arg(col_name)).alias(col_name)\n    for col_name in scores_df.columns\n    if col_name != \"loss\"\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Even if the score distributions overlap due to the variance in the dataset,\nit is true that the average RMSE is lower when `loss=\"squared_error\"`, whereas\nthe average MAPE is lower when `loss=\"absolute_error\"` as expected. That is\nalso the case for the Mean Pinball Loss with the quantiles 5 and 95. The score\ncorresponding to the 50 quantile loss is overlapping with the score obtained\nby minimizing other loss functions, which is also the case for the MAE.\n\n## A qualitative look at the predictions\nWe can now visualize the performance of the model with regards\nto the 5th percentile, median and the 95th percentile:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "all_splits = list(ts_cv.split(X, y))\ntrain_idx, test_idx = all_splits[0]\n\nX_train, X_test = X[train_idx, :], X[test_idx, :]\ny_train, y_test = y[train_idx], y[test_idx]\n\nmax_iter = 50\ngbrt_mean_poisson = HistGradientBoostingRegressor(loss=\"poisson\", max_iter=max_iter)\ngbrt_mean_poisson.fit(X_train, y_train)\nmean_predictions = gbrt_mean_poisson.predict(X_test)\n\ngbrt_median = HistGradientBoostingRegressor(\n    loss=\"quantile\", quantile=0.5, max_iter=max_iter\n)\ngbrt_median.fit(X_train, y_train)\nmedian_predictions = gbrt_median.predict(X_test)\n\ngbrt_percentile_5 = HistGradientBoostingRegressor(\n    loss=\"quantile\", quantile=0.05, max_iter=max_iter\n)\ngbrt_percentile_5.fit(X_train, y_train)\npercentile_5_predictions = gbrt_percentile_5.predict(X_test)\n\ngbrt_percentile_95 = HistGradientBoostingRegressor(\n    loss=\"quantile\", quantile=0.95, max_iter=max_iter\n)\ngbrt_percentile_95.fit(X_train, y_train)\npercentile_95_predictions = gbrt_percentile_95.predict(X_test)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can now take a look at the predictions made by the regression models:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "last_hours = slice(-96, None)\nfig, ax = plt.subplots(figsize=(15, 7))\nplt.title(\"Predictions by regression models\")\nax.plot(\n    y_test[last_hours],\n    \"x-\",\n    alpha=0.2,\n    label=\"Actual demand\",\n    color=\"black\",\n)\nax.plot(\n    median_predictions[last_hours],\n    \"^-\",\n    label=\"GBRT median\",\n)\nax.plot(\n    mean_predictions[last_hours],\n    \"x-\",\n    label=\"GBRT mean (Poisson)\",\n)\nax.fill_between(\n    np.arange(96),\n    percentile_5_predictions[last_hours],\n    percentile_95_predictions[last_hours],\n    alpha=0.3,\n    label=\"GBRT 90% interval\",\n)\n_ = ax.legend()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Here it's interesting to notice that the blue area between the 5% and 95%\npercentile estimators has a width that varies with the time of the day:\n\n- At night, the blue band is much narrower: the pair of models is quite\n  certain that there will be a small number of bike rentals. And furthermore\n  these seem correct in the sense that the actual demand stays in that blue\n  band.\n- During the day, the blue band is much wider: the uncertainty grows, probably\n  because of the variability of the weather that can have a very large impact,\n  especially on week-ends.\n- We can also see that during week-days, the commute pattern is still visible in\n  the 5% and 95% estimations.\n- Finally, it is expected that 10% of the time, the actual demand does not lie\n  between the 5% and 95% percentile estimates. On this test span, the actual\n  demand seems to be higher, especially during the rush hours. It might reveal that\n  our 95% percentile estimator underestimates the demand peaks. This could be be\n  quantitatively confirmed by computing empirical coverage numbers as done in\n  the `calibration of confidence intervals <calibration-section>`.\n\nLooking at the performance of non-linear regression models vs\nthe best models:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from sklearn.metrics import PredictionErrorDisplay\n\nfig, axes = plt.subplots(ncols=3, figsize=(15, 6), sharey=True)\nfig.suptitle(\"Non-linear regression models\")\npredictions = [\n    median_predictions,\n    percentile_5_predictions,\n    percentile_95_predictions,\n]\nlabels = [\n    \"Median\",\n    \"5th percentile\",\n    \"95th percentile\",\n]\nfor ax, pred, label in zip(axes, predictions, labels):\n    PredictionErrorDisplay.from_predictions(\n        y_true=y_test,\n        y_pred=pred,\n        kind=\"residual_vs_predicted\",\n        scatter_kwargs={\"alpha\": 0.3},\n        ax=ax,\n    )\n    ax.set(xlabel=\"Predicted demand\", ylabel=\"True demand\")\n    ax.legend([\"Best model\", label])\n\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Conclusion\nThrough this example we explored time series forecasting using lagged\nfeatures. We compared a naive regression (using the standardized\n:class:`~sklearn.model_selection.train_test_split`) with a proper time\nseries evaluation strategy using\n:class:`~sklearn.model_selection.TimeSeriesSplit`. We observed that the\nmodel trained using :class:`~sklearn.model_selection.train_test_split`,\nhaving a default value of `shuffle` set to `True` produced an overly\noptimistic Mean Average Percentage Error (MAPE). The results\nproduced from the time-based split better represent the performance\nof our time-series regression model. We also analyzed the predictive uncertainty\nof our model via Quantile Regression. Predictions based on the 5th and\n95th percentile using `loss=\"quantile\"` provide us with a quantitative estimate\nof the uncertainty of the forecasts made by our time series regression model.\nUncertainty estimation can also be performed\nusing [MAPIE](https://mapie.readthedocs.io/en/latest/index.html),\nthat provides an implementation based on recent work on conformal prediction\nmethods and estimates both aleatoric and epistemic uncertainty at the same time.\nFurthermore, functionalities provided\nby [sktime](https://www.sktime.net/en/latest/users.html)\ncan be used to extend scikit-learn estimators by making use of recursive time\nseries forecasting, that enables dynamic predictions of future values.\n\n"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.11.14"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}