标准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 方法 (Assimakopoulos & Nikolopoulos, 2000，以下简称 A&N) 应用于非季节性或季节调整后的时间序列，其中季节调整通常通过乘法经典分解来进行。该方法通过所谓的 theta 系数将原始时间序列分解为两条新线，记作 \({\theta}_1\) 和 \({\theta}_2\)，其中 \({\theta}_1, {\theta}_2 \in \mathbb{R}\)，这些系数应用于数据的第二差分。当 \({\theta}<1\) 时，第二差分会减少，从而更好地逼近序列的长期行为（Assimakopoulos, 1995）。如果 \({\theta}\) 等于零，新线是一条直线。当 \({\theta}>1\) 时，局部曲率增加，放大了时间序列的短期波动（A&N）。产生的新线称为 theta 线，记作 \(\text{Z}(\theta_1)\) 和 \(\text{Z}(\theta_2)\)。这些线具有与原始数据相同的均值和斜率，但局部曲率根据 \(\theta\) 系数的值被过滤或增强。

换句话说，分解过程的优点在于利用了数据中通常无法通过原始时间序列的外推完全捕捉和建模的信息。Theta 线可以视为新的时间序列，并使用适当的预测方法单独外推。一旦每条 theta 线的外推完成，就会通过组合方案进行重组，以计算原始时间序列的点预测。组合在预测文献中长期被认为是一种有用的实践（例如，Clemen, 1989, Makridakis and Winkler, 1983, Petropoulos et al., 2014），因此将其应用于 Theta 方法预计会带来更准确和更稳健的预测。

Theta方法在选择theta线的数量、theta系数和外推方法方面非常灵活，并结合这些方法来获得稳健的预测。然而，A&N提出了一种简化版本，仅使用两个theta线，并使用预设的\(\theta\)系数，通过线性回归（LR）模型对theta线进行时间外推，其中\({\theta}_1 =0\)，而对于\({\theta}_2 =2\)的theta线，则使用简单指数平滑（SES）。最终的预测是通过将两个theta线的预测以相等权重结合而产生的。

Theta方法的性能也得到了其他实证研究的确认（例如Nikolopoulos等，2012年，Petropoulos和Nikolopoulos，2013年）。此外，Hyndman和Billah（2003年，以下简称H&B）展示了带漂移的简单指数平滑模型（SES-d）是Theta方法简化版本的统计模型。最近，Thomakos和Nikolopoulos（2014年）提供了更多的理论见解，而Thomakos和Nikolopoulos（2015年）为该方法在多元时间序列中的应用推导了新的理论公式，并研究了二元Theta方法比单元Theta方法更优秀的预测条件。尽管取得了这些进展，我们认为Theta方法应当受到预测社区更多的关注，因为它简单且具有优越的预测性能。

Theta方法的一个关键方面是，它在定义上是动态的。可以选择不同的theta线，并使用相等或不等的权重结合产生的预测。然而，AN通过将theta系数固定为预定义值，限制了这一重要属性。

标准Theta模型

Assimakopoulos和Nikolopoulo为标准Theta模型提出了Theta线作为方程的解

\[ \begin{equation} D^2 \zeta_t(\theta) = \theta D^2 Y_t, t = 1,\cdots,T \tag 1 \end{equation} \]

其中 \(Y_1, \cdots , Y_T\) 代表原始时间序列数据，而 \(DX_t = (X_t − X_{t−1})\)。初始值 \(\zeta_1\) 和 \(\zeta_2\) 是通过最小化 \(\sum_{i=1}^{T} [Y_t - \zeta_t (\theta) ]^2\) 来获得的。然而，方程(1)的解析解为

\[ \begin{equation} \zeta_t(\theta)=\theta Y_t +(1−\theta)(A_T +B_T t),\ t=1, \cdots, T, \tag 2 \end{equation} \]

其中 \(A_T\) 和 \(B_T\) 是对 \(Y_1, \cdots,Y_T\) 与 \(1, \cdots , T\) 进行简单线性回归的最小平方系数，这些系数仅依赖于原始数据，具体如下

\[ \begin{equation} A_T=\frac{1}{T} \sum_{i=1}^{T} Y_t - \frac{T+1}{2} B_T \tag 3 \end{equation} \]

\[ \begin{equation} B_T=\frac{6}{T^2 - 1} (\frac{2}{T} \sum_{t=1}^{T} tY_t - \frac{T+1}{T} \sum_{t=1}^{T} Y_t \tag 4 \end{equation} \]

从这个角度来看，Theta线可以理解为直接应用于数据的线性回归模型的函数。事实上，Theta方法在h步预测时，是对\(\zeta(0)\)和\(\zeta(2)\)线性外推的一个临时组合（50% - 50%）。

  • 当\(\theta < 1\) 应用于数据的二阶差分时，分解过程由Theta系数定义，该系数减少了二阶差分并改善了序列行为的近似。

  • 如果\(\theta = 0\)，则分解后的线变为一条常数直线。 （见图）

  • 如果\(\theta > 1\)，则被分析系列的短期波动显示出更多的局部曲率（见图）

图

我们将上述设定称为标准Theta方法。构建Theta方法的步骤如下：

  1. 去季节化： 首先，对时间序列数据进行统计显著性季节行为的测试。如果时间序列是季节性的，则满足以下条件：

\[|\rho_m| > q_{1- \frac{\alpha}{2} } \sqrt{\frac{1+2 \sum_{i=1}^{m-1} \rho_{i}^{2} }{T} }\]

其中ρk表示滞后\(k\)自相关函数，\(m\)是一个季节周期内的周期数（例如，月度数据为12），\(T\)是样本容量，\(q\)是标准正态分布的分位数函数，\((1 − a)\%\)是置信水平。Assimakopoulos和Nikolopoulo [标准Theta模型]选择了90%的置信水平。如果时间序列被识别为季节性，则通过经典分解法去季节化，假设季节成分具有乘法关系。

  1. 分解： 第二步是将经季节调整的时间序列分解为两个Theta线：线性回归线\(\zeta(0)\)和Theta线\(\zeta(2)\)。

  2. 外推： 使用简单指数平滑（SES）对\(\zeta(2)\)进行外推，而\(\zeta(0)\)则作为正常的线性回归线进行外推。

  3. 组合： 最终预测是两个\(\theta\)线预测的组合，使用相等的权重。

  4. 季节化： 如果在第一步中存在季节性，则最终预测值将乘以相应的季节指数。

加载库和数据

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方法绘制一些系列。此方法从数据集中打印一个随机系列，对于基本的探索性数据分析（EDA）非常有用。

from statsforecast import StatsForecast

StatsForecast.plot(df, engine="matplotlib")

自相关图

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 模型的数据 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()

StandardTheta的实现与StatsForecast

为了更好地了解StandardTheta Model函数的参数，它们在下面列出。有关更多信息，请访问文档。

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

加载库

from statsforecast import StatsForecast
from statsforecast.models import Theta

实例化模型

导入并实例化模型。设置参数有时比较困难。大师 Rob Hyndmann 的这篇关于 季节周期 的文章对 season_length 会很有帮助。

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

models = [Theta(season_length=season_length, 
                decomposition_type="additive")] # 乘法  加法

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

模型：一个模型列表。从模型中选择你想要的模型并导入。

  • 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=[Theta])

让我们看看我们的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([225.82002697,   0.76015625]), fn=10.638733596938769, nit=19, simplex=array([[241.83142594,   0.76274414],
       [225.82002697,   0.76015625],
       [212.41789302,   0.76391602]]))

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

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

residual=pd.DataFrame(result.get("residuals"), columns=["residual Model"])
residual
residual Model
0 -17.596375
1 -46.997192
2 23.093933
... ...
153 -59.003235
154 -91.150085
155 -42.749451

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 预测值列，以及不确定性区间的列。根据您的计算机，这一步应该大约需要1分钟。

# 预测
Y_hat = sf.forecast(horizon, fitted=True)
Y_hat
ds Theta
unique_id
1 1975-01-01 838.559814
1 1975-02-01 800.188232
1 1975-03-01 893.472900
... ... ...
1 1975-10-01 816.166931
1 1975-11-01 786.962036
1 1975-12-01 823.826538

12 rows × 2 columns

values=sf.forecast_fitted_values()
values.head()
ds y Theta
unique_id
1 1962-01-01 589.0 606.596375
1 1962-02-01 561.0 607.997192
1 1962-03-01 640.0 616.906067
1 1962-04-01 656.0 608.873047
1 1962-05-01 727.0 607.395142 
StatsForecast.plot(values)

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

sf.forecast(h=horizon, level=[95])
ds Theta Theta-lo-95 Theta-hi-95
unique_id
1 1975-01-01 838.559814 741.324280 954.365540
1 1975-02-01 800.188232 640.785645 944.996887
1 1975-03-01 893.472900 705.123901 1064.757324
... ... ... ... ...
1 1975-10-01 816.166931 539.706848 1083.791870
1 1975-11-01 786.962036 487.946075 1032.028931
1 1975-12-01 823.826538 512.674866 1101.965942

12 rows × 4 columns

Y_hat=Y_hat.reset_index()
Y_hat
unique_id ds Theta
0 1 1975-01-01 838.559814
1 1 1975-02-01 800.188232
2 1 1975-03-01 893.472900
... ... ... ...
9 1 1975-10-01 816.166931
10 1 1975-11-01 786.962036
11 1 1975-12-01 823.826538

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 Theta
0 1975-01-01 834 1 838.559814
1 1975-02-01 782 1 800.188232
2 1975-03-01 892 1 893.472900
... ... ... ... ...
9 1975-10-01 827 1 816.166931
10 1975-11-01 797 1 786.962036
11 1975-12-01 843 1 823.826538

12 rows × 4 columns

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

带置信区间的预测方法

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

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

  • h（整数）：表示预测未来h步。在此情况下，指的是向前12个月。

  • level（浮点数列表）：此可选参数用于概率预测。设置预测区间的水平（或置信百分位）。例如，level=[95]意味着模型预计真实值将在该区间内95%的时间。

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

此步骤应耗时不到1秒。

sf.predict(h=horizon)
ds Theta
unique_id
1 1975-01-01 838.559814
1 1975-02-01 800.188232
1 1975-03-01 893.472900
... ... ...
1 1975-10-01 816.166931
1 1975-11-01 786.962036
1 1975-12-01 823.826538

12 rows × 2 columns

forecast_df = sf.predict(h=horizon, level=[80,95]) 

forecast_df
ds Theta Theta-lo-80 Theta-hi-80 Theta-lo-95 Theta-hi-95
unique_id
1 1975-01-01 838.559814 765.496155 927.260071 741.324280 954.365540
1 1975-02-01 800.188232 701.729797 898.807434 640.785645 944.996887
1 1975-03-01 893.472900 758.481018 1006.847656 705.123901 1064.757324
... ... ... ... ... ... ...
1 1975-10-01 816.166931 611.404541 991.667419 539.706848 1083.791870
1 1975-11-01 786.962036 561.990845 969.637634 487.946075 1032.028931
1 1975-12-01 823.826538 591.283691 1029.491577 512.674866 1101.965942

12 rows × 6 columns

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

pd.concat([df, forecast_df]).set_index('ds')
y unique_id Theta Theta-lo-80 Theta-hi-80 Theta-lo-95 Theta-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 816.166931 611.404541 991.667419 539.706848 1083.791870
1975-11-01 NaN NaN 786.962036 561.990845 969.637634 487.946075 1032.028931
1975-12-01 NaN NaN 823.826538 591.283691 1029.491577 512.674866 1101.965942

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=2 )
    colors = ['green']
    ax.fill_between(df_plot.index, 
                df_plot['Theta-lo-80'], 
                df_plot['Theta-hi-80'],
                alpha=.20,
                color='lime',
                label='Theta_level_80')
    ax.fill_between(df_plot.index, 
                df_plot['Theta-lo-95'], 
                df_plot['Theta-hi-95'],
                alpha=.2,
                color='white',
                label='Theta_level_95')
    ax.set_title('', fontsize=22)
    ax.set_ylabel("Return", fontsize=20)
    ax.set_xlabel('Month-Days', fontsize=20)
    ax.legend(prop={'size': 15})
    ax.grid(True)
    plt.show()
plot_forecasts(train, test, forecast_df, models=['Theta'])

让我们使用Statsforecast中的plot函数绘制相同的图，如下所示。

sf.plot(df, forecast_df, level=[95])

交叉验证

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

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

下图展示了这样的交叉验证策略：

执行时间序列交叉验证

时间序列模型的交叉验证被视为最佳实践，但大多数实现速度非常慢。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": 包含模型名称和拟合值的列。

评估模型

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

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

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

  1. 实际值。
  2. 预测值，在本例中为Theta。
rmse = rmse(crossvalidation_df['y'], crossvalidation_df["Theta"])
print("RMSE using cross-validation: ", rmse)
RMSE using cross-validation:  12.643162

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

现在我们将使用预测结果来评估我们的模型，我们将使用不同类型的指标 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="Theta")
mae mape mase rmse smape
Theta 8.111287 0.964855 0.36478 9.730347 0.965874

感谢

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

参考文献

  1. Jose A. Fiorucci, Tiago R. Pellegrini, Francisco Louzada, Fotios Petropoulos, Anne B. Koehler (2016). “优化theta方法的模型及其与状态空间模型的关系”。国际预测期刊。
  2. Nixtla参数。
  3. Pandas可用频率。
  4. Rob J. Hyndman和George Athanasopoulos (2018). “预测原理与实践，时间序列交叉验证”。。
  5. 季节性周期 - Rob J Hyndman。

Give us a ⭐ on Github