ARIMA时间序列（自学使用，侵权即删）

引子：什么是时间序列？

在我们学习如何处理时间序列数据之前，首先应理解什么是时间序列，以及它与其他类型的数据的区别。时间序列的正式定义如下：它是一系列在相同时间间隔内测量到的数据点。

简言之，时间序列是指以固定的时间间隔记录下的特定的值，时间间隔可以是小时、每天、每周、每10天等等。时间序列的特殊性是：该序列中的每个数据点都与先前的数据点相关。

        EX1:假设你从某公司获得了一个贷款人员的数据集(如下表所示)。你认为每一行都与前面的行相关吗？当然不是！一个人的贷款金额取决于他的经济状况和需要(可能还有其他因素，如家庭规模等，但为了简单起见，我们只考虑收入和贷款类型)。此外，这些数据不是在特定时间间隔内收集的，它仅与公司何时收到贷款申请相关。

        EX2:假设你有一个数据集，其中包含每天空气中的二氧化碳水平(下面是截图)。那么可以通过过去几天的数值来预测第二天的二氧化碳水平吗？当然可以。如果你观察到的数据是每天记录下来的，那么，时间间隔便是恒定的(24小时)。

现在我们已经对时间序列的定义有了比较清楚的认知，下面我们来介绍一种流行且广泛使用的时间序列预测统计方法：ARIMA模型。

ARIMA时间序列

 ARIMA 是代表autoRegressive I integrated Moving a average，自回归综合移动平均线的首字母缩写词，它是一类在时间序列数据中捕获一组不同标准时间结构的模型。预测方程中平稳序列的滞后称为“自回归”项，预测误差的滞后称为“移动平均”项，需要差分才能使其平稳的时间序列被称为平稳序列的“综合”版本。

随机游走和随机趋势模型、自回归模型和指数平滑模型都是 ARIMA 模型的特例。

ARIMA 模型可以被视为一个“过滤器”，它试图将信号与噪声分开，然后将信号外推到未来以获得预测。ARIMA模型特别适合于拟合显示非平稳性的数据。

为了能够使用ARIMA，我们还需要补充一些概念。

1.平稳性

时间序列模型的数据必须要满足平稳性！

平稳性就是要求经由样本时间序列所得到的拟合曲线在未来的一段时间内仍然能够按照现有的形态延续下去。

平稳性要求序列的期望和方差不发生明显变化，根据平稳性，可将数据分为严平稳和宽平稳

        严平稳：也称为强平稳或严格平稳，它要求时间序列的所有统计特性（如均值、方差、协方差等）在时间上完全不变。
        宽平稳：也称为弱平稳，它要求时间序列的一阶和二阶矩（即均值和方差）在时间上保持不变，但不需要协方差在时间上保持不变。

实际数据大致上都是宽平稳，如果一个时间序列不是平稳的，通常需要通过差分的方式将其转化为平稳时间序列。

2.差分

从统计学的角度来看，数据差分是指将非平稳数据转换为平稳的过程，去除其非恒定的趋势。“差分“消除了时间序列水平的变化，消除了趋势和季节性，从而稳定了时间序列的平均值。”  季节性差分应用于季节性时间序列以去除季节性成分。

3.ARIMA模型拆解

剖析ARIMA的各个部分，以便更好地理解它如何帮助我们时间序列建模，并对其进行预测。

1）AR - 自回归

自回归模型，顾名思义，就是及时地“回顾”过去，分析数据中先前的值，并对它们做出假设。这些先前的值称为“滞后”。一个例子是显示每月铅笔销售的数据。每个月的销售总额将被认为是数据集中的一个“进化变量”。这个模型是作为“利益的演化变量根据其自身的滞后值(即先验值)进行回归”而建立的。

自回归模型用于描述当前值和历史值之间的关系，用变量自身的历史数据对自身进行预测，数据必须要满足平稳性要求，只适用于预测与自身前期相关的现象（时间序列的自相关性）

2）I - 表示差分(数据的”综合“)

与“ARMA”模型的区别在于，ARIMA中的“I”指的是它的综合方面。当应用差分步骤时，数据是“综合”的，以消除非平稳性。表示原始观测值的差异，以允许时间序列变得平稳，即数据值被数据值和以前的值之间的差异替换。

3）MA - 移动平均线

移动平均模型关注的是自回归模型中误差项的累计，是将观测值与应用于滞后观测值的移动平均模型的残差之间的相关性合并。

ARIMA用于使模型尽可能地符合时间序列数据的特殊形式。

 

4.ARIMA模型建立

一般步骤

1. 加载数据：构建模型的第一步当然是加载数据集。
2. 预处理：根据数据集定义预处理步骤。包括创建时间戳、日期/时间列转换为d类型、序列单变量化等。
3. 序列平稳化：为了满足假设，应确保序列平稳。这包括检查序列的平稳性和执行所需的转换。
4. 确定d值：为了使序列平稳，执行差分操作的次数将确定为d值。
5. 创建ACF和PACF图：这是ARIMA实现中最重要的一步。用ACF PACF图来确定ARIMA模型的输入参数。
6. 确定p值和q值：从上一步的ACF和PACF图中读取p和q的值。
7. 拟合ARIMA模型：利用我们从前面步骤中计算出来的数据和参数值，拟合ARIMA模型。
8. 在验证集上进行预测：预测未来的值。
9. 计算RMSE：通过检查RMSE值来检查模型的性能，用验证集上的预测值和实际值检查RMSE值。

5.详细操作步骤

Step1：导入必要的库

import os 
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import datetime as dt
import math
 
from pandas.plotting import autocorrelation_plot
from statsmodels.tsa.statespace.sarimax import SARIMAX
from sklearn.preprocessing import MinMaxScaler
from common.utils import load_data, mape
from IPython.display import Image
 
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf  
from statsmodels.tsa.stattools import adfuller  # adf检验库
from statsmodels.stats.diagnostic import acorr_ljungbox  # 随机性检验库
from statsmodels.tsa.arima_model import ARMA 
 
%matplotlib inline
plt.rcParams['figure.figsize'] = (12,6)
pd.options.display.float_format = '{:,.2f}'.format
np.set_printoptions(precision=2)
warnings.filterwarnings("ignore") # 忽略一些可能的红字警告

Step2：稳定性检验

adfuller(Augmented Dickey-Fuller)测试可用于在存在串行相关的情况下在单变量过程中测试单位根。

statsmodels.tsa.stattools.adfuller(x，maxlag = None，regression ='c'，autolag ='AIC'，store = False，regresults = False ) 

adfuller中可进行adf校验，一般传入一个data就行，包括 list, numpy array 和 pandas series都可以作为输入，其他参数可以保留默认。

返回值：

adf（float）
测试统计

pvalue（float）
MacKinnon基于MacKinnon的近似p值（1994年，2010年）

usedlag（int）
使用的滞后数量

nobs（int）
用于ADF回归的观察数和临界值的计算

critical values（dict）
测试统计数据的临界值为1％，5％和10％。基于MacKinnon（2010）

icbest（float）
如果autolag不是None，则最大化信息标准。

resstore （ResultStore，可选）
一个虚拟类，其结果作为属性附加

如何确定该序列能否平稳呢？主要看：

1%、5%、10%不同程度拒绝原假设的统计值和ADF Test result的比较，ADF Test result同时小于1%、5%、10%即说明非常好地拒绝该假设。另外，P-value是否非常接近0，接近0，则是平稳的，否则，不平稳。

若不平稳，则需要进行差分，差分后再进行检测。

Step3：白噪声检测

白噪声检验也称为纯随机性检验，当数据是纯随机数据时，再对数据进行分析就没有任何意义了，所以拿到数据后最好对数据进行一个纯随机性检验。acorr_ljungbox(x, lags=None, boxpierce=False,         model_df=0, period=None, return_df=True, auto_lag=False)

主要参数

lags为延迟期数，如果为整数，则是包含在内的延迟期数，如果是一个列表或数组，那么所有时滞都包含在列表中最大的时滞中。

boxpierce为True时表示除开返回LB统计量还会返回Box和Pierce的Q统计量

返回值

lbvalue: (float or array)
测试的统计量

pvalue: (float or array)
基于卡方分布的p统计量

bpvalue: ((optionsal), float or array)
基于 Box-Pierce 的检验的p统计量

bppvalue: ((optional), float or array)
基于卡方分布下的Box-Pierce检验的p统计量

若p值远小于0.01，因此我们拒绝原假设，认为该时间序列是平稳的。（这里原假设是存在单位根，即时间序列为非平稳的。)

Step4：建立ARIMA模型

ACF 和 PACF 图： 通过差分对时间序列进行平稳化后，拟合 ARIMA 模型的下一步是确定是否需要 AR 或 MA 项来校正差分序列中剩余的任何自相关。结合自相关图和偏自相关图共同进行判断时间序列模型。

关于ARMA通用判断标准说明如下表格：

模型	自相关图	偏自相关图
AR(p)	拖尾	p阶截尾
MA(q)	q阶截尾	拖尾
ARMA(p,q)	拖尾	拖尾
模型不适合	截尾	截尾

拖尾和截尾说明如下：

拖尾： 始终有非零取值，不会在大于某阶后就快速趋近于0（而是在0附近波动），可简单理解为无论如何都不会为0，而是在某阶之后在0附近随机变化。

截尾： 在大于某阶(k)后快速趋于0为k阶截尾，可简单理解为从某阶之后直接就变为0。

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
fig, axes = plt.subplots(2, 1, figsize=(12, 4*2))
# 自相关
plot_acf(data['co2'],lags=12,title='raw_acf', ax=axes[0])  
# 偏自相关
plot_pacf(data['co2'],lags=12,title='raw_pacf', ax=axes[1])  
plt.show()
# 其中lags 表示滞后的阶数，以上分别得到acf 图和pacf 图。

Step5：选择超参数

为ARIMA模型的参数选择最佳值可能是一个挑战，因为这有点主观，也有点耗时。

可以考虑使用'pyramid'库中的 'auto_arima()' 函数。

本文提供一种网格搜索方法来自动选择超参数。将网格搜索定义为一个函数evaluate_arima_model()，该函数以时间序列数据集作为输入，以及元组(p,d,q)作为参数用于评估模型。

数据集分为两部分：初始训练数据集为 66％，测试数据集为剩余的 34％。

迭代测试集的每个时间步。一次迭代就可以训练一个模型，然后使用该模型对新数据进行预测。每次迭代都进行预测并存储在列表中。最后用测试集将所有预测值与预期值列表进行比较，并计算并返回均方误差分数。

# 为 ARIMA (p,d,q)模型选取合理的超参数

def evaluate_arima_model(X, arima_order):

         # 数据集划分

         train_size = int(len(X) * 0.66)

         train, test = X[0:train_size], X[train_size:]

         history = [x for x in train]

         # 预测

         predictions = list()

         for t in range(len(test)):

          model = ARIMA(history, order=arima_order)

          model_fit = model.fit(disp=0)

          yhat = model_fit.forecast()[0]

          predictions.append(yhat)

          history.append(test[t])

         # 计算样本误差

         error = mean_squared_error(test, predictions)

         return error

# evaluate combinations of p, d and q values for an ARIMA model
def evaluate_models(dataset, p_values, d_values, q_values):
  # 确保输入数据是浮点值（而不是整数或字符串）
         dataset = dataset.astype('float32')
         best_score, best_cfg = float("inf"), None

#暴力迭代搜寻
         for p in p_values:
                  for d in d_values:
                           for q in q_values:
    order = (p,d,q)
    try:
             mse = evaluate_arima_model(dataset, order)
             if mse < best_score:
                      best_score, best_cfg = mse, order
                     print('ARIMA%s MSE=%.3f' % (order,mse))
            except:
                     continue
 print('Best ARIMA%s MSE=%.3f' % (best_cfg, best_score))

Step6：滑动窗口进行模型检验

为了评估模型，可以使用walk forward验证。在实践中，每次有新的数据可用时，时间序列模型都要重新训练。这使得模型可以在每个时间步骤中做出最好的预测。

从使用该技术的时间序列的开始，在训练数据集上训练模型。然后对下一个时间步骤进行预测。根据已知值对预测进行评估。然后将训练集扩展到包含已知值，并重复该过程。

Note: 为了更有效的训练，应该保持训练集窗口固定，以便每次向训练集添加新的观察值时，并将该观察值从集合的开始处删除。

滑动窗口过程为模型在实践执行提供了更可靠的估计。然而，这是以创建众多模型的计算成本为代价的。如果数据较小或模型简单的话，这是可以接受的，但如果数据量大，或者模型规模大可能是一个问题。

training_window = 720 # 投入30天(720小时)进行训练

train_ts = train['load']

test_ts = test_shifted

history = [x for x in train_ts]

history = history[(-training_window):]

predictions = list()

order = (2, 1, 0)

seasonal_order = (1, 1, 0, 24)

for t in range(test_ts.shape[0]):

    model = SARIMAX(endog=history, order=order, seasonal_order=seasonal_order)

    model_fit = model.fit()

    yhat = model_fit.forecast(steps = HORIZON)

    predictions.append(yhat)

    obs = list(test_ts.iloc[t])

    # move the training window

    history.append(obs[0])

    history.pop(0)

    print(test_ts.index[t])

    print(t+1, ': predicted =', yhat, 'expected =', obs)

Walk-forward validation 是时间序列模型评估的黄金标准，可以考虑用于你自己的项目。

PS:

注意一点，ARIMA适用于短期 单变量预测，长期的预测值都会用均值填充，效果十分不好

原文链接：理论加实践，终于把时间序列预测ARIMA模型讲明白了-CSDN博客独家 | 利用Auto ARIMA构建高性能时间序列模型（附Python和R代码） - 知乎 (zhihu.com)python使用ARIMA进行时间序列的预测（基础教程）_python arima-CSDN博客