Multiple seasonalities
In this example we will show how to forecast data with multiple seasonalities using an MSTL.
Tip
For this task, StatsForecast’s MSTL is 68% more accurate and 600% faster than Prophet and NeuralProphet. (Reproduce experiments here)
Multiple seasonal data refers to time series that have more than one clear seasonality. Multiple seasonality is traditionally present in data that is sampled at a low frequency. For example, hourly electricity data exhibits daily and weekly seasonality. That means that there are clear patterns of electricity consumption for specific hours of the day like 6:00pm vs 3:00am or for specific days like Sunday vs Friday.
Traditional statistical models are not able to model more than one
seasonal length. In this example, we will show how to model the two
seasonalities efficiently using Multiple SeasonalTrend decompositions
with LOESS
(MSTL
).
For this example, we will use hourly electricity load data from Pennsylvania, New Jersey, and Maryland (PJM). The original data can be found here. (Click here for info on PJM)
First, we will load the data, then we will use the StatsForecast.fit
and StatsForecast.predict
methods to predict the next 24 hours. We
will then decompose the different elements of the time series into
trends and its multiple seasonalities. At the end, you will use the
StatsForecast.forecast
for productionready forecasting.
Outline
 Install libraries
 Load and explore the data
 Fit a multipleseasonality model
 Decompose the series in trend and seasonality
 Predict the next 24 hours
 Optional: Forecast in production
Tip
You can use Colab to run this Notebook interactively
Install libraries
We assume you have StatsForecast already installed. Check this guide for instructions on how to install StatsForecast.
Install the necessary packages using pip install statsforecast
“
Load Data
The input to StatsForecast is always a data frame in long
format with
three columns: unique_id
, ds
and y
:

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 YYYYMMDD for a date or YYYYMMDD HH:MM:SS for a timestamp. 
The
y
(numeric) represents the measurement we wish to forecast. We will rename the
You will read the data with pandas and change the necessary names. This step should take around 2s.
import pandas as pd
df = pd.read_csv('https://raw.githubusercontent.com/jnagura/Energyconsumptionpredictionanalysis/master/PJM_Load_hourly.csv')
df.columns = ['ds', 'y']
df.insert(0, 'unique_id', 'PJM_Load_hourly')
df.tail()
unique_id  ds  y  

32891  PJM_Load_hourly  20010101 20:00:00  35209.0 
32892  PJM_Load_hourly  20010101 21:00:00  34791.0 
32893  PJM_Load_hourly  20010101 22:00:00  33669.0 
32894  PJM_Load_hourly  20010101 23:00:00  31809.0 
32895  PJM_Load_hourly  20010102 00:00:00  29506.0 
StatsForecast can handle unsorted data, however, for plotting purposes, it is convenient to sort the data frame.
df = df.sort_values(['unique_id', 'ds']).reset_index(drop=True)
Plot the series using the plot
method from the
StatsForecast
class. This method prints up to 8 random series from the dataset and is
useful for basic EDA. In this case, it will print just one series given
that we have just one unique_id.
Note
The
StatsForecast.plot
method uses Plotly as a default engine. You can change to MatPlotLib by settingengine="matplotlib"
.
from statsforecast import StatsForecast
StatsForecast.plot(df)
/Users/max.mergenthaler/Nixtla/statsforecast/statsforecast/core.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
from tqdm.autonotebook import tqdm
The time series exhibits seasonal patterns. Moreover, the time series
contains 32,896
observations, so it is necessary to use very
computationally efficient methods.
Fit an MSTL model
The
MSTL
(Multiple SeasonalTrend decompositions using LOESS) model, originally
developed by Kasun Bandara, Rob J Hyndman and Christoph
Bergmeir, decomposes the time series
in multiple seasonalities using a Local Polynomial Regression (LOESS).
Then it forecasts the trend using a nonseasonal model and each
seasonality using a
SeasonalNaive
model. You can choose the nonseasonal model you want to use to forecast
the trend component of the MSTL model. In this example, we will use an
AutoARIMA.
Import the
StatsForecast
class and the models you need.
from statsforecast import StatsForecast
from statsforecast.models import MSTL, AutoARIMA
First, we must define the model parameters. As mentioned before, the
electricity load presents seasonalities every 24 hours (Hourly) and
every 24 * 7 (Daily) hours. Therefore, we will use [24, 24 * 7]
for
season length. The trend component will be forecasted with an
AutoARIMA
model. (You can also try with:
AutoTheta
,
AutoCES
,
and
AutoETS
)
# 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
)]
We fit the models by instantiating a new
StatsForecast
object with the following required parameters:

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.)
Any settings are passed into the constructor. Then you call its fit method and pass in the historical data frame.
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)
Use the fit
method to fit each model to each time series. In this
case, we are just fitting one model to one series. Check this guide to
learn how to fit many models to many
series.
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
Once the model is fitted, access the decomposition using the fitted_
attribute of
StatsForecast
.
This attribute stores all relevant information of the fitted models for
each of the time series.
In this case, we are fitting a single model for a single time series, so
by accessing the fitted_ location [0, 0] we will find the relevant
information of our model. The
MSTL
class generates a model_
attribute that contains the way the series
was decomposed.
sf.fitted_[0, 0].model_
data  trend  seasonal24  seasonal168  remainder  

0  22259.0  26183.898892  5215.124554  609.000432  681.225229 
1  21244.0  26181.599305  6255.673234  603.823918  714.250011 
2  20651.0  26179.294886  6905.329895  636.820423  740.214587 
3  20421.0  26176.985472  7073.420118  615.825999  701.608647 
4  20713.0  26174.670877  7062.395760  991.521912  609.202971 
…  …  …  …  …  … 
32891  36392.0  33123.552727  4387.149171  488.177882  630.524015 
32892  35082.0  33148.242575  3479.852929  682.928737  863.166767 
32893  33890.0  33172.926165  2307.808829  650.566775  940.168219 
32894  32590.0  33197.603322  748.587723  555.177849  801.013195 
32895  31569.0  33222.273902  967.124123  265.895357  420.254422 
We will use matplotlib, to visualize the different components of the series.
import matplotlib.pyplot as plt
sf.fitted_[0, 0].model_.tail(24 * 28).plot(subplots=True, grid=True)
plt.tight_layout()
plt.show()
We observe a clear upward trend (orange line) and seasonality repeating every day (24H) and every week (168H).
Predict the next 24 hours
Probabilistic forecasting with levels
To generate forecasts use the predict
method.
The predict
method takes two arguments: forecasts the next h
(for
horizon) and level
.

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 thelevel
(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.
The forecast object here is a new data frame that includes a column with the name of the model and the y hat values, as well as columns for the uncertainty intervals.
This step should take less than 1 second.
forecasts = sf.predict(h=24, level=[90])
forecasts.head()
ds  MSTL  MSTLlo90  MSTLhi90  

unique_id  
PJM_Load_hourly  20020101 01:00:00  29956.744141  29585.187500  30328.298828 
PJM_Load_hourly  20020101 02:00:00  29057.691406  28407.498047  29707.884766 
PJM_Load_hourly  20020101 03:00:00  28654.699219  27767.101562  29542.298828 
PJM_Load_hourly  20020101 04:00:00  28499.009766  27407.640625  29590.378906 
PJM_Load_hourly  20020101 05:00:00  28821.716797  27552.236328  30091.197266 
You can plot the forecast by calling the StatsForecast.plot
method and
passing in your forecast dataframe.
sf.plot(df, forecasts, max_insample_length=24 * 7)
Forecast in production
If you want to gain speed in productive settings where you have multiple
series or models we recommend using the
StatsForecast.forecast
method instead of .fit
and .predict
.
The main difference is that the .forecast
doest not store the fitted
values and is highly scalable in distributed environments.
The
forecast
method takes two arguments: forecasts next h
(horizon) and level
.

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 thelevel
(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.
The forecast object here is a new data frame that includes a column with the name of the model and the y hat values, as well as columns for the uncertainty intervals. Depending on your computer, this step should take around 1min. (If you want to speed things up to a couple of seconds, remove the AutoModels like ARIMA and Theta)
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(h=24, level=[90])
forecasts_df.head()
ds  MSTL  MSTLlo90  MSTLhi90  

unique_id  
PJM_Load_hourly  20020101 01:00:00  29956.744141  29585.187500  30328.298828 
PJM_Load_hourly  20020101 02:00:00  29057.691406  28407.498047  29707.884766 
PJM_Load_hourly  20020101 03:00:00  28654.699219  27767.101562  29542.298828 
PJM_Load_hourly  20020101 04:00:00  28499.009766  27407.640625  29590.378906 
PJM_Load_hourly  20020101 05:00:00  28821.716797  27552.236328  30091.197266 