目录

介绍

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

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

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

简单指数平滑模型

指数平滑方法中最简单的方法自然称为简单指数平滑 (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,预测值等于朴素预测值。

优化

应用每种指数平滑方法都需要选择平滑参数和初始值。特别是对于简单指数平滑,我们需要选择 α\alpha0\ell_0 的值。一旦知道这些值,就可以从数据计算所有预测。对于后续方法,通常需要选择多个平滑参数和多个初始分量。

在某些情况下,平滑参数可以主观选择——预测者根据先前的经验指定平滑参数的值。然而,一种更可靠和客观的方法来获取未知参数的值是从观测数据中估计它们。

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

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

加载库和数据

提示

需要 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.head()
dsyunique_id
02017-09-13T00:00:00801151
12017-09-13T01:00:00798851
22017-09-13T02:00:00893251
32017-09-13T03:00:001019301
42017-09-13T04:00:001216301
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. 用于训练我们的 Simple Exponential Smoothing Optimized Model 的数据
  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 实现 SimpleExponentialSmoothingOptimized

加载库

from statsforecast import StatsForecast
from statsforecast.models import SimpleExponentialSmoothingOptimized

实例化模型

horizon = len(test) # number of predictions

models = [SimpleExponentialSmoothingOptimized()] # multiplicative   additive

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

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

  • freq: 表示数据频率的字符串。(参见 pandas 可用频率。)

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

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

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

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

拟合模型

sf.fit(df=train)
StatsForecast(models=[SESOpt])

让我们看看 Simple Exponential Smoothing Optimized model 的结果。我们可以使用以下指令观察它

result=sf.fitted_[0,0].model_
result
{'mean': array([139526.04792941]),
 '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  ], dtype=float32)}

现在让我们可视化模型的残差。

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

fitted=pd.DataFrame(result.get("fitted"), columns=["fitted"])
fitted["ds"]=df["ds"]
fitted
fittedds
0NaN2017-09-13 00:00:00
180115.0000002017-09-13 01:00:00
279887.2968752017-09-13 02:00:00
183141519.4062502017-09-20 15:00:00
184141589.2968752017-09-20 16:00:00
185140619.7968752017-09-20 17:00:00
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 不存储拟合值,并且在分布式环境中具有很高的可伸缩性。

forecast 方法接受两个参数:预测接下来的 h(预测范围)和 level

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

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

# Prediction
Y_hat = sf.forecast(df=train, h=horizon, fitted=True)
Y_hat
unique_iddsSESOpt
012017-09-20 18:00:00139526.046875
112017-09-20 19:00:00139526.046875
212017-09-20 20:00:00139526.046875
2712017-09-21 21:00:00139526.046875
2812017-09-21 22:00:00139526.046875
2912017-09-21 23:00:00139526.046875

让我们可视化拟合值

values=sf.forecast_fitted_values()
values.head()
unique_iddsySESOpt
012017-09-13 00:00:0080115.0NaN
112017-09-13 01:00:0079885.080115.000000
212017-09-13 02:00:0089325.079887.296875
312017-09-13 03:00:00101930.089230.625000
412017-09-13 04:00:00121630.0101803.007812

带置信区间的预测方法

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

predict 方法接受两个参数:预测接下来的 h(表示预测范围)和 level

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

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

此步骤应花费不到 1 秒。

forecast_df = sf.predict(h=horizon) 
forecast_df
unique_iddsSESOpt
012017-09-20 18:00:00139526.046875
112017-09-20 19:00:00139526.046875
212017-09-20 20:00:00139526.046875
2712017-09-21 21:00:00139526.046875
2812017-09-21 22:00:00139526.046875
2912017-09-21 23:00:00139526.046875
sf.plot(train, 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)

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

  • unique_id: 索引。如果您不喜欢使用索引,只需运行 crossvalidation_df.reset_index()
  • ds: 时间戳或时间索引
  • cutoff: n_windows 的最后一个时间戳或时间索引。
  • y: 真实值
  • model: 包含模型名称和拟合值的列。
crossvalidation_df
unique_iddscutoffySESOpt
012017-09-18 06:00:002017-09-18 05:00:0099440.0111447.140625
112017-09-18 07:00:002017-09-18 05:00:0097655.0111447.140625
212017-09-18 08:00:002017-09-18 05:00:0097655.0111447.140625
8712017-09-21 21:00:002017-09-20 17:00:00103080.0139526.046875
8812017-09-21 22:00:002017-09-20 17:00:0095155.0139526.046875
8912017-09-21 23:00:002017-09-20 17:00:0080285.0139526.046875

模型评估

现在我们将使用预测结果来评估我们的模型,我们将使用不同类型的指标 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_idmetricSESOpt
01mae29230.182292
11mape0.314203
21mase3.611444
31rmse35866.963426
41smape0.124271

参考

  1. Changquan Huang • Alla Petukhina. Springer 系列 (2022)。Applied Time Series Analysis and Forecasting with Python。
  2. Ivan Svetunkov. 使用增强动态自适应模型 (ADAM) 进行预测和分析
  3. James D. Hamilton. 时间序列分析 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.