Multiple Seasonal Trend (MSTL)

​

Table of Contents

• Introduction
• IMAPA Model
• Loading libraries and data
• Explore data with the plot method
• Split the data into training and testing
• Implementation of IMAPA with StatsForecast
• Cross-validation
• Model evaluation
• References

​

Introduction

1. Trend estimation: LOESS is used to estimate the trend component of the time series. LOESS is a non-parametric smoothing method that locally fits data and allows complex trend patterns to be captured.
2. Estimation of seasonal components: Seasonal decomposition techniques are applied to identify and model the different seasonal patterns present in the data. This involves extracting and modeling seasonal components, such as daily, weekly, or yearly patterns.
3. Estimation of the residuals: The residuals are calculated as the difference between the original time series and the sum of the estimates of trend and seasonal components. Residuals represent variation not explained by trend and seasonal patterns and may contain additional information or noise.

​

MSTL

1. a long-term tendency or secular trend;
2. cyclical movements superimposed upon the long-term trend. These cycles appear to reach their peaks during periods of industrial prosperity and their troughs during periods of depressions, their rise and fall constituting the business cycle;
3. a seasonal movement within each year, the shape of which depends on the nature of the series;
4. residual variations due to changes impacting individual variables or other major events, such as wars and national catastrophes affecting a number of variables.

​

LOESS (Local Regression Smoothing)

• For each point t in the time series, a nearest neighbor window is selected.
• Weights are assigned to neighbors based on their proximity to t, using a weighting function, such as the Gaussian kernel.
• A weighted regression is performed using the neighbors and their assigned weights.
• The fitted value for point t is obtained based on local regression.
• The process is repeated for all points in the time series, thus obtaining a smoothed estimate of the trend.

​

MSTL General Properties

1. Decomposition of multiple seasonal components: The MSTL model is capable of handling time series that exhibit multiple seasonal patterns simultaneously. You can effectively identify and model different seasonal components present in the data.
2. Flexibility in detecting complex trends: Thanks to the use of LOESS, the MSTL model can capture complex trend patterns in the data. This includes non-linear trends and abrupt changes in the time series.
3. Adaptability to different seasonal frequencies: The MSTL model is capable of handling data with different seasonal frequencies, such as daily, weekly, monthly, or even yearly patterns. You can identify and model seasonal patterns of different cycle lengths. (see) Seasonal periods

1. Ability to smooth noise and outliers: The smoothing process used in LOESS allows to reduce the impact of noise and outliers in the time series. This can improve detection of underlying patterns and make it easier to analyze trend and seasonality.
2. Improved forecasting: By decomposing the time series into seasonal, trend, and residual components, the MSTL model can provide more accurate forecasts. Forecasts can be generated by extrapolating trend and seasonal patterns into the future, and adding the stochastic residuals.
3. More detailed interpretation and analysis: The MSTL decomposition allows you to analyze and understand the different components of the time series in a more detailed way. This facilitates the identification of seasonal patterns, changes in trend, and the evaluation of residual variability.
4. Efficient Implementation: Although the specific implementation may vary, the MSTL model can be calculated efficiently, especially when LOESS is used in combination with optimized calculation algorithms.

​

Loading libraries and data

Tip Statsforecast will be needed. To install, see instructions.

import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
plt.style.use('grayscale') # fivethirtyeight  grayscale  classic
plt.rcParams['lines.linewidth'] = 1.5
dark_style = {
    'figure.facecolor': '#008080',  # #212946
    'axes.facecolor': '#008080',
    'savefig.facecolor': '#008080',
    'axes.grid': True,
    'axes.grid.which': 'both',
    'axes.spines.left': False,
    'axes.spines.right': False,
    'axes.spines.top': False,
    'axes.spines.bottom': False,
    'grid.color': '#000000',  #2A3459
    'grid.linewidth': '1',
    'text.color': '0.9',
    'axes.labelcolor': '0.9',
    'xtick.color': '0.9',
    'ytick.color': '0.9',
    'font.size': 12 }
plt.rcParams.update(dark_style)


from pylab import rcParams
rcParams['figure.figsize'] = (18,7)

import pandas as pd
df=pd.read_csv("https://raw.githubusercontent.com/Naren8520/Serie-de-tiempo-con-Machine-Learning/main/Data/ads.csv")

df.head()

• The unique_id (string, int or category) represents an identifier for the series.
• The ds (datestamp) column should be of a format expected by Pandas, ideally 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.

df["unique_id"]="1"
df.columns=["ds", "y", "unique_id"]
df.head()

print(df.dtypes)

ds           object
y             int64
unique_id    object
dtype: object

df["ds"] = pd.to_datetime(df["ds"])

​

Explore Data with the plot method

from statsforecast import StatsForecast

StatsForecast.plot(df)

​

Autocorrelation plots

• Look for patterns in charts. Common patterns include repeating peaks and valleys, sawtooth patterns, and plateau patterns.
• Compare ACF and PACF charts. The PACF chart generally has fewer spikes than the ACF chart.
• Consider the length of the time series. ACF and PACF charts for longer time series will have more spikes.
• Use a confidence interval. The ACF and PACF plots also show confidence intervals for the autocorrelation values. If an autocorrelation value is outside the confidence interval, it is likely to be significant.

fig, axs = plt.subplots(nrows=1, ncols=2)

plot_acf(df["y"],  lags=30, ax=axs[0],color="fuchsia")
axs[0].set_title("Autocorrelation");

# Grafico
plot_pacf(df["y"],  lags=30, ax=axs[1],color="lime")
axs[1].set_title('Partial Autocorrelation')

plt.show();

​

Decomposition of the time series

• Level: This is the primary value that averages over time.
• Trend: The trend is the value that causes increasing or decreasing patterns in a time series.
• Seasonality: This is a cyclical event that occurs in a time series for a short time and causes short-term increasing or decreasing patterns in a time series.
• Residual/Noise: These are the random variations in the time series.

​

Additive time series

​

Multiplicative time series

from statsmodels.tsa.seasonal import seasonal_decompose
from plotly.subplots import make_subplots
import plotly.graph_objects as go

def plotSeasonalDecompose(
    x,
    model='additive',
    filt=None,
    period=None,
    two_sided=True,
    extrapolate_trend=0,
    title="Seasonal Decomposition"):

    result = seasonal_decompose(
            x, model=model, filt=filt, period=period,
            two_sided=two_sided, extrapolate_trend=extrapolate_trend)
    fig = make_subplots(
            rows=4, cols=1,
            subplot_titles=["Observed", "Trend", "Seasonal", "Residuals"])
    for idx, col in enumerate(['observed', 'trend', 'seasonal', 'resid']):
        fig.add_trace(
            go.Scatter(x=result.observed.index, y=getattr(result, col), mode='lines'),
                row=idx+1, col=1,
            )
    return fig

plotSeasonalDecompose(
    df["y"],
    model="additive",
    period=24,
    title="Seasonal Decomposition")

​

Split the data into training and testing

train = df[df.ds<='2017-09-20 17:00:00'] 
test = df[df.ds>'2017-09-20 17:00:00']

train.shape, test.shape

((186, 3), (30, 3))

sns.lineplot(train,x="ds", y="y", label="Train", linestyle="--",linewidth=2)
sns.lineplot(test, x="ds", y="y", label="Test", linewidth=2, color="yellow")
plt.title("Ads watched (hourly data)");
plt.xlabel("Hours")
plt.show()

​

Implementation of MSTL Method with StatsForecast

​

Load libraries

from statsforecast import StatsForecast
from statsforecast.models import MSTL, AutoARIMA

​

Instantiating Model

from statsforecast.utils import ConformalIntervals
horizon = len(test) # number of predictions

models = [MSTL(season_length=[24, 168], # seasonalities of the time series 
trend_forecaster=AutoARIMA(prediction_intervals=ConformalIntervals(n_windows=3, h=horizon)))]

• freq: a string indicating the frequency of the data. (See pandas’ available frequencies.)
• n_jobs: n_jobs: int, number of jobs used in the parallel processing, use -1 for all cores.
• fallback_model: a model to be used if a model fails.

sf = StatsForecast(models=models, freq='h')

​

Fit Model

sf.fit(df=train)

StatsForecast(models=[MSTL])

result=sf.fitted_[0,0].model_
result

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

​

Forecast Method

• h (int): represents the forecast h steps into the future. In this case, 30 hours 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.

Y_hat = sf.forecast(df=train, h=horizon, fitted=True)
Y_hat

values=sf.forecast_fitted_values()
values.head()

StatsForecast.plot(values)

sf.forecast(df=train, h=horizon, level=[95])

sf.plot(train, Y_hat)

​

Predict method with confidence interval

• h (int): represents the forecast h steps into the future. In this case, 30 hours 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=[95] means that the model expects the real value to be inside that interval 95% of the times.

sf.predict(h=horizon)

forecast_df = sf.predict(h=horizon, level=[80,95]) 
forecast_df

sf.plot(train, forecast_df, level=[80, 95])

​

Cross-validation

​

Perform time series cross-validation

• df: training data frame
• h (int): represents h steps into the future that are being forecasted. In this case, 500 hours ahead.
• step_size (int): step size between each window. In other words: how often do you want to run the forecasting processes.
• n_windows(int): number of windows used for cross validation. In other words: what number of forecasting processes in the past do you want to evaluate.

crossvalidation_df = sf.cross_validation(df=df,
                                         h=horizon,
                                         step_size=30,
                                         n_windows=5)

• unique_id: series identifier
• ds: datestamp or temporal index
• cutoff: the last datestamp or temporal index for the n_windows.
• y: true value
• model: columns with the model’s name and fitted value.

crossvalidation_df

from IPython.display import display

cross_validation=crossvalidation_df.copy()
cross_validation.rename(columns = {'y' : 'actual'}, inplace = True) # rename actual values 

cutoff = cross_validation['cutoff'].unique()

for k in range(len(cutoff)): 
    cv = cross_validation[cross_validation['cutoff'] == cutoff[k]]
    display(StatsForecast.plot(df, cv.loc[:, cv.columns != 'cutoff']))

​

Model Evaluation

from functools import partial

import utilsforecast.losses as ufl
from utilsforecast.evaluation import evaluate

evaluate(
    test.merge(Y_hat),
    metrics=[ufl.mae, ufl.mape, partial(ufl.mase, seasonality=24), ufl.rmse, ufl.smape],
    train_df=train,
)

​

References

1. Changquan Huang • Alla Petukhina. Springer series (2022). Applied Time Series Analysis and Forecasting with Python.
2. Ivan Svetunkov. Forecasting and Analytics with the Augmented Dynamic Adaptive Model (ADAM)
3. James D. Hamilton. Time Series Analysis Princeton University Press, Princeton, New Jersey, 1st Edition, 1994.
4. Nixtla Parameters.
5. Pandas available frequencies.
6. Rob J. Hyndman and George Athanasopoulos (2018). “Forecasting Principles and Practice (3rd ed)”.
7. Seasonal periods- Rob J Hyndman.