# 数据的处理与加工
# ==============================================================================
import numpy as np
import pandas as pd
import scipy.stats as stats
# 日期（时间）数据的处理与加工
# ==============================================================================
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
预测区间
以下文章的目的是获得关于使用mlforecast构建预测模型中的预测区间的逐步指南。
在此过程中，我们将熟悉主要的MlForecast类以及一些相关的方法，如MLForecast.fit、MLForecast.predict和MLForecast.cross_validation等。
让我们开始吧！！！
目录
介绍
我们预测的目标是一些未知的东西（否则我们就不会进行预测），因此我们可以将其视为一个随机变量。例如，下个月的总销售额可能有不同的可能值，而我们直到月底获得实际销售额之前都不知道确切的值。在下个月的销售额未确定之前，这就是一个随机数。
当下个月临近时，我们通常对可能的销售值有一个相当好的了解。然而，如果我们是在预测明年同月的销售额，可能的值就会变化得更多。在大多数预测情况下，随着我们距离事件的临近，预测相关的变动性会减少。换句话说，我们在时间上越往回进行预测，不确定性就越大。
我们可以想象许多可能的未来情景，每个情景都为我们试图预测的内容带来不同的值。
当我们获得预测时，我们是在估计随机变量可能取值范围的中间值。通常，预测伴随着一个预测区间，给出了随机变量可能以相对较高的概率取值的范围。例如，95%的预测区间包含的值范围应当以95%的概率包括实际未来值。
我们通常不绘制单个可能的未来，而是显示这些预测区间。
当我们生成预测时，通常会产生一个被称为点预测的单一值。然而，这个值并没有告诉我们与预测相关的不确定性。为了衡量这种不确定性，我们需要预测区间。
预测区间是一个值范围，该范围是在给定概率下预测可以取的值。因此，95%的预测区间应该包含一个值范围，以95%的概率包括实际未来值。概率预测旨在生成完整的预测分布。另一方面，点预测通常返回该分布的均值或中位数。然而，在现实场景中，最好不仅预测最可能的未来结果，还预测许多替代结果。
问题是，一些时间序列模型提供预测分布，而其他模型仅提供点预测。那么我们该如何估计预测的不确定性呢？
预测和预测区间
在使用时间序列模型进行预测时，至少有四个不确定性来源：
- 随机误差项；
- 参数估计；
- 对历史数据的模型选择；
- 历史数据生成过程向未来的延续。
当我们为时间序列模型生成预测区间时，通常只考虑以上不确定性来源中的第一个。可以通过模拟来考虑来源2和3，但几乎从不这样做，因为计算会耗费太多时间。随着计算速度的提高，在未来这可能成为一种可行的方法。
即使我们忽略模型不确定性和数据生成过程的不确定性（来源3和4），仅仅考虑参数不确定性以及随机误差项（来源1和2），除了某些简单的特例外，也没有封闭形式的解决方案。更多信息请参见完整文章 Rob J Hyndman
预测分布
我们使用预测分布来表达我们预测中的不确定性。这些概率分布描述了使用拟合模型观察到不同未来值的概率。点预测对应于这个分布的均值。大多数时间序列模型生成遵循正态分布的预测，这意味着我们假设可能的未来值遵循正态分布。然而，在本节后面，我们将查看一些正态分布的替代方案。
时间序列中置信区间预测的重要性：
不确定性估计：置信区间提供了与时间序列预测相关的不确定性度量。它使变异性和可能的未来值范围可以量化，这对于做出明智的决策至关重要。
精度评估：通过拥有置信区间，可以评估预测的精度。如果区间狭窄，表明预测更准确和可靠。另一方面，如果区间宽广，则表示预测的不确定性和精度较低。
风险管理：置信区间通过提供有关可能未来场景的信息来帮助风险管理。它允许识别真实值可能位于的范围，并基于这些可能场景做出决策。
有效沟通：置信区间是清晰准确地传达预测的有用工具。它允许将与预测相关的变异性和不确定性传达给利益相关者，避免对结果的误解或过于乐观的解读。
因此，在时间序列中，置信区间预测对于理解和管理不确定性、评估预测的准确性以及基于可能的未来场景做出明智的决策至关重要。
预测区间
预测区间为我们提供了一个范围，我们期望 \(y_t\) 以指定的概率落在这个范围内。例如，如果我们假设未来观测值的分布遵循正态分布，则步长 h 的 95% 预测区间可以表示为以下范围
\[\hat{y}_{T+h|T} \pm 1.96 \hat\sigma_h,\]
其中 \(\hat\sigma_h\) 是 h 步预测分布标准差的估计值。
更一般地，预测区间可以写成
\[\hat{y}_{T+h|T} \pm c \hat\sigma_h\]
在这个上下文中，“乘数 c” 与覆盖概率相关。在本文中，通常计算 80% 和 95% 的区间，但可以使用任何其他百分比。下表显示了对应于不同覆盖概率的 c 值，假设预测分布呈正态分布。
| 百分比 | 乘数 |
|---|---|
| 50 | 0.67 |
| 55 | 0.76 |
| 60 | 0.84 |
| 65 | 0.93 |
| 70 | 1.04 |
| 75 | 1.15 |
| 80 | 1.28 |
| 85 | 1.44 |
| 90 | 1.64 |
| 95 | 1.96 |
| 96 | 2.05 |
| 97 | 2.17 |
| 98 | 2.33 |
| 99 | 2.58 |
预测区间是有价值的，因为它们反映了预测中的不确定性。如果我们只生成点预测，就无法评估这些预测的准确性。然而，通过提供预测区间，每个预测相关的不确定性变得明显。因此，点预测在没有对应预测区间的情况下可能缺乏重要价值。
一步预测区间
在对未来步骤进行预测时，可以使用残差的标准差来估算预测分布的标准差，具体计算方式为
\[\begin{equation} \hat{\sigma} = \sqrt{\frac{1}{T-K-M}\sum_{t=1}^T e_t^2}, \tag{1} \end{equation}\]
其中 \(K\) 是在预测方法中估计的参数数量，\(M\) 是残差中缺失值的数量。（例如，对于简单预测，\(M=1\)，因为我们无法预测第一个观察值。）
多步预测区间
预测区间的一个典型特征是，当预测范围延长时，其长度往往会增加。随着时间的推移，预测的不确定性增大，导致预测区间变得更宽。一般而言，σh 趋向于随着 h 的增加而增加（尽管有一些非线性预测方法不遵循这一特性）。
要生成一个预测区间，需要有 σh 的估计。如上所述，对于一步预测（h=1），公式(1)提供了预测标准差 σ1 的良好估计。然而，对于多步预测，则需要更复杂的计算方法。这些计算假设残差之间没有相关性。
基准方法
对于四种基准方法，可以在假设残差不相关的条件下，数学推导出预测标准差。如果 \(\hat{\sigma}_h\) 表示 h 步预测分布的标准差，而 \(\hat{\sigma}\) 是由 (1) 给出的残差标准差，则我们可以使用下一表中所示的表达式。注意，当 \(h=1\) 且 \(T\) 较大时，这些都给出了相同的近似值 \(\hat{\sigma}\)。
| 方法 | h步预测标准差 |
|---|---|
| 均值预测 | \(\hat\sigma_h = \hat\sigma\sqrt{1 + 1/T}\) |
| 简单预测 | \(\hat\sigma_h = \hat\sigma\sqrt{h}\) |
| 季节性简单预测 | \(\hat\sigma_h = \hat\sigma\sqrt{k+1}\) |
| 漂移预测 | \(\hat\sigma_h = \hat\sigma\sqrt{h(1+h/T)}\) |
注意，当 \(h=1\) 且 \(T\) 较大时，这些都给出了相同的近似值 \(\hat{\sigma}\)。
基于自助抽样残差的预测区间
当对残差的正态分布假设不合理时，一种替代方法是使用自助抽样，这只假设残差是不相关的并且具有恒定的方差。我们将使用一种简单的预测方法来说明这一过程。
一步预测误差定义为 \(e_t = y_t - \hat{y}_{t|t-1}\)。对于简单的预测方法， \(\hat{y}_{t|t-1} = y_{t-1}\)，因此我们可以将其重写为 \[y_t = y_{t-1} + e_t.\]
假设未来的误差将与过去的误差相似，当 \(t>T\) 时，我们可以通过从我们过去看到的误差集合中抽样来替换 \(e_{t}\)（即残差）。因此，我们可以使用以下公式模拟时间序列的下一个观测值：
\[y^*_{T+1} = y_{T} + e^*_{T+1}\]
其中 \(e^*_{T+1}\) 是从过去随机抽取的误差，\(y^*_{T+1}\) 是如果发生该特定误差值时可能出现的未来值。我们使用 “*” 来表示这不是观察到的 \(y_{T+1}\) 值，而是可能发生的一个未来值。将新的模拟观测值添加到我们的数据集中，我们可以重复该过程以获得
\[y^*_{T+2} = y_{T+1}^* + e^*_{T+2},\]
其中 \(e^*_{T+2}\) 是从残差集合中再次抽取的另一个值。以这种方式继续，我们可以模拟出时间序列的一整组未来值。
拟合预测
多分位损失和统计模型可以提供预测区间，但问题在于这些区间是未经校准的，这意味着实际落在区间内的观察频率与其相关的置信水平并不一致。例如，一个经过校准的95%预测区间应该在重复抽样中95%的时间包含真实值。另一方面，一个未经校准的95%预测区间可能只在80%的时间内包含真实值，或者可能在99%的时间内包含。在第一种情况下，区间过窄并低估了不确定性，而在第二种情况下，区间过宽并高估了不确定性。
统计方法通常假设正态分布。在这里，我们讨论另一种不需要任何分布假设的方法，称为拟合预测。
拟合预测区间使用交叉验证对点预测模型进行生成。这意味着不需要先验概率，输出是良好校准的。不需要额外的训练，模型被视为一个黑箱。该方法与任何模型兼容。
mlforecast现在支持所有可用模型上的拟合预测。
安装 mlforecast
使用 pip：
pip install mlforecast
使用 conda：
conda install -c conda-forge mlforecast
加载库和数据
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# 情节
# ==============================================================================
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
from statsmodels.graphics.tsaplots import plot_acf, plot_pacfplt.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)
# 定义地块大小
# ==============================================================================
plt.rcParams['figure.figsize'] = (18,7)import warnings# 隐藏警告
warnings.filterwarnings("ignore")读取数据
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()| timestamp | value | |
|---|---|---|
| 0 | 2014-07-01 00:00:00 | 10844 |
| 1 | 2014-07-01 00:30:00 | 8127 |
| 2 | 2014-07-01 01:00:00 | 6210 |
| 3 | 2014-07-01 01:30:00 | 4656 |
| 4 | 2014-07-01 02:00:00 | 3820 |
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()| ds | y | unique_id | |
|---|---|---|---|
| 0 | 2014-07-01 00:00:00 | 10844 | 1 |
| 1 | 2014-07-01 00:30:00 | 8127 | 1 |
| 2 | 2014-07-01 01:00:00 | 6210 | 1 |
| 3 | 2014-07-01 01:30:00 | 4656 | 1 |
| 4 | 2014-07-01 02:00:00 | 3820 | 1 |
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)fig.savefig('../../figs/prediction_intervals_in_forecasting_models__eda.png', bbox_inches='tight')
增强型迪基-福勒检验
增强型迪基-福勒（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. 设 \(\{x_t;1 ≤ t ≤ n\}\) 是来自 \(\{X_t\}\) 的大小为 n 的时间序列样本。 1. \(\bar x = \sum_{t=1}^n \frac{x_t}{n}\) 被称为 \(\{X_t\}\) 的样本均值。 2. \(c_k =\sum_{t=1}^{n−k} (x_{t+k}- \bar x)(x_t−\bar x)/n\) 被称为 \(\{X_t\}\) 的样本自协方差函数。 3. \(r_k = c_k /c_0\) 被称为 \(\{X_t\}\) 的样本自相关函数。
关于这个定义注意以下几点：
像大多数文献一样，本指南使用 ACF 来表示样本自相关函数以及自相关函数。在上下文中可以很容易识别出 ACF 代表的内容。
显然，c0 是 \(\{X_t\}\) 的样本方差。此外，\(r_0 = c_0/c_0 = 1\)，对于任何整数 k，\(|r_k| ≤ 1\)。
当我们计算长度为 n 的任何样本序列的 ACF 时，对于大的 k，我们无法过于自信于 \(r_k\) 的值，因为当 k 很大时，可用于计算 \(r_k\) 的 \((x_{t +k }, x_t )\) 的对数较少。一条经验法则是不对 \(r_k\) 进行估计当 \(k > n/3\)，另一条是 \(n ≥ 50, k ≤ n/4\)。无论如何，保持谨慎总是一个好主意。
我们还可以通过定义 1 来计算非平稳时间序列样本的 ACF。然而，在这种情况下，随着 k 的增加，ACF 或 \(r_k\) 的衰减非常缓慢或几乎没有。
将 ACF \((r_k)\) 绘制为滞后 k 的图形很简单，但在分析时间序列样本方面非常有帮助。这样的 ACF 图称为自相关图。
如果 \(\{X_t\}\) 是平稳的，且 \(E(X_t)=0\) 并且对于所有 \(k \neq 0\) 有 \(\rho_k =0\)，也就是说它是一个白噪声序列，则 \(r_k\) 的采样分布渐近于正态分布，均值为 0，方差为 \(1/n\)。因此，有大约 95% 的机会使得 \(r_k\) 落在区间 \([−1.96/√n, 1.96/√n]\) 内。
现在我们可以总结一下： (1) 如果时间序列图清晰地显示出趋势或/和季节性，那么它肯定是非平稳的；(2) 如果 ACF \(r_k\) 随着滞后 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");
# 格拉菲科
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();
时间序列的分解
如何分解时间序列以及为什么要这样做？
在时间序列分析中，为了预测新值，了解过去的数据非常重要。更正式地说，了解数值随时间变化所遵循的模式是非常重要的。导致我们的预测值出现错误方向的原因可能有很多。基本上，时间序列由四个组件组成。这些组件的变化导致时间序列模式的改变。这些组件是：
- 水平：这是随时间平均的主要值。
- 趋势：趋势是导致时间序列中出现增加或减少模式的值。
- 季节性：这是在短时间内在时间序列中发生的周期性事件，导致时间序列中出现短期增加或减少模式。
- 残差/噪声：这些是时间序列中的随机变化。
随着时间的推移，组合这些组件形成时间序列。大多数时间序列由水平和噪声/残差组成，而趋势或季节性是可选值。
如果季节性和趋势是时间序列的一部分，那么将会对预测值产生影响。因为预测时间序列的模式可能与先前的时间序列不同。
时间序列中组件的组合可以有两种类型： * 加法型 * 乘法型
加法时间序列
如果时间序列的组件是相加形成时间序列的，那么这个时间序列被称为加法时间序列。通过可视化，可以说如果时间序列的增加或减少模式在整个系列中是相似的，那么该时间序列就是加法的。任何加法时间序列的数学函数可以表示为： \[y(t) = 水平 + 趋势 + 季节性 + 噪声\]
乘法时间序列
如果时间序列的组件是相乘形成的，那么这个时间序列被称为乘法时间序列。通过可视化，如果时间序列随着时间的推移呈现指数增长或下降，那么这个时间序列可以被视为乘法时间序列。乘法时间序列的数学函数可以表示为： \[y(t) = 水平 * 趋势 * 季节性 * 噪声\]
加性
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)fig.savefig('../../figs/prediction_intervals_in_forecasting_models__train_test.png')
使用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| ds | y | unique_id | |
|---|---|---|---|
| 1 | 2014-07-01 00:30:00 | -2717.0 | 1 |
| 2 | 2014-07-01 01:00:00 | -1917.0 | 1 |
| 3 | 2014-07-01 01:30:00 | -1554.0 | 1 |
| 4 | 2014-07-01 02:00:00 | -836.0 | 1 |
| 5 | 2014-07-01 02:30:00 | -947.0 | 1 |
| ... | ... | ... | ... |
| 10315 | 2015-01-31 21:30:00 | 951.0 | 1 |
| 10316 | 2015-01-31 22:00:00 | 1051.0 | 1 |
| 10317 | 2015-01-31 22:30:00 | 1588.0 | 1 |
| 10318 | 2015-01-31 23:00:00 | -718.0 | 1 |
| 10319 | 2015-01-31 23:30:00 | -303.0 | 1 |
10319 rows × 3 columns
这将从每个值中减去滞后1，我们可以看看我们的系列现在的样子。
fig = plot_series(prep)fig.savefig('../../figs/prediction_intervals_in_forecasting_models__plot_values.png')
添加特征
延迟
看起来季节性消失了，我们现在可以尝试添加一些延迟特征。
mlf = MLForecast(models=model1,
freq='30min',
lags=[1,24],
target_transforms=[Differences([1])],
) prep = mlf.preprocess(df)
prep| ds | y | unique_id | lag1 | lag24 | |
|---|---|---|---|---|---|
| 25 | 2014-07-01 12:30:00 | -22.0 | 1 | 445.0 | -2717.0 |
| 26 | 2014-07-01 13:00:00 | -708.0 | 1 | -22.0 | -1917.0 |
| 27 | 2014-07-01 13:30:00 | 1281.0 | 1 | -708.0 | -1554.0 |
| 28 | 2014-07-01 14:00:00 | 87.0 | 1 | 1281.0 | -836.0 |
| 29 | 2014-07-01 14:30:00 | 1045.0 | 1 | 87.0 | -947.0 |
| ... | ... | ... | ... | ... | ... |
| 10315 | 2015-01-31 21:30:00 | 951.0 | 1 | 428.0 | 4642.0 |
| 10316 | 2015-01-31 22:00:00 | 1051.0 | 1 | 951.0 | -519.0 |
| 10317 | 2015-01-31 22:30:00 | 1588.0 | 1 | 1051.0 | 2411.0 |
| 10318 | 2015-01-31 23:00:00 | -718.0 | 1 | 1588.0 | 214.0 |
| 10319 | 2015-01-31 23:30:00 | -303.0 | 1 | -718.0 | 2595.0 |
10295 rows × 5 columns
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| ds | y | unique_id | lag1 | lag24 | expanding_mean_lag1 | rolling_mean_lag24_window_size7 | |
|---|---|---|---|---|---|---|---|
| 31 | 2014-07-01 15:30:00 | -836.0 | 1 | -1211.0 | -305.0 | 284.533325 | -1254.285767 |
| 32 | 2014-07-01 16:00:00 | -2316.0 | 1 | -836.0 | 157.0 | 248.387100 | -843.714294 |
| 33 | 2014-07-01 16:30:00 | -1215.0 | 1 | -2316.0 | -63.0 | 168.250000 | -578.857117 |
| 34 | 2014-07-01 17:00:00 | 2190.0 | 1 | -1215.0 | 357.0 | 126.333336 | -305.857147 |
| 35 | 2014-07-01 17:30:00 | 2322.0 | 1 | 2190.0 | 1849.0 | 187.029419 | 77.714287 |
| ... | ... | ... | ... | ... | ... | ... | ... |
| 10315 | 2015-01-31 21:30:00 | 951.0 | 1 | 428.0 | 4642.0 | 1.248303 | 2064.285645 |
| 10316 | 2015-01-31 22:00:00 | 1051.0 | 1 | 951.0 | -519.0 | 1.340378 | 1873.428589 |
| 10317 | 2015-01-31 22:30:00 | 1588.0 | 1 | 1051.0 | 2411.0 | 1.442129 | 2179.000000 |
| 10318 | 2015-01-31 23:00:00 | -718.0 | 1 | 1588.0 | 214.0 | 1.595910 | 1888.714233 |
| 10319 | 2015-01-31 23:30:00 | -303.0 | 1 | -718.0 | 2595.0 | 1.526168 | 2071.714355 |
10289 rows × 7 columns
您可以看到这两种方法得出了相同的结果，您可以使用您觉得最舒服的方式。
日期特征
如果你的时间列由时间戳组成，那么提取诸如周、星期几、季度等特征可能是有意义的。你可以通过传递一个包含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"]) # 季节性数据prep = mlf.preprocess(df)
prep| ds | y | unique_id | lag1 | lag24 | expanding_mean_lag1 | rolling_mean_lag24_window_size7 | year | month | day | hour | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 31 | 2014-07-01 15:30:00 | -836.0 | 1 | -1211.0 | -305.0 | 284.533325 | -1254.285767 | 2014 | 7 | 1 | 15 |
| 32 | 2014-07-01 16:00:00 | -2316.0 | 1 | -836.0 | 157.0 | 248.387100 | -843.714294 | 2014 | 7 | 1 | 16 |
| 33 | 2014-07-01 16:30:00 | -1215.0 | 1 | -2316.0 | -63.0 | 168.250000 | -578.857117 | 2014 | 7 | 1 | 16 |
| 34 | 2014-07-01 17:00:00 | 2190.0 | 1 | -1215.0 | 357.0 | 126.333336 | -305.857147 | 2014 | 7 | 1 | 17 |
| 35 | 2014-07-01 17:30:00 | 2322.0 | 1 | 2190.0 | 1849.0 | 187.029419 | 77.714287 | 2014 | 7 | 1 | 17 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 10315 | 2015-01-31 21:30:00 | 951.0 | 1 | 428.0 | 4642.0 | 1.248303 | 2064.285645 | 2015 | 1 | 31 | 21 |
| 10316 | 2015-01-31 22:00:00 | 1051.0 | 1 | 951.0 | -519.0 | 1.340378 | 1873.428589 | 2015 | 1 | 31 | 22 |
| 10317 | 2015-01-31 22:30:00 | 1588.0 | 1 | 1051.0 | 2411.0 | 1.442129 | 2179.000000 | 2015 | 1 | 31 | 22 |
| 10318 | 2015-01-31 23:00:00 | -718.0 | 1 | 1588.0 | 214.0 | 1.595910 | 1888.714233 | 2015 | 1 | 31 | 23 |
| 10319 | 2015-01-31 23:30:00 | -303.0 | 1 | -718.0 | 2595.0 | 1.526168 | 2071.714355 | 2015 | 1 | 31 | 23 |
10289 rows × 11 columns
拟合模型
# 拟合模型
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| ds | y | XGBRegressor | |
|---|---|---|---|
| unique_id | |||
| 1 | 2014-07-01 15:30:00 | 18544.0 | 18441.443359 |
| 1 | 2014-07-01 16:00:00 | 16228.0 | 16391.152344 |
| 1 | 2014-07-01 16:30:00 | 15013.0 | 15260.714844 |
| 1 | 2014-07-01 17:00:00 | 17203.0 | 17066.148438 |
| 1 | 2014-07-01 17:30:00 | 19525.0 | 19714.404297 |
| ... | ... | ... | ... |
| 1 | 2015-01-31 21:30:00 | 24670.0 | 24488.646484 |
| 1 | 2015-01-31 22:00:00 | 25721.0 | 25868.865234 |
| 1 | 2015-01-31 22:30:00 | 27309.0 | 27290.125000 |
| 1 | 2015-01-31 23:00:00 | 26591.0 | 27123.226562 |
| 1 | 2015-01-31 23:30:00 | 26288.0 | 26241.205078 |
10289 rows × 3 columns
from statsmodels.stats.diagnostic import normal_ad
from scipy import statssw_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图（概率-概率图）将数据样本与正态分布图进行比较，如果符合正态分布，则数据点将形成一条直线。
这三种正态性测试通过p值来确定数据样本来自正态分布总体的可能性。每个测试的原假设是“样本来自正态分布人口”。这意味着，如果结果的p值低于选择的显著性水平（alpha）值，则拒绝原假设。因此，有证据表明数据来自非正态分布。对于本文，我们将使用0.01的显著性水平。
result=mlf.forecast_fitted_values()
fig, axs = plt.subplots(nrows=2, ncols=2)
# 情节[1,1]
result["XGBRegressor"].plot(ax=axs[0,0])
axs[0,0].set_title("Residuals model");
# 情节
axs[0,1].hist(result["XGBRegressor"], density=True,bins=50, alpha=0.5 )
axs[0,1].set_title("Density plot - Residual");
# 情节
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_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();
带预测区间的预测方法
要生成预测，可以使用预测方法。
forecast_df = mlf.predict(h=30, level=[80,95])
forecast_df.head()| unique_id | ds | XGBRegressor | XGBRegressor-lo-95 | XGBRegressor-lo-80 | XGBRegressor-hi-80 | XGBRegressor-hi-95 | |
|---|---|---|---|---|---|---|---|
| 0 | 1 | 2015-02-01 00:00:00 | 26320.298828 | 25559.884241 | 25680.228369 | 26960.369287 | 27080.713416 |
| 1 | 1 | 2015-02-01 00:30:00 | 26446.472656 | 24130.429614 | 25195.461621 | 27697.483691 | 28762.515698 |
| 2 | 1 | 2015-02-01 01:00:00 | 24909.970703 | 23094.950537 | 23579.583398 | 26240.358008 | 26724.990869 |
| 3 | 1 | 2015-02-01 01:30:00 | 24405.402344 | 21548.628296 | 22006.662598 | 26804.142090 | 27262.176392 |
| 4 | 1 | 2015-02-01 02:00:00 | 22292.390625 | 20666.736963 | 21130.215430 | 23454.565820 | 23918.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')
置信区间是一个值的范围，具有高概率包含变量的真实值。在机器学习时间序列模型中，置信区间用于估计预测的不确定性。
使用置信区间的主要好处之一是，它允许用户理解预测的准确性。例如，如果置信区间非常宽，这意味着预测的准确性较低。相反，如果置信区间非常窄，这意味着预测的准确性较高。
置信区间的另一个好处是，它帮助用户做出明智的决策。例如，如果一个预测在置信区间内，这意味着它很可能会成真。相反，如果一个预测在置信区间外，这意味着它不太可能成真。
总体而言，置信区间是机器学习时间序列模型的重要工具。它帮助用户理解预测的准确性并做出明智的决策。
参考文献
- Changquan Huang • Alla Petukhina. Springer系列（2022年）。使用Python进行应用时间序列分析与预测。
- Ivan Svetunkov. 增强动态自适应模型（ADAM）的预测与分析。
- James D. Hamilton. 时间序列分析，普林斯顿大学出版社，新泽西州普林斯顿，第一版，1994年。
- Nixtla的Mlforecast参数。
- Pandas可用频率。
- Rob J. Hyndman和George Athanasopoulos（2018年）。“预测原则与实践，时间序列交叉验证”。。
- 季节周期- Rob J Hyndman。
Give us a ⭐ on Github