季节性指数平滑优化模型


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

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 的 SeasonalExponentialSmoothingOptimized Model 的分步指南。

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

引言

季节性指数平滑优化（SESO）模型是一种预测技术，用于预测表现出季节性模式的时间序列的未来值。它是指数平滑方法的一种变体，该方法利用过去值和预测值的组合来生成预测。

SESO算法采用优化方法来寻找季节性指数平滑参数的最佳值。这些参数包括时间序列的水平、趋势和季节成分的平滑系数。

SESO模型特别适用于具有明显季节模式的时间序列的预测，例如季节性产品销售或季节性温度，以及许多其他领域。通过使用SESO，可以为商业规划和决策制定生成准确且有用的预测。

季节指数平滑模型

SESO模型基于指数平滑方法，该方法使用过去值和预测值的组合来生成预测。SESO模型的数学公式如下：

\[\hat{y}{t+1,s} = \alpha y_t + (1-\alpha) \hat{y}{t-1,s}\]

其中： - \(\hat{y}{t+1,s}\) 是季节 \(s\) 下一期的预测值。 - \(\alpha\) 是通过最小化平方误差来优化的平滑参数。 - \(y_t\) 是在时刻 \(t\) 站点 \(s\) 的当前观测值。 - \(\hat{y}{t-1,s}\) 是季节 \(s\) 前一期的预测值。

该方程表明，下一季节周期 \(s\) 的预测值是当前观测值和相同站点的前一个预测值的加权组合。平滑参数 \(\alpha\) 控制这两个项对最终预测的相对影响。较高的 α 值使当前观测值的权重更大，而前一个预测值的权重更小，使模型对时间序列中的近期变化更加敏感。另一方面，较低的 α 值则使前一个预测值的权重更大，而当前观测值的权重更小，使模型更加稳定和平滑。

平滑参数 \(\alpha\) 的最佳值是通过最小化模型生成的预测值与时间序列实际值之间的平方误差来确定的。

模型选择

在 SESO 模型的背景下，模型选择是指选择平滑参数和季节成分的最优值的过程。这些参数的最优值是指对于给定数据集产生最佳预测性能的值。

ETS 统计框架的一个很大优势是可以使用信息准则进行模型选择。\(AIC, AIC_c\) 和 \(BIC\) 也可以在此处用来确定哪些 ETS 模型最适合给定的时间序列。

对于 ETS 模型，赤池信息量准则（AIC）定义为 \[\text{AIC} = -2\log(L) + 2k,\]

其中 \(L\) 是模型的似然，\(k\) 是已估计的参数总数和初始状态数（包括残差方差）。

经过小样本偏差修正的 AIC (\(AIC_c\)) 定义为 \[\text{AIC}_{\text{c}} = \text{AIC} + \frac{2k(k+1)}{T-k-1},\]

而贝叶斯信息准则（BIC）是 \[\text{BIC} = \text{AIC} + k[\log(T)-2].\]

这些准则在模型的拟合优度和复杂性之间进行平衡，并提供了一种选择模型的方法，以最大化数据的似然性，同时最小化参数的数量。

除了这些技术外，专家判断和领域知识也可以用来选择最优的 SESO 模型。这涉及考虑时间序列的基本动态、季节性模式和任何可能影响模型选择的相关因素。

总体而言，SESO 模型的模型选择过程涉及统计技术、信息准则和专家判断的结合，以识别平滑参数和季节成分的最优值，从而为给定数据集产生最佳预测性能。

加载库和数据

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)

增强型迪基-福勒检验

增强型迪基-福勒（ADF）检验是一种统计检验，旨在确定时间序列数据中是否存在单位根。单位根可能导致时间序列分析中的结果不可预测。在单位根检验中构建一个原假设，以确定时间序列数据受趋势影响的程度。接受原假设意味着我们接受时间序列数据不平稳的证据。拒绝原假设或接受备择假设则意味着我们接受时间序列数据由一个平稳过程生成的证据。这个过程也被称为平稳趋势。ADF检验统计量的值为负数。较低的ADF值表示对原假设的更强拒绝。

增强型迪基-福勒检验是一种常用的统计检验，用于测试给定的时间序列是否平稳。我们可以通过定义原假设和备择假设来实现这一点。

原假设：时间序列是非平稳的。它给出一个依赖于时间的趋势。备择假设：时间序列是平稳的。换句话说，序列不依赖于时间。

当ADF或t统计量 < 临界值时：拒绝原假设，时间序列是平稳的。当ADF或t统计量 > 临界值时：未能拒绝原假设，时间序列是非平稳的。

from statsmodels.tsa.stattools import adfuller

def Augmented_Dickey_Fuller_Test_func(series , column_name):
    print (f'Dickey-Fuller test results for columns: {column_name}')
    dftest = adfuller(series, autolag='AIC')
    dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','No Lags Used','Number of observations used'])
    for key,value in dftest[4].items():
       dfoutput['Critical Value (%s)'%key] = value
    print (dfoutput)
    if dftest[1] <= 0.05:
        print("Conclusion:====>")
        print("Reject the null hypothesis")
        print("The data is stationary")
    else:
        print("Conclusion:====>")
        print("The null hypothesis cannot be rejected")
        print("The data is not stationary")

Augmented_Dickey_Fuller_Test_func(df["y"],'Ads')

Dickey-Fuller test results for columns: Ads
Test Statistic         -7.089634e+00
p-value                 4.444804e-10
No Lags Used            9.000000e+00
                            ...     
Critical Value (1%)    -3.462499e+00
Critical Value (5%)    -2.875675e+00
Critical Value (10%)   -2.574304e+00
Length: 7, dtype: float64
Conclusion:====>
Reject the null hypothesis
The data is stationary

自相关图

自相关（ACF）和偏自相关（PACF）的重要特征如下：

自相关（ACF）： 1. 识别时间依赖的模式：ACF显示观察值与其滞后值在不同时间间隔的相关性。帮助识别时间序列中的时间依赖模式，例如趋势或季节性的存在。

表示序列的“记忆”：ACF使我们能够确定过去的观察值对未来观察值的影响程度。如果ACF在几个滞后期显示显著的自相关性，这表明序列具有长期记忆，过去的观察值对预测未来的观察值是相关的。
有助于识别MA（移动平均）模型：ACF的形状可以揭示时间序列中移动平均成分的存在。在ACF显示显著相关性的滞后期可能指示MA模型的阶数。

偏自相关（PACF）： 1. 识别直接依赖关系：与ACF不同，PACF消除了中间滞后的间接影响，测量观察值与其滞后值之间的直接相关性。它有助于识别观察值与其滞后值之间的直接依赖关系，而不受中间滞后的影响。

有助于识别AR（自回归）模型：PACF的形状可以揭示时间序列中自回归成分的存在。在PACF显示显著相关性的滞后期可能指示AR模型的阶数。
与ACF结合使用：PACF与ACF结合使用以确定AR或MA模型的阶数。通过分析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.savefig("Gráfico de Densidad y qq")
plt.show();

时间序列的分解

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

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

水平（Level）： 这是一个随时间平均的主要值。
趋势（Trend）： 趋势是导致时间序列中出现增加或减少模式的值。
季节性（Seasonality）： 这是时间序列中发生的周期性事件，持续时间较短，导致时间序列中出现短期的增加或减少模式。
残差/噪音（Residual/Noise）： 这些是时间序列中的随机变化。

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

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

时间序列中组件的组合可以分为两种类型： * 加法（Additive） * 乘法（Multiplicative）

加法时间序列

如果时间序列的组件是相加构成时间序列的，那么这个时间序列被称为加法时间序列。通过可视化，我们可以说，如果时间序列的增加或减少模式在整个序列中是相似的，则时间序列是加法的。任何加法时间序列的数学函数可以表示为： \[y(t) = level + Trend + seasonality + noise\]

乘法时间序列

如果时间序列的组件是相乘在一起的，那么这个时间序列被称为乘法时间序列。通过可视化，如果时间序列随着时间的推移呈指数增长或下降，那么可以将其视为乘法时间序列。乘法时间序列的数学函数可以表示为：

\[y(t) = Level * Trend * seasonality * Noise\]

加法性质

from statsmodels.tsa.seasonal import seasonal_decompose 
a = seasonal_decompose(df["y"], model = "additive", period=12)
a.plot();

乘法性质

from statsmodels.tsa.seasonal import seasonal_decompose 
a = seasonal_decompose(df["y"], model = "Multiplicative", period=12)
a.plot();

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

让我们将数据分成几个部分

用于训练我们的 季节性指数平滑优化模型 的数据。
用于测试我们模型的数据。

对于测试数据，我们将使用最后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="--")
sns.lineplot(test, x="ds", y="y", label="Test")
plt.title("Ads watched (hourly data)");
plt.show()

使用StatsForecast实现季节性指数平滑优化

要了解SeasonalExponentialSmoothingOptimized Model函数的参数，以下是相关列表。更多信息，请访问文档。

alpha : float
    平滑参数。
season_length : int
    每单位时间的观察次数。例如：24小时数据。
alias : str
    模型的自定义名称。
prediction_intervals : Optional[ConformalIntervals]
    用于计算符合预测区间的信息。
    默认情况下，模型将计算本地预测
    区间。

导入库

from statsforecast import StatsForecast
from statsforecast.models import SeasonalExponentialSmoothingOptimized

构建模型

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

season_length = 24 # 每小时数据 
horizon = len(test) # 预测数量

models = [SeasonalExponentialSmoothingOptimized(season_length=season_length)]

我们通过实例化一个新的 StatsForecast 对象来拟合模型，所需的参数如下：

models: 模型的列表。从模型中选择所需的模型并导入它们。

freq: 一个字符串，表示数据的频率。（请查看 pandas 可用的频率。）
n_jobs: n_jobs: int，进行并行处理时使用的工作数量，使用 -1 以利用所有核心。
fallback_model: 如果某个模型失败，将使用的模型。

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

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

拟合模型

sf.fit()

StatsForecast(models=[SeasESOpt])

让我们看看我们的季节性指数平滑优化模型的结果。我们可以通过以下指令进行观察：

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

{'mean': array([ 80116.2  ,  79812.98 ,  84372.74 ,  99043.75 , 121514.94 ,
        116611.32 , 108388.39 , 103424.516, 108289.664, 116023.664,
        121505.336, 139466.77 , 146281.33 , 149037.22 , 149215.86 ,
        146367.08 , 148183.7  , 150684.8  , 157380.17 , 156142.53 ,
        127853.48 , 103063.23 ,  94190.15 ,  82520.19 ], dtype=float32),
 'fitted': array([       nan,        nan,        nan,        nan,        nan,
               nan,        nan,        nan,        nan,        nan,
               nan,        nan,        nan,        nan,        nan,
               nan,        nan,        nan,        nan,        nan,
               nan,        nan,        nan,        nan,  80115.   ,
         79885.   ,  89325.   , 101930.   , 121630.   , 116475.   ,
        106495.   , 102795.   , 108055.   , 116125.   , 131030.   ,
        149020.   , 157590.   , 150715.   , 149295.   , 150100.   ,
        144780.   , 150690.   , 163840.   , 166235.   , 139520.   ,
        105895.   ,  96780.   ,  82520.   ,  80115.1  ,  79839.5  ,
         87987.92 , 101954.1  , 121665.05 , 116490.55 , 106601.53 ,
        102791.8  , 107996.85 , 116122.15 , 130866.836, 147694.23 ,
        154525.55 , 150612.58 , 149282.45 , 149723.33 , 145698.03 ,
        150688.8  , 160717.08 , 162397.19 , 135501.27 , 103835.8  ,
         95715.3  ,  82522.45 ,  80077.1  ,  79792.41 ,  86293.46 ,
         99839.95 , 121632.7  , 116477.55 , 106770.836, 102752.484,
        107958.734, 116047.58 , 129459.34 , 145644.39 , 153794.8  ,
        150328.69 , 149269.83 , 149142.73 , 145707.48 , 150674.77 ,
        160501.92 , 162076.73 , 135508.52 , 112853.91 ,  96752.19 ,
         82573.375,  80154.68 ,  79882.93 ,  88212.44 , 100583.016,
        121575.77 , 116602.266, 108121.59 , 103169.36 , 108311.64 ,
        116219.   , 130052.28 , 144750.84 , 155067.58 , 150470.8  ,
        149314.48 , 149742.   , 146605.06 , 150642.36 , 158771.97 ,
        158364.27 , 131538.7  , 117874.29 ,  96740.12 ,  82683.74 ,
         80243.734,  79977.555,  88961.   , 100214.62 , 121485.71 ,
        116730.945, 109420.42 , 103663.266, 108754.33 , 116468.516,
        135878.83 , 149370.3  , 159073.2  , 151538.19 , 149452.73 ,
        151950.38 , 148868.33 , 150736.14 , 160848.06 , 161181.45 ,
        135859.64 , 113004.195,  96879.98 ,  82673.66 ,  80236.39 ,
         79961.27 ,  88670.77 , 100146.4  , 121508.66 , 116676.89 ,
        109030.95 , 103603.18 , 108643.336, 116329.48 , 130568.05 ,
        145525.64 , 152335.23 , 150896.27 , 149380.4  , 150026.   ,
        148226.   , 150732.88 , 160993.8  , 159284.78 , 135418.84 ,
        107124.39 ,  96455.72 ,  82642.07 ,  80217.38 ,  79908.37 ,
         86554.016,  99793.52 , 121487.02 , 116641.266, 108634.83 ,
        103507.15 , 108493.5  , 116208.03 , 126965.76 , 142833.   ,
        150244.78 , 150128.48 , 149358.84 , 148539.44 , 148797.55 ,
        150786.34 , 161078.64 , 160682.95 , 134904.86 , 105600.39 ,
         95623.21 ,  82608.34 ,  80215.01 ,  79890.38 ,  86310.35 ,
         99828.305, 121510.95 , 116638.2  , 108465.28 , 103486.48 ,
        108384.914, 116128.6  , 125062.48 , 142273.05 , 146089.02 ,
        149530.4  , 149280.52 , 146510.22 , 147309.14 , 150673.64 ,
        157855.16 , 156224.12 , 130665.7  , 101402.41 ,  93899.984,
         82542.766], dtype=float32)}

现在让我们可视化我们模型的拟合值。

正如我们所看到的，上述结果以字典的形式输出，接下来我们将使用 .get() 函数从字典中提取每个元素，然后将其保存到 pd.DataFrame() 中。

fitted=pd.DataFrame(result.get("fitted"), columns=["fitted"])
fitted["ds"]=df["ds"]
fitted

	fitted	ds
0	NaN	2017-09-13 00:00:00
1	NaN	2017-09-13 01:00:00
2	NaN	2017-09-13 02:00:00
...	...	...
213	101402.406250	2017-09-21 21:00:00
214	93899.984375	2017-09-21 22:00:00
215	82542.765625	2017-09-21 23:00:00

216 rows × 2 columns

sns.lineplot(df, x="ds", y="y", label="Actual", linewidth=2)
sns.lineplot(fitted,x="ds", y="fitted", label="Fitted", linestyle="--" )

plt.title("Ads watched (hourly data)");
plt.show()

预测方法

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

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

预测方法接受两个参数：预测下一个 h（时间范围）和 level。

h (int): 表示预测未来的 h 步。在这种情况下，是指向前 12 个月。

此处的预测对象是一个新的数据框，包含一个列，该列具有模型名称和 y hat 值，以及用于不确定性区间的列。根据您的计算机，此步骤应大约需要 1 分钟。

# 预测
Y_hat = sf.forecast(horizon, fitted=True)
Y_hat

	ds	SeasESOpt
unique_id
1	2017-09-22 00:00:00	80116.203125
1	2017-09-22 01:00:00	79812.976562
1	2017-09-22 02:00:00	84372.742188
...	...	...
1	2017-09-23 03:00:00	99043.750000
1	2017-09-23 04:00:00	121514.937500
1	2017-09-23 05:00:00	116611.320312

30 rows × 2 columns

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

	ds	y	SeasESOpt
unique_id
1	2017-09-13 00:00:00	80115.0	NaN
1	2017-09-13 01:00:00	79885.0	NaN
1	2017-09-13 02:00:00	89325.0	NaN
1	2017-09-13 03:00:00	101930.0	NaN
1	2017-09-13 04:00:00	121630.0	NaN

StatsForecast.plot(values)

Y_hat=Y_hat.reset_index()
Y_hat

	unique_id	ds	SeasESOpt
0	1	2017-09-22 00:00:00	80116.203125
1	1	2017-09-22 01:00:00	79812.976562
2	1	2017-09-22 02:00:00	84372.742188
...	...	...	...
27	1	2017-09-23 03:00:00	99043.750000
28	1	2017-09-23 04:00:00	121514.937500
29	1	2017-09-23 05:00:00	116611.320312

30 rows × 3 columns

Y_hat1= pd.concat([df,Y_hat])
Y_hat1

	ds	y	unique_id	SeasESOpt
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
...	...	...	...	...
27	2017-09-23 03:00:00	NaN	1	99043.750000
28	2017-09-23 04:00:00	NaN	1	121514.937500
29	2017-09-23 05:00:00	NaN	1	116611.320312

246 rows × 4 columns

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

带置信区间的预测方法

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

预测方法接受两个参数：预测下一个 h（表示时间范围）和 level。

h (int): 表示预测未来 h 步。在本例中，是 30 小时。

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

这一步应该在 1 秒钟内完成。

forecast_df = sf.predict(h=horizon) 
forecast_df

	ds	SeasESOpt
unique_id
1	2017-09-22 00:00:00	80116.203125
1	2017-09-22 01:00:00	79812.976562
1	2017-09-22 02:00:00	84372.742188
...	...	...
1	2017-09-23 03:00:00	99043.750000
1	2017-09-23 04:00:00	121514.937500
1	2017-09-23 05:00:00	116611.320312

30 rows × 2 columns

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

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

	y	unique_id	SeasESOpt
ds
2017-09-13 00:00:00	80115.0	1	NaN
2017-09-13 01:00:00	79885.0	1	NaN
2017-09-13 02:00:00	89325.0	1	NaN
...	...	...	...
2017-09-23 03:00:00	NaN	NaN	99043.750000
2017-09-23 04:00:00	NaN	NaN	121514.937500
2017-09-23 05:00:00	NaN	NaN	116611.320312

246 rows × 3 columns

现在让我们可视化我们预测的结果和时间序列的历史数据。

sns.lineplot(df_plot,x="ds", y="y", label="Actual")
sns.lineplot(df_plot, x="ds", y="SeasESOpt", label="SeasESOpt")
plt.show()

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

sf.plot(df, forecast_df)

交叉验证

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

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

以下图表描绘了这种交叉验证策略：

执行时间序列交叉验证

时间序列模型的交叉验证被认为是一种最佳实践，但大多数实现都非常缓慢。statsforecast库将交叉验证实现为分布式操作，使得这个过程的执行不那么耗时。如果您有大数据集，您还可以使用Ray、Dask或Spark在分布式集群中执行交叉验证。

在这种情况下，我们希望评估每个模型在过去5个月的表现(n_windows=)，每隔两个月进行一次预测(step_size=12)。根据您的计算机，这一步骤大约需要1分钟。

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

df: 训练数据框
h (int): 表示预测的未来h步。在本例中，预测12个月的未来。
step_size (int): 每个窗口之间的步长。换句话说：您希望多长时间运行一次预测过程。
n_windows (int): 用于交叉验证的窗口数量。换句话说：您希望评估过去多少个预测过程。

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

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

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

crossvalidation_df

	ds	cutoff	y	SeasESOpt
unique_id
1	2017-09-18 06:00:00	2017-09-18 05:00:00	99440.0	141401.750000
1	2017-09-18 07:00:00	2017-09-18 05:00:00	97655.0	152474.250000
1	2017-09-18 08:00:00	2017-09-18 05:00:00	97655.0	152482.796875
...	...	...	...	...
1	2017-09-21 21:00:00	2017-09-20 17:00:00	103080.0	105600.390625
1	2017-09-21 22:00:00	2017-09-20 17:00:00	95155.0	96717.390625
1	2017-09-21 23:00:00	2017-09-20 17:00:00	80285.0	82608.343750

90 rows × 4 columns

评估模型

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

%%capture
!pip install datasetsforecast

from datasetsforecast.losses import rmse

计算RMSE的函数接受两个参数：

实际值。
预测值，在本例中为“季节性指数平滑优化模型”。

rmse = rmse(crossvalidation_df['y'], crossvalidation_df["SeasESOpt"])
print("RMSE using cross-validation: ", rmse)

RMSE using cross-validation:  16560.705

致谢

我们要感谢Naren Castellon为编写此教程所做的贡献。

参考文献

Changquan Huang • Alla Petukhina. Springer 系列 (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

目录

引言

季节指数平滑模型

模型选择

加载库和数据

读取数据

使用plot方法探索数据

增强型迪基-福勒检验

自相关图

时间序列的分解

加法时间序列

乘法时间序列

加法性质

乘法性质

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

使用StatsForecast实现季节性指数平滑优化

导入库

构建模型

拟合模型

预测方法

带置信区间的预测方法

交叉验证

执行时间序列交叉验证

评估模型

致谢

参考文献