多季节性趋势 (MSTL)


# 此单元格将不会被渲染，但用于隐藏警告并限制显示的行数。

import warnings
warnings.filterwarnings("ignore")

import logging
logging.getLogger('statsforecast').setLevel(logging.ERROR)

import pandas as pd
pd.set_option('display.max_rows', 6)

使用 Statsforecast 的 MSTL 模型 的逐步指南。

在本指南中，我们将熟悉主要的 StatsForecast 类以及一些相关的方法，例如 StatsForecast.plot、StatsForecast.forecast 和 StatsForecast.cross_validation 等。

介绍

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

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

MSTL分解过程分为几个阶段：

趋势估计：使用LOESS来估计时间序列的趋势分量。LOESS是一种非参数平滑方法，可以局部拟合数据，从而捕捉复杂的趋势模式。
季节分量估计：应用季节分解技术以识别和建模数据中存在的不同季节模式。这涉及提取和建模季节分量，例如每日、每周或每年的模式。
残差估计：残差被计算为原始时间序列与趋势和季节分量估计值之和之间的差。残差表示未被趋势和季节模式解释的变异，可能包含额外的信息或噪声。

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

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

MSTL

时间序列分析中的一个重要目标是将序列分解为一组不可观察（潜在）成分，这些成分能够与不同类型的时间变化相关联。时间序列分解的思想非常古老，17 世纪的天文学家就已使用它来计算行星轨道。Persons 是首位明确阐述未观察成分假设的人。根据 Persons 的观点，时间序列由四种波动组成：

长期趋势或世俗趋势；
在长期趋势上叠加的周期性运动。这些周期在工业繁荣时期似乎达到高峰，而在经济萧条时期则达到谷底，其上涨和下跌构成了商业周期；
每年内的季节性运动，其形状取决于系列的性质；
由于影响单个变量或其他重大事件（如战争和影响多个变量的国家灾难）而产生的残差波动。

传统上，四种变异被假定互相独立，并通过加法分解模型指定：

\[ \begin{equation} y_t= T_t +C_t +S_t +I_t, t=1,\ \cdots, n \tag 1 \end{equation} \]

其中 \(y_t\) 表示时间 \(t\) 的观察序列，\(T_t\) 为长期趋势，\(C_t\) 为商业周期，\(S_t\) 是季节性，\(I_t\) 为不规则成分。

如果潜在成分之间存在依赖关系，这种关系通过乘法模型指定：

\[ \begin{equation} y_t= T_t \times C_t \times S_t \times I_t, t=1,\ \cdots, n \tag 2 \end{equation} \]

在这里，\(S_t\) 和 \(I_t\) 以趋势-周期 \(T_t \times C_t\) 的比例表达。在某些情况下，会使用混合的加法-乘法模型。

LOESS（局部回归平滑）

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

LOESS 计算涉及以下步骤：

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

MSTL的一般特性

MSTL模型（使用LOESS的多季节-趋势分解）具有多个使其在时间序列分析中非常有用的特性。以下是其一些特性的列表：

多季节组件的分解：MSTL模型能够处理同时表现出多种季节模式的时间序列。您可以有效识别和建模数据中存在的不同季节组件。
检测复杂趋势的灵活性：由于使用了LOESS，MSTL模型可以捕捉数据中的复杂趋势模式。这包括非线性趋势和时间序列中的突变变化。
适应不同季节频率的能力：MSTL模型能够处理具有不同季节频率的数据，例如每日、每周、每月，甚至每年的模式。您可以识别和建模不同周期长度的季节模式。（参见）季节周期

频率

数据	分钟	小时	天	周	年
每日				7	365.25
每小时			24	168	8766
每半小时			48	336	17532
分钟		60	1440	10080	525960
秒	60	3600	86400	604800	31557600

平滑噪声和异常值的能力：LOESS中使用的平滑过程可以减少时间序列中噪声和异常值的影响。这可以改善对潜在模式的检测，使得趋势和季节性分析更加容易。
改进的预测：通过将时间序列分解成季节、趋势和残差组件，MSTL模型可以提供更准确的预测。通过将趋势和季节模式外推到未来，并添加随机残差，可以生成预测。
更详细的解释和分析：MSTL分解允许您以更详细的方式分析和理解时间序列的不同组件。这有助于识别季节模式、趋势变化和评估残差的变异性。
高效实现：尽管具体实现可能有所不同，但MSTL模型的计算可以高效完成，特别是当LOESS与优化计算算法结合使用时。

这些特性使得MSTL模型成为探索性时间序列分析、数据预测和在多季节组件和复杂趋势中进行模式检测的有用工具。

加载库和数据

Tip

需要使用Statsforecast。安装方法请参见说明。

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

import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
plt.style.use('grayscale') # 五三八 灰度 经典
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()

	Time	Ads
0	2017-09-13T00:00:00	80115
1	2017-09-13T01:00:00	79885
2	2017-09-13T02:00:00	89325
3	2017-09-13T03:00:00	101930
4	2017-09-13T04:00:00	121630

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

	ds	y	unique_id
0	2017-09-13T00:00:00	80115	1
1	2017-09-13T01:00:00	79885	1
2	2017-09-13T02:00:00	89325	1
3	2017-09-13T03:00:00	101930	1
4	2017-09-13T04:00:00	121630	1

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

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

plt.show();

时间序列的分解

如何分解时间序列，以及为什么要这样做？

在时间序列分析中，为了预测新的值，了解过去的数据是非常重要的。更正式地说，了解值随时间变化所遵循的模式是非常重要的。可能有许多原因导致我们的预测值朝错误的方向发展。基本上，时间序列由四个组成部分构成。这些组成部分的变化导致了时间序列模式的变化。这些组成部分包括：

水平： 这是随时间平均的主要值。
趋势： 趋势是导致时间序列中出现上升或下降模式的值。
季节性： 这是时间序列中短时间内发生的周期性事件，导致时间序列中出现短期的上升或下降模式。
残差/噪声： 这些是时间序列中的随机变化。

随着时间的推移，将这些组成部分结合起来会形成一个时间序列。大多数时间序列由水平和噪声/残差组成，趋势或季节性是可选值。

如果季节性和趋势是时间序列的一部分，那么它们将对预测值产生影响，因为预测时间序列的模式可能与之前的时间序列不同。

时间序列中组成部分的组合可以分为两种类型： * 加法型 * 乘法型

加法时间序列

如果时间序列的组成部分是相加而形成时间序列的，那么这个时间序列称为加法时间序列。从可视化上看，如果时间序列的上升或下降模式在整个序列中类似，我们可以说这个时间序列是加法的。任何加法时间序列的数学函数可以表示为： \[y(t) = 水平 + 趋势 + 季节性 + 噪声\]

乘法时间序列

如果时间序列的组成部分是相乘而形成的，那么这个时间序列称为乘法时间序列。从可视化上看，如果时间序列随着时间的推移呈指数增长或下降，那么这个时间序列可以被认为是乘法时间序列。乘法时间序列的数学函数可以表示为：

\[y(t) = 水平 * 趋势 * 季节性 * 噪声\]

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

Unable to display output for mime type(s): application/vnd.plotly.v1+json

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

让我们将数据划分为以下几组： 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()

MSTL方法的实现与StatsForecast

要了解更多关于MSTL Model函数的参数，列表如下。有关更多信息，请访问[文档][这里](https://nixtla.github.io/statsforecast/src/core/models.html#multiple-seasonalities)。

season_length : Union[int, List[int]
    每单位时间的观测值数量。对于多个季节性，请使用列表。
trend_forecaster : model
    用于预测趋势成分的StatsForecast模型。
alias : str
    模型的自定义名称。

加载库

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) # 预测数量

models = [MSTL(season_length=[24, 168], # 时间序列的季节性 
trend_forecaster=AutoARIMA(prediction_intervals=ConformalIntervals(n_windows=3, h=horizon)))]

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

模型：一个模型的列表。从模型中选择你想要的模型并导入它们。

freq: 一个字符串，表示数据的频率。（请参阅 pandas 的可用频率。）
n_jobs: n_jobs: 整数，用于并行处理的作业数量，使用 -1 表示所有核心。
fallback_model: 如果一个模型失败，使用的备用模型。

任何设置都将传递给构造函数。然后你调用它的 fit 方法并传入历史数据框。

sf = StatsForecast(df=df,
                   models=models,
                   freq='H', 
                   n_jobs=-1)

拟合模型

sf.fit()

StatsForecast(models=[MSTL])

让我们来看一下我们的 MSTL 模型 的结果。我们可以通过以下指令查看它：

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

	data	trend	seasonal24	seasonal168	remainder
0	80115.0	126671.223975	-42015.126409	-1882.157178	-2658.940388
1	79885.0	126633.864017	-43562.628320	-1744.429759	-1441.805937
2	89325.0	126596.519614	-36535.489263	104.986307	-841.016659
...	...	...	...	...	...
213	103080.0	119433.722676	-13488.006561	-4840.338623	1974.622508
214	95155.0	119403.419736	-27159.429895	-928.487998	3839.498158
215	80285.0	119373.145455	-39617.228552	-929.451981	1458.535078

216 rows × 5 columns

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

预测方法

如果您希望在生产环境中实现更高的效率，特别是在处理多个序列或模型时，我们建议使用 StatsForecast.forecast 方法，而不是使用 .fit 和 .predict。

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

预测方法需要两个参数：预测未来 h 步（时间范围）和 level。

h (int): 表示在未来预测的 h 步。在此案例中，是指提前30小时。
level (list of floats): 此可选参数用于概率预测。设定预测区间的水平（或置信百分位数）。例如，level=[90] 意味着模型期望实际值在该区间内的概率为90%。

预测对象是一个新的数据框，其中包含一个列，显示模型名称和 y hat 值，以及不确定性区间的列。根据您的计算机，这一步骤大约需要1分钟。（如果您希望将时间缩短到几秒钟，请删除诸如 ARIMA 和 Theta 的自动模型）

Y_hat = sf.forecast(horizon, fitted=True)
Y_hat

	ds	MSTL
unique_id
1	2017-09-22 00:00:00	73093.007812
1	2017-09-22 01:00:00	73268.078125
1	2017-09-22 02:00:00	80297.703125
...	...	...
1	2017-09-23 03:00:00	102918.906250
1	2017-09-23 04:00:00	115512.109375
1	2017-09-23 05:00:00	127950.929688

30 rows × 2 columns

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

	ds	y	MSTL
unique_id
1	2017-09-13 00:00:00	80115.0	79990.984375
1	2017-09-13 01:00:00	79885.0	78820.171875
1	2017-09-13 02:00:00	89325.0	88571.578125
1	2017-09-13 03:00:00	101930.0	102300.695312
1	2017-09-13 04:00:00	121630.0	122994.265625

StatsForecast.plot(values)

添加95%的置信区间与预测方法

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

	ds	MSTL	MSTL-lo-95	MSTL-hi-95
unique_id
1	2017-09-22 00:00:00	73093.007812	69844.109375	76341.914062
1	2017-09-22 01:00:00	73268.078125	67407.664062	79128.492188
1	2017-09-22 02:00:00	80297.703125	78214.578125	82380.820312
...	...	...	...	...
1	2017-09-23 03:00:00	102918.906250	100072.039062	105765.781250
1	2017-09-23 04:00:00	115512.109375	111221.421875	119802.789062
1	2017-09-23 05:00:00	127950.929688	124300.039062	131601.828125

30 rows × 4 columns

Y_hat=Y_hat.reset_index()
Y_hat

	unique_id	ds	MSTL
0	1	2017-09-22 00:00:00	73093.007812
1	1	2017-09-22 01:00:00	73268.078125
2	1	2017-09-22 02:00:00	80297.703125
...	...	...	...
27	1	2017-09-23 03:00:00	102918.906250
28	1	2017-09-23 04:00:00	115512.109375
29	1	2017-09-23 05:00:00	127950.929688

30 rows × 3 columns

# 将预测值与真实值连接起来
Y_hat1 = pd.concat([df,Y_hat],  keys=['unique_id', 'ds'])
Y_hat1

		ds	y	unique_id	MSTL
unique_id	0	2017-09-13 00:00:00	80115.0	1	NaN
	1	2017-09-13 01:00:00	79885.0	1	NaN
	2	2017-09-13 02:00:00	89325.0	1	NaN
...	...	...	...	...	...
ds	27	2017-09-23 03:00:00	NaN	1	102918.906250
	28	2017-09-23 04:00:00	NaN	1	115512.109375
	29	2017-09-23 05:00:00	NaN	1	127950.929688

246 rows × 4 columns

fig, ax = plt.subplots(1, 1)
plot_df = pd.concat([df, Y_hat1]).set_index('ds')
plot_df['y'].plot(ax=ax, linewidth=2)
plot_df["MSTL"].plot(ax=ax, linewidth=2, color="yellow")
ax.set_title(' Forecast', fontsize=22)
ax.set_ylabel('Ads watched (hourly data)', fontsize=20)
ax.set_xlabel('Monthly [t]', fontsize=20)
ax.legend(prop={'size': 15})
ax.grid(True)

预测方法及置信区间

要生成预测，请使用预测方法。

预测方法需要两个参数：预测接下来的 h（代表时间范围）和 level。

h (int): 表示预测向未来延伸的步数。在这种情况下，向前预测30个小时。
level (list of floats): 这个可选参数用于概率预测。设置预测区间的水平（或置信百分位）。例如，level=[95] 意味着模型期望实际值在该区间内的概率为95%。

这里的预测对象是一个新的数据框，其中包含模型名称列和y hat值，以及不确定性区间的列。

这一步应该花费少于1秒的时间。

sf.predict(h=horizon)

	ds	MSTL
unique_id
1	2017-09-22 00:00:00	73093.007812
1	2017-09-22 01:00:00	73268.078125
1	2017-09-22 02:00:00	80297.703125
...	...	...
1	2017-09-23 03:00:00	102918.906250
1	2017-09-23 04:00:00	115512.109375
1	2017-09-23 05:00:00	127950.929688

30 rows × 2 columns

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

	ds	MSTL	MSTL-lo-95	MSTL-lo-80	MSTL-hi-80	MSTL-hi-95
unique_id
1	2017-09-22 00:00:00	73093.007812	69844.109375	71204.914062	74981.101562	76341.914062
1	2017-09-22 01:00:00	73268.078125	67407.664062	69866.078125	76670.070312	79128.492188
1	2017-09-22 02:00:00	80297.703125	78214.578125	79040.671875	81554.726562	82380.820312
...	...	...	...	...	...	...
1	2017-09-23 03:00:00	102918.906250	100072.039062	100317.101562	105520.718750	105765.781250
1	2017-09-23 04:00:00	115512.109375	111221.421875	112638.765625	118385.445312	119802.789062
1	2017-09-23 05:00:00	127950.929688	124300.039062	124856.750000	131045.117188	131601.828125

30 rows × 6 columns

我们可以使用pandas函数pd.concat()将预测结果与历史数据进行合并，然后就可以利用这个结果进行图表绘制。

pd.concat([df, forecast_df]).set_index('ds')

	y	unique_id	MSTL	MSTL-lo-95	MSTL-lo-80	MSTL-hi-80	MSTL-hi-95
ds
2017-09-13 00:00:00	80115.0	1	NaN	NaN	NaN	NaN	NaN
2017-09-13 01:00:00	79885.0	1	NaN	NaN	NaN	NaN	NaN
2017-09-13 02:00:00	89325.0	1	NaN	NaN	NaN	NaN	NaN
...	...	...	...	...	...	...	...
2017-09-23 03:00:00	NaN	NaN	102918.906250	100072.039062	100317.101562	105520.718750	105765.781250
2017-09-23 04:00:00	NaN	NaN	115512.109375	111221.421875	112638.765625	118385.445312	119802.789062
2017-09-23 05:00:00	NaN	NaN	127950.929688	124300.039062	124856.750000	131045.117188	131601.828125

246 rows × 7 columns

df_plot= pd.concat([df, forecast_df]).set_index('ds')
df_plot

	y	unique_id	MSTL	MSTL-lo-95	MSTL-lo-80	MSTL-hi-80	MSTL-hi-95
ds
2017-09-13 00:00:00	80115.0	1	NaN	NaN	NaN	NaN	NaN
2017-09-13 01:00:00	79885.0	1	NaN	NaN	NaN	NaN	NaN
2017-09-13 02:00:00	89325.0	1	NaN	NaN	NaN	NaN	NaN
...	...	...	...	...	...	...	...
2017-09-23 03:00:00	NaN	NaN	102918.906250	100072.039062	100317.101562	105520.718750	105765.781250
2017-09-23 04:00:00	NaN	NaN	115512.109375	111221.421875	112638.765625	118385.445312	119802.789062
2017-09-23 05:00:00	NaN	NaN	127950.929688	124300.039062	124856.750000	131045.117188	131601.828125

246 rows × 7 columns

现在，让我们可视化我们预测的结果和我们的时间序列的历史数据，并绘制我们在以95%置信度进行预测时获得的置信区间。

def plot_forecasts(y_hist, y_pred, models):
  _, ax = plt.subplots(1, 1, figsize = (20, 7))
  df_plot = pd.concat([y_hist, y_pred]).set_index('ds').tail(12*10)
  df_plot[['y'] + models].plot(ax=ax, linewidth=3)
  colors = [['blue', "red"]]
  for model, color in zip(models, colors):
    ax.fill_between(df_plot.index, 
                        df_plot[f'{model}-lo-95'], 
                        df_plot[f'{model}-hi-95'],
                        alpha=.35,
                        color="blue",
                        label=f'{model}-level-95')
  for model, color in zip(models, colors):
      ax.fill_between(df_plot.index, 
                        df_plot[f'{model}-lo-80'], 
                        df_plot[f'{model}-hi-80'],
                        alpha=.20,
                        color="white",
                        label=f'{model}-level-80')
  ax.set_title('', fontsize=22)
  ax.set_ylabel("Ads watched (hourly data)", fontsize=20)
  ax.set_xlabel('Hourly', fontsize=20)
  ax.legend(prop={'size': 20})
  ax.grid(True)
  plt.show()

plot_forecasts(df, forecast_df, models=["MSTL"])

让我们使用Statsforecast中提供的plot函数绘制相同的图形，如下所示。

交叉验证

在之前的步骤中，我们利用历史数据来预测未来。然而，为了评估其准确性，我们还希望了解模型在过去的表现。为了评估模型在数据上的准确性和稳健性，执行交叉验证。

对于时间序列数据，交叉验证是通过在历史数据上定义滑动窗口并预测其后的时间段来进行的。这种形式的交叉验证使我们能够在更广泛的时间实例中更好地估算模型的预测能力，同时保持训练集中的数据是连续的，这正是我们模型所需的。

以下图展示了这种交叉验证策略：

执行时间序列交叉验证

时间序列模型的交叉验证被认为是最佳实践，但大多数实现速度很慢。statsforecast库将交叉验证实现为分布式操作，使得进行该过程更加节省时间。如果您有大数据集，还可以使用Ray、Dask或Spark在分布式集群中执行交叉验证。

在这种情况下，我们希望评估每个模型在过去5个月的性能(n_windows=)，每隔两个月进行一次预测(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)

crossvalidation_df 对象是一个新的数据框，包括以下列：

unique_id: 索引。如果您不喜欢使用索引，只需运行 crossvalidation_df.resetindex()。
ds: 日期戳或时间索引
cutoff: n_windows 的最后日期戳或时间索引。
y: 真实值
model: 包含模型名称和拟合值的列。

crossvalidation_df

	ds	cutoff	y	MSTL
unique_id
1	2017-09-15 18:00:00	2017-09-15 17:00:00	159725.0	158384.250000
1	2017-09-15 19:00:00	2017-09-15 17:00:00	161085.0	162015.171875
1	2017-09-15 20:00:00	2017-09-15 17:00:00	135520.0	138495.093750
...	...	...	...	...
1	2017-09-21 21:00:00	2017-09-20 17:00:00	103080.0	98109.875000
1	2017-09-21 22:00:00	2017-09-20 17:00:00	95155.0	86342.015625
1	2017-09-21 23:00:00	2017-09-20 17:00:00	80285.0	76815.976562

150 rows × 4 columns

现在我们将为每个截止期绘制预测图。为了使图表更加清晰，我们将重新命名每个期间的实际值。

cross_validation=crossvalidation_df.copy()
cross_validation.rename(columns = {'y' : 'actual'}, inplace = True) # 重命名实际值 

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

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

模型评估

我们现在可以使用适当的准确性指标来计算预测的准确性。在这里，我们将使用均方根误差（RMSE）。为了做到这一点，我们首先需要install datasetsforecast，这是一个由Nixtla开发的Python库，其中包含计算RMSE的函数。

%%capture
!pip install datasetsforecast

from datasetsforecast.losses import rmse

计算RMSE的函数需要两个参数：

实际值。
预测值，在这种情况下为MSTL模型。

rmse = rmse(cross_validation['actual'], cross_validation["MSTL"])
print("RMSE using cross-validation: ", rmse)

RMSE using cross-validation:  12548.12

致谢

我们要感谢Naren Castellon编写本教程。

参考文献

Changquan Huang • Alla Petukhina. 施普林格系列 (2022). 《使用Python进行应用时间序列分析和预测》。
Ivan Svetunkov. 使用增强动态自适应模型 (ADAM) 进行预测和分析。
James D. Hamilton. 《时间序列分析》普林斯顿大学出版社, 新泽西州普林斯顿, 第1版, 1994.
Nixtla 参数。
Pandas 可用频率。
Rob J. Hyndman 和 George Athanasopoulos (2018). “预测原则与实践，时间序列交叉验证”。。
季节性周期 - Rob J Hyndman。

Give us a ⭐ on Github

目录

介绍

MSTL

LOESS（局部回归平滑）

MSTL的一般特性

加载库和数据

使用 plot 方法探索数据

自相关图

时间序列的分解

加法时间序列

乘法时间序列

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

MSTL方法的实现与StatsForecast

加载库

实例化模型

拟合模型

预测方法

预测方法及置信区间

交叉验证

执行时间序列交叉验证

模型评估

致谢

参考文献