目录

引言

指数平滑于 20 世纪 50 年代末被提出(Brown, 1959; Holt, 1957; Winters, 1960),并催生了一些最成功的预测方法。使用指数平滑方法产生的预测是过去观测值的加权平均,权重随观测值的时效性呈指数衰减。换句话说,越近期的观测值具有越高的关联权重。这个框架能够快速为广泛的时间序列生成可靠的预测,这在行业应用中具有巨大优势和重要意义。

简单指数平滑模型是一种用于时间序列分析的方法,用于基于历史观测值预测未来值。该模型基于这样一个理念:时间序列的未来值将受到过去值的影响,且过去值的影响力随时间回溯呈指数衰减。

简单指数平滑模型使用一个平滑因子,这是一个介于 0 和 1 之间的数字,表示在预测未来值时过去观测值的相对重要性。值为 1 表示所有过去的观测值都具有同等重要性,而值为 0 表示只考虑最新的观测值。

简单指数平滑模型可以用数学形式表示为

y^T+1T=αyT+α(1α)yT1+α(1α)2yT2+,\hat{y}_{T+1|T} = \alpha y_T + \alpha(1-\alpha) y_{T-1} + \alpha(1-\alpha)^2 y_{T-2}+ \cdots,

其中 yTy_T 是周期 tt 的观测值,y^T+1T\hat{y}_{T+1|T} 是下一周期的预测值,y (t1)(t-1) 是前一周期的观测值,α\alpha 是平滑因子。

简单指数平滑模型因其简单性和易用性而被广泛使用。然而,它也有其局限性,因为它无法捕捉数据中复杂的模式,不适用于具有趋势或季节性模式的时间序列。

构建简单指数平滑模型

指数平滑方法中最简单的一种自然称为简单指数平滑(SES)。这种方法适用于预测没有明显趋势或季节性模式的数据。

使用朴素方法,未来的所有预测值都等于序列的最后一个观测值,y^T+hT=yT,\hat{y}_{T+h|T} = y_{T},

其中 h=1,2,h=1,2,\dots。因此,朴素方法假设最新观测值是唯一重要的,所有之前的观测值都不提供未来信息。这可以被视为一种加权平均,其中所有权重都赋给了最后一个观测值。

使用平均方法,所有未来预测值都等于观测数据的简单平均值,y^T+hT=1Tt=1Tyt,\hat{y}_{T+h|T} = \frac1T \sum_{t=1}^T y_t,

其中 h=1,2,h=1,2,\dots。因此,平均方法假设所有观测值都具有同等重要性,并在生成预测时赋予它们相等的权重。

我们通常希望介于这两个极端之间。例如,对于较近期的观测值,赋予比遥远过去的观测值更大的权重可能是合理的。这正是简单指数平滑背后的概念。预测是使用加权平均计算的,其中权重随观测值时效性呈指数衰减——最小的权重与最旧的观测值相关。

其中 0α10 \le \alpha \le 1 是平滑参数。时间 T+1T+1 的一步超前预测是序列 y1,,yTy_1,\dots,y_T 中所有观测值的加权平均值。权重的衰减速度由参数 α\alpha 控制。

对于介于 0 和 1 之间的任何 α\alpha,附加在观测值上的权重随时间回溯呈指数衰减,因此得名“指数平滑”。如果 α\alpha 较小(即接近 0),则更多权重赋给了更遥远过去的观测值。如果 α\alpha 较大(即接近 1),则更多权重赋给了较近期的观测值。在极端情况下,当 α=1\alpha=1 时,y^T+1T=yT\hat{y}_{T+1|T}=y_T 并且预测值等于朴素预测值。

我们给出简单指数平滑的两种等价形式,每种形式都可以得到预测方程 (1)。

加权平均形式

时间 T+1T+1 的预测值等于最新观测值 yTy_T 与前一次预测值 y^TT1\hat{y}_{T|T-1} 的加权平均值

y^T+1T=αyT+(1α)y^TT1,\hat{y}_{T+1|T} = \alpha y_T + (1-\alpha) \hat{y}_{T|T-1},

其中 0α10 \le \alpha \le 1 是平滑参数。类似地,我们可以将拟合值写为 y^t+1t=αyt+(1α)y^tt1,\hat{y}_{t+1|t} = \alpha y_t + (1-\alpha) \hat{y}_{t|t-1},

对于 t=1,,Tt=1,\dots,T。(请注意,拟合值只是训练数据的一步超前预测。)

该过程必须从某个地方开始,所以我们将时间 1 的第一个拟合值记为 0\ell_{0}(我们将对其进行估计)。然后

将每个方程代入后续方程,我们得到

对于大的 TT,最后一项变得非常小。因此,加权平均形式得到了与预测方程 (1) 相同的结论。

分量形式

另一种表示形式是分量形式。对于简单指数平滑,包含的唯一分量是水平 t\ell_{t}。指数平滑方法的分量形式表示包含一个预测方程和每个包含分量的平滑方程。简单指数平滑的分量形式由下式给出

其中 t\ell_{t} 是时间 tt 的序列水平(或平滑值)。设置 h=1h=1 得到拟合值,而设置 t=Tt=T 得到训练数据之外的真实预测值。

预测方程表明,时间 t+1t+1 的预测值是时间 tt 的估计水平。水平的平滑方程(通常称为水平方程)给出每个周期 tt 的序列估计水平。

如果在平滑方程中用 y^t+1t\hat{y}_{t+1|t} 替换 t\ell_{t},用 y^tt1\hat{y}_{t|t-1} 替换 t1\ell_{t-1},我们将恢复简单指数平滑的加权平均形式。

简单指数平滑的分量形式本身并不是特别有用,但当我们开始添加其他分量时,它将成为最容易使用的形式。

平稳预测

简单指数平滑具有“平稳”预测函数

y^T+hT=y^T+1T=T,h=2,3,.\hat{y}_{T+h|T} = \hat{y}_{T+1|T}=\ell_T, \qquad h=2,3,\dots.

也就是说,所有预测值都取相同的值,等于最后一个水平分量。请记住,这些预测仅适用于时间序列没有趋势或季节性分量的情况。

加载库和数据

提示

需要 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()
时间广告
02017-09-13T00:00:0080115
12017-09-13T01:00:0079885
22017-09-13T02:00:0089325
32017-09-13T03:00:00101930
42017-09-13T04:00:00121630

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
dsyunique_id
02017-09-13T00:00:00801151
12017-09-13T01:00:00798851
22017-09-13T02:00:00893251
2132017-09-21T21:00:001030801
2142017-09-21T22:00:00951551
2152017-09-21T23:00:00802851
print(df.dtypes)
ds           object
y             int64
unique_id    object
dtype: object
df["ds"] = pd.to_datetime(df["ds"])

使用 plot 方法探索数据

使用 StatsForecast 类中的 plot 方法绘制一些时间序列。此方法会打印数据集中随机选择的时间序列,对于基础 EDA 很有用。

from statsforecast import StatsForecast

StatsForecast.plot(df)

自相关图

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.savefig("Gráfico de Densidad y qq")
plt.show();

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

我们将数据拆分为数据集

  1. 用于训练我们的 简单指数平滑 (SES) 的数据。
  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))

在 StatsForecast 中实现 SimpleExponentialSmoothing

加载库

from statsforecast import StatsForecast
from statsforecast.models import SimpleExponentialSmoothing

实例化模型

我们将为不同的 alpha 值构建不同的模型。

horizon = len(test)
# We call the model that we are going to use
models = [SimpleExponentialSmoothing(alpha=0.1, alias="SES01"),
          SimpleExponentialSmoothing(alpha=0.5,alias="SES05"),
          SimpleExponentialSmoothing(alpha=0.8,alias="SES08")
          ]

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

models:模型列表。从模型中选择您想要的模型并导入。

  • freq:指示数据频率的字符串。(请参阅 pandas 可用的频率。)

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

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

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

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

拟合模型

sf.fit(df=train)
StatsForecast(models=[SES01,SES05,SES08])

让我们看看我们的简单 简单指数平滑模型 (SES) 的结果。我们可以通过以下指令查看

result01=sf.fitted_[0,0].model_
result05=sf.fitted_[0,1].model_
result08=sf.fitted_[0,2].model_
result01
{'mean': array([126112.90072589]),
 'fitted': array([       nan,  80115.   ,  80092.   ,  81015.3  ,  83106.77 ,
         86959.09 ,  89910.69 ,  91569.12 ,  92691.7  ,  94228.03 ,
         96417.73 ,  99878.96 , 104793.06 , 110072.76 , 114136.98 ,
        117652.78 , 120897.5  , 123285.75 , 126026.18 , 129807.56 ,
        133450.3  , 134057.28 , 131241.05 , 127794.945, 123267.445,
        118953.2  , 114591.38 , 111642.74 , 110686.47 , 112131.32 ,
        112721.19 , 112371.57 , 111381.914, 110467.73 , 111004.95 ,
        112958.45 , 116095.11 , 119382.6  , 122359.336, 124927.41 ,
        127315.664, 129567.1  , 131667.39 , 133444.66 , 135152.19 ,
        134549.97 , 131476.47 , 127546.32 , 123068.19 , 118392.875,
        114066.586, 110923.92 , 108711.03 , 109682.93 , 110233.64 ,
        110304.27 , 109159.84 , 108662.36 , 108662.625, 110460.36 ,
        113457.83 , 117359.05 , 120250.64 , 123027.58 , 125498.32 ,
        127523.484, 129699.64 , 132702.17 , 135540.45 , 135538.4  ,
        133279.06 , 129971.164, 125735.55 , 121945.49 , 118635.445,
        116006.9  , 114852.71 , 114961.44 , 116360.3  , 118862.766,
        121420.484, 123603.44 , 124562.09 , 125229.88 , 126954.9  ,
        129996.91 , 132247.22 , 134396.   , 136075.89 , 137532.81 ,
        138523.03 , 139923.22 , 140618.4  , 139081.06 , 136965.45 ,
        132938.9  , 129006.016, 125011.414, 121444.77 , 118357.8  ,
        116351.016, 115972.914, 117322.625, 119730.86 , 123013.77 ,
        125970.4  , 127490.36 , 129496.32 , 132657.69 , 136025.42 ,
        139100.88 , 141504.8  , 143084.81 , 144681.83 , 146215.64 ,
        148428.58 , 150575.72 , 149789.16 , 146105.73 , 141229.66 ,
        135274.2  , 129697.77 , 124563.   , 120911.2  , 118799.08 ,
        119297.17 , 118499.95 , 116593.96 , 114700.06 , 112995.555,
        111952.5  , 112750.25 , 115050.73 , 117557.66 , 119974.89 ,
        122199.4  , 124515.46 , 126597.414, 128978.67 , 132232.81 ,
        134351.03 , 134387.92 , 131655.62 , 127994.57 , 123146.61 ,
        118665.445, 114265.91 , 111038.31 , 109729.484, 110691.03 ,
        110933.43 , 109728.086, 108155.28 , 106705.75 , 106453.68 ,
        107783.305, 110603.98 , 114189.08 , 116686.67 , 119740.51 ,
        122259.95 , 125170.96 , 128261.86 , 131574.17 , 134917.77 ,
        134834.98 , 131909.98 , 128004.484, 123131.04 , 118815.94 ,
        114745.34 , 111849.305, 110665.375, 111986.836, 112421.66 ,
        111608.49 , 110591.64 , 109295.98 , 109192.875, 110398.59 ,
        113443.734, 115954.86 , 118458.375, 120765.04 , 122847.53 ,
        124623.78 ], dtype=float32)}

如我们所见,上面获得的结果是一个字典输出,为了从字典中提取每个元素,我们将使用 .get() 函数提取元素,然后将其保存在 pd.DataFrame() 中。

fitted=pd.DataFrame(result01.get("fitted"), columns=["fitted01"])
fitted["fitted05"]=result05.get("fitted")
fitted["fitted08"]=result08.get("fitted")
fitted["ds"]=df["ds"]
fitted
fitted01fitted05fitted08ds
0NaNNaNNaN2017-09-13 00:00:00
180115.00000080115.0080115.0000002017-09-13 01:00:00
280092.00000080000.0079931.0000002017-09-13 02:00:00
183120765.039062139195.00141302.8281252017-09-20 15:00:00
184122847.531250140392.50141532.5625002017-09-20 16:00:00
185124623.781250140501.25140794.5156252017-09-20 17:00:00
sns.lineplot(df, x="ds", y="y", label="Actual", linewidth=2)
sns.lineplot(fitted,x="ds", y="fitted01", label="Fitted01", linestyle="--", )
sns.lineplot(fitted, x="ds", y="fitted05", label="Fitted05", color="lime")
sns.lineplot(fitted, x="ds", y="fitted08", label="Fitted08")
plt.title("Ads watched (hourly data)");
plt.show()

预测方法

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

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

forecast 方法接受两个参数:预测未来 h(范围)和 level(置信水平)。

  • h (整数):表示预测未来 h 个步长。在本例中,预测未来 30 小时。

这里的 forecast 对象是一个新的数据框,其中包含模型名称和 y hat 值的列,以及用于不确定性区间的列。根据您的计算机性能,此步骤大约需要 1 分钟。

# Prediction
Y_hat = sf.forecast(df=train, h=horizon, fitted=True)
Y_hat.head()
unique_iddsSES01SES05SES08
012017-09-20 18:00:00126112.898438140008.125139770.90625
112017-09-20 19:00:00126112.898438140008.125139770.90625
212017-09-20 20:00:00126112.898438140008.125139770.90625
312017-09-20 21:00:00126112.898438140008.125139770.90625
412017-09-20 22:00:00126112.898438140008.125139770.90625
values=sf.forecast_fitted_values()
values.head()
unique_iddsySES01SES05SES08
012017-09-13 00:00:0080115.0NaNNaNNaN
112017-09-13 01:00:0079885.080115.00000080115.0080115.000000
212017-09-13 02:00:0089325.080092.00000080000.0079931.000000
312017-09-13 03:00:00101930.081015.29687584662.5087446.203125
412017-09-13 04:00:00121630.083106.77343893296.2599033.242188

predict 方法

要生成预测,请使用 predict 方法。

predict 方法接受两个参数:预测未来 h(范围)。* h (整数):表示预测未来 hh 个步长。在本例中,预测未来 30 小时。

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

此步骤应不到 1 秒。

forecast_df = sf.predict(h=horizon) 
forecast_df
unique_iddsSES01SES05SES08
012017-09-20 18:00:00126112.898438140008.125139770.90625
112017-09-20 19:00:00126112.898438140008.125139770.90625
212017-09-20 20:00:00126112.898438140008.125139770.90625
2712017-09-21 21:00:00126112.898438140008.125139770.90625
2812017-09-21 22:00:00126112.898438140008.125139770.90625
2912017-09-21 23:00:00126112.898438140008.125139770.90625
sf.plot(train, forecast_df)

交叉验证

在前面的步骤中,我们使用历史数据来预测未来。然而,为了评估其准确性,我们还想了解模型在过去可能如何表现。为了评估模型在数据上的准确性和稳健性,请执行交叉验证。

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

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

执行时间序列交叉验证

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

在本例中,我们希望评估每个模型在最近 30 小时 (n_windows=) 的性能,每隔一个月预测一次 (step_size=30)。根据您的计算机性能,此步骤大约需要 1 分钟。

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

  • df:训练数据框

  • h (整数):表示预测未来 h 个步长。在本例中,预测未来 30 小时。

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

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

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

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

  • unique_id:时间序列标识符
  • ds:日期戳或时间索引
  • cutoffn_windows 的最后一个日期戳或时间索引。
  • y:真实值
  • model:包含模型名称和拟合值的列。
crossvalidation_df
unique_iddscutoffySES01SES05SES08
012017-09-18 06:00:002017-09-18 05:00:0099440.0118499.953125109816.250112747.695312
112017-09-18 07:00:002017-09-18 05:00:0097655.0118499.953125109816.250112747.695312
212017-09-18 08:00:002017-09-18 05:00:0097655.0118499.953125109816.250112747.695312
8712017-09-21 21:00:002017-09-20 17:00:00103080.0126112.898438140008.125139770.906250
8812017-09-21 22:00:002017-09-20 17:00:0095155.0126112.898438140008.125139770.906250
8912017-09-21 23:00:002017-09-20 17:00:0080285.0126112.898438140008.125139770.906250

模型评估

现在我们将使用预测结果来评估我们的模型,我们将使用不同类型的指标 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,
)
unique_id指标SES01SES05SES08
01mae25173.93958329390.87500029311.802083
11mape0.2550880.3164400.315339
21mase3.1102883.6312983.621528
31rmse28923.39538136184.34086936027.710540
41smape0.1099720.1248030.124542

参考

  1. Changquan Huang • Alla Petukhina. Springer series (2022). Applied Time Series Analysis and Forecasting with Python.
  2. Ivan Svetunkov. 使用增强动态自适应模型 (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.