动态优化Theta模型


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

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的动态优化Theta模型的逐步指南。

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

目录

介绍

动态优化Theta模型（DOTM）是一种用于预测时间序列未来值的预测技术。它是Theta方法的一个变体，该方法结合了指数平滑和线性趋势来预测未来值。

DOTM通过引入动态优化过程扩展了Theta方法，该过程基于历史数据选择指数平滑组件的最佳平滑参数以及线性趋势组件的最佳权重。这个优化过程是通过遗传算法迭代进行的，该算法搜索能够最小化预测误差的参数组合。

DOTM旨在处理时间序列数据，这些数据随时间表现出非线性和非平稳行为。它特别适用于预测具有复杂模式的时间序列，如季节性、趋势和周期波动。

DOTM相对于其他预测方法有几个优点。首先，它是一种简单易行的方法，不需要广泛的统计知识。其次，它是一种高度适应的方法，可以根据多种时间序列数据进行定制。第三，它是一种健壮的方法，可以处理时间序列中的缺失数据、异常值和其他异常情况。

总体而言，动态优化Theta模型是一种强大的预测技术，可以用于生成广泛时间序列数据的准确和可靠的预测。

动态优化Theta模型（DOTM）

到目前为止，我们将\(A_n\)和\(B_n\)设为所有\(t\)的固定系数。现在，我们将这些系数视为动态函数；即，在将状态\(t\)更新为\(t+1\)时，我们将仅考虑先前的信息\(Y_1, \cdots, Y_t\)来计算\(A_t\)和\(B_t\)。因此，我们在优化Theta模型的笔记本的方程(3)和(4)中用\(A_t\)和\(B_t\)替代\(A_n\)和\(B_n\)。然后，在将新的方程(4)应用于新的方程(3)并在时间\(t\)以\(h=1\)重写结果后，我们得到

\[ \begin{equation} \hat Y_{t+1|t}=\ell_{t} + \left(1 - \frac{1}{\theta} \right) \left( (1-\alpha)^t A_t + \left[ \frac{1 - ( 1 - \alpha)^{t+1}}{\alpha} \right] B_t \tag 1 \right) \end{equation} \]

然后，假设添加一阶前向误差并重写方程(3)（见AutoTheta模型）、(1)，我们得到

\[ \begin{equation} Y_t=\mu_t +\varepsilon_t \tag 2 \end{equation} \]

\[ \begin{equation} \mu_t=\ell_{t-1}+ \left(1-\frac{1}{\theta}\right) \left( \left(1-\alpha\right)^{t-1} A_{t-1} + \left(\frac{1-(1-\alpha)^{t}}{\alpha}\right) B_{t-1} \tag 3 \right) \end{equation} \]

\[ \begin{equation} \ell_{t}=\alpha Y_t+ (1-\alpha) \ell_{t-1} \tag 4 \end{equation} \]

\[ \begin{equation} A_t=\bar Y_t - \frac{t+1}{2} B_t \tag 5 \end{equation} \]

\[ \begin{equation} B_t=\frac{1}{t+1} \left((t-2) B_{t-1} +\frac{6}{t} (Y_t - \bar Y_{t-1}) \right) \tag 6 \end{equation} \]

\[ \begin{equation} \bar Y_t=\frac{1}{t} \left((t-1) \bar Y_{t-1} + Y_t \right) \tag 7 \end{equation} \]

对于\(t=1, \cdots ,n\)。方程(2), (3), (4), (5), (6), (7)配置了一个状态空间模型，参数为\(\ell_{0} \in \mathbb{R}, \alpha \in (0,1)\)，和\(\theta \in [1,\infty )\)。状态的初始化假设为\(A_0 =B_0=B_1=\bar Y_0 =0\)。此后，我们将该模型称为动态优化Theta模型（DOTM）。

DOTM的一个重要属性是，当\(\theta=1\)时，这意味着\(Z_t(1)=Y_t\)，由方程(3)给出的预测向量（见OTM）将等于 \[\hat Y_{t+h|t} = \tilde Z_{t+h|t}\]

因此，当\(\theta=1\)时，DOTM就是SES方法。当\(\theta>1\)时，DOTM（作为SES-d）作为SES的扩展，添加了一个长期组件。

DOTM在原点产生的样本外一阶前向预测为

\[ \begin{equation} \hat Y_{n+1|n}=E[Y_{n+1|Y_1, \cdots, Y_n} ]=\ell_{n} + \left(1-\frac{1}{\theta}\right) \left( (1-\alpha)^n A_n + \left(\frac{1-(1-\alpha)^{n+1}}{\alpha}\right) B_n \right) \tag 8 \end{equation} \]

对于\(h \geq 2\)的预测，\(\hat Y_{n+2|n}, \cdots , \hat Y_{n+h|n}\)是通过使用方程(3), (4), (5), (6), (7), (8)递归计算的，将未观察的值\(Y_{n+1}, \cdots , Y_{n+h-1}\)替换为它们的期望值\(\hat Y_{n+1|n}, \cdots , \hat Y_{n+h-1|n}\)。条件方差\(Var[Y_{n+h}|Y_{1}, \cdots, Y_n ]\)很难以分析的方式写出。然而，\(Y_{n+h}\)的方差和预测区间可以使用自助法技术进行估计，其中从估计模型中模拟出一个（通常很大）可能的\(Y_{n+h}\)值的样本。

请注意，与STheta、STM和OTM不同，DSTM和DOTM产生的预测不一定是线性的。这也是DSTM/DOTM与SES-d之间的一个根本区别：在SES-d中，长期趋势是常量，而对于DSTM/DOTM，这一点在样本内拟合或样本外预测中都不是这样。

加载库和数据

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/milk_production.csv", usecols=[1,2])
df.head()
month production
0 1962-01-01 589
1 1962-02-01 561
2 1962-03-01 640
3 1962-04-01 656
4 1962-05-01 727

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 1962-01-01 589 1
1 1962-02-01 561 1
2 1962-03-01 640 1
3 1962-04-01 656 1
4 1962-05-01 727 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, )

自相关图

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

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

让我们将数据划分为不同的集合

  1. 用于训练我们的 动态优化Theta模型(DOTM) 的数据。
  2. 用于测试我们模型的数据。

对于测试数据，我们将使用过去12个月的数据来测试和评估我们模型的表现。

train = df[df.ds<='1974-12-01'] 
test = df[df.ds>'1974-12-01'] 
train.shape, test.shape
((156, 3), (12, 3))

现在让我们绘制训练数据和测试数据。

sns.lineplot(train,x="ds", y="y", label="Train", linestyle="--")
sns.lineplot(test, x="ds", y="y", label="Test")
plt.title("Monthly Milk Production");
plt.show()

DynamicOptimizedTheta的实现与StatsForecast

要了解DynamicOptimizedTheta模型的函数参数，下面列出了相关信息。有关更多信息，请访问文档。

season_length : int
    每单位时间的观测次数。例如：24小时数据。
decomposition_type : str
    季节分解类型，'multiplicative'（默认）或'additive'。
alias : str
    模型的自定义名称。
prediction_intervals : Optional[ConformalIntervals]
    计算保留预测区间的信息。
    默认情况下，模型将计算原生预测
    区间。

加载库

from statsforecast import StatsForecast
from statsforecast.models import DynamicOptimizedTheta

实例化模型

导入并实例化模型。设置参数有时比较棘手。Rob Hyndmann的大师文章 季节周期 对于 season_length 可能会很有用。

season_length = 12 # 月度数据 
horizon = len(test) # 预测数量

# 我们将要使用的模型称为
models = [DynamicOptimizedTheta(season_length=season_length, 
                decomposition_type="additive")] # 乘法  加法

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

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

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

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

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

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

sf = StatsForecast(df=train,
                   models=models,
                   freq='MS', 
                   n_jobs=-1)

拟合模型

sf.fit()
StatsForecast(models=[DynamicOptimizedTheta])

让我们看看我们动态优化Theta模型的结果。我们可以通过以下指令来观察它：

result=sf.fitted_[0,0].model_
print(result.keys())
print(result['fit'])
dict_keys(['mse', 'amse', 'fit', 'residuals', 'm', 'states', 'par', 'n', 'modeltype', 'mean_y', 'decompose', 'decomposition_type', 'seas_forecast', 'fitted'])
results(x=array([250.83207246,   0.75624902,   4.67964777]), fn=10.697567725248804, nit=55, simplex=array([[237.42075735,   0.75306547,   4.46023813],
       [250.83207246,   0.75624902,   4.67964777],
       [257.164453  ,   0.75229688,   4.42377059],
       [256.90854919,   0.75757957,   4.43171897]]))

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

如我们所见，上述结果以字典的形式输出，为了从字典中提取每个元素，我们将使用 .get() 函数来提取元素，随后将其保存到 pd.DataFrame() 中。

residual=pd.DataFrame(result.get("residuals"), columns=["residual Model"])
residual
residual Model
0 -18.247131
1 -88.625732
2 2.864929
... ...
153 -59.747070
154 -91.901550
155 -43.503296

156 rows × 1 columns

import scipy.stats as stats

fig, axs = plt.subplots(nrows=2, ncols=2)

residual.plot(ax=axs[0,0])
axs[0,0].set_title("Residuals");

sns.distplot(residual, ax=axs[0,1]);
axs[0,1].set_title("Density plot - Residual");

stats.probplot(residual["residual Model"], dist="norm", plot=axs[1,0])
axs[1,0].set_title('Plot Q-Q')

plot_acf(residual,  lags=35, ax=axs[1,1],color="fuchsia")
axs[1,1].set_title("Autocorrelation");

plt.show();

预测方法

如果您想在生产环境中提高速度，特别是在处理多个系列或模型时，我们建议使用 StatsForecast.forecast 方法，而不是 .fit 和 .predict。

主要区别在于 .forecast 不会存储拟合值，并且在分布式环境中具有很高的可扩展性。

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

  • h (int): 表示预测向未来 h 步的值。在这种情况下，是指 12 个月之后。

  • level (list of floats): 这个可选参数用于概率预测。设置您的预测区间的水平（或置信分位数）。例如，level=[90] 意味着模型期望真实值在该区间内的概率为 90%。

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

# 预测
Y_hat = sf.forecast(horizon, fitted=True)

Y_hat
ds DynamicOptimizedTheta
unique_id
1 1975-01-01 839.259705
1 1975-02-01 801.399170
1 1975-03-01 895.189148
... ... ...
1 1975-10-01 821.271179
1 1975-11-01 792.530396
1 1975-12-01 829.854492

12 rows × 2 columns

values=sf.forecast_fitted_values()
values.head()
ds y DynamicOptimizedTheta
unique_id
1 1962-01-01 589.0 607.247131
1 1962-02-01 561.0 649.625732
1 1962-03-01 640.0 637.135071
1 1962-04-01 656.0 609.225830
1 1962-05-01 727.0 604.995300 
StatsForecast.plot(values)

添加95%置信区间与预测方法

sf.forecast(h=horizon, level=[95])
ds DynamicOptimizedTheta DynamicOptimizedTheta-lo-95 DynamicOptimizedTheta-hi-95
unique_id
1 1975-01-01 839.259705 741.963562 955.137634
1 1975-02-01 801.399170 641.886230 946.029053
1 1975-03-01 895.189148 707.210693 1066.337036
... ... ... ... ...
1 1975-10-01 821.271179 546.113647 1088.162476
1 1975-11-01 792.530396 494.657928 1037.432007
1 1975-12-01 829.854492 519.697144 1108.181885

12 rows × 4 columns

Y_hat=Y_hat.reset_index()
Y_hat
unique_id ds DynamicOptimizedTheta
0 1 1975-01-01 839.259705
1 1 1975-02-01 801.399170
2 1 1975-03-01 895.189148
... ... ... ...
9 1 1975-10-01 821.271179
10 1 1975-11-01 792.530396
11 1 1975-12-01 829.854492

12 rows × 3 columns

# 将预测结果与真实值合并
test['unique_id'] = test['unique_id'].astype(int)
Y_hat1 = test.merge(Y_hat, how='left', on=['unique_id', 'ds'])
Y_hat1
ds y unique_id DynamicOptimizedTheta
0 1975-01-01 834 1 839.259705
1 1975-02-01 782 1 801.399170
2 1975-03-01 892 1 895.189148
... ... ... ... ...
9 1975-10-01 827 1 821.271179
10 1975-11-01 797 1 792.530396
11 1975-12-01 843 1 829.854492

12 rows × 4 columns

fig, ax = plt.subplots(1, 1)
plot_df = pd.concat([train, Y_hat1]).set_index('ds')
plot_df[['y', "DynamicOptimizedTheta"]].plot(ax=ax, linewidth=2)
ax.set_title(' Forecast', fontsize=22)
ax.set_ylabel('Monthly Milk Production', fontsize=20)
ax.set_xlabel('Timestamp [t]', fontsize=20)
ax.legend(prop={'size': 15})
ax.grid(True)

带置信区间的预测方法

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

预测方法接受两个参数：预测下一个 h（代表时间范围）和 level。

  • h (int): 表示预测未来 h 步。在这种情况下，表示 12 个月。

  • level (list of floats): 此可选参数用于概率预测。设置您预测区间的水平（或置信百分位数）。例如，level=[95] 表示模型期望实际值在该区间内的概率为 95%。

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

此步骤应在 1 秒内完成。

sf.predict(h=horizon)
ds DynamicOptimizedTheta
unique_id
1 1975-01-01 839.259705
1 1975-02-01 801.399170
1 1975-03-01 895.189148
... ... ...
1 1975-10-01 821.271179
1 1975-11-01 792.530396
1 1975-12-01 829.854492

12 rows × 2 columns

forecast_df = sf.predict(h=horizon, level=[80,95]) 
forecast_df
ds DynamicOptimizedTheta DynamicOptimizedTheta-lo-80 DynamicOptimizedTheta-hi-80 DynamicOptimizedTheta-lo-95 DynamicOptimizedTheta-hi-95
unique_id
1 1975-01-01 839.259705 766.150513 928.015320 741.963562 955.137634
1 1975-02-01 801.399170 702.992554 899.872864 641.886230 946.029053
1 1975-03-01 895.189148 760.141418 1008.321960 707.210693 1066.337036
... ... ... ... ... ... ...
1 1975-10-01 821.271179 617.415344 996.678284 546.113647 1088.162476
1 1975-11-01 792.530396 568.329102 975.049255 494.657928 1037.432007
1 1975-12-01 829.854492 598.124878 1035.452515 519.697144 1108.181885

12 rows × 6 columns

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

pd.concat([df, forecast_df]).set_index('ds')
y unique_id DynamicOptimizedTheta DynamicOptimizedTheta-lo-80 DynamicOptimizedTheta-hi-80 DynamicOptimizedTheta-lo-95 DynamicOptimizedTheta-hi-95
ds
1962-01-01 589.0 1 NaN NaN NaN NaN NaN
1962-02-01 561.0 1 NaN NaN NaN NaN NaN
1962-03-01 640.0 1 NaN NaN NaN NaN NaN
... ... ... ... ... ... ... ...
1975-10-01 NaN NaN 821.271179 617.415344 996.678284 546.113647 1088.162476
1975-11-01 NaN NaN 792.530396 568.329102 975.049255 494.657928 1037.432007
1975-12-01 NaN NaN 829.854492 598.124878 1035.452515 519.697144 1108.181885

180 rows × 7 columns

现在让我们可视化我们的预测结果和时间序列的历史数据，同时绘制出我们在以95%的置信度进行预测时获得的置信区间。

def plot_forecasts(y_hist, y_true, y_pred, models):
    _, ax = plt.subplots(1, 1, figsize = (20, 7))
    y_true = y_true.merge(y_pred, how='left', on=['unique_id', 'ds'])
    df_plot = pd.concat([y_hist, y_true]).set_index('ds').tail(12*10)
    df_plot[['y'] + models].plot(ax=ax, linewidth=3 , )
    colors = ['green', "lime"]
    ax.fill_between(df_plot.index, 
                df_plot['DynamicOptimizedTheta-lo-80'], 
                df_plot['DynamicOptimizedTheta-hi-80'],
                alpha=.20,
                color='orange',
                label='DynamicOptimizedTheta_level_80')
    ax.fill_between(df_plot.index, 
                df_plot['DynamicOptimizedTheta-lo-95'], 
                df_plot['DynamicOptimizedTheta-hi-95'],
                alpha=.3,
                color='lime',
                label='DynamicOptimizedTheta_level_95')
    ax.set_title('', fontsize=22)
    ax.set_ylabel("Montly Mil Production", fontsize=20)
    ax.set_xlabel('Month-Days', fontsize=20)
    ax.legend(prop={'size': 20})
    ax.grid(True)
    plt.show()
plot_forecasts(train, test, forecast_df, models=['DynamicOptimizedTheta'])

交叉验证

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

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

以下图表描述了这种交叉验证策略：

执行时间序列交叉验证

时间序列模型的交叉验证被认为是最佳实践，但大多数实现速度非常慢。statsforecast 库将交叉验证实现为分布式操作，使得这一过程不再耗时。如果您有大数据集，还可以使用 Ray、Dask 或 Spark 在分布式集群中执行交叉验证。

在这种情况下，我们想要评估每个模型在过去 5 个月的表现 (n_windows=5)，每隔一个月进行一次预测 (step_size=12)。根据您的计算机，这一步大约需要 1 分钟。

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

  • df: 训练数据框

  • h (int): 表示预测的未来 h 步。这种情况下是 12 个月。

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

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

crossvalidation_df = sf.cross_validation(df=train,
                                         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 DynamicOptimizedTheta
unique_id
1 1972-01-01 1971-12-01 826.0 828.692017
1 1972-02-01 1971-12-01 799.0 792.444092
1 1972-03-01 1971-12-01 890.0 883.122620
... ... ... ... ...
1 1974-10-01 1973-12-01 812.0 810.342834
1 1974-11-01 1973-12-01 773.0 781.845703
1 1974-12-01 1973-12-01 813.0 818.855103

36 rows × 4 columns

模型评估

现在我们可以使用合适的准确性指标来计算预测的准确性。在这里，我们将使用均方根误差（RMSE）。为此，我们首先需要install datasetsforecast，这是一个由Nixtla开发的Python库，其中包含计算RMSE的函数。

%%capture
!pip install datasetsforecast
from datasetsforecast.losses import rmse

计算RMSE的函数接受两个参数：

  1. 实际值。
  2. 预测值，在这种情况下为动态优化Theta模型。
rmse = rmse(crossvalidation_df['y'], crossvalidation_df["DynamicOptimizedTheta"])
print("RMSE using cross-validation: ", rmse)
RMSE using cross-validation:  14.320624

正如您所注意到的，我们已经使用交叉验证的结果来评估我们的模型。

现在我们将使用预测结果来评估我们的模型，我们将使用不同类型的指标 MAE, MAPE, MASE, RMSE, SMAPE 来评估 准确性。

from datasetsforecast.losses import (mae, mape, mase, rmse, smape)
def evaluate_performace(y_hist, y_true, y_pred, model):
    y_true = y_true.merge(y_pred, how='left', on=['unique_id', 'ds'])
    evaluation = {}
    evaluation[model] = {}
    for metric in [mase, mae, mape, rmse, smape]:
        metric_name = metric.__name__
        if metric_name == 'mase':
            evaluation[model][metric_name] = metric(y_true['y'].values, 
                                                y_true[model].values, 
                                                y_hist['y'].values, seasonality=12)
        else:
            evaluation[model][metric_name] = metric(y_true['y'].values, y_true[model].values)
    return pd.DataFrame(evaluation).T
evaluate_performace(train, test, Y_hat, model="DynamicOptimizedTheta")
mae mape mase rmse smape
DynamicOptimizedTheta 6.861959 0.804529 0.308595 8.64746 0.801941

致谢

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

参考文献

  1. Kostas I. Nikolopoulos, Dimitrios D. Thomakos. 使用Theta方法进行预测——理论与应用。2019年，约翰·威利与儿子有限公司。
  2. Jose A. Fiorucci, Tiago R. Pellegrini, Francisco Louzada, Fotios Petropoulos, Anne B. Koehler (2016). “优化Theta方法的模型及其与状态空间模型的关系”。 国际预测学期刊。
  3. Nixtla参数。
  4. Pandas可用频率。
  5. Rob J. Hyndman 和 George Athanasopoulos (2018). “预测的原则与实践，时间序列交叉验证”。。
  6. 季节性周期 - Rob J Hyndman。

Give us a ⭐ on Github