简单指数平滑优化模型

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

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 的 SimpleExponentialSmoothingOptimized Model 的逐步指南。

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

目录

• 介绍
• 简单指数平滑优化模型
• 加载库和数据
• 使用绘图方法探索数据
• 将数据分为训练集和测试集
• 使用StatsForecast实现简单指数平滑优化
• 交叉验证
• 模型评估
• 参考文献

介绍

简单指数平滑优化（SES 优化）是一种用于预测单变量时间序列未来值的预测模型。它是简单指数平滑（SES）方法的一种变体，采用优化方法更准确地估计模型参数。

SES 优化方法使用单一平滑参数来估计时间序列数据中的趋势和季节性。该模型试图通过优化算法最小化预测值与训练样本中实际值之间的均方误差（MSE）。

SES 优化方法特别适用于具有强趋势和季节性模式的时间序列，或对于数据噪声较大的时间序列。然而，重要的是要注意该模型假设时间序列是平稳的，并且数据的变化是随机的，且数据中不存在非随机模式。如果这些假设不成立，SES 优化模型可能表现不佳，可能需要其他预测方法。

简单指数平滑模型

最简单的指数平滑方法自然被称为简单指数平滑（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\) 的一步ahead预测是系列 \(y_1,\dots,y_T\) 中所有观测值的加权平均。权重下降的速率由参数 \(\alpha\) 控制。

对于任何在 0 和 1 之间的 \(\alpha\)，附加到观测值的权重随着时间的推移而指数下降，因此得名“指数平滑”。如果 \(\alpha\) 很小（即接近 0），则给来自更远过去的观测值更多的权重。如果 \(\alpha\) 很大（即接近 1），则给最近的观测值更多的权重。在极端情况下，当 \(\alpha=1\) 时，\(\hat{y}_{T+1|T}=y_T\)，预测等于天真的预测。

优化

每种指数平滑方法的应用都需要选择平滑参数和初始值。特别是对于简单指数平滑，我们需要选择 \(\alpha\) 和 \(\ell_0\) 的值。一旦知道这些值，就可以根据数据计算所有的预测。对于接下来的方法，通常有多个平滑参数和多个初始组件需要选择。

在某些情况下，平滑参数可以以主观方式选择——预测者根据以往经验指定平滑参数的值。然而，获取未知参数值的更可靠和客观的方式是从观察到的数据中进行估计。

从回归模型中，我们通过最小化平方残差的和（通常称为 SSE 或“平方误差之和”）来估计回归模型的系数。类似地，可以通过最小化 SSE 来估计任何指数平滑方法的未知参数和初始值。残差被定义为 \(e_t=y_t - \hat{y}_{t|t-1}\)，其中 \(t=1,\dots,T\)。因此，我们寻找未知参数和初始值的值，以最小化

\[\begin{equation} \text{SSE}=\sum_{t=1}^T(y_t - \hat{y}_{t|t-1})^2=\sum_{t=1}^Te_t^2. \tag{2} \end{equation}\]

与回归情况不同（我们有公式返回最小化 SSE 的回归系数值），这涉及非线性最小化问题，我们需要使用优化工具来解决它。

加载库和数据

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)

自相关图

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

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

让我们将数据分为几组

1. 用于训练我们的简单指数平滑优化模型的数据
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="--")
sns.lineplot(test, x="ds", y="y", label="Test")
plt.title("Ads watched (hourly data)");
plt.show()

使用StatsForecast实现SimpleExponentialSmoothingOptimized

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

alpha : float
    平滑参数。
alias : str
    模型的自定义名称。
prediction_intervals : Optional[ConformalIntervals]
    用于计算符合预测区间的信息。
    默认情况下，模型将计算原生预测
    区间。

加载库

from statsforecast import StatsForecast
from statsforecast.models import SimpleExponentialSmoothingOptimized

实例化模型

horizon = len(test) # 预测数量

models = [SimpleExponentialSmoothingOptimized()] # 乘法的  加法的

我们通过实例化一个新的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=[SESOpt])

让我们看看我们的 简单指数平滑优化模型 的结果。我们可以通过以下指令来观察它：

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

{'mean': array([80434.516], dtype=float32),
 'fitted': array([       nan,  80115.   ,  79887.3  ,  89230.625, 101803.01 ,
        121431.73 , 116524.57 , 106595.3  , 102833.   , 108002.78 ,
        116043.78 , 130880.14 , 148838.6  , 157502.48 , 150782.88 ,
        149309.88 , 150092.1  , 144833.12 , 150631.44 , 163707.92 ,
        166209.73 , 139786.89 , 106233.92 ,  96874.54 ,  82663.55 ,
         80150.38 ,  75383.16 ,  85007.78 , 101909.28 , 124902.74 ,
        118098.73 , 109313.734, 102543.39 , 102243.03 , 115704.03 ,
        130391.64 , 144185.67 , 148922.16 , 149147.72 , 148051.08 ,
        148802.4  , 149819.72 , 150562.5  , 149451.22 , 150509.31 ,
        129343.8  , 104070.29 ,  92293.95 ,  82860.29 ,  76380.45 ,
         75142.51 ,  82565.02 ,  88732.7  , 118133.02 , 115219.43 ,
        110982.8  ,  98981.23 , 104132.96 , 108619.68 , 126459.8  ,
        140295.25 , 152348.25 , 146335.73 , 148003.16 , 147737.69 ,
        145769.88 , 149249.84 , 159620.25 , 161070.36 , 135775.5  ,
        113173.305, 100329.734,  87742.15 ,  87834.07 ,  88834.89 ,
         92314.85 , 104343.5  , 115824.03 , 128818.74 , 141259.34 ,
        144408.19 , 143261.58 , 133290.72 , 131260.5  , 142367.81 ,
        157224.92 , 152547.25 , 153723.12 , 151220.28 , 150650.75 ,
        147467.16 , 152474.42 , 146931.   , 125461.86 , 118000.37 ,
         96913.   ,  93643.03 ,  89105.83 ,  89342.61 ,  90562.68 ,
         98212.73 , 112426.43 , 129299.56 , 141283.95 , 152447.23 ,
        152578.67 , 141284.1  , 147487.34 , 160973.77 , 166281.39 ,
        166775.02 , 163176.34 , 157363.72 , 159038.1  , 160010.19 ,
        168261.66 , 169883.61 , 142981.73 , 113255.266,  97504.1  ,
         81833.29 ,  79533.234,  78361.836,  87948.17 ,  99671.58 ,
        123538.914, 111447.14 ,  99560.07 ,  97674.05 ,  97655.19 ,
        102515.9  , 119755.86 , 135595.02 , 140074.75 , 141713.45 ,
        142214.94 , 145328.55 , 145334.94 , 150359.25 , 161408.39 ,
        153494.94 , 134907.75 , 107343.43 ,  95167.984,  79671.53 ,
         78348.37 ,  74706.78 ,  81917.164,  97789.67 , 119129.445,
        113175.14 ,  99022.95 ,  94050.23 ,  93663.9  , 104079.79 ,
        119593.3  , 135826.03 , 146348.7  , 139236.84 , 147145.12 ,
        144957.1  , 151305.88 , 156032.27 , 161331.47 , 164973.22 ,
        134398.83 , 105873.14 ,  92985.18 ,  79407.15 ,  79974.27 ,
         78128.64 ,  85708.44 ,  99866.984, 123639.87 , 116408.05 ,
        104411.18 , 101469.71 ,  97673.34 , 108159.086, 121119.09 ,
        140652.69 , 138575.98 , 140965.86 , 141519.4  , 141589.3  ,
        140619.8  , 139526.05 , 146148.11 , 142462.23 , 124130.17 ,
        101587.7  ,  88304.18 ,  76172.54 ,  70393.375,  72132.44 ,
         80114.375,  94796.695, 121638.87 , 114026.89 , 106570.32 ,
         97382.805,  98845.23 , 105567.1  , 114291.87 , 132154.56 ,
        146485.25 , 142039.9  , 142807.25 , 145987.88 , 152058.67 ,
        151792.69 , 155626.28 , 155887.36 , 123719.92 , 103286.4  ,
         95236.31 ], 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	80115.000000	2017-09-13 01:00:00
2	79887.296875	2017-09-13 02:00:00
...	...	...
213	123719.921875	2017-09-21 21:00:00
214	103286.398438	2017-09-21 22:00:00
215	95236.312500	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	SESOpt
unique_id
1	2017-09-22 00:00:00	80434.515625
1	2017-09-22 01:00:00	80434.515625
1	2017-09-22 02:00:00	80434.515625
...	...	...
1	2017-09-23 03:00:00	80434.515625
1	2017-09-23 04:00:00	80434.515625
1	2017-09-23 05:00:00	80434.515625

30 rows × 2 columns

让我们可视化拟合值

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

	ds	y	SESOpt
unique_id
1	2017-09-13 00:00:00	80115.0	NaN
1	2017-09-13 01:00:00	79885.0	80115.000000
1	2017-09-13 02:00:00	89325.0	79887.296875
1	2017-09-13 03:00:00	101930.0	89230.625000
1	2017-09-13 04:00:00	121630.0	101803.007812

StatsForecast.plot(values)

sf.forecast(h=horizon)

	ds	SESOpt
unique_id
1	2017-09-22 00:00:00	80434.515625
1	2017-09-22 01:00:00	80434.515625
1	2017-09-22 02:00:00	80434.515625
...	...	...
1	2017-09-23 03:00:00	80434.515625
1	2017-09-23 04:00:00	80434.515625
1	2017-09-23 05:00:00	80434.515625

30 rows × 2 columns

Y_hat=Y_hat.reset_index()
Y_hat.head(10)

	unique_id	ds	SESOpt
0	1	2017-09-22 00:00:00	80434.515625
1	1	2017-09-22 01:00:00	80434.515625
2	1	2017-09-22 02:00:00	80434.515625
...	...	...	...
7	1	2017-09-22 07:00:00	80434.515625
8	1	2017-09-22 08:00:00	80434.515625
9	1	2017-09-22 09:00:00	80434.515625

10 rows × 3 columns

预测方法与置信区间

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

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

• h (int): 表示预测的未来 h 步。在这种情况下，是 30 小时前。

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

此步骤应少于 1 秒。

forecast_df = sf.predict(h=horizon) 

forecast_df

	ds	SESOpt
unique_id
1	2017-09-22 00:00:00	80434.515625
1	2017-09-22 01:00:00	80434.515625
1	2017-09-22 02:00:00	80434.515625
...	...	...
1	2017-09-23 03:00:00	80434.515625
1	2017-09-23 04:00:00	80434.515625
1	2017-09-23 05:00:00	80434.515625

30 rows × 2 columns

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

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

	y	unique_id	SESOpt
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	80434.515625
2017-09-23 04:00:00	NaN	NaN	80434.515625
2017-09-23 05:00:00	NaN	NaN	80434.515625

246 rows × 3 columns

# 绘制数据及其指数平滑后的数据
plt.plot(df_plot['y'],label="Actual", linewidth=2.5)
plt.plot(df_plot['SESOpt'], label="SESOpt") # '-', '--', '-.', ':',
plt.title("Ads watched (hourly data)");
plt.ylabel("")
plt.legend()
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步。这种情况下是30小时。

• 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	SESOpt
unique_id
1	2017-09-18 06:00:00	2017-09-18 05:00:00	99440.0	111447.140625
1	2017-09-18 07:00:00	2017-09-18 05:00:00	97655.0	111447.140625
1	2017-09-18 08:00:00	2017-09-18 05:00:00	97655.0	111447.140625
...	...	...	...	...
1	2017-09-21 21:00:00	2017-09-20 17:00:00	103080.0	139526.046875
1	2017-09-21 22:00:00	2017-09-20 17:00:00	95155.0	139526.046875
1	2017-09-21 23:00:00	2017-09-20 17:00:00	80285.0	139526.046875

90 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的函数接受两个参数：

1. 实际值。
2. 预测值，在此情况下，为“简单指数平滑优化模型”。

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

RMSE using cross-validation:  30098.377

感谢

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

参考文献

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

Give us a ⭐ on Github