Multiple seasonalities

Tip For this task, StatsForecast’s MSTL is 68% more accurate and 600% faster than Prophet and NeuralProphet. (Reproduce experiments here)

1. Install libraries
2. Load and explore the data
3. Fit a multiple-seasonality model
4. Decompose the series in trend and seasonality
5. Predict the next 24 hours
6. Optional: Forecast in production

Tip You can use Colab to run this Notebook interactively

​

Install libraries

​

Load Data

• The unique_id (string, int or category) represents an identifier for the series.
• The ds (datestamp or int) column should be either an integer indexing time or a datestamp ideally like YYYY-MM-DD for a date or YYYY-MM-DD HH:MM:SS for a timestamp.
• The y (numeric) represents the measurement we wish to forecast. We will rename the

import pandas as pd

df = pd.read_csv('https://raw.githubusercontent.com/panambY/Hourly_Energy_Consumption/master/data/PJM_Load_hourly.csv')
df.columns = ['ds', 'y']
df.insert(0, 'unique_id', 'PJM_Load_hourly')
df = df.sort_values(['unique_id', 'ds']).reset_index(drop=True)
df.tail()

Note The StatsForecast.plot method uses matplotlib as a default engine. You can change to plotly by setting engine="plotly".

from statsforecast import StatsForecast

StatsForecast.plot(df)

​

Fit an MSTL model

from statsforecast.models import MSTL, AutoARIMA

# Create a list of models and instantiation parameters

models = [MSTL(
    season_length=[24, 24 * 7], # seasonalities of the time series 
    trend_forecaster=AutoARIMA() # model used to forecast trend
)]

• models: a list of models. Select the models you want from models and import them.
• freq: a string indicating the frequency of the data. (See panda’s available frequencies.)

sf = StatsForecast(
    models=models, # model used to fit each time series 
    freq='h', # frequency of the data
)

Tip StatsForecast also supports this optional parameter.

• n_jobs: n_jobs: int, number of jobs used in the parallel processing, use -1 for all cores. (Default: 1)
• fallback_model: a model to be used if a model fails. (Default: none)

Note StatsForecast achieves its blazing speed using JIT compiling through Numba. The first time you call the statsforecast class, the fit method should take around 10 seconds. The second time -once Numba compiled your settings- it should take less than 5s.

sf = sf.fit(df=df)

​

Decompose the series

sf.fitted_[0, 0].model_

import matplotlib.pyplot as plt

sf.fitted_[0, 0].model_.tail(24 * 28).plot(subplots=True, grid=True)
plt.tight_layout()
plt.show()

​

Predict the next 24 hours

Probabilistic forecasting with levels

• h (int): represents the forecast h steps into the future. In this case, 12 months ahead.
• level (list of floats): this optional parameter is used for probabilistic forecasting. Set the level (or confidence percentile) of your prediction interval. For example, level=[90] means that the model expects the real value to be inside that interval 90% of the times.

forecasts = sf.predict(h=24, level=[90])
forecasts.head()

sf.plot(df, forecasts, max_insample_length=24 * 7)

​

Forecast in production

• h (int): represents the forecast h steps into the future. In this case, 12 months ahead.
• level (list of floats): this optional parameter is used for probabilistic forecasting. Set the level (or confidence percentile) of your prediction interval. For example, level=[90] means that the model expects the real value to be inside that interval 90% of the times.

Note StatsForecast achieves its blazing speed using JIT compiling through Numba. The first time you call the statsforecast class, the fit method should take around 10 seconds. The second time -once Numba compiled your settings- it should take less than 5s.

forecasts_df = sf.forecast(df=df, h=24, level=[90])
forecasts_df.head()

​

References

• Bandara, Kasun & Hyndman, Rob & Bergmeir, Christoph. (2021). “MSTL: A Seasonal-Trend Decomposition Algorithm for Time Series with Multiple Seasonal Patterns”.

​

Next Steps

• Learn how to use cross-validation to assess the robustness of your model.