Where Is the Seasonal Naive Benchmark?

First Look at the Bike Sharing Data

The data used in the notebook are bike sharing demand observed as hourly time series for a period of about two years in 2011 and 2012. In contrast to other datasets of its kind, this dataset consists of a single series. The tutorial loads it from the OpenML repository:

from sklearn.datasets import fetch_openml

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

To take a sneak peek at what the time series looks like, the tutorial computes the average hourly demand over observed across the entire dataset:

Unsurprisingly, the bike sharing demand is highly seasonal. There is a within-day pattern that repeats on weekdays corresponding to demand from commuters. Also, weekends are clearly distinct from weekdays but likely similar across weekends (in theory this is not visible from this plot, but…). Absent any drastic trends, and ignoring the slower yearly seasonality given that we will focus on short-term forecasts, most of the signal from one week to the next will come from the daily and weekly seasonality. Any model we use should be designed to capture it. And that’s what the tutorial does, it demonstrates different ways to capture the seasonal patterns via feature engineering.

While most of the feature engineering is done to reduce the error of linear regression-based models, the best model of the tutorial turns out to be the gradient boosting model that is evaluated first and uses none of the fancy transformations. That’s the model we focus on now.

Setup of the Original Gradient Boosting Pipeline

First, the target column count is standardized by dividing by the largest observation. Thus we don’t predict the count of bike rentals directly, but the share of the largest possible hourly demand. This means that a mean absolute error of 0.01 corresponds to 1pp difference from the actually observed percent of maximum demand. The transformed column is assigned to y, then dropped from the dataframe to create the feature matrix X.

y = df["count"] / df["count"].max()
X = df.drop("count", axis="columns")

The initial model matrix consists of the following features (the range shown here is a bit misleading for the categorical columns). Note the weather-related columns such as temp, humidity, and windspeed. These are the problematic features mentioned in the introduction.

for col in X.columns:
    print(f"'{col}' with values in [{min(X[col])}, {max(X[col])}]")

## 'season' with values in [fall, winter]
## 'year' with values in [0.0, 1.0]
## 'month' with values in [1.0, 12.0]
## 'hour' with values in [0.0, 23.0]
## 'holiday' with values in [False, True]
## 'weekday' with values in [0.0, 6.0]
## 'workingday' with values in [False, True]
## 'weather' with values in [clear, rain]
## 'temp' with values in [0.8200000000000001, 41.0]
## 'feel_temp' with values in [0.0, 50.0]
## 'humidity' with values in [0.0, 1.0]
## 'windspeed' with values in [0.0, 56.996900000000004]

As heavy_rain was a rare occurence, it’s grouped together with rain.

X["weather"].replace(to_replace="heavy_rain", value="rain", inplace=True)

The tutorial uses a time series cross validation with five splits and a gap between train and validation of 48 hours. This gap of 48 hours is what we’ll use as the definition of the necessary lag between the moment of creating the prediction, and the moment that we predict. For example, if we need to schedule 48 hours in advance the number of bikes brought into circulation from repairs, then we can’t wait to observe the demand at 3pm today to plan for tomorrow’s demand. Instead, the latest information we can take into account is the demand observed yesterday, 48 hours in advance.

Thus, we won’t consider a daily seasonal naive model that projects the observation from 24 hours ago into the future. And because of this lag, the weather columns are somewhat unrealistic. But more on this in a moment.

We define the splits, each of which evaluates on a test set of about 42 days.

from sklearn.model_selection import TimeSeriesSplit

ts_cv = TimeSeriesSplit(
    n_splits=5,
    gap=48,
    max_train_size=10000,
    test_size=1000,
)

all_splits = list(ts_cv.split(X, y))

# this split is used to plot predictions from the models below
train_0, test_0 = all_splits[0]

The splits will be used to evaluate the model pipelines. For that, an evaluation helper is defined that uses both MAE and RMSE as metrics.

def print_metrics(mae, rmse):
  print(
        f"Mean Absolute Error:     {mae.mean():.3f} +/- {mae.std():.3f}\n"
        f"Root Mean Squared Error: {rmse.mean():.3f} +/- {rmse.std():.3f}"
    )

def evaluate(model, X, y, cv):
    cv_results = cross_validate(
        model,
        X,
        y,
        cv=cv,
        scoring=["neg_mean_absolute_error", "neg_root_mean_squared_error"],
    )
    mae = -cv_results["test_neg_mean_absolute_error"]
    rmse = -cv_results["test_neg_root_mean_squared_error"]
    print_metrics(mae, rmse)

Next, the tutorial defines the pipeline object (note the inclusion of the weather column).

import numpy as np
import matplotlib.pyplot as plt
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import OrdinalEncoder
from sklearn.compose import ColumnTransformer
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.model_selection import cross_validate

categorical_columns = [
    "weather",
    "season",
    "holiday",
    "workingday",
]
categories = [
    ["clear", "misty", "rain"],
    ["spring", "summer", "fall", "winter"],
    ["False", "True"],
    ["False", "True"],
]
ordinal_encoder = OrdinalEncoder(categories=categories)

gbrt_pipeline = make_pipeline(
    ColumnTransformer(
        transformers=[
            ("categorical", ordinal_encoder, categorical_columns),
        ],
        remainder="passthrough",
    ),
    HistGradientBoostingRegressor(
        categorical_features=range(4),
    ),
)

And finally the cross validation performance is being evaluated. Note that we pass the entire model matrix X here—including columns describing the hourly temperature, humidity, and wind speed!

evaluate(gbrt_pipeline, X, y, cv=ts_cv)

## Mean Absolute Error:     0.044 +/- 0.003
## Root Mean Squared Error: 0.068 +/- 0.005

This reproduces the performance shown in the tutorial notebook exactly, perfect. Given the standardization that was applied on the series, the mean absolute error corresponds to an error of about 4.4% of the highest hourly demand observed across the time period. Given the maximum of 977 hourly bike rentals, that’s an error of about 42 bike rentals.

Let’s drop the weather-related features and evaluate again.

Gradient Boosting Without Weather Features

We define the pipeline as before, dropping weather from the categorical columns and then explicitly selecting columns from X to not pass through the temperature, humidity, and wind speed. We only keep the columns for which we can know the values in advance.

categorical_columns_no_weather = [
    "season",
    "holiday",
    "workingday",
]
categories_no_weather = [
    ["spring", "summer", "fall", "winter"],
    ["False", "True"],
    ["False", "True"],
]
ordinal_encoder_no_weather = OrdinalEncoder(categories=categories_no_weather)

gbrt_no_weather_pipeline = make_pipeline(
    ColumnTransformer(
        transformers=[
            ("categorical", ordinal_encoder_no_weather, categorical_columns_no_weather),
        ],
        remainder="passthrough",
    ),
    HistGradientBoostingRegressor(
        categorical_features=range(3),
    ),
)

X_no_weather = X[
  ['season', 'year', 'month', 'hour', 'holiday', 'weekday', 'workingday']
]

Evaluate the new pipeline on the same time series splits:

evaluate(gbrt_no_weather_pipeline, X_no_weather, y, cv=ts_cv)

## Mean Absolute Error:     0.062 +/- 0.011
## Root Mean Squared Error: 0.094 +/- 0.015

The MAE increased by 1.8pp. Note also the increase in the standard deviation across splits. The performance clearly worsened, but still isn’t too shabby given that we didn’t put much effort into building this model.

Can we get away with even less work?

Comparison with Seasonal Naive Method

Given that the data has several seasonal effects (daily, weekly, yearly), there is no one way to set up a seasonal naive prediction: We could use the value from 24 hours ago, the value from a week ago (168 hours), or even a mix of the two.¹ As desribed above, we’ll focus on the weekly seasonal naive prediction due to the gap of 48 hours between train and test sets. A lag of one week is anyway more realistic for operational decisions.

Given the simplicity of the seasonal naive method, we just need to shift the label series. Since the cross validation splits use a gap of 48 hours that is smaller than our lag, we don’t need to do this for each cross validation split. We can do it once across the entire data set and still don’t introduce target leakage.²

snaive = y.shift(168)
print(snaive[165:172])

## 165         NaN
## 166         NaN
## 167         NaN
## 168    0.016377
## 169    0.040942
## 170    0.032753
## 171    0.013306
## Name: count, dtype: float64

A custom cross validation loop gives us the performance of the seasonal naive method (i.e., of the new column) compared to the target, on the same splits that we used for the gradient boosting models.

mae = np.empty(5)
rmse = np.empty(5)

for i in range(5):
    _, test_idx = all_splits[i]
    mae[i] = np.mean(np.abs(
      y.iloc[test_idx].values - snaive.iloc[test_idx].values
    ))
    rmse[i] = np.sqrt(np.mean(np.square(
      y.iloc[test_idx].values - snaive.iloc[test_idx].values
    )))

print_metrics(mae, rmse)

## Mean Absolute Error:     0.073 +/- 0.022
## Root Mean Squared Error: 0.123 +/- 0.036

How is that for a one-liner: 7.3% mean absolute error compared to the 6.2% of the gradient boosted trees model without weather features. Given that the seasonal naive method requires neither code nor infrastructure, and is inherently interpretable, one has to check the details of the business case whether any other method is worth the effort. A large chunk of the time series’ signal is already captured by this method.

Then again, the gradient boosted trees pipeline we compare against is the simplest version possible and only served as baseline in the tutorial. It could likely be improved further.

But the dataset at hand is also quite simple. Real use cases would likely have multiple series and require uncertainty quantification, thus increasing the effort going into building a model. This makes it even more important to compare against the simplest of benchmarks.

---

We can round this analysis off by plotting the predictions and errors of the models, similar to what was done in the scikit-learn notebook.

Graph of Predictions Over Time

We generate predictions for the first test set as defined by the cross validation split.

gbrt_pipeline.fit(X.iloc[train_0], y.iloc[train_0])
gbrt_predictions = gbrt_pipeline.predict(X.iloc[test_0])

gbrt_no_weather_pipeline.fit(X_no_weather.iloc[train_0], y.iloc[train_0])
gbrt_no_weather_predictions = gbrt_no_weather_pipeline.predict(
  X_no_weather.iloc[test_0]
)

snaive_predictions = snaive.iloc[test_0].values

Those we can then use to visualize the last week of the test set:³

last_hours = slice(-168, None)
fig, ax = plt.subplots(figsize=(12, 4))
fig.suptitle("Predictions by models and naive method")
ax.plot(
    gbrt_predictions[last_hours],
    "x-",
    label="Gradient Boosted Trees",
)
ax.plot(
    gbrt_no_weather_predictions[last_hours],
    "x-",
    label="Gradient Boosted Trees without Weather",
)
ax.plot(
    snaive_predictions[last_hours],
    "x-",
    label="Weekly Seasonal Naive",
)
ax.plot(
    y.iloc[test_0].values[last_hours],
    "x-",
    alpha=0.5,
    label="Actual demand",
    color="black",
)
_ = ax.legend()

While the graph is a bit busy, the seasonal patterns are clear. While the gradient bossted trees model without the weather information underestimates the demand on most of the weekdays and isn’t able to predict the rush-hour peaks, the seasonal naive does a bit better in these respects. But it also makes larger errors on this weekend in particular, presumably due to a very different previous weekend.

To compare more than just a week of predictions against actuals, we can use a scatterplot.

Graph of Predictions Against Actuals

The graph of predictions over time is great to see the actual behavior of the time series. But consistent problems can become apparent in a scatterplot of all predictions against realized values for the current test set. The predictions from a perfect model would fall on the diagonal line, any deviation from it is the error of the prediction.

fig, axes = plt.subplots(ncols=3, figsize=(12, 4), sharey=True)
predictions = [
    gbrt_predictions,
    gbrt_no_weather_predictions,
    snaive_predictions,
]
labels = [
    "GBT",
    "GBT without Weather",
    "Seasonal Naive",
]
for ax, pred, label in zip(axes, predictions, labels):
    ax.scatter(y.iloc[test_0].values, pred, alpha=0.3, label=label)
    ax.plot([0, 1], [0, 1], "--", label="Perfect model")
    ax.set(
        xlim=(0, 1),
        ylim=(0, 1),
        xlabel="True demand",
        ylabel="Predicted demand",
    )
    ax.legend()

plt.show()

As discussed in the tutorial notebook, the gradient boosted trees models systematically underestimate the peak demand. For the model without weather information, this also appears to hold for the range of lower demand—similar to what is visible in the graph of demand over time above.

In comparison, the seasonal naive method appears to make some errors that are larger in absolute size, but it doesn’t systematically under- or overestimate. This makes sense as its predictions adjust entirely from week to week: A too large prediction in one week tends to imply a too low prediction the week thereafter if the data is noisy around some average level. Then the seasonal naive method is low bias, high variance. A bias would easily be introduced by a strong trend that persisted over time in the series, though.