目录

引言

季节性指数平滑优化 (SESO) 模型是一种用于预测具有季节性模式的时间序列未来值的预测技术。它是指数平滑方法的一种变体,该方法结合使用过去值和预测值来生成预测。

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

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

季节性指数平滑模型

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

y^t+1,s=αyt+(1α)y^t1,s\hat{y}{t+1,s} = \alpha y_t + (1-\alpha) \hat{y}{t-1,s}

其中: - y^t+1,s\hat{y}{t+1,s} 是季节 ss 下一期的预测值。 - α\alpha 是通过最小化平方误差来优化的平滑参数。 - yty_t 是季节 ss 在时期 tt 的当前观测值。 - y^t1,s\hat{y}{t-1,s} 是季节 ss 前一期的预测值。

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

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

模型选择

SESO 模型中的模型选择是指为模型选择平滑参数和季节性分量的最优值的过程。这些参数的最优值是能使给定数据集获得最佳预测性能的值。

ETS 统计框架的一个巨大优势在于可以使用信息准则进行模型选择。AIC,AICcAIC, AIC_cBICBIC 也可以用于确定哪种 ETS 模型最适合给定的时间序列。

对于 ETS 模型,Akaike 信息准则 (AIC) 定义为 AIC=2log(L)+2k,\text{AIC} = -2\log(L) + 2k,

其中 LL 是模型的似然性,kk 是估计的总参数和初始状态数(包括残差方差)。

对小样本偏差进行校正的 AIC (AICcAIC_c) 定义为 AICc=AIC+2k(k+1)Tk1,\text{AIC}_{\text{c}} = \text{AIC} + \frac{2k(k+1)}{T-k-1},

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

这些准则平衡了模型的拟合优度和复杂度,提供了一种方法来选择最大化数据似然性同时最小化参数数量的模型。

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

总的来说,SESO 模型的模型选择过程结合了统计技术、信息准则和专家判断,以确定平滑参数和季节性分量的最优值,从而使给定数据集获得最佳预测性能。

加载库和数据

提示

需要 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)

增广迪基-福勒检验

增广迪基-福勒 (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 显示了观测值与其在不同时间间隔的滞后值之间的相关性。有助于识别时间序列中的时间依赖模式,例如趋势或季节性的存在。

  1. 指示序列的“记忆”:ACF 使我们能够确定过去的观测值对未来观测值的影响程度。如果 ACF 在多个滞后处显示出显著的自相关,则表明该序列具有长期记忆,并且过去的观测值与预测未来观测值相关。

  2. 有助于识别 MA(移动平均)模型:ACF 的形状可以揭示时间序列中移动平均分量的存在。ACF 显示显著相关性的滞后可能指示 MA 模型的阶数。

偏自相关 (PACF):1. 识别直接依赖性:与 ACF 不同,PACF 消除了中间滞后的间接影响,并测量了观测值与其滞后值之间的直接相关性。它有助于识别观测值与其滞后值之间的直接依赖性,而没有中间滞后的影响。

  1. 有助于识别 AR(自回归)模型:PACF 的形状可以揭示时间序列中自回归分量的存在。PACF 显示显著相关性的滞后可能指示 AR 模型的阶数。

  2. 与 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");

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

时间序列分解

如何以及为何分解时间序列?

在时间序列分析中预测新值时,了解过去数据非常重要。更正式地说,了解值随时间变化的模式非常重要。有很多原因可能导致我们的预测值方向错误。基本上,时间序列由四个分量组成。这些分量的变化导致时间序列模式的变化。这些分量是

  • 水平 (Level):这是随时间平均的主要值。
  • 趋势 (Trend):趋势是导致时间序列中出现增长或下降模式的值。
  • 季节性 (Seasonality):这是一种在时间序列中短期发生的周期性事件,导致时间序列中出现短期增长或下降模式。
  • 残差/噪声 (Residual/Noise):这是时间序列中的随机变化。

将这些分量随时间结合起来就形成了时间序列。大多数时间序列包含水平和噪声/残差,而趋势或季节性是可选分量。

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

时间序列中分量的组合可以是两种类型: * 加法型 * 乘法型

加法时间序列

如果时间序列的分量相加以构成时间序列。则该时间序列称为加法时间序列。通过可视化,如果时间序列的增长或下降模式在整个序列中相似,则可以说该时间序列是加法型的。任何加法时间序列的数学函数可以表示为:y(t)=level+Trend+seasonality+noisey(t) = level + Trend + seasonality + noise

乘法时间序列

如果时间序列的分量相乘构成时间序列,则该时间序列称为乘法时间序列。从可视化角度来看,如果时间序列随时间呈指数增长或下降,则可以认为是乘法时间序列。乘法时间序列的数学函数可以表示为。

y(t)=LevelTrendseasonalityNoisey(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();

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

让我们将数据分成几个集合

  1. 用于训练我们的 Seasonal 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 实现 SeasonalExponentialSmoothingOptimized 模型

加载库

from statsforecast import StatsForecast
from statsforecast.models import SeasonalExponentialSmoothingOptimized

构建模型

导入并实例化模型。设置参数有时很棘手。Rob Hyndmann 大师撰写的关于季节周期的文章对于设置 season_length 参数可能很有帮助。

season_length = 24 # Hourly data 
horizon = len(test) # number of predictions

models = [SeasonalExponentialSmoothingOptimized(season_length=season_length)]

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

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

  • freq:一个字符串,表示数据的频率。(请参阅panda 的可用频率。)

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

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

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

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

拟合模型

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

让我们看看 Seasonal Exponential Smoothing Optimized Model 的结果。我们可以通过以下指令查看

result=sf.fitted_[0,0].model_
result
{'mean': array([161532.05 , 161051.69 , 135531.64 , 105600.39 ,  96717.39 ,
         82608.34 ,  80224.33 ,  78075.98 ,  85233.23 , 100179.336,
        122245.62 , 118087.57 , 109614.81 , 104729.91 , 104895.02 ,
        115862.96 , 130370.98 , 144231.89 , 149036.73 , 149072.73 ,
        148110.77 , 148760.73 , 149767.53 , 150561.8  ], 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,        nan,
               nan,        nan,        nan,        nan,        nan,
               nan,        nan,        nan,        nan,        nan,
               nan,        nan,        nan,        nan,        nan,
               nan,        nan, 163840.   , 166235.   , 139520.   ,
        105895.   ,  96780.   ,  82520.   ,  80125.   ,  75335.   ,
         85105.   , 102080.   , 125135.   , 118030.   , 109225.   ,
        102475.   , 102240.   , 115840.   , 130540.   , 144325.   ,
        148970.   , 149150.   , 148040.   , 148810.   , 149830.   ,
        150570.   , 162030.27 , 163222.1  , 137347.33 , 103835.8  ,
         96733.95 ,  82522.45 ,  80086.9  ,  75132.05 ,  85074.36 ,
        100452.66 , 121044.03 , 118001.6  , 109242.15 , 102349.03 ,
        102321.49 , 115768.25 , 130501.   , 144286.1  , 149005.   ,
        149121.25 , 148039.8  , 148799.25 , 149789.2  , 150557.16 ,
        161740.55 , 162812.36 , 136965.22 , 112853.91 ,  96768.61 ,
         82573.375,  80164.38 ,  88707.87 ,  85164.8  , 100944.15 ,
        117929.875, 118111.086, 109563.58 , 103815.7  , 104036.375,
        115942.47 , 130508.39 , 144268.03 , 149088.7  , 149155.03 ,
        148096.75 , 148823.2  , 149797.77 , 150525.92 , 160582.38 ,
        159756.83 , 134514.39 , 117874.29 ,  96767.92 ,  82683.74 ,
         80253.336,  89338.625,  85232.055, 100619.03 , 114659.62 ,
        118224.67 , 109881.99 , 105514.21 , 106070.33 , 116194.74 ,
        130678.805, 144436.45 , 149261.16 , 149331.28 , 148247.19 ,
        148908.03 , 149890.33 , 150620.88 , 161557.95 , 161701.48 ,
        136228.19 , 113004.195,  96773.695,  82673.66 ,  80245.91 ,
         78459.88 ,  85267.016, 100517.48 , 120224.3  , 118155.68 ,
        109777.57 , 105240.35 , 105717.734, 116058.445, 130571.32 ,
        144349.64 , 149169.75 , 149255.28 , 148186.9  , 148872.55 ,
        149844.78 , 150618.77 , 161553.17 , 160112.8  , 135912.81 ,
        107124.39 ,  96756.41 ,  82642.07 ,  80226.8  ,  74707.9  ,
         85226.28 , 100202.98 , 119687.805, 118105.27 , 109668.59 ,
        104848.68 , 105212.516, 115939.71 , 130463.1  , 144266.05 ,
        149142.61 , 149154.38 , 148177.3  , 148833.17 , 149860.03 ,
        150673.38 ], dtype=float32)}

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

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

fitted=pd.DataFrame(result.get("fitted"), columns=["fitted"])
fitted["ds"]=df["ds"]
fitted
拟合值ds
0NaN2017-09-13 00:00:00
1NaN2017-09-13 01:00:00
2NaN2017-09-13 02:00:00
183148833.1718752017-09-20 15:00:00
184149860.0312502017-09-20 16:00:00
185150673.3750002017-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 步。在此示例中,提前 12 个月。

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

# Prediction
Y_hat = sf.forecast(df=train, h=horizon, fitted=True)
Y_hat
unique_iddsSeasESOpt
012017-09-20 18:00:00161532.046875
112017-09-20 19:00:00161051.687500
212017-09-20 20:00:00135531.640625
2712017-09-21 21:00:00105600.390625
2812017-09-21 22:00:0096717.390625
2912017-09-21 23:00:0082608.343750
values=sf.forecast_fitted_values()
values.head()
unique_iddsySeasESOpt
012017-09-13 00:00:0080115.0NaN
112017-09-13 01:00:0079885.0NaN
212017-09-13 02:00:0089325.0NaN
312017-09-13 03:00:00101930.0NaN
412017-09-13 04:00:00121630.0NaN
sf.plot(train, Y_hat)

使用置信区间的预测方法

使用 predict 方法生成预测。

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

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

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

此步骤应少于 1 秒。

forecast_df = sf.predict(h=horizon) 
forecast_df
unique_iddsSeasESOpt
012017-09-20 18:00:00161532.046875
112017-09-20 19:00:00161051.687500
212017-09-20 20:00:00135531.640625
2712017-09-21 21:00:00105600.390625
2812017-09-21 22:00:0096717.390625
2912017-09-21 23:00:0082608.343750

交叉验证

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

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

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

执行时间序列交叉验证

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

在此示例中,我们希望评估每个模型在过去 5 个窗口 (n_windows=) 的性能,每次预测间隔 12 步 (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:序列标识符。
  • ds:日期戳或时间索引
  • cutoff: n_windows 中每个窗口的最后一个日期戳或时间索引。
  • y:真实值
  • model:包含模型名称和拟合值的列。
crossvalidation_df
unique_idds截止点ySeasESOpt
012017-09-18 06:00:002017-09-18 05:00:0099440.0141401.750000
112017-09-18 07:00:002017-09-18 05:00:0097655.0152474.250000
212017-09-18 08:00:002017-09-18 05:00:0097655.0152482.796875
8712017-09-21 21:00:002017-09-20 17:00:00103080.0105600.390625
8812017-09-21 22:00:002017-09-20 17:00:0095155.096717.390625
8912017-09-21 23:00:002017-09-20 17:00:0080285.082608.343750

模型评估

现在我们将使用预测结果评估我们的模型,我们将使用不同类型的指标 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=season_length), ufl.rmse, ufl.smape],
    train_df=train,
)
unique_id指标SeasESOpt
01mae6694.042188
11mape0.060392
21mase0.827062
31rmse8118.297509
41smape0.028961

参考文献

  1. Changquan Huang • Alla Petukhina. Springer 系列 (2022). Applied Time Series Analysis and Forecasting with Python.
  2. Ivan Svetunkov. 使用增广动态自适应模型 (ADAM) 进行预测和分析
  3. James D. Hamilton. Time Series Analysis Princeton University Press, 普林斯顿, 新泽西, 第 1 版, 1994.
  4. Nixtla 参数.
  5. Pandas 可用频率.
  6. Rob J. Hyndman 和 George Athanasopoulos (2018). “Forecasting Principles and Practice (第 3 版)”.
  7. 季节周期 - Rob J Hyndman.