多重季节趋势 (MSTL)

​

目录

• 引言
• IMAPA 模型
• 加载库和数据
• 使用 plot 方法探索数据
• 将数据分为训练集和测试集
• 使用 StatsForecast 实现 MSTL
• 交叉验证
• 模型评估
• 参考资料

​

引言

MSTL 模型（使用 LOESS 进行多重季节趋势分解）是一种用于将时间序列分解为其季节、趋势和残差分量的方法。该方法基于使用 LOESS（局部回归平滑）来估计时间序列的各个分量。

MSTL 分解是经典季节趋势分解方法（也称为 Holt-Winters 分解）的扩展，旨在处理数据中存在多种季节模式的情况。例如，当时间序列同时表现出每日、每周和每年模式时，就可能出现这种情况。

MSTL 分解过程分几个阶段进行

1. 趋势估计：使用 LOESS 来估计时间序列的趋势分量。LOESS 是一种非参数平滑方法，它在局部拟合数据并允许捕获复杂的趋势模式。

2. 季节分量估计：应用季节分解技术来识别和建模数据中存在的不同季节模式。这包括提取和建模季节分量，例如每日、每周或每年的模式。

3. 残差估计：残差计算为原始时间序列与趋势和季节分量估计值之和之间的差值。残差表示趋势和季节模式无法解释的变动，并可能包含额外信息或噪声。

MSTL 分解使您能够更详细地分析和理解时间序列的不同分量，从而更容易进行预测和检测模式或异常。此外，LOESS 的使用提供了适应数据中不同趋势和季节模式的灵活性。

需要注意的是，MSTL 模型只是时间序列分解的可用方法之一，其选择将取决于数据的具体特征和应用背景。

​

MSTL

时间序列分析的一个重要目标是将序列分解为一组不可观测的（潜在的）分量，这些分量可以与不同类型的时间变动相关联。时间序列分解的思想非常古老，并被十七世纪的天文学家用于计算行星轨道。Persons 是第一个明确提出未观测分量假设的人。在 Persons 看来，时间序列由四种类型的波动组成

1. 长期趋势或长期倾向；
2. 叠加在长期趋势之上的周期性变动。这些周期似乎在工业繁荣时期达到顶峰，在经济萧条时期跌入谷底，其上升和下降构成了商业周期；
3. 每年内部的季节性变动，其形状取决于序列的性质；
4. 由于影响单个变量的变化或其他重大事件（例如影响许多变量的战争和国家灾难）而产生的残差变动。

传统上，这四种变动被认为是相互独立的，并通过加法分解模型来指定

其中 $y_t$ 表示时间 $t$ 的观测序列，$T_t$ 是长期趋势，$C_t$ 是商业周期，$S_t$ 是季节性，$I_t$ 是不规则变动。

如果潜在分量之间存在依赖关系，这种关系可以通过乘法模型来指定

其中 $S_t$ 和 $I_t$ 现在以趋势-周期 $T_t \times C_t$ 的比例表示。在某些情况下，会使用混合加法-乘法模型。

​

LOESS（局部回归平滑）

LOESS 是一种非参数平滑方法，用于估计局部拟合数据的光滑函数。对于时间序列中的每个点，LOESS 使用最近邻点执行加权回归。

LOESS 计算包括以下步骤

• 对于时间序列中的每个点 t，选择一个最近邻窗口。
• 根据邻点与 t 的接近程度，使用加权函数（例如高斯核）为邻点分配权重。
• 使用邻点及其分配的权重执行加权回归。
• 基于局部回归获得点 t 的拟合值。
• 对时间序列中的所有点重复该过程，从而获得趋势的平滑估计值。

​

MSTL 一般属性

MSTL 模型（使用 LOESS 的多重季节趋势分解）具有使其在时间序列分析中非常有用的几个属性。以下是其部分属性列表

1. 多重季节分量分解：MSTL 模型能够处理同时表现出多种季节模式的时间序列。您可以有效地识别和建模数据中存在的不同季节分量。

2. 检测复杂趋势的灵活性：由于使用了 LOESS，MSTL 模型可以捕获数据中的复杂趋势模式。这包括非线性趋势和时间序列中的突然变化。

3. 适应不同季节频率：MSTL 模型能够处理具有不同季节频率的数据，例如每日、每周、每月甚至每年的模式。您可以识别和建模不同周期长度的季节模式。（参见）季节周期

1. 平滑噪声和异常值的能力：LOESS 中使用的平滑过程可以减少噪声和异常值对时间序列的影响。这可以改进对潜在模式的检测，并使分析趋势和季节性更容易。

2. 改进预测：通过将时间序列分解为季节性、趋势和残差分量，MSTL 模型可以提供更准确的预测。预测可以通过将趋势和季节模式外推到未来，并添加随机残差来生成。

3. 更详细的解释和分析：MSTL 分解使您能够更详细地分析和理解时间序列的不同分量。这有助于识别季节模式、趋势变化以及评估残差变异性。

4. 高效实现：尽管具体实现可能有所不同，但 MSTL 模型可以高效计算，特别是在将 LOESS 与优化计算算法结合使用时。

这些属性使 MSTL 模型成为在存在多种季节分量和复杂趋势的情况下进行探索性时间序列分析、数据预测和模式检测的有用工具。

​

加载库和数据

提示

需要 Statsforecast。要安装，请参阅说明。

接下来，我们导入绘图库并配置绘图样式。

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()

StatsForecast 的输入始终是具有三列（unique_id、ds 和 y）的长格式数据框。

• unique_id（字符串、整数或类别）表示序列的标识符。

• ds（日期戳）列应采用 Pandas 期望的格式，理想情况下日期为 YYYY-MM-DD，时间戳为 YYYY-MM-DD HH:MM:SS。

• y（数字）表示我们希望预测的测量值。

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

print(df.dtypes)

ds           object
y             int64
unique_id    object
dtype: object

我们可以看到我们的时间变量 (ds) 是对象格式，我们需要将其转换为日期格式。

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

​

使用 plot 方法探索数据

使用 StatsForecast 类中的 plot 方法绘制一些序列。此方法从数据集中打印一个随机序列，对基本 EDA 很有用。

from statsforecast import StatsForecast

StatsForecast.plot(df)

​

自相关图

自相关 (ACF) 和偏自相关 (PACF) 图是用于分析时间序列的统计工具。ACF 图显示时间序列值与其滞后值之间的相关性，而 PACF 图显示时间序列值与其滞后值之间的相关性，但已去除先前滞后值的影响。

ACF 和 PACF 图可用于识别时间序列的结构，这有助于选择合适的时间序列模型。例如，如果 ACF 图显示重复的波峰和波谷模式，则表明时间序列是平稳的，这意味着它随时间具有相同的统计属性。如果 PACF 图显示快速下降的尖峰模式，则表明时间序列是可逆的，这意味着它可以反转以获得平稳时间序列。

ACF 和 PACF 图的重要性在于它们可以帮助分析师更好地理解时间序列的结构。这种理解有助于选择合适的时间序列模型，从而提高预测时间序列未来值的能力。

分析 ACF 和 PACF 图的方法

• 在图中寻找模式。常见模式包括重复的波峰和波谷、锯齿模式和平台模式。
• 比较 ACF 和 PACF 图。PACF 图通常比 ACF 图的尖峰少。
• 考虑时间序列的长度。较长时间序列的 ACF 和 PACF 图会有更多尖峰。
• 使用置信区间。ACF 和 PACF 图还显示了自相关值的置信区间。如果自相关值超出置信区间，则可能具有显著性。

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();

​

时间序列分解

如何以及为何分解时间序列？

在时间序列分析中预测新值时，了解过去数据非常重要。更正式地说，了解数值随时间遵循的模式非常重要。导致预测值方向错误的原因可能有很多。基本上，时间序列由四个分量组成。这些分量的变化导致时间序列模式的变化。这些分量是

• 水平项： 这是随时间平均的主要值。
• 趋势项： 趋势是导致时间序列中出现增加或减少模式的值。
• 季节项： 这是一种周期性事件，在时间序列中短时间发生，导致时间序列中出现短期增加或减少的模式。
• 残差/噪声项： 这是时间序列中的随机变动。

这些分量随时间组合形成时间序列。大多数时间序列包含水平项和噪声/残差项，而趋势项或季节项是可选值。

如果季节项和趋势项是时间序列的一部分，那么预测值会受到影响。因为预测时间序列的模式可能与之前的时间序列不同。

时间序列中分量的组合可以是两种类型：* 加法模型 * 乘法模型

​

加法时间序列

如果时间序列的分量相加构成时间序列。那么该时间序列称为加法时间序列。从可视化角度来看，如果时间序列的增加或减少模式在整个序列中相似，则可以说该时间序列是加法的。任何加法时间序列的数学函数可以表示为：$y(t) = level + Trend + seasonality + noise$

​

乘法时间序列

如果时间序列的分量相乘构成时间序列，则该时间序列称为乘法时间序列。从可视化角度来看，如果时间序列随时间呈指数增长或下降，则该时间序列可以视为乘法时间序列。乘法时间序列的数学函数可以表示为。

$y(t) = Level * Trend * seasonality * Noise$

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")

​

将数据分为训练集和测试集

让我们将数据分为两组：1. 用于训练 `MSTL 模型` 的数据。2. 用于测试模型的数据。

对于测试数据，我们将使用最后 30 小时的数据来测试和评估模型的性能。

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()

​

使用 StatsForecast 实现 MSTL 方法

​

加载库

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

​

实例化模型

导入并实例化模型。设置参数有时很棘手。大师 Rob Hyndmann 关于季节周期的这篇文章对于 `season_length` 可能会有帮助。

首先，我们必须定义模型参数。如前所述，糖果生产负荷每 24 小时（每时）和每 24 * 7 小时（每日）呈现季节性。因此，我们将使用 `[24, 24 * 7]` 作为季节长度。趋势分量将使用 AutoARIMA 模型进行预测。（您也可以尝试使用：AutoTheta、AutoCES 和 AutoETS）

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)))]

我们通过实例化一个新的 StatsForecast 对象并使用以下参数来拟合模型

models：模型列表。从 models 中选择您需要的模型并导入它们。

• freq: 一个字符串，表示数据的频率。（参见 pandas 的可用频率。）

• n_jobs: n_jobs: int，并行处理中使用的作业数，使用 -1 表示所有核心。

• fallback_model: 如果某个模型失败，则使用的备用模型。

任何设置都传递到构造函数中。然后调用其 fit 方法并传入历史数据框。

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

​

拟合模型

sf.fit(df=train)

StatsForecast(models=[MSTL])

让我们看看我们的 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 方法

如果您想在拥有多个序列或模型的生产环境中提高速度，我们建议使用 StatsForecast.forecast 方法而不是 .fit 和 .predict。

主要区别在于 .forecast 不存储拟合值，并且在分布式环境中具有高度可扩展性。

forecast 方法接受两个参数：预测接下来的 h（预测范围）和 level。

• h (int): 表示向未来预测 h 步。在此示例中，向前预测 30 小时。

• level (list of floats): 这个可选参数用于概率预测。设置预测区间的置信水平（或置信百分位数）。例如，level=[90] 意味着模型期望真实值有 90% 的几率落在此区间内。

这里的 forecast 对象是一个新的数据框，包含模型名称和 y_hat 值列，以及不确定性区间列。根据您的计算机性能，此步骤大约需要 1 分钟。（如果您想将速度加快到几秒钟，请删除像 ARIMA 和 Theta 这样的 AutoModels）。

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

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

StatsForecast.plot(values)

使用 forecast 方法添加 95% 置信区间

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

sf.plot(train, Y_hat)

​

带置信区间的 predict 方法

使用 predict 方法生成预测。

predict 方法接受两个参数：预测接下来的 h（预测范围）和 level。

• h (int): 表示向未来预测 h 步。在此示例中，向前预测 30 小时。

• level (list of floats): 这个可选参数用于概率预测。设置预测区间的置信水平（或置信百分位数）。例如，level=[95] 意味着模型期望真实值有 95% 的几率落在此区间内。

这里的 forecast 对象是一个新的数据框，包含模型名称和 y_hat 值列，以及不确定性区间列。

此步骤应少于 1 秒。

sf.predict(h=horizon)

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

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

​

交叉验证

在前面的步骤中，我们使用历史数据来预测未来。然而，为了评估其准确性，我们也想知道模型在过去会如何表现。为了评估模型在您的数据上的准确性和鲁棒性，请执行交叉验证。

对于时间序列数据，交叉验证通过在历史数据上定义一个滑动窗口并预测其后的时期来完成。这种形式的交叉验证使我们能够在更广泛的时间实例范围内更好地估计模型的预测能力，同时保持训练集中的数据连续，这是我们的模型所要求的。

下图描述了这种交叉验证策略

​

执行时间序列交叉验证

时间序列模型的交叉验证被认为是最佳实践，但大多数实现非常慢。statsforecast 库将交叉验证实现为分布式操作，从而减少了执行时间。如果您有大数据集，还可以使用 Ray、Dask 或 Spark 在分布式集群中执行交叉验证。

在本例中，我们希望评估每个模型在过去 5 个窗口（n_windows=）的性能，每隔 50 个步长（step_size=50）进行预测。根据您的计算机性能，此步骤大约需要 1 分钟。

StatsForecast 类中的 cross_validation 方法接受以下参数。

• df: 训练数据框

• h (int): 表示正在预测的未来 h 步。在本例中，向前预测 500 小时。

• step_size (int): 每个窗口之间的步长。换句话说：您希望多久运行一次预测过程。

• n_windows(int): 用于交叉验证的窗口数。换句话说：您希望评估过去多少个预测过程。

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

crossvaldation_df 对象是一个新的数据框，包含以下列

• unique_id: 序列标识符
• ds: 日期戳或时间索引
• cutoff: n_windows 的最后一个日期戳或时间索引。
• y: 真实值
• model: 包含模型名称和拟合值的列。

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']))

​

模型评估

现在我们将使用预测结果评估我们的模型，我们将使用不同类型的指标 MAE、MAPE、MASE、RMSE、SMAPE 来评估准确性。

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,
)

​

参考资料

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 参数.
5. Pandas 可用频率.
6. Rob J. Hyndman and George Athanasopoulos (2018). “Forecasting Principles and Practice (3rd ed)”.
7. 季节周期 - Rob J Hyndman.