.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/applications/plot_time_series_lagged_features.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_examples_applications_plot_time_series_lagged_features.py>`
        to download the full example code or to run this example in your browser via JupyterLite or Binder

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_applications_plot_time_series_lagged_features.py:


===========================================
Lagged features for time series forecasting
===========================================

This example demonstrates how pandas-engineered lagged features can be used
for time series forecasting with
:class:`~sklearn.ensemble.HistGradientBoostingRegressor` on the Bike Sharing
Demand dataset.

See the example on
:ref:`sphx_glr_auto_examples_applications_plot_cyclical_feature_engineering.py`
for some data exploration on this dataset and a demo on periodic feature
engineering.

.. GENERATED FROM PYTHON SOURCE LINES 19-23

Analyzing the Bike Sharing Demand dataset
-----------------------------------------

We start by loading the data from the OpenML repository.

.. GENERATED FROM PYTHON SOURCE LINES 23-33

.. code-block:: Python

    import numpy as np
    import pandas as pd

    from sklearn.datasets import fetch_openml

    bike_sharing = fetch_openml(
        "Bike_Sharing_Demand", version=2, as_frame=True, parser="pandas"
    )
    df = bike_sharing.frame








.. GENERATED FROM PYTHON SOURCE LINES 34-36

Next, we take a look at the statistical summary of the dataset
so that we can better understand the data that we are working with.

.. GENERATED FROM PYTHON SOURCE LINES 36-45

.. code-block:: Python

    summary = pd.DataFrame(df.describe())
    summary = (
        summary.style.background_gradient()
        .set_table_attributes("style = 'display: inline'")
        .set_caption("Statistics of the Dataset")
        .set_table_styles([{"selector": "caption", "props": [("font-size", "16px")]}])
    )
    summary






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <style type="text/css">
    #T_2318b caption {
      font-size: 16px;
    }
    #T_2318b_row0_col0, #T_2318b_row0_col1, #T_2318b_row0_col2, #T_2318b_row0_col3, #T_2318b_row0_col4, #T_2318b_row0_col5, #T_2318b_row0_col6, #T_2318b_row0_col7, #T_2318b_row0_col8 {
      background-color: #023858;
      color: #f1f1f1;
    }
    #T_2318b_row1_col0, #T_2318b_row1_col1, #T_2318b_row1_col2, #T_2318b_row1_col3, #T_2318b_row1_col4, #T_2318b_row1_col5, #T_2318b_row1_col6, #T_2318b_row1_col7, #T_2318b_row2_col0, #T_2318b_row2_col1, #T_2318b_row2_col2, #T_2318b_row2_col3, #T_2318b_row2_col4, #T_2318b_row2_col5, #T_2318b_row2_col6, #T_2318b_row2_col7, #T_2318b_row3_col0, #T_2318b_row3_col1, #T_2318b_row3_col2, #T_2318b_row3_col3, #T_2318b_row3_col4, #T_2318b_row3_col5, #T_2318b_row3_col6, #T_2318b_row3_col7, #T_2318b_row3_col8, #T_2318b_row4_col0, #T_2318b_row4_col1, #T_2318b_row4_col2, #T_2318b_row4_col3, #T_2318b_row4_col4, #T_2318b_row4_col5, #T_2318b_row4_col6, #T_2318b_row4_col7, #T_2318b_row4_col8, #T_2318b_row5_col0, #T_2318b_row5_col1, #T_2318b_row5_col2, #T_2318b_row5_col3, #T_2318b_row5_col4, #T_2318b_row5_col5, #T_2318b_row5_col6, #T_2318b_row5_col7, #T_2318b_row6_col0, #T_2318b_row6_col1, #T_2318b_row6_col2, #T_2318b_row6_col3, #T_2318b_row6_col4, #T_2318b_row6_col5, #T_2318b_row6_col6, #T_2318b_row6_col7, #T_2318b_row7_col0, #T_2318b_row7_col1, #T_2318b_row7_col2, #T_2318b_row7_col3, #T_2318b_row7_col4, #T_2318b_row7_col5, #T_2318b_row7_col6, #T_2318b_row7_col7 {
      background-color: #fff7fb;
      color: #000000;
    }
    #T_2318b_row1_col8, #T_2318b_row2_col8, #T_2318b_row5_col8 {
      background-color: #fef6fa;
      color: #000000;
    }
    #T_2318b_row6_col8 {
      background-color: #fdf5fa;
      color: #000000;
    }
    #T_2318b_row7_col8 {
      background-color: #f7f0f7;
      color: #000000;
    }
    </style>
    <table id="T_2318b" style = 'display: inline'>
      <caption>Statistics of the Dataset</caption>
      <thead>
        <tr>
          <th class="blank level0" >&nbsp;</th>
          <th id="T_2318b_level0_col0" class="col_heading level0 col0" >year</th>
          <th id="T_2318b_level0_col1" class="col_heading level0 col1" >month</th>
          <th id="T_2318b_level0_col2" class="col_heading level0 col2" >hour</th>
          <th id="T_2318b_level0_col3" class="col_heading level0 col3" >weekday</th>
          <th id="T_2318b_level0_col4" class="col_heading level0 col4" >temp</th>
          <th id="T_2318b_level0_col5" class="col_heading level0 col5" >feel_temp</th>
          <th id="T_2318b_level0_col6" class="col_heading level0 col6" >humidity</th>
          <th id="T_2318b_level0_col7" class="col_heading level0 col7" >windspeed</th>
          <th id="T_2318b_level0_col8" class="col_heading level0 col8" >count</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th id="T_2318b_level0_row0" class="row_heading level0 row0" >count</th>
          <td id="T_2318b_row0_col0" class="data row0 col0" >17379.000000</td>
          <td id="T_2318b_row0_col1" class="data row0 col1" >17379.000000</td>
          <td id="T_2318b_row0_col2" class="data row0 col2" >17379.000000</td>
          <td id="T_2318b_row0_col3" class="data row0 col3" >17379.000000</td>
          <td id="T_2318b_row0_col4" class="data row0 col4" >17379.000000</td>
          <td id="T_2318b_row0_col5" class="data row0 col5" >17379.000000</td>
          <td id="T_2318b_row0_col6" class="data row0 col6" >17379.000000</td>
          <td id="T_2318b_row0_col7" class="data row0 col7" >17379.000000</td>
          <td id="T_2318b_row0_col8" class="data row0 col8" >17379.000000</td>
        </tr>
        <tr>
          <th id="T_2318b_level0_row1" class="row_heading level0 row1" >mean</th>
          <td id="T_2318b_row1_col0" class="data row1 col0" >0.502561</td>
          <td id="T_2318b_row1_col1" class="data row1 col1" >6.537775</td>
          <td id="T_2318b_row1_col2" class="data row1 col2" >11.546752</td>
          <td id="T_2318b_row1_col3" class="data row1 col3" >3.003683</td>
          <td id="T_2318b_row1_col4" class="data row1 col4" >20.376474</td>
          <td id="T_2318b_row1_col5" class="data row1 col5" >23.788755</td>
          <td id="T_2318b_row1_col6" class="data row1 col6" >0.627229</td>
          <td id="T_2318b_row1_col7" class="data row1 col7" >12.736540</td>
          <td id="T_2318b_row1_col8" class="data row1 col8" >189.463088</td>
        </tr>
        <tr>
          <th id="T_2318b_level0_row2" class="row_heading level0 row2" >std</th>
          <td id="T_2318b_row2_col0" class="data row2 col0" >0.500008</td>
          <td id="T_2318b_row2_col1" class="data row2 col1" >3.438776</td>
          <td id="T_2318b_row2_col2" class="data row2 col2" >6.914405</td>
          <td id="T_2318b_row2_col3" class="data row2 col3" >2.005771</td>
          <td id="T_2318b_row2_col4" class="data row2 col4" >7.894801</td>
          <td id="T_2318b_row2_col5" class="data row2 col5" >8.592511</td>
          <td id="T_2318b_row2_col6" class="data row2 col6" >0.192930</td>
          <td id="T_2318b_row2_col7" class="data row2 col7" >8.196795</td>
          <td id="T_2318b_row2_col8" class="data row2 col8" >181.387599</td>
        </tr>
        <tr>
          <th id="T_2318b_level0_row3" class="row_heading level0 row3" >min</th>
          <td id="T_2318b_row3_col0" class="data row3 col0" >0.000000</td>
          <td id="T_2318b_row3_col1" class="data row3 col1" >1.000000</td>
          <td id="T_2318b_row3_col2" class="data row3 col2" >0.000000</td>
          <td id="T_2318b_row3_col3" class="data row3 col3" >0.000000</td>
          <td id="T_2318b_row3_col4" class="data row3 col4" >0.820000</td>
          <td id="T_2318b_row3_col5" class="data row3 col5" >0.000000</td>
          <td id="T_2318b_row3_col6" class="data row3 col6" >0.000000</td>
          <td id="T_2318b_row3_col7" class="data row3 col7" >0.000000</td>
          <td id="T_2318b_row3_col8" class="data row3 col8" >1.000000</td>
        </tr>
        <tr>
          <th id="T_2318b_level0_row4" class="row_heading level0 row4" >25%</th>
          <td id="T_2318b_row4_col0" class="data row4 col0" >0.000000</td>
          <td id="T_2318b_row4_col1" class="data row4 col1" >4.000000</td>
          <td id="T_2318b_row4_col2" class="data row4 col2" >6.000000</td>
          <td id="T_2318b_row4_col3" class="data row4 col3" >1.000000</td>
          <td id="T_2318b_row4_col4" class="data row4 col4" >13.940000</td>
          <td id="T_2318b_row4_col5" class="data row4 col5" >16.665000</td>
          <td id="T_2318b_row4_col6" class="data row4 col6" >0.480000</td>
          <td id="T_2318b_row4_col7" class="data row4 col7" >7.001500</td>
          <td id="T_2318b_row4_col8" class="data row4 col8" >40.000000</td>
        </tr>
        <tr>
          <th id="T_2318b_level0_row5" class="row_heading level0 row5" >50%</th>
          <td id="T_2318b_row5_col0" class="data row5 col0" >1.000000</td>
          <td id="T_2318b_row5_col1" class="data row5 col1" >7.000000</td>
          <td id="T_2318b_row5_col2" class="data row5 col2" >12.000000</td>
          <td id="T_2318b_row5_col3" class="data row5 col3" >3.000000</td>
          <td id="T_2318b_row5_col4" class="data row5 col4" >20.500000</td>
          <td id="T_2318b_row5_col5" class="data row5 col5" >24.240000</td>
          <td id="T_2318b_row5_col6" class="data row5 col6" >0.630000</td>
          <td id="T_2318b_row5_col7" class="data row5 col7" >12.998000</td>
          <td id="T_2318b_row5_col8" class="data row5 col8" >142.000000</td>
        </tr>
        <tr>
          <th id="T_2318b_level0_row6" class="row_heading level0 row6" >75%</th>
          <td id="T_2318b_row6_col0" class="data row6 col0" >1.000000</td>
          <td id="T_2318b_row6_col1" class="data row6 col1" >10.000000</td>
          <td id="T_2318b_row6_col2" class="data row6 col2" >18.000000</td>
          <td id="T_2318b_row6_col3" class="data row6 col3" >5.000000</td>
          <td id="T_2318b_row6_col4" class="data row6 col4" >27.060000</td>
          <td id="T_2318b_row6_col5" class="data row6 col5" >31.060000</td>
          <td id="T_2318b_row6_col6" class="data row6 col6" >0.780000</td>
          <td id="T_2318b_row6_col7" class="data row6 col7" >16.997900</td>
          <td id="T_2318b_row6_col8" class="data row6 col8" >281.000000</td>
        </tr>
        <tr>
          <th id="T_2318b_level0_row7" class="row_heading level0 row7" >max</th>
          <td id="T_2318b_row7_col0" class="data row7 col0" >1.000000</td>
          <td id="T_2318b_row7_col1" class="data row7 col1" >12.000000</td>
          <td id="T_2318b_row7_col2" class="data row7 col2" >23.000000</td>
          <td id="T_2318b_row7_col3" class="data row7 col3" >6.000000</td>
          <td id="T_2318b_row7_col4" class="data row7 col4" >41.000000</td>
          <td id="T_2318b_row7_col5" class="data row7 col5" >50.000000</td>
          <td id="T_2318b_row7_col6" class="data row7 col6" >1.000000</td>
          <td id="T_2318b_row7_col7" class="data row7 col7" >56.996900</td>
          <td id="T_2318b_row7_col8" class="data row7 col8" >977.000000</td>
        </tr>
      </tbody>
    </table>

    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 46-48

Let us look at the count of the seasons `"fall"`, `"spring"`, `"summer"`
and `"winter"` present in the dataset to confirm they are balanced.

.. GENERATED FROM PYTHON SOURCE LINES 48-54

.. code-block:: Python


    import matplotlib.pyplot as plt

    df["season"].value_counts()






.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    season
    fall      4496
    summer    4409
    spring    4242
    winter    4232
    Name: count, dtype: int64



.. GENERATED FROM PYTHON SOURCE LINES 55-62

Generating pandas-engineered lagged features
--------------------------------------------
Let's consider the problem of predicting the demand at the
next hour given past demands. Since the demand is a continuous
variable, one could intuitively use any regression model. However, we do
not have the usual `(X_train, y_train)` dataset. Instead, we just have
the `y_train` demand data sequentially organized by time.

.. GENERATED FROM PYTHON SOURCE LINES 62-84

.. code-block:: Python

    count = df["count"]
    lagged_df = pd.concat(
        [
            count,
            count.shift(1).rename("lagged_count_1h"),
            count.shift(2).rename("lagged_count_2h"),
            count.shift(3).rename("lagged_count_3h"),
            count.shift(24).rename("lagged_count_1d"),
            count.shift(24 + 1).rename("lagged_count_1d_1h"),
            count.shift(7 * 24).rename("lagged_count_7d"),
            count.shift(7 * 24 + 1).rename("lagged_count_7d_1h"),
            count.shift(1).rolling(24).mean().rename("lagged_mean_24h"),
            count.shift(1).rolling(24).max().rename("lagged_max_24h"),
            count.shift(1).rolling(24).min().rename("lagged_min_24h"),
            count.shift(1).rolling(7 * 24).mean().rename("lagged_mean_7d"),
            count.shift(1).rolling(7 * 24).max().rename("lagged_max_7d"),
            count.shift(1).rolling(7 * 24).min().rename("lagged_min_7d"),
        ],
        axis="columns",
    )
    lagged_df.tail(10)






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>count</th>
          <th>lagged_count_1h</th>
          <th>lagged_count_2h</th>
          <th>lagged_count_3h</th>
          <th>lagged_count_1d</th>
          <th>lagged_count_1d_1h</th>
          <th>lagged_count_7d</th>
          <th>lagged_count_7d_1h</th>
          <th>lagged_mean_24h</th>
          <th>lagged_max_24h</th>
          <th>lagged_min_24h</th>
          <th>lagged_mean_7d</th>
          <th>lagged_max_7d</th>
          <th>lagged_min_7d</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>17369</th>
          <td>247</td>
          <td>203.0</td>
          <td>224.0</td>
          <td>157.0</td>
          <td>160.0</td>
          <td>169.0</td>
          <td>70.0</td>
          <td>135.0</td>
          <td>93.500000</td>
          <td>224.0</td>
          <td>1.0</td>
          <td>67.732143</td>
          <td>271.0</td>
          <td>1.0</td>
        </tr>
        <tr>
          <th>17370</th>
          <td>315</td>
          <td>247.0</td>
          <td>203.0</td>
          <td>224.0</td>
          <td>138.0</td>
          <td>160.0</td>
          <td>46.0</td>
          <td>70.0</td>
          <td>97.125000</td>
          <td>247.0</td>
          <td>1.0</td>
          <td>68.785714</td>
          <td>271.0</td>
          <td>1.0</td>
        </tr>
        <tr>
          <th>17371</th>
          <td>214</td>
          <td>315.0</td>
          <td>247.0</td>
          <td>203.0</td>
          <td>133.0</td>
          <td>138.0</td>
          <td>33.0</td>
          <td>46.0</td>
          <td>104.500000</td>
          <td>315.0</td>
          <td>1.0</td>
          <td>70.386905</td>
          <td>315.0</td>
          <td>1.0</td>
        </tr>
        <tr>
          <th>17372</th>
          <td>164</td>
          <td>214.0</td>
          <td>315.0</td>
          <td>247.0</td>
          <td>123.0</td>
          <td>133.0</td>
          <td>33.0</td>
          <td>33.0</td>
          <td>107.875000</td>
          <td>315.0</td>
          <td>1.0</td>
          <td>71.464286</td>
          <td>315.0</td>
          <td>1.0</td>
        </tr>
        <tr>
          <th>17373</th>
          <td>122</td>
          <td>164.0</td>
          <td>214.0</td>
          <td>315.0</td>
          <td>125.0</td>
          <td>123.0</td>
          <td>26.0</td>
          <td>33.0</td>
          <td>109.583333</td>
          <td>315.0</td>
          <td>1.0</td>
          <td>72.244048</td>
          <td>315.0</td>
          <td>1.0</td>
        </tr>
        <tr>
          <th>17374</th>
          <td>119</td>
          <td>122.0</td>
          <td>164.0</td>
          <td>214.0</td>
          <td>102.0</td>
          <td>125.0</td>
          <td>26.0</td>
          <td>26.0</td>
          <td>109.458333</td>
          <td>315.0</td>
          <td>1.0</td>
          <td>72.815476</td>
          <td>315.0</td>
          <td>1.0</td>
        </tr>
        <tr>
          <th>17375</th>
          <td>89</td>
          <td>119.0</td>
          <td>122.0</td>
          <td>164.0</td>
          <td>72.0</td>
          <td>102.0</td>
          <td>18.0</td>
          <td>26.0</td>
          <td>110.166667</td>
          <td>315.0</td>
          <td>1.0</td>
          <td>73.369048</td>
          <td>315.0</td>
          <td>1.0</td>
        </tr>
        <tr>
          <th>17376</th>
          <td>90</td>
          <td>89.0</td>
          <td>119.0</td>
          <td>122.0</td>
          <td>47.0</td>
          <td>72.0</td>
          <td>23.0</td>
          <td>18.0</td>
          <td>110.875000</td>
          <td>315.0</td>
          <td>1.0</td>
          <td>73.791667</td>
          <td>315.0</td>
          <td>1.0</td>
        </tr>
        <tr>
          <th>17377</th>
          <td>61</td>
          <td>90.0</td>
          <td>89.0</td>
          <td>119.0</td>
          <td>36.0</td>
          <td>47.0</td>
          <td>22.0</td>
          <td>23.0</td>
          <td>112.666667</td>
          <td>315.0</td>
          <td>1.0</td>
          <td>74.190476</td>
          <td>315.0</td>
          <td>1.0</td>
        </tr>
        <tr>
          <th>17378</th>
          <td>49</td>
          <td>61.0</td>
          <td>90.0</td>
          <td>89.0</td>
          <td>49.0</td>
          <td>36.0</td>
          <td>12.0</td>
          <td>22.0</td>
          <td>113.708333</td>
          <td>315.0</td>
          <td>1.0</td>
          <td>74.422619</td>
          <td>315.0</td>
          <td>1.0</td>
        </tr>
      </tbody>
    </table>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 85-87

Watch out however, the first lines have undefined values because their own
past is unknown. This depends on how much lag we used:

.. GENERATED FROM PYTHON SOURCE LINES 87-89

.. code-block:: Python

    lagged_df.head(10)






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>count</th>
          <th>lagged_count_1h</th>
          <th>lagged_count_2h</th>
          <th>lagged_count_3h</th>
          <th>lagged_count_1d</th>
          <th>lagged_count_1d_1h</th>
          <th>lagged_count_7d</th>
          <th>lagged_count_7d_1h</th>
          <th>lagged_mean_24h</th>
          <th>lagged_max_24h</th>
          <th>lagged_min_24h</th>
          <th>lagged_mean_7d</th>
          <th>lagged_max_7d</th>
          <th>lagged_min_7d</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>16</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>1</th>
          <td>40</td>
          <td>16.0</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>2</th>
          <td>32</td>
          <td>40.0</td>
          <td>16.0</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>3</th>
          <td>13</td>
          <td>32.0</td>
          <td>40.0</td>
          <td>16.0</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>4</th>
          <td>1</td>
          <td>13.0</td>
          <td>32.0</td>
          <td>40.0</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>5</th>
          <td>1</td>
          <td>1.0</td>
          <td>13.0</td>
          <td>32.0</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>6</th>
          <td>2</td>
          <td>1.0</td>
          <td>1.0</td>
          <td>13.0</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>7</th>
          <td>3</td>
          <td>2.0</td>
          <td>1.0</td>
          <td>1.0</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>8</th>
          <td>8</td>
          <td>3.0</td>
          <td>2.0</td>
          <td>1.0</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>9</th>
          <td>14</td>
          <td>8.0</td>
          <td>3.0</td>
          <td>2.0</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
        </tr>
      </tbody>
    </table>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 90-92

We can now separate the lagged features in a matrix `X` and the target variable
(the counts to predict) in an array of the same first dimension `y`.

.. GENERATED FROM PYTHON SOURCE LINES 92-97

.. code-block:: Python

    lagged_df = lagged_df.dropna()
    X = lagged_df.drop("count", axis="columns")
    y = lagged_df["count"]
    print("X shape: {}\ny shape: {}".format(X.shape, y.shape))





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    X shape: (17210, 13)
    y shape: (17210,)




.. GENERATED FROM PYTHON SOURCE LINES 98-108

Naive evaluation of the next hour bike demand regression
--------------------------------------------------------
Let's randomly split our tabularized dataset to train a gradient
boosting regression tree (GBRT) model and evaluate it using Mean
Absolute Percentage Error (MAPE). If our model is aimed at forecasting
(i.e., predicting future data from past data), we should not use training
data that are ulterior to the testing data. In time series machine learning
the "i.i.d" (independent and identically distributed) assumption does not
hold true as the data points are not independent and have a temporal
relationship.

.. GENERATED FROM PYTHON SOURCE LINES 108-117

.. code-block:: Python

    from sklearn.ensemble import HistGradientBoostingRegressor
    from sklearn.model_selection import train_test_split

    X_train, X_test, y_train, y_test = train_test_split(
        X, y, test_size=0.2, random_state=42
    )

    model = HistGradientBoostingRegressor().fit(X_train, y_train)








.. GENERATED FROM PYTHON SOURCE LINES 118-119

Taking a look at the performance of the model.

.. GENERATED FROM PYTHON SOURCE LINES 119-124

.. code-block:: Python

    from sklearn.metrics import mean_absolute_percentage_error

    y_pred = model.predict(X_test)
    mean_absolute_percentage_error(y_test, y_pred)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    0.3889873516666431



.. GENERATED FROM PYTHON SOURCE LINES 125-131

Proper next hour forecasting evaluation
---------------------------------------
Let's use a proper evaluation splitting strategies that takes into account
the temporal structure of the dataset to evaluate our model's ability to
predict data points in the future (to avoid cheating by reading values from
the lagged features in the training set).

.. GENERATED FROM PYTHON SOURCE LINES 131-141

.. code-block:: Python

    from sklearn.model_selection import TimeSeriesSplit

    ts_cv = TimeSeriesSplit(
        n_splits=3,  # to keep the notebook fast enough on common laptops
        gap=48,  # 2 days data gap between train and test
        max_train_size=10000,  # keep train sets of comparable sizes
        test_size=3000,  # for 2 or 3 digits of precision in scores
    )
    all_splits = list(ts_cv.split(X, y))








.. GENERATED FROM PYTHON SOURCE LINES 142-143

Training the model and evaluating its performance based on MAPE.

.. GENERATED FROM PYTHON SOURCE LINES 143-151

.. code-block:: Python

    train_idx, test_idx = all_splits[0]
    X_train, X_test = X.iloc[train_idx], X.iloc[test_idx]
    y_train, y_test = y.iloc[train_idx], y.iloc[test_idx]

    model = HistGradientBoostingRegressor().fit(X_train, y_train)
    y_pred = model.predict(X_test)
    mean_absolute_percentage_error(y_test, y_pred)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    0.44300751539296973



.. GENERATED FROM PYTHON SOURCE LINES 152-157

The generalization error measured via a shuffled trained test split
is too optimistic. The generalization via a time-based split is likely to
be more representative of the true performance of the regression model.
Let's assess this variability of our error evaluation with proper
cross-validation:

.. GENERATED FROM PYTHON SOURCE LINES 157-164

.. code-block:: Python

    from sklearn.model_selection import cross_val_score

    cv_mape_scores = -cross_val_score(
        model, X, y, cv=ts_cv, scoring="neg_mean_absolute_percentage_error"
    )
    cv_mape_scores





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    array([0.44300752, 0.27772182, 0.3697178 ])



.. GENERATED FROM PYTHON SOURCE LINES 165-168

The variability across splits is quite large! In a real life setting
it would be advised to use more splits to better assess the variability.
Let's report the mean CV scores and their standard deviation from now on.

.. GENERATED FROM PYTHON SOURCE LINES 168-170

.. code-block:: Python

    print(f"CV MAPE: {cv_mape_scores.mean():.3f} ± {cv_mape_scores.std():.3f}")





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    CV MAPE: 0.363 ± 0.068




.. GENERATED FROM PYTHON SOURCE LINES 171-173

We can compute several combinations of evaluation metrics and loss functions,
which are reported a bit below.

.. GENERATED FROM PYTHON SOURCE LINES 173-223

.. code-block:: Python

    from collections import defaultdict

    from sklearn.metrics import (
        make_scorer,
        mean_absolute_error,
        mean_pinball_loss,
        root_mean_squared_error,
    )
    from sklearn.model_selection import cross_validate


    def consolidate_scores(cv_results, scores, metric):
        if metric == "MAPE":
            scores[metric].append(f"{value.mean():.2f} ± {value.std():.2f}")
        else:
            scores[metric].append(f"{value.mean():.1f} ± {value.std():.1f}")

        return scores


    scoring = {
        "MAPE": make_scorer(mean_absolute_percentage_error),
        "RMSE": make_scorer(root_mean_squared_error),
        "MAE": make_scorer(mean_absolute_error),
        "pinball_loss_05": make_scorer(mean_pinball_loss, alpha=0.05),
        "pinball_loss_50": make_scorer(mean_pinball_loss, alpha=0.50),
        "pinball_loss_95": make_scorer(mean_pinball_loss, alpha=0.95),
    }
    loss_functions = ["squared_error", "poisson", "absolute_error"]
    scores = defaultdict(list)
    for loss_func in loss_functions:
        model = HistGradientBoostingRegressor(loss=loss_func)
        cv_results = cross_validate(
            model,
            X,
            y,
            cv=ts_cv,
            scoring=scoring,
            n_jobs=2,
        )
        time = cv_results["fit_time"]
        scores["loss"].append(loss_func)
        scores["fit_time"].append(f"{time.mean():.2f} ± {time.std():.2f} s")

        for key, value in cv_results.items():
            if key.startswith("test_"):
                metric = key.split("test_")[1]
                scores = consolidate_scores(cv_results, scores, metric)









.. GENERATED FROM PYTHON SOURCE LINES 224-240

Modeling predictive uncertainty via quantile regression
-------------------------------------------------------
Instead of modeling the expected value of the distribution of
:math:`Y|X` like the least squares and Poisson losses do, one could try to
estimate quantiles of the conditional distribution.

:math:`Y|X=x_i` is expected to be a random variable for a given data point
:math:`x_i` because we expect that the number of rentals cannot be 100%
accurately predicted from the features. It can be influenced by other
variables not properly captured by the existing lagged features. For
instance whether or not it will rain in the next hour cannot be fully
anticipated from the past hours bike rental data. This is what we
call aleatoric uncertainty.

Quantile regression makes it possible to give a finer description of that
distribution without making strong assumptions on its shape.

.. GENERATED FROM PYTHON SOURCE LINES 240-292

.. code-block:: Python

    quantile_list = [0.05, 0.5, 0.95]

    for quantile in quantile_list:
        model = HistGradientBoostingRegressor(loss="quantile", quantile=quantile)
        cv_results = cross_validate(
            model,
            X,
            y,
            cv=ts_cv,
            scoring=scoring,
            n_jobs=2,
        )
        time = cv_results["fit_time"]
        scores["fit_time"].append(f"{time.mean():.2f} ± {time.std():.2f} s")

        scores["loss"].append(f"quantile {int(quantile*100)}")
        for key, value in cv_results.items():
            if key.startswith("test_"):
                metric = key.split("test_")[1]
                scores = consolidate_scores(cv_results, scores, metric)

    df = pd.DataFrame(scores)

    styled_df_copy = df.copy()


    def extract_numeric(value):
        parts = value.split("±")
        mean_value = float(parts[0])
        std_value = float(parts[1].split()[0])

        return mean_value, std_value


    # Convert columns containing "±" to tuples of numerical values
    cols_to_convert = df.columns[1:]  # Exclude the "loss" column
    for col in cols_to_convert:
        df[col] = df[col].apply(extract_numeric)

    min_values = df.min()

    # Create a mask for highlighting minimum values
    mask = pd.DataFrame("", index=df.index, columns=df.columns)
    for col in cols_to_convert:
        mask[col] = df[col].apply(
            lambda x: "font-weight: bold" if x == min_values[col] else ""
        )

    styled_df_copy = styled_df_copy.style.apply(lambda x: mask, axis=None)
    styled_df_copy







.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <style type="text/css">
    #T_556da_row0_col1, #T_556da_row0_col3, #T_556da_row0_col4, #T_556da_row0_col6, #T_556da_row2_col2, #T_556da_row3_col5, #T_556da_row4_col2, #T_556da_row5_col7 {
      font-weight: bold;
    }
    </style>
    <table id="T_556da">
      <thead>
        <tr>
          <th class="blank level0" >&nbsp;</th>
          <th id="T_556da_level0_col0" class="col_heading level0 col0" >loss</th>
          <th id="T_556da_level0_col1" class="col_heading level0 col1" >fit_time</th>
          <th id="T_556da_level0_col2" class="col_heading level0 col2" >MAPE</th>
          <th id="T_556da_level0_col3" class="col_heading level0 col3" >RMSE</th>
          <th id="T_556da_level0_col4" class="col_heading level0 col4" >MAE</th>
          <th id="T_556da_level0_col5" class="col_heading level0 col5" >pinball_loss_05</th>
          <th id="T_556da_level0_col6" class="col_heading level0 col6" >pinball_loss_50</th>
          <th id="T_556da_level0_col7" class="col_heading level0 col7" >pinball_loss_95</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th id="T_556da_level0_row0" class="row_heading level0 row0" >0</th>
          <td id="T_556da_row0_col0" class="data row0 col0" >squared_error</td>
          <td id="T_556da_row0_col1" class="data row0 col1" >0.31 ± 0.01 s</td>
          <td id="T_556da_row0_col2" class="data row0 col2" >0.36 ± 0.07</td>
          <td id="T_556da_row0_col3" class="data row0 col3" >62.3 ± 3.5</td>
          <td id="T_556da_row0_col4" class="data row0 col4" >39.1 ± 2.3</td>
          <td id="T_556da_row0_col5" class="data row0 col5" >17.7 ± 1.3</td>
          <td id="T_556da_row0_col6" class="data row0 col6" >19.5 ± 1.1</td>
          <td id="T_556da_row0_col7" class="data row0 col7" >21.4 ± 2.4</td>
        </tr>
        <tr>
          <th id="T_556da_level0_row1" class="row_heading level0 row1" >1</th>
          <td id="T_556da_row1_col0" class="data row1 col0" >poisson</td>
          <td id="T_556da_row1_col1" class="data row1 col1" >0.33 ± 0.01 s</td>
          <td id="T_556da_row1_col2" class="data row1 col2" >0.32 ± 0.07</td>
          <td id="T_556da_row1_col3" class="data row1 col3" >64.2 ± 4.0</td>
          <td id="T_556da_row1_col4" class="data row1 col4" >39.3 ± 2.8</td>
          <td id="T_556da_row1_col5" class="data row1 col5" >16.7 ± 1.5</td>
          <td id="T_556da_row1_col6" class="data row1 col6" >19.7 ± 1.4</td>
          <td id="T_556da_row1_col7" class="data row1 col7" >22.6 ± 3.0</td>
        </tr>
        <tr>
          <th id="T_556da_level0_row2" class="row_heading level0 row2" >2</th>
          <td id="T_556da_row2_col0" class="data row2 col0" >absolute_error</td>
          <td id="T_556da_row2_col1" class="data row2 col1" >0.42 ± 0.01 s</td>
          <td id="T_556da_row2_col2" class="data row2 col2" >0.32 ± 0.06</td>
          <td id="T_556da_row2_col3" class="data row2 col3" >64.6 ± 3.8</td>
          <td id="T_556da_row2_col4" class="data row2 col4" >39.9 ± 3.2</td>
          <td id="T_556da_row2_col5" class="data row2 col5" >17.1 ± 1.1</td>
          <td id="T_556da_row2_col6" class="data row2 col6" >19.9 ± 1.6</td>
          <td id="T_556da_row2_col7" class="data row2 col7" >22.7 ± 3.1</td>
        </tr>
        <tr>
          <th id="T_556da_level0_row3" class="row_heading level0 row3" >3</th>
          <td id="T_556da_row3_col0" class="data row3 col0" >quantile 5</td>
          <td id="T_556da_row3_col1" class="data row3 col1" >0.56 ± 0.01 s</td>
          <td id="T_556da_row3_col2" class="data row3 col2" >0.41 ± 0.01</td>
          <td id="T_556da_row3_col3" class="data row3 col3" >145.6 ± 20.9</td>
          <td id="T_556da_row3_col4" class="data row3 col4" >92.5 ± 16.2</td>
          <td id="T_556da_row3_col5" class="data row3 col5" >5.9 ± 0.9</td>
          <td id="T_556da_row3_col6" class="data row3 col6" >46.2 ± 8.1</td>
          <td id="T_556da_row3_col7" class="data row3 col7" >86.6 ± 15.3</td>
        </tr>
        <tr>
          <th id="T_556da_level0_row4" class="row_heading level0 row4" >4</th>
          <td id="T_556da_row4_col0" class="data row4 col0" >quantile 50</td>
          <td id="T_556da_row4_col1" class="data row4 col1" >0.60 ± 0.00 s</td>
          <td id="T_556da_row4_col2" class="data row4 col2" >0.32 ± 0.06</td>
          <td id="T_556da_row4_col3" class="data row4 col3" >64.6 ± 3.8</td>
          <td id="T_556da_row4_col4" class="data row4 col4" >39.9 ± 3.2</td>
          <td id="T_556da_row4_col5" class="data row4 col5" >17.1 ± 1.1</td>
          <td id="T_556da_row4_col6" class="data row4 col6" >19.9 ± 1.6</td>
          <td id="T_556da_row4_col7" class="data row4 col7" >22.7 ± 3.1</td>
        </tr>
        <tr>
          <th id="T_556da_level0_row5" class="row_heading level0 row5" >5</th>
          <td id="T_556da_row5_col0" class="data row5 col0" >quantile 95</td>
          <td id="T_556da_row5_col1" class="data row5 col1" >0.68 ± 0.07 s</td>
          <td id="T_556da_row5_col2" class="data row5 col2" >1.07 ± 0.27</td>
          <td id="T_556da_row5_col3" class="data row5 col3" >99.6 ± 8.7</td>
          <td id="T_556da_row5_col4" class="data row5 col4" >72.0 ± 6.1</td>
          <td id="T_556da_row5_col5" class="data row5 col5" >62.9 ± 7.4</td>
          <td id="T_556da_row5_col6" class="data row5 col6" >36.0 ± 3.1</td>
          <td id="T_556da_row5_col7" class="data row5 col7" >9.1 ± 1.3</td>
        </tr>
      </tbody>
    </table>

    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 293-304

Even if the score distributions overlap due to the variance in the dataset, it
is true that the average RMSE is lower when `loss="squared_error"`, whereas
the average MAPE is lower when `loss="absolute_error"` as expected. That is
also the case for the Mean Pinball Loss with the quantiles 5 and 95. The score
corresponding to the 50 quantile loss is overlapping with the score obtained
by minimizing other loss functions, which is also the case for the MAE.

A qualitative look at the predictions
-------------------------------------
We can now visualize the performance of the model with regards
to the 5th percentile, median and the 95th percentile:

.. GENERATED FROM PYTHON SOURCE LINES 304-333

.. code-block:: Python

    all_splits = list(ts_cv.split(X, y))
    train_idx, test_idx = all_splits[0]

    X_train, X_test = X.iloc[train_idx], X.iloc[test_idx]
    y_train, y_test = y.iloc[train_idx], y.iloc[test_idx]

    max_iter = 50
    gbrt_mean_poisson = HistGradientBoostingRegressor(loss="poisson", max_iter=max_iter)
    gbrt_mean_poisson.fit(X_train, y_train)
    mean_predictions = gbrt_mean_poisson.predict(X_test)

    gbrt_median = HistGradientBoostingRegressor(
        loss="quantile", quantile=0.5, max_iter=max_iter
    )
    gbrt_median.fit(X_train, y_train)
    median_predictions = gbrt_median.predict(X_test)

    gbrt_percentile_5 = HistGradientBoostingRegressor(
        loss="quantile", quantile=0.05, max_iter=max_iter
    )
    gbrt_percentile_5.fit(X_train, y_train)
    percentile_5_predictions = gbrt_percentile_5.predict(X_test)

    gbrt_percentile_95 = HistGradientBoostingRegressor(
        loss="quantile", quantile=0.95, max_iter=max_iter
    )
    gbrt_percentile_95.fit(X_train, y_train)
    percentile_95_predictions = gbrt_percentile_95.predict(X_test)








.. GENERATED FROM PYTHON SOURCE LINES 334-335

We can now take a look at the predictions made by the regression models:

.. GENERATED FROM PYTHON SOURCE LINES 335-364

.. code-block:: Python

    last_hours = slice(-96, None)
    fig, ax = plt.subplots(figsize=(15, 7))
    plt.title("Predictions by regression models")
    ax.plot(
        y_test.values[last_hours],
        "x-",
        alpha=0.2,
        label="Actual demand",
        color="black",
    )
    ax.plot(
        median_predictions[last_hours],
        "^-",
        label="GBRT median",
    )
    ax.plot(
        mean_predictions[last_hours],
        "x-",
        label="GBRT mean (Poisson)",
    )
    ax.fill_between(
        np.arange(96),
        percentile_5_predictions[last_hours],
        percentile_95_predictions[last_hours],
        alpha=0.3,
        label="GBRT 90% interval",
    )
    _ = ax.legend()




.. image-sg:: /auto_examples/applications/images/sphx_glr_plot_time_series_lagged_features_001.png
   :alt: Predictions by regression models
   :srcset: /auto_examples/applications/images/sphx_glr_plot_time_series_lagged_features_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 365-386

Here it's interesting to notice that the blue area between the 5% and 95%
percentile estimators has a width that varies with the time of the day:

- At night, the blue band is much narrower: the pair of models is quite
  certain that there will be a small number of bike rentals. And furthermore
  these seem correct in the sense that the actual demand stays in that blue
  band.
- During the day, the blue band is much wider: the uncertainty grows, probably
  because of the variability of the weather that can have a very large impact,
  especially on week-ends.
- We can also see that during week-days, the commute pattern is still visible in
  the 5% and 95% estimations.
- Finally, it is expected that 10% of the time, the actual demand does not lie
  between the 5% and 95% percentile estimates. On this test span, the actual
  demand seems to be higher, especially during the rush hours. It might reveal that
  our 95% percentile estimator underestimates the demand peaks. This could be be
  quantitatively confirmed by computing empirical coverage numbers as done in
  the :ref:`calibration of confidence intervals <calibration-section>`.

Looking at the performance of non-linear regression models vs
the best models:

.. GENERATED FROM PYTHON SOURCE LINES 386-413

.. code-block:: Python

    from sklearn.metrics import PredictionErrorDisplay

    fig, axes = plt.subplots(ncols=3, figsize=(15, 6), sharey=True)
    fig.suptitle("Non-linear regression models")
    predictions = [
        median_predictions,
        percentile_5_predictions,
        percentile_95_predictions,
    ]
    labels = [
        "Median",
        "5th percentile",
        "95th percentile",
    ]
    for ax, pred, label in zip(axes, predictions, labels):
        PredictionErrorDisplay.from_predictions(
            y_true=y_test.values,
            y_pred=pred,
            kind="residual_vs_predicted",
            scatter_kwargs={"alpha": 0.3},
            ax=ax,
        )
        ax.set(xlabel="Predicted demand", ylabel="True demand")
        ax.legend(["Best model", label])

    plt.show()




.. image-sg:: /auto_examples/applications/images/sphx_glr_plot_time_series_lagged_features_002.png
   :alt: Non-linear regression models
   :srcset: /auto_examples/applications/images/sphx_glr_plot_time_series_lagged_features_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 414-437

Conclusion
----------
Through this example we explored time series forecasting using lagged
features. We compared a naive regression (using the standardized
:class:`~sklearn.model_selection.train_test_split`) with a proper time
series evaluation strategy using
:class:`~sklearn.model_selection.TimeSeriesSplit`. We observed that the
model trained using :class:`~sklearn.model_selection.train_test_split`,
having a default value of `shuffle` set to `True` produced an overly
optimistic Mean Average Percentage Error (MAPE). The results
produced from the time-based split better represent the performance
of our time-series regression model. We also analyzed the predictive uncertainty
of our model via Quantile Regression. Predictions based on the 5th and
95th percentile using `loss="quantile"` provide us with a quantitative estimate
of the uncertainty of the forecasts made by our time series regression model.
Uncertainty estimation can also be performed
using `MAPIE <https://mapie.readthedocs.io/en/latest/index.html>`_,
that provides an implementation based on recent work on conformal prediction
methods and estimates both aleatoric and epistemic uncertainty at the same time.
Furthermore, functionalities provided
by `sktime <https://www.sktime.net/en/latest/users.html>`_
can be used to extend scikit-learn estimators by making use of recursive time
series forecasting, that enables dynamic predictions of future values.


.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 14.788 seconds)


.. _sphx_glr_download_auto_examples_applications_plot_time_series_lagged_features.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: binder-badge

      .. image:: images/binder_badge_logo.svg
        :target: https://mybinder.org/v2/gh/scikit-learn/scikit-learn/1.4.X?urlpath=lab/tree/notebooks/auto_examples/applications/plot_time_series_lagged_features.ipynb
        :alt: Launch binder
        :width: 150 px

    .. container:: lite-badge

      .. image:: images/jupyterlite_badge_logo.svg
        :target: ../../lite/lab/?path=auto_examples/applications/plot_time_series_lagged_features.ipynb
        :alt: Launch JupyterLite
        :width: 150 px

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_time_series_lagged_features.ipynb <plot_time_series_lagged_features.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_time_series_lagged_features.py <plot_time_series_lagged_features.py>`


.. include:: plot_time_series_lagged_features.recommendations


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_