本文的目的是提供使用 mlforecast 构建 预测区间在预测模型中 的分步指南。

在本演练中,我们将熟悉主要的 MlForecast 类以及一些相关方法,例如 MLForecast.fitMLForecast.predictMLForecast.cross_validation

让我们开始吧!!!

目录

  1. 引言
  2. 预测与预测区间
  3. 安装 mlforecast
  4. 加载库和数据
  5. 使用 plot 方法探索数据
  6. 将数据分割为训练集和测试集
  7. 使用 mlforecast 进行建模
  8. 参考文献

引言

我们的预测目标是未知的事物(否则我们就不会进行预测),因此我们可以将其视为一个随机变量。例如,下个月的总销售额可能有不同的可能值,在月底获得实际销售额之前,我们不会知道确切的值是多少。在下个月的销售额已知之前,这是一个随机量。

当下个月临近时,我们通常对可能的销售值有一个很好的概念。然而,如果我们预测明年同月的销售额,可能的值变化会大得多。在大多数预测情况下,与我们预测目标相关的变异性会随着事件的临近而减少。换句话说,我们预测的时间越早,不确定性就越大。

我们可以想象许多可能的未来场景,每个场景都会产生我们试图预测的不同值。

当我们获得预测时,我们正在估计随机变量可能取的可能值范围的中间值。通常,预测会伴随一个预测区间,该区间给出了随机变量以相对较高的概率可能取值的范围。例如,一个 95% 的预测区间包含一个值范围,该范围应以 95% 的概率包含实际的未来值。

我们通常不绘制单个可能的未来,而是显示这些预测区间。

当我们生成预测时,通常会产生一个称为点预测的单个值。然而,这个值并没有告诉我们与预测相关的不确定性。要衡量这种不确定性,我们需要预测区间。

预测区间是指预测在给定概率下可以取值的范围。因此,95% 的预测区间应包含一个以 95% 概率包含实际未来值的范围。概率预测旨在生成完整的预测分布。另一方面,点预测通常返回该分布的均值或中位数。然而,在实际场景中,最好不仅预测最可能的未来结果,还预测许多替代结果。

问题在于,有些时间序列模型提供预测分布,而另一些则只提供点预测。那么我们如何估计预测的不确定性呢?

预测与预测区间

使用时间序列模型进行预测至少有四个不确定性来源

  1. 随机误差项;
  2. 参数估计;
  3. 历史数据模型的选择;
  4. 历史数据生成过程在未来的延续。

当我们为时间序列模型生成预测区间时,通常只考虑第一个不确定性来源。可以使用模拟来考虑来源 2 和 3,但这几乎从未做过,因为它计算时间太长。随着计算速度的提高,这在未来可能成为一种可行的方法。

即使我们忽略模型不确定性和 DGP 不确定性(来源 3 和 4),而只尝试考虑参数不确定性以及随机误差项(来源 1 和 2),除了某些简单的特殊情况外,也没有封闭形式的解。参见 Rob J Hyndman 的完整文章

预测分布

我们使用预测分布来表达预测中的不确定性。这些概率分布描述了使用拟合模型观察不同未来值的概率。点预测对应于该分布的均值。大多数时间序列模型生成的预测服从正态分布,这意味着我们假设可能的未来值服从正态分布。但是,在本节后面,我们将探讨正态分布的一些替代方案。

时间序列中置信区间预测的重要性

  1. 不确定性估计:置信区间提供了衡量时间序列预测相关不确定性的指标。它能够量化变异性和可能的未来值范围,这对于做出明智决策至关重要。

  2. 精度评估:通过拥有置信区间,可以评估预测的精度。如果区间较窄,则表明预测更准确可靠。另一方面,如果区间较宽,则表明预测的不确定性更大,精度较低。

  3. 风险管理:置信区间通过提供有关可能未来场景的信息来帮助进行风险管理。它允许确定实际值可能所在的范围,并根据这些可能的场景做出决策。

  4. 有效沟通:置信区间是清晰准确地传达预测的有用工具。它可以将与预测相关的变异性和不确定性传达给利益相关者,避免对结果产生错误或过于乐观的解读。

因此,时间序列中的置信区间预测对于理解和管理不确定性、评估预测的准确性以及根据可能的未来场景做出明智决策至关重要。

预测区间

预测区间提供了一个范围,我们期望 yty_t 以指定的概率落入其中。例如,如果我们假设未来观测值的分布服从正态分布,则步骤 h 的预测的 95% 预测区间将表示为该范围

y^T+hT±1.96σ^h,\hat{y}_{T+h|T} \pm 1.96 \hat\sigma_h,

其中 σ^h\hat\sigma_h 是 h 步预测分布标准差的估计值。

更一般地,预测区间可以写成

y^T+hT±cσ^h\hat{y}_{T+h|T} \pm c \hat\sigma_h

在此背景下,“乘数 c” 一词与覆盖概率相关。在本文中,通常计算 80% 和 95% 的区间,但可以使用任何其他百分比。下表显示了对应于不同覆盖概率的 c 值,假设预测分布为正态分布。

百分比乘数
500.67
550.76
600.84
650.93
701.04
751.15
801.28
851.44
901.64
951.96
962.05
972.17
982.33
992.58

预测区间很有价值,因为它们反映了预测中的不确定性。如果只生成点预测,我们就无法评估这些预测的准确性。然而,通过提供预测区间,与每个预测相关的不确定性变得显而易见。因此,如果不包含相应的预测区间,点预测可能缺乏重要价值。

一步预测区间

在对未来步骤进行预测时,可以使用残差的标准差来估计预测分布的标准差,计算方法为

其中 KK 是预测方法中估计的参数数量,MM 是残差中缺失值的数量。(例如,对于朴素预测,M=1M=1,因为我们无法预测第一个观测值。)

多步预测区间

预测区间的一个典型特征是,随着预测范围的延长,它们的长度会增加。随着我们在时间上进一步向前推进,与预测相关的不确定性会增加,从而导致预测区间变宽。一般来说,σh 随着 h 的增加而增加(尽管有些非线性预测方法不遵循此属性)。

要生成预测区间,需要估计 σh。如上所述,对于一步预测 (h=1),公式 (1) 提供了预测标准差 σ1 的良好估计。然而,对于多步预测,需要更复杂的计算方法。这些计算假定残差彼此不相关。

基准方法

对于四种基准方法,可以在残差不相关的假设下数学推导出预测标准差。如果 σ^h\hat{\sigma}_h 表示 h 步预测分布的标准差,σ^\hat{\sigma} 是由 (1) 给出的残差标准差,那么我们可以使用下表所示的表达式。请注意,当 h=1h=1TT 较大时,这些都给出相同的近似值 σ^\hat{\sigma}

方法h 步预测标准差
均值预测σ^h=σ^1+1/T\hat\sigma_h = \hat\sigma\sqrt{1 + 1/T}
朴素预测σ^h=σ^h\hat\sigma_h = \hat\sigma\sqrt{h}
季节性朴素预测σ^h=σ^k+1\hat\sigma_h = \hat\sigma\sqrt{k+1}
漂移预测σ^h=σ^h(1+h/T)\hat\sigma_h = \hat\sigma\sqrt{h(1+h/T)}

请注意,当 h=1h=1TT 较大时,这些都给出相同的近似值 σ^\hat{\sigma}

来自自举残差的预测区间

当残差的正态分布是一个不合理的假设时,一个替代方法是使用自举法,它只假设残差不相关且方差恒定。我们将使用朴素预测方法来说明该过程。

一步预测误差定义为 et=yty^tt1e_t = y_t - \hat{y}_{t|t-1}。对于朴素预测方法,y^tt1=yt1\hat{y}_{t|t-1} = y_{t-1},因此我们可以将其重写为 yt=yt1+et.y_t = y_{t-1} + e_t.

假设未来的误差与过去的误差相似,当 t>Tt>T 时,我们可以通过从过去观察到的误差集合(即残差)中采样来替换 ete_{t}。因此,我们可以使用以下方法模拟时间序列的下一个观测值

yT+1=yT+eT+1y^*_{T+1} = y_{T} + e^*_{T+1}

其中 eT+1e^*_{T+1} 是从过去随机采样的误差,yT+1y^*_{T+1} 是如果发生该特定误差值可能出现的未来值。我们使用 * 来表示这并非观测到的 yT+1y_{T+1} 值,而是可能发生的未来情况之一。将新的模拟观测值添加到我们的数据集中,我们可以重复该过程以获得

yT+2=yT+1+eT+2,y^*_{T+2} = y_{T+1}^* + e^*_{T+2},

其中 eT+2e^*_{T+2} 是从残差集合中进行的另一次抽取。以此类推,我们可以模拟时间序列的一整套未来值。

保形预测

多分位数损失和统计模型可以提供预测区间,但问题在于这些区间未经校准,这意味着落在区间内的观测值的实际频率与与之相关的置信水平不一致。例如,一个经过校准的 95% 预测区间在重复抽样中应包含真实值 95% 的时间。另一方面,未经校准的 95% 预测区间可能仅包含真实值 80% 的时间,或者可能是 99% 的时间。在第一种情况下,区间太窄,低估了不确定性;而在第二种情况下,区间太宽,高估了不确定性。

统计方法也假设正态性。在这里,我们讨论另一种称为保形预测的方法,它不需要任何分布假设。

保形预测区间使用点预测器模型的交叉验证来生成区间。这意味着无需先验概率,输出经过良好校准。无需额外的训练,模型被视为黑盒。该方法与任何模型兼容

mlforecast 现在支持所有可用模型的保形预测。

安装 mlforecast

  • 使用 pip

    • pip install mlforecast
  • 使用 conda

    • conda install -c conda-forge mlforecast

加载库和数据

# Handling and processing of Data
# ==============================================================================
import numpy as np
import pandas as pd

import scipy.stats as stats

# Handling and processing of Data for Date (time)
# ==============================================================================
import datetime
import time
from datetime import datetime, timedelta

# 
# ==============================================================================
from statsmodels.tsa.stattools import adfuller
import statsmodels.api as sm
import statsmodels.tsa.api as smt
from statsmodels.tsa.seasonal import seasonal_decompose 
# 
# ==============================================================================
from utilsforecast.plotting import plot_series
import xgboost as xgb

from mlforecast import MLForecast
from mlforecast.lag_transforms import ExpandingMean, ExponentiallyWeightedMean, RollingMean
from mlforecast.target_transforms import Differences
from mlforecast.utils import PredictionIntervals
# Plot
# ==============================================================================
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

读取数据

data_url = "https://raw.githubusercontent.com/Naren8520/Serie-de-tiempo-con-Machine-Learning/main/Data/nyc_taxi.csv"
df = pd.read_csv(data_url, parse_dates=["timestamp"])
df.head()
时间戳
02014-07-01 00:00:0010844
12014-07-01 00:30:008127
22014-07-01 01:00:006210
32014-07-01 01:30:004656
42014-07-01 02:00:003820

MlForecast 的输入始终是长格式的数据框,包含三列: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
02014-07-01 00:00:00108441
12014-07-01 00:30:0081271
22014-07-01 01:00:0062101
32014-07-01 01:30:0046561
42014-07-01 02:00:0038201
df.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 10320 entries, 0 to 10319
Data columns (total 3 columns):
 #   Column     Non-Null Count  Dtype         
---  ------     --------------  -----         
 0   ds         10320 non-null  datetime64[ns]
 1   y          10320 non-null  int64         
 2   unique_id  10320 non-null  object        
dtypes: datetime64[ns](1), int64(1), object(1)
memory usage: 242.0+ KB

使用 plot 方法探索数据

使用 StatsForecast 类中的 plot 方法绘制一些序列。此方法会打印数据集中的 8 个随机序列,对于基本 EDA 很有用。

fig = plot_series(df)

增强迪基-富勒检验

增强迪基-富勒 (ADF) 检验是一种统计检验,用于确定时间序列数据中是否存在单位根。单位根可能导致时间序列分析结果不可预测。在单位根检验中形成一个零假设,以确定时间序列数据受趋势影响的强度。通过接受零假设,我们接受时间序列数据非平稳的证据。通过拒绝零假设或接受备择假设,我们接受时间序列数据是由平稳过程生成的证据。这个过程也称为平稳趋势。ADF 检验统计量的值是负的。较低的 ADF 值表示对零假设的更强拒绝。

增强迪基-富勒检验 (ADF) 是一种常用的统计检验,用于检验给定的时间序列是否平稳。我们可以通过定义零假设和备择假设来实现这一点。

  • 零假设:时间序列非平稳。它呈现时间依赖的趋势。

  • 备择假设:时间序列平稳。换句话说,该序列不依赖于时间。

  • ADF 或 t 统计量 < 临界值:拒绝零假设,时间序列平稳。

  • ADF 或 t 统计量 > 临界值:未能拒绝零假设,时间序列非平稳。

def augmented_dickey_fuller_test(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(df["y"],'Ads')
Dickey-Fuller test results for columns: Ads
Test Statistic                -1.076452e+01
p-value                        2.472132e-19
No Lags Used                   3.900000e+01
Number of observations used    1.028000e+04
Critical Value (1%)           -3.430986e+00
Critical Value (5%)           -2.861821e+00
Critical Value (10%)          -2.566920e+00
dtype: float64
Conclusion:====>
Reject the null hypothesis
The data is stationary

自相关图

自相关函数

定义 1.{xt;1tn}\{x_t;1 ≤ t ≤ n\} 是来自 {Xt}\{X_t\} 的大小为 n 的时间序列样本。 1. xˉ=t=1nxtn\bar x = \sum_{t=1}^n \frac{x_t}{n} 称为 {Xt}\{X_t\} 的样本均值。 2. ck=t=1nk(xt+kxˉ)(xtxˉ)/nc_k =\sum_{t=1}^{n−k} (x_{t+k}- \bar x)(x_t−\bar x)/n 称为 {Xt}\{X_t\} 的样本自协方差函数。 3. rk=ck/c0r_k = c_k /c_0 称为 {Xt}\{X_t\} 的样本自相关函数。

关于此定义的说明如下

  • 与大多数文献类似,本指南使用 ACF 表示样本自相关函数和自相关函数。通过上下文很容易识别 ACF 表示的内容。

  • 显然 c0 是 {Xt}\{X_t\} 的样本方差。此外,r0=c0/c0=1r_0 = c_0/c_0 = 1,且对于任何整数 k,rk1k, |r_k| ≤ 1

  • 当计算任何固定长度为 n 的样本序列的 ACF 时,对于较大的 k 值,我们不能对其 rkr_k 值过于自信,因为当 kk 较大时,可用于计算 rkr_k(xt+k,xt)(x_{t +k }, x_t ) 对会更少。一个经验法则是不要估计 k>n/3k > n/3 时的 rkr_k,另一个经验法则是当 n50n ≥ 50 时,kn/4k ≤ n/4。无论如何,小心总是好的。

  • 我们也可以通过定义 1 计算非平稳时间序列样本的 ACF。然而,在这种情况下,随着 k 的增加,ACF 或 rkr_k 衰减得非常慢或几乎不衰减。

  • 将 ACF (rk)(r_k) 对滞后 kk 作图很容易,但在分析时间序列样本时非常有用。这种 ACF 图被称为相关图。

  • 如果 {Xt}\{X_t\} 是平稳的,且对于所有 k0k \neq 0,有 E(Xt)=0E(X_t)=0ρk=0\rho_k =0,即它是白噪声序列,那么 rkr_k 的抽样分布渐近正态,均值为 0,方差为 1/n1/n。因此,rkr_k 落入区间 [1.96/n,1.96/n][−1.96/√n, 1.96/√n] 的概率约为 95%。

现在我们可以总结如下:(1) 如果时间序列图清楚地显示出趋势或/和季节性,则它肯定是非平稳的;(2) 如果随着滞后 kk 的增加,ACF rkr_k 衰减得非常慢或几乎不衰减,则该时间序列也应是非平稳的。

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("../../figs/prediction_intervals_in_forecasting_models__autocorrelation.png", bbox_inches='tight')
plt.close();

时间序列分解

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

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

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

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

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

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

加法时间序列模型

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

乘法时间序列模型

如果时间序列的分量是相乘的,则该时间序列称为乘法时间序列模型。通过可视化,如果时间序列随时间呈现指数增长或下降,则该时间序列可以被认为是乘法时间序列模型。乘法时间序列模型的数学函数可以表示为:

y(t)=LevelTrendseasonalityNoisey(t) = Level * Trend * seasonality * Noise

加法

a = seasonal_decompose(df["y"], model = "additive", period=24).plot()
a.savefig('../../figs/prediction_intervals_in_forecasting_models__seasonal_decompose_aditive.png', bbox_inches='tight')
plt.close()

乘法

b = seasonal_decompose(df["y"], model = "Multiplicative", period=24).plot()
b.savefig('../../figs/prediction_intervals_in_forecasting_models__seasonal_decompose_multiplicative.png', bbox_inches='tight')
plt.close();

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

我们将数据分为以下几个集合:1. 用于训练模型的数据。2. 用于测试模型的数据。

对于测试数据,我们将使用最后 500 小时的数据来测试和评估模型的性能。

train = df[df.ds<='2015-01-21 13:30:00'] 
test = df[df.ds>'2015-01-21 13:30:00']
train.shape, test.shape
((9820, 3), (500, 3))

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

fig = plot_series(train,test)

使用 mlforecast 进行建模

构建模型

我们定义要使用的模型,在本例中我们将使用 XGBoost 模型。

model1 = [xgb.XGBRegressor()]

我们可以使用 MLForecast.preprocess 方法来探索不同的变换。

如果确实我们正在处理的序列是平稳序列(参见 Dickey-Fuller 检验),然而,为了本指南的练习和说明目的,我们将对序列应用差分。我们将使用 target_transforms 参数并调用 diff 函数来实现,例如:mlforecast.target_transforms.Differences

mlf = MLForecast(models=model1,
                 freq='30min', 
                 target_transforms=[Differences([1])],
                 )

在使用参数 target_transforms=[Differences([1])] 时,需要考虑的是如果序列是平稳的,我们可以使用一阶差分;如果序列是非平稳的,我们可以使用多个差分,以使序列随时间保持恒定,即均值和方差保持恒定。

prep = mlf.preprocess(df)
prep
dsyunique_id
12014-07-01 00:30:00-2717.01
22014-07-01 01:00:00-1917.01
32014-07-01 01:30:00-1554.01
42014-07-01 02:00:00-836.01
52014-07-01 02:30:00-947.01
103152015-01-31 21:30:00951.01
103162015-01-31 22:00:001051.01
103172015-01-31 22:30:001588.01
103182015-01-31 23:00:00-718.01
103192015-01-31 23:30:00-303.01

这已经从每个值中减去了滞后 1 的值,现在我们可以看看我们的序列是什么样子的。

fig = plot_series(prep)

添加特征

滞后特征

季节性似乎已经消失了,我们现在可以尝试添加一些滞后特征。

mlf = MLForecast(models=model1,
                 freq='30min',  
                 lags=[1,24],
                 target_transforms=[Differences([1])],
                 )
prep = mlf.preprocess(df)
prep
dsyunique_id滞后1滞后24
252014-07-01 12:30:00-22.01445.0-2717.0
262014-07-01 13:00:00-708.01-22.0-1917.0
272014-07-01 13:30:001281.01-708.0-1554.0
282014-07-01 14:00:0087.011281.0-836.0
292014-07-01 14:30:001045.0187.0-947.0
103152015-01-31 21:30:00951.01428.04642.0
103162015-01-31 22:00:001051.01951.0-519.0
103172015-01-31 22:30:001588.011051.02411.0
103182015-01-31 23:00:00-718.011588.0214.0
103192015-01-31 23:30:00-303.01-718.02595.0
prep.drop(columns=['unique_id', 'ds']).corr()['y']
y        1.000000
lag1     0.663082
lag24    0.155366
Name: y, dtype: float64

滞后变换

滞后变换被定义为一个字典,其中键是滞后值,值是我们要应用于该滞后的变换列表。更多详细信息请参阅滞后变换指南

mlf = MLForecast(models=model1,
                 freq='30min',  
                 lags=[1,24],
                 lag_transforms={1: [ExpandingMean()], 24: [RollingMean(window_size=7)]},
                 target_transforms=[Differences([1])],
                 )
prep = mlf.preprocess(df)
prep
dsyunique_id滞后1滞后24expanding_mean_lag1rolling_mean_lag24_window_size7
312014-07-01 15:30:00-836.01-1211.0-305.0284.533325-1254.285767
322014-07-01 16:00:00-2316.01-836.0157.0248.387100-843.714294
332014-07-01 16:30:00-1215.01-2316.0-63.0168.250000-578.857117
342014-07-01 17:00:002190.01-1215.0357.0126.333336-305.857147
352014-07-01 17:30:002322.012190.01849.0187.02941977.714287
103152015-01-31 21:30:00951.01428.04642.01.2483032064.285645
103162015-01-31 22:00:001051.01951.0-519.01.3403781873.428589
103172015-01-31 22:30:001588.011051.02411.01.4421292179.000000
103182015-01-31 23:00:00-718.011588.0214.01.5959101888.714233
103192015-01-31 23:30:00-303.01-718.02595.01.5261682071.714355

您可以看到这两种方法得到相同的结果,您可以使用您觉得最舒服的一种。

日期特征

如果您的时间列由时间戳组成,那么提取诸如 week、dayofweek、quarter 等特征可能是有意义的。您可以通过传入包含 pandas 时间/日期组件的字符串列表来实现这一点。您也可以传入将时间列作为输入的函数,正如我们将在这里展示的那样。

mlf = MLForecast(models=model1,
                 freq='30min', 
                 lags=[1,24],
                 lag_transforms={1: [ExpandingMean()], 24: [RollingMean(window_size=7)]},
                 target_transforms=[Differences([1])],
                 date_features=["year", "month", "day", "hour"]) # Seasonal data
prep = mlf.preprocess(df)
prep
dsyunique_id滞后1滞后24expanding_mean_lag1rolling_mean_lag24_window_size7年 (year)月 (month)日 (day)小时 (hour)
312014-07-01 15:30:00-836.01-1211.0-305.0284.533325-1254.28576720147115
322014-07-01 16:00:00-2316.01-836.0157.0248.387100-843.71429420147116
332014-07-01 16:30:00-1215.01-2316.0-63.0168.250000-578.85711720147116
342014-07-01 17:00:002190.01-1215.0357.0126.333336-305.85714720147117
352014-07-01 17:30:002322.012190.01849.0187.02941977.71428720147117
103152015-01-31 21:30:00951.01428.04642.01.2483032064.285645201513121
103162015-01-31 22:00:001051.01951.0-519.01.3403781873.428589201513122
103172015-01-31 22:30:001588.011051.02411.01.4421292179.000000201513122
103182015-01-31 23:00:00-718.011588.0214.01.5959101888.714233201513123
103192015-01-31 23:30:00-303.01-718.02595.01.5261682071.714355201513123

拟合模型

# fit the models
mlf.fit(df,  
 fitted=True, 
prediction_intervals=PredictionIntervals(n_windows=5, h=30, method="conformal_distribution" )  )
MLForecast(models=[XGBRegressor], freq=30min, lag_features=['lag1', 'lag24', 'expanding_mean_lag1', 'rolling_mean_lag24_window_size7'], date_features=['year', 'month', 'day', 'hour'], num_threads=1)

让我们看看我们模型的结果,在本例中是 XGBoost 模型。我们可以通过以下指令来观察它

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

result = mlf.forecast_fitted_values()
result = result.set_index("unique_id")
result
dsyXGBRegressor
unique_id
12014-07-01 15:30:0018544.018441.443359
12014-07-01 16:00:0016228.016391.152344
12014-07-01 16:30:0015013.015260.714844
12014-07-01 17:00:0017203.017066.148438
12014-07-01 17:30:0019525.019714.404297
12015-01-31 21:30:0024670.024488.646484
12015-01-31 22:00:0025721.025868.865234
12015-01-31 22:30:0027309.027290.125000
12015-01-31 23:00:0026591.027123.226562
12015-01-31 23:30:0026288.026241.205078
from statsmodels.stats.diagnostic import normal_ad
from scipy import stats
sw_result = stats.shapiro(result["XGBRegressor"])
ad_result = normal_ad(np.array(result["XGBRegressor"]), axis=0)
dag_result = stats.normaltest(result["XGBRegressor"], axis=0, nan_policy='propagate')

需要注意的是,只有当我们假设验证预测的残差呈正态分布时,才能使用此方法。为了查看是否如此,我们将使用 PP 图,并通过 Anderson-Darling、Kolmogorov-Smirnov 和 D’Agostino K^2 检验来测试其正态性。

PP 图(Probability-to-Probability)将数据样本与正态分布图进行对比绘制,使得如果数据呈正态分布,数据点将形成一条直线。

这三个正态性检验使用 p 值来确定数据样本来自正态分布总体的可能性。每个检验的零假设是“样本来自正态分布总体”。这意味着如果得到的 p 值低于所选的 alpha 值,则拒绝零假设。因此,有证据表明数据来自非正态分布。对于本文,我们将使用 0.01 的 Alpha 值。

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

# plot[1,1]
result["XGBRegressor"].plot(ax=axs[0,0])
axs[0,0].set_title("Residuals model");

# plot
axs[0,1].hist(result["XGBRegressor"], density=True,bins=50, alpha=0.5 )
axs[0,1].set_title("Density plot - Residual");

# plot
stats.probplot(result["XGBRegressor"], dist="norm", plot=axs[1,0])
axs[1,0].set_title('Plot Q-Q')
axs[1,0].annotate("SW p-val: {:.4f}".format(sw_result[1]), xy=(0.05,0.9), xycoords='axes fraction', fontsize=15,
            bbox=dict(boxstyle="round", fc="none", ec="gray", pad=0.6))

axs[1,0].annotate("AD p-val: {:.4f}".format(ad_result[1]), xy=(0.05,0.8), xycoords='axes fraction', fontsize=15,
            bbox=dict(boxstyle="round", fc="none", ec="gray", pad=0.6))

axs[1,0].annotate("DAG p-val: {:.4f}".format(dag_result[1]), xy=(0.05,0.7), xycoords='axes fraction', fontsize=15,
            bbox=dict(boxstyle="round", fc="none", ec="gray", pad=0.6))
# plot
plot_acf(result["XGBRegressor"],  lags=35, ax=axs[1,1],color="fuchsia")
axs[1,1].set_title("Autocorrelation");

plt.savefig("../../figs/prediction_intervals_in_forecasting_models__plot_residual_model.png", bbox_inches='tight')
plt.close();

带预测区间的预测方法

使用 predict 方法生成预测。

forecast_df = mlf.predict(h=30, level=[80,95])
forecast_df.head()
unique_iddsXGBRegressorXGBRegressor-lo-95XGBRegressor-lo-80XGBRegressor-hi-80XGBRegressor-hi-95
012015-02-01 00:00:0026320.29882825559.88424125680.22836926960.36928727080.713416
112015-02-01 00:30:0026446.47265624130.42961425195.46162127697.48369128762.515698
212015-02-01 01:00:0024909.97070323094.95053723579.58339826240.35800826724.990869
312015-02-01 01:30:0024405.40234421548.62829622006.66259826804.14209027262.176392
412015-02-01 02:00:0022292.39062520666.73696321130.21543023454.56582023918.044287

绘制预测区间

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

fig = plot_series(df, forecast_df, level=[80,95], max_insample_length=200,engine="matplotlib")
fig.get_axes()[0].set_title("Prediction intervals")
fig.savefig('../../figs/prediction_intervals_in_forecasting_models__plot_forecasting_intervals.png', bbox_inches='tight')

置信区间是一个值范围,具有高概率包含变量的真实值。在机器学习时间序列模型中,置信区间用于估计预测中的不确定性。

使用置信区间的主要好处之一是它允许用户了解预测的准确性。例如,如果置信区间很宽,则意味着预测不太准确。相反,如果置信区间很窄,则意味着预测更准确。

置信区间的另一个好处是它有助于用户做出明智的决策。例如,如果预测在置信区间内,则意味着它很可能实现。相反,如果预测超出置信区间,则意味着它不太可能实现。

总的来说,置信区间是机器学习时间序列模型的重要工具。它帮助用户了解预测的准确性并做出明智的决策。

参考文献

  1. Changquan Huang • Alla Petukhina. Springer series (2022). 应用时间序列分析和 Python 预测。
  2. Ivan Svetunkov. 使用增强动态自适应模型 (ADAM) 进行预测和分析
  3. James D. Hamilton. 时间序列分析 Princeton University Press, Princeton, New Jersey, 1st Edition, 1994.
  4. 用于 Mlforecast 的 Nixtla 参数.
  5. Pandas 可用的频率.
  6. Rob J. Hyndman and George Athanasopoulos (2018). “预测原理与实践,时间序列交叉验证”。.
  7. 季节性周期 - Rob J Hyndman.