# 此单元格将不会被渲染，但用于隐藏警告并限制显示的行数。
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的SeasonalExponentialSmoothing Model的逐步指南。
在这个演练中，我们将熟悉主要的 StatsForecast 类以及一些相关的方法，如 StatsForecast.plot、StatsForecast.forecast 和 StatsForecast.cross_validation 等。
目录
介绍
简单指数平滑 (SES) 是一种预测方法，它使用历史值的加权平均来预测下一个值。权重分配给最近的值，而最旧的值则得到较低的权重。这是因为SES假设最近的值在预测未来时比旧值更相关。
SES通过一个简单的公式实现：
\[\hat{y}_{T+1|T} = \alpha y_T + \alpha(1-\alpha) y_{T-1} + \alpha(1-\alpha)^2 y_{T-2}+ \cdots,\]
平滑因子控制分配给最近值的权重。较高的α值意味着赋予较新的值更大的权重，而较低的α值则意味着更大的权重赋予旧值。
时间序列中的季节性是指在特定时间段内，时间序列变化的规律性重复模式。
季节性在时间序列分析中可能是一大挑战，因为它可能掩盖数据中的潜在趋势。
季节性是分析时间序列数据时一个重要的因素。通过理解数据中的季节模式，可以做出更准确的预测和更好的决策。
季节性指数平滑模型
最简单的指数平滑方法被称为简单指数平滑（SES）。该方法适用于没有明显趋势或季节模式的数据预测。
使用朴素方法，未来的所有预测都等于系列的最后一个观察值， \[\hat{y}_{T+h|T} = y_{T},\]
对于\(h=1,2,\dots\)。因此，朴素方法假设最近的观察值是唯一重要的，所有以前的观察值对未来没有信息。这可以被看作是一个加权平均，其中所有的权重都赋给了最后的观察值。
使用平均方法，所有未来的预测等于观测数据的简单平均值， \[\hat{y}_{T+h|T} = \frac1T \sum_{t=1}^T y_t,\]
对于\(h=1,2,\dots\)。因此，平均方法假设所有观察值的重要性相等，并在生成预测时给予它们相等的权重。
我们通常希望在这两个极端之间找到一些东西。例如，将更大的权重附加到更近期的观察值上，而不是遥远过去的观察值，这可能是合理的。这正是简单指数平滑背后的概念。预测是使用加权平均进行计算的，其中权重随着观察值来自更遥远的过去而指数减少——最小的权重与最古老的观察值相关联：
\[\begin{equation} \hat{y}_{T+1|T} = \alpha y_T + \alpha(1-\alpha) y_{T-1} + \alpha(1-\alpha)^2 y_{T-2}+ \cdots, \tag{1} \end{equation}\]
其中\(0 \le \alpha \le 1\)是平滑参数。时间\(T+1\)的一步预测是系列\(y_1,\dots,y_T\)中所有观察值的加权平均。权重减少的速率由参数\(\alpha\)控制。
对于任何在0和1之间的\(\alpha\)，附加到观察值的权重会随着时间的推移而指数下降，从而得名“指数平滑”。如果\(\alpha\)很小（即接近0），则更远的过去的观察值将获得更多的权重。如果\(\alpha\)很大（即接近1），则最近的观察值将获得更多的权重。对于极端情况，当\(\alpha=1\)时，\(\hat{y}_{T+1|T}=y_T\)，预测等于朴素预测。
如何确定季节性参数的值？
为了确定季节性调整的简单指数平滑（SES季节性调整）模型中的季节性参数s的值，可以根据数据的性质和分析的目标采用不同的方法。
以下是一些常见的方法来确定季节性参数\(s\)的值：
视觉分析： 可以对时间序列数据进行视觉分析，以识别任何季节性模式。如果在数据中观察到明显的季节性模式，则可以将季节性周期的长度用作\(s\)的值。
统计方法： 可以使用统计技术，如自相关，来识别数据中的季节性模式。\(s\)的值可以是自相关函数中观察到显著峰值的周期数。
频率分析： 可以对数据进行频率分析以识别季节性模式。\(s\)的值可以是频谱中观察到显著峰值的周期数。见
试错法： 你可以尝试不同的\(s\)值，并选择使模型与数据最佳拟合的值。
需要注意的是，\(s\)值的选择会显著影响季节性调整SES模型预测的准确性。因此，建议在选择\(s\)的最终值之前，测试不同的\(s\)值，并使用适当的评估指标评估模型的性能。
如何验证季节调整的简单指数平滑模型？
为了验证季节调整的简单指数平滑（SES季节调整）模型，可以根据分析的目标和数据的性质使用不同的理论和评估指标。
以下是一些用于验证季节调整SES模型的常见理论：
高斯-马尔可夫定理：该定理指出，如果满足某些条件，最小二乘估计量是最佳的线性无偏估计量。在季节调整SES的情况下，模型参数是通过最小二乘法估计的，因此可以使用高斯-马尔可夫定理来评估模型拟合的质量。
单根定理：该定理用于确定时间序列是否平稳。如果时间序列是非平稳的，季节调整的SES模型就不合适，因为它假设时间序列是平稳的。因此，单根定理用于评估时间序列的平稳性，并确定季节调整的SES模型是否适用。
Ljung-Box定理：该定理用于评估模型的拟合优度，并确定模型残差是否为白噪声。如果残差是白噪声，模型很好地拟合了数据，模型预测是准确的。Ljung-Box定理用于检验模型残差是否独立且不相关。
除了这些理论外，还可以使用各种评估指标，如均方根误差（MSE）、平均绝对误差（MAE）和决定系数（R²），来评估季节调整SES模型的性能，并与其他预测模型进行比较。
加载库和数据
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 方法绘制一些系列。该方法打印数据集中随机的一系列数据，对于基本的探索性数据分析非常有用。
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.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的季节性指数平滑的实现
若要了解季节性指数平滑模型函数的参数，以下列出了相关信息。有关更多信息，请访问文档。
alpha : float
平滑参数。
season_length : int
每个时间单位的观察次数。例如：24个小时数据。
alias : str
模型的自定义名称。
prediction_intervals : Optional[ConformalIntervals]
用于计算保合预测区间的信息。
默认情况下，该模型将计算原生预测
区间。加载库
from statsforecast import StatsForecast
from statsforecast.models import SeasonalExponentialSmoothing实例化模型
导入并实例化模型。设置参数有时会很复杂。罗伯·海因德曼（Rob Hyndmann）关于季节周期的文章可能对season_length有帮助。
season_length = 24 # 每小时数据
horizon = len(test) # 预测数量
models = [SeasonalExponentialSmoothing(alpha=0.8, season_length=season_length)]我们通过实例化一个新的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=[SeasonalES])让我们看看我们的季节性指数平滑模型的结果。我们可以通过以下指令进行观察：
result=sf.fitted_[0,0].model_
result{'mean': array([ 72221.8 , 73250.41 , 81213.68 , 95889.06 , 122141.46 ,
114311.9 , 105907.24 , 97934.02 , 98570.34 , 106042.58 ,
115734.17 , 133876.25 , 145319.06 , 141776.67 , 142775.06 ,
145288.16 , 150219.11 , 149963.72 , 154389.11 , 154027.47 ,
123923.64 , 102927.42 , 93966.52 , 79575.586], 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. , 80123. , 76245. ,
85949. , 102050. , 124434. , 117719. , 108679. ,
102539. , 103403. , 115897. , 130638. , 145264. ,
150694. , 149463. , 148291. , 149068. , 148820. ,
150594. , 152320. , 153663. , 131208. , 104231. ,
93096. , 82716. , 77076.6 , 75353. , 83301.8 ,
91446. , 119630.8 , 115695.8 , 110487.8 , 99595.8 ,
104028.6 , 110111.4 , 127439.6 , 141400.8 , 152114.8 ,
146912.6 , 148074.2 , 148001.6 , 146364. , 149546.8 ,
158244. , 159600.6 , 134657.6 , 111202.2 , 98779.2 ,
86635.2 , 85683.32 , 86146.6 , 90540.36 , 101861.2 ,
116678.16 , 126299.16 , 135205.56 , 135471.16 , 135405.72 ,
128574.28 , 130479.92 , 142264.16 , 156322.95 , 151382.52 ,
152602.84 , 150556.31 , 149788.8 , 147857.36 , 153668.8 ,
149420.12 , 127127.52 , 116580.44 , 97115.84 , 92215.04 ,
88384.664, 88705.32 , 90568.07 , 99004.24 , 113391.63 ,
128835.83 , 140165.11 , 149142.23 , 149145.14 , 138650.86 ,
144135.98 , 157340.83 , 164332.6 , 163700.5 , 161032.56 ,
155955.27 , 157201.77 , 157587.47 , 165409.77 , 165804.03 ,
139593.5 , 113680.086, 97299.17 , 83783.01 , 81284.93 ,
80421.06 , 88549.62 , 99632.85 , 121702.33 , 114827.164,
107585.02 , 107952.445, 107953.03 , 109782.17 , 124771.195,
140072.17 , 144962.52 , 146124.1 , 145982.52 , 147479.05 ,
147708.36 , 151845.5 , 162297.95 , 155892.81 , 135694.7 ,
108388.016, 95495.836, 80368.6 , 78924.984, 75820.21 ,
83301.92 , 98286.57 , 119816.47 , 113457.43 , 100621.01 ,
96790.49 , 96518.61 , 105304.44 , 120754.24 , 136806.44 ,
146156.5 , 140556.81 , 146976.5 , 145443.81 , 150637.67 ,
155233.1 , 161567.6 , 163186.56 , 134410.94 , 106145.6 ,
93383.164, 79489.72 , 79769. , 77652.04 , 85288.38 ,
99665.31 , 123067.3 , 115759.484, 103556.2 , 100510.1 ,
97411.72 , 107672.89 , 121150.85 , 140041.28 , 140075.3 ,
140903.36 , 142615.3 , 142360.77 , 142615.53 , 142658.62 ,
149285.52 , 146577.31 , 126038.19 , 102317.12 , 89212.63 ,
76737.945], 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 | 102317.117188 | 2017-09-21 21:00:00 |
| 214 | 89212.632812 | 2017-09-21 22:00:00 |
| 215 | 76737.945312 | 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 步。在此情况下，表示向前 30 小时。
此处的预测对象是一个新的数据框，包含一个名为模型的列和 y hat 值，以及不确定性区间的列。根据您的计算机性能，此步骤大约需要 1 分钟。
Y_hat = sf.forecast(horizon, fitted=True)
Y_hat| ds | SeasonalES | |
|---|---|---|
| unique_id | ||
| 1 | 2017-09-22 00:00:00 | 72221.796875 |
| 1 | 2017-09-22 01:00:00 | 73250.406250 |
| 1 | 2017-09-22 02:00:00 | 81213.679688 |
| ... | ... | ... |
| 1 | 2017-09-23 03:00:00 | 95889.062500 |
| 1 | 2017-09-23 04:00:00 | 122141.460938 |
| 1 | 2017-09-23 05:00:00 | 114311.898438 |
30 rows × 2 columns
values=sf.forecast_fitted_values()
values.head()| ds | y | SeasonalES | |
|---|---|---|---|
| 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 | SeasonalES | |
|---|---|---|---|
| 0 | 1 | 2017-09-22 00:00:00 | 72221.796875 |
| 1 | 1 | 2017-09-22 01:00:00 | 73250.406250 |
| 2 | 1 | 2017-09-22 02:00:00 | 81213.679688 |
| ... | ... | ... | ... |
| 27 | 1 | 2017-09-23 03:00:00 | 95889.062500 |
| 28 | 1 | 2017-09-23 04:00:00 | 122141.460938 |
| 29 | 1 | 2017-09-23 05:00:00 | 114311.898438 |
30 rows × 3 columns
Y_hat1= pd.concat([df,Y_hat])
Y_hat1| ds | y | unique_id | SeasonalES | |
|---|---|---|---|---|
| 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 | 95889.062500 |
| 28 | 2017-09-23 04:00:00 | NaN | 1 | 122141.460938 |
| 29 | 2017-09-23 05:00:00 | NaN | 1 | 114311.898438 |
246 rows × 4 columns
fig, ax = plt.subplots(1, 1)
plot_df = pd.concat([df,Y_hat]).set_index('ds')
plot_df[['y', "SeasonalES"]].plot(ax=ax, linewidth=2)
ax.set_title(' Forecast', fontsize=22)
ax.set_ylabel("Ads watched (hourly data)", fontsize=20)
ax.set_xlabel('Hourly', fontsize=20)
ax.legend(prop={'size': 15})
ax.grid(True)
预测方法与置信区间
要生成预测，请使用预测方法。
预测方法接受两个参数：预测下一个 h（代表时间范围）和 level。
h (int):表示未来的预测步数。在此情况下，预测30个小时后。
这里的预测对象是一个新的数据框，其中包含一个名为模型的列和y hat值，以及用于不确定性区间的列。
此步骤应该在1秒钟内完成。
forecast_df = sf.predict(h=horizon)
forecast_df| ds | SeasonalES | |
|---|---|---|
| unique_id | ||
| 1 | 2017-09-22 00:00:00 | 72221.796875 |
| 1 | 2017-09-22 01:00:00 | 73250.406250 |
| 1 | 2017-09-22 02:00:00 | 81213.679688 |
| ... | ... | ... |
| 1 | 2017-09-23 03:00:00 | 95889.062500 |
| 1 | 2017-09-23 04:00:00 | 122141.460938 |
| 1 | 2017-09-23 05:00:00 | 114311.898438 |
30 rows × 2 columns
我们可以使用pandas函数pd.concat()将预测结果与历史数据连接在一起，然后能够使用这个结果进行绘图。
df_plot= pd.concat([df, forecast_df]).set_index('ds')
df_plot| y | unique_id | SeasonalES | |
|---|---|---|---|
| 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 | 95889.062500 |
| 2017-09-23 04:00:00 | NaN | NaN | 122141.460938 |
| 2017-09-23 05:00:00 | NaN | NaN | 114311.898438 |
246 rows × 3 columns
现在让我们可视化我们的预测结果和时间序列的历史数据。
sns.lineplot(df_plot,x="ds", y="y", label="Actual")
sns.lineplot(df_plot, x="ds", y="SeasonalES", label="SeasonalES")
plt.title("Ads watched (hourly data)")
plt.show()
交叉验证
在之前的步骤中，我们使用历史数据预测未来。然而，为了评估其准确性，我们还希望知道模型在过去的表现。为了评估您的模型在数据上的准确性和稳健性，进行交叉验证。
对于时间序列数据，交叉验证是通过在历史数据上定义滑动窗口，并预测其后的一段时间来完成的。这种形式的交叉验证使我们能够在更广泛的时间实例中更好地估计模型的预测能力，同时保持训练集中数据的连续性，这正是我们的模型所要求的。
以下图展示了这种交叉验证策略：
执行时间序列交叉验证
时间序列模型的交叉验证被认为是一种最佳实践，但大多数实现速度都很慢。statsforecast库将交叉验证实现为分布式操作，使得执行过程更加高效。如果您有大型数据集，还可以使用Ray、Dask或Spark在分布式集群中执行交叉验证。
在本例中，我们希望评估每个模型在过去5个月的表现(n_windows=5)，每隔两个月进行一次预测(step_size=12)。根据您的计算机，此步骤大约需要1分钟。
StatsForecast类的cross_validation方法接受以下参数。
df:训练数据框h (int):代表向未来预测的h步。在本例中，预测30小时。step_size (int):每个窗口之间的步长。换句话说：您希望多频繁地运行预测过程。n_windows(int):用于交叉验证的窗口数。换句话说：您希望评估过去多少次预测过程。
crossvalidation_df = sf.cross_validation(df=df,
h=horizon,
step_size=12,
n_windows=3)crossvaldation_df对象是一个新的数据框，包括以下列：
unique_id:索引。如果你不喜欢使用索引，直接运行 crossvalidation_df.resetindex()。ds:日期戳或时间索引cutoff:n_windows的最后日期戳或时间索引。y:真实值"model":包含模型名称和拟合值的列。
crossvalidation_df| ds | cutoff | y | SeasonalES | |
|---|---|---|---|---|
| unique_id | ||||
| 1 | 2017-09-19 18:00:00 | 2017-09-19 17:00:00 | 161385.0 | 162297.953125 |
| 1 | 2017-09-19 19:00:00 | 2017-09-19 17:00:00 | 165010.0 | 155892.812500 |
| 1 | 2017-09-19 20:00:00 | 2017-09-19 17:00:00 | 134090.0 | 135694.703125 |
| ... | ... | ... | ... | ... |
| 1 | 2017-09-21 21:00:00 | 2017-09-20 17:00:00 | 103080.0 | 106145.601562 |
| 1 | 2017-09-21 22:00:00 | 2017-09-20 17:00:00 | 95155.0 | 93383.164062 |
| 1 | 2017-09-21 23:00:00 | 2017-09-20 17:00:00 | 80285.0 | 79489.718750 |
90 rows × 4 columns
评估模型
现在我们可以使用适当的准确性指标来计算预测的准确性。在这里，我们将使用均方根误差（RMSE）。为此，我们首先需要安装datasetsforecast，这是由Nixtla开发的一个Python库，包含一个计算RMSE的函数。
%%capture
!pip install datasetsforecastfrom datasetsforecast.losses import rmse计算RMSE的函数需要两个参数：
- 实际值。
- 预测值，在这个例子中是
季节性指数平滑模型。
rmse = rmse(crossvalidation_df['y'], crossvalidation_df["SeasonalES"])
print("RMSE using cross-validation: ", rmse)RMSE using cross-validation: 7127.6504感谢
我们要感谢Naren Castellon撰写本教程。
参考文献
- Changquan Huang • Alla Petukhina. Springer系列（2022）。使用Python进行应用时间序列分析与预测。
- Ivan Svetunkov. 使用增强型动态自适应模型（ADAM）进行预测与分析。
- James D. Hamilton。时间序列分析，普林斯顿大学出版社，普林斯顿，新泽西，第一版，1994年。
- Nixtla参数。
- Pandas可用频率。
- Rob J. Hyndman 和 George Athanasopoulos（2018）。 “预测原则与实践，时间序列交叉验证”。。
- 季节性周期 - Rob J Hyndman。
Give us a ⭐ on Github