时间序列预测 — — ARIMA模型（理论分析与代码详解）

时间序列预测 — — ARIMA模型

一、理论分析

ARIMA模型的全称为：差分整合移动平均自回归模型

时间序列是一种数据类型，它记录了在连续时间点上观测到的数值。这些数值可以是任何可以量化的度量，比如经济指标、股票价格、温度、销售额等。

时间序列数据是一种按照时间顺序排列的数据集，其中每个数据点都与一个特定的时间点相关联。这种数据类型通常用于记录和分析随时间变化的数值信息。

时间序列根据时间和数值的性质不同，可以分为时期时间序列和时点时间序列。

• 时期时间序列：数值要素反映现象在一定时期内发展的结果。
• 时点时间序列：数值要素反应现象在一定时点上的瞬时水平。

性质：时期序列可加，时点序列不可加

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

平稳性

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

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

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

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

差分法

时间序列在$t$和$t-1$时刻的差值，将非平稳序列变平稳。
$$\Delta y(x)=y(x+1)-y(x),(x=0,1,2,\cdots )$$
ARIMA模型可以分为三部分：AR表示自回归模型，MA表示移动平均模型，I表示差分法。

模型详解

一、自回归模型：AR( p )

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

$p$阶自回归过程的公式：
$$y_t=\mu +\sum _{i=1}^p{\gamma }_i{y}_{t-i}+{ϵ}_t$$
$p$表示使用几期的历史值来预测，其中$y_t$是当前值，$\mu$是常数项，p是阶数，${\gamma }_i$是自相关系数

二、移动平均模型：MA( q )

移动平均模型关注的是自回归模型中误差项的累计

$q$阶自回归过程的公式定义：
$$y_t=\mu +\sum _{i=1}^q{\theta }_t{ϵ}_{t-i}+{ϵ}_t$$

公式表明：时间序列的当前值与历史值没有关系，而只依赖于历史白噪声的线性组合，移动平均法能有效地消除预测中的随机波动

三、自回归移动平均模型：ARMA(p, q)

自回归移动平均模型是自回归模型与移动平均模型的结合，公式：
$$y_t=\mu +\sum _{i=1}^p{\gamma }_i{y}_{t-1}+\sum _{i-1}^q{\theta }_i{ϵ}_{t-1}+{ϵ}_t$$
公式含义：

1. 一个随机时间序列可以通过一个自回归移动平均模型来表示，即该序列可以由其自身的过去或滞后值以及随机扰动项来解释。
2. 如果该序列是平稳的，即它的行为不会随着时间的推移而变化，那么我们就可以通过该序列过去的行为来预测未来。

四、差分自回归移动平均模型：ARIMA(p, d, q)

差分自回归移动平均模型是将自回归模型AR，移动平均模型MA，和差分法进行结合，这样就得到了差分自回归移动平均模型ARIMA(p, d, q)

其中$p$是自回归项数，$q$是移动平均项数，$d$为时间序列成为平稳时所做的差分次数。

ARIMA模型的建模步骤

1. 对序列绘图进行平稳性检验，观察序列是否为平稳，对于非平稳时间序列要先进行$d$阶差分，转化为平稳序列。
2. 经过第一步处理，已经得到平稳时间序列。要对平稳时间序列分别求得其自身相关系数ACF，和偏自相关系数PACF，通过对自相关图和偏自相关图的分析或通过 AIC/BIC搜索，得到最佳的阶数$p$、$q$
3. 由以上得到$p$、$q$、$d$的值，得到ARIMA模型，然后开始对得到的模型进行模型检验。

ACF与PACF

ACF（自相关函数Autocorrelation Function）

自相关函数测量的是时间序列在不同时间滞后下的值之间的相关性，对于时间序列$y_t$，$y_t$与$y_{t-k}$的相关系数称为$y_t$间隔k的自相关系数，公式：
$$ACF(k)={\rho }_k=\frac{Cov(y_t,y_{t-k})}{Var(y_t)}$$

PACF（偏自相关函数，Partial Autocorrelation Function）

对于一个平稳AR§模型，求出滞后k自相关系数$\rho (k)$时，实际上得到的并不是$x(t)$与$x(t-k)$之间单纯的相关关系，还会受到中间k - 1个随机变量$x(t-1),x(t-2),\cdots ,x(t-k+1)$的影响，为了能单纯测度$x(t-k)$对$x(t)$的影响，引进偏自相关系数PACF的概念，PACF在滞后k时的值表示在已知时间序列在$t-1,t-2,\cdots ,t-k+1$的值的情况下，时间点t和t−k的值之间的相关性。

公式：
$$PACF(k)=\frac{Cov[(z_t-\overline{z_t}),(z_{t-k}-\overline{z_{t-k}})]}{\sqrt{Var(z_t-\overline{z_t})}\sqrt{Var(z_{t-k}-\overline{z_{t-k}})}}$$

ADF检验

ADF检验，全称为Augmented Dickey-Fuller检验，是一种用于测试时间序列数据是否具有单位根（unit root）的统计检验方法。单位根是非平稳时间序列的一个特征，它意味着时间序列的均值和方差会随着时间的推移而变化，从而使得时间序列的长期行为变得不可预测。

对于一个时间序列，如何确定它是否满足平稳性要求？

1. 通过图检验法（通过时间序列趋势图或者自相关函数图判断）。
2. ADF检验。

ADF大致的思想就是基于随机游走（不平稳的一个特殊序列），对其进行回归，如果发现p = 1，说明序列满足随机游走，就是非平稳的。

ADF的参数

ADF检验的结果共有五个参数：

• 第一个值：表示Test Statistic，即T检验，表T统计量，假设检验值。
• 第二个值：p-value，即$p$值，表T统计量对应的概率值。
• 第三和第四个值：Lags Used，表示延迟和测试次数。。
• 第五个参数${'1%'：xxx, '5%'：xxx, '10%'：xxx}：表示不同程度拒接原假设的统计值。

如何确定是否平稳？

• 1%，5%，10%不同程度拒绝原假设的统计值和ADF假设检验值比较，ADF假设检验值同时小于1%，5%，10%即说明非常好地拒绝该假设。
• p-value的值是否非常接近0。

拖尾和截尾

拖尾指序列以指数率单调递减或震荡衰减，而截尾指序列从某个时点变得非常小。

• 截尾
  1. 在最初d阶明显大于2倍标准差范围
  2. 之后几乎95%的（偏）自相关系数都落在了2倍标准差范围以内。
  3. 并且由零自相关系数衰减为在零附近小值波动的过程非常突然。
• 拖尾
  1. 如果有超过5%的样本（偏）自相关系数都落入2倍标准差范围之外。
  2. 或者是由显著非0的（偏）自相关系数衰减为小值波动的过程比较缓慢或非常连续。

确定p，q

	ACF拖尾	ACF截尾
PACF拖尾	ARMA(p, q)****	MA(q)
PACF截尾	AR§	序列本身不存在明显的自相关性，ARMA类模型可能不适用。

MA与AR模型，p， q的值都是通过截尾的阶数来判定的，MA相当于ARMA(0, q)， AR相当于ARMA(p, 0)

通过拖尾和截尾来定阶，具有很强的主观性，我们可以根据信息准则函数法，来确定模型阶数。

• AIC准则（全称：最小化信息量 [Akaike Information Criterion]），公式为：
  $$AIC=-z\ln (L)+2K$$
  其中L表示模型的极大似然函数，K表示模型参数个数。AIC准则存在一定的不足，当样本容量很大时，AIC准则模型不收敛于真实模型。

• BIC（全称：贝叶斯信息准则[Bayesian Information Criterion]），BIC弥补了AIC的不足，公式为：
  $$BIC=-2\ln (L)+L\ln (n)$$
  其中n表样本容量。显然，这两个评价（AIC与BIC）指标越小越好。我们通过网格搜索法来确定AIC、BIC的最优模型**(p, q)**。

模型检验

检验残差序列的随机性，即残差之间是独立的$e_i=y_i-\widehat{y_i}$，残差序列的随机性可以通过自相关函数法来检验，即做残差的自相关函数图。

二、代码实验

导入相关的库

import pandas as pd
import matplotlib.pyplot as plt

加载数据

ChinaBank = pd.read_csv('ChinaBank.csv',index_col = 'Date',parse_dates=['Date'])

ChinaBank.head()

提取Close列

ChinaBank.index = pd.to_datetime(ChinaBank.index)
sub = ChinaBank.loc['2014-01':'2014-06','Close']

sub.head()

划分数据集

train = sub.loc['2014-01':'2014-03']
test = sub.loc['2014-04':'2014-06']

#查看训练集的时间序列与数据(只包含训练集)
plt.figure(figsize=(12,6))
plt.plot(train)
plt.xticks(rotation=45) #旋转45度
plt.show()

差分法

#.diff(1)做一个时间间隔
ChinaBank['diff_1'] = ChinaBank['Close'].diff(1) #1阶差分

#对一阶差分数据在划分时间间隔
ChinaBank['diff_2'] = ChinaBank['diff_1'].diff(1) #2阶差分

fig = plt.figure(figsize=(12,10))
#原数据
ax1 = fig.add_subplot(311)
ax1.plot(ChinaBank['Close'])
#1阶差分
ax2 = fig.add_subplot(312)
ax2.plot(ChinaBank['diff_1'])
#2阶差分
ax3 = fig.add_subplot(313)
ax3.plot(ChinaBank['diff_2'])
plt.show()

ADF检验

import matplotlib.pyplot as plt
import pandas as pd
import statsmodels.api as sm
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.stattools import adfuller as ADF

# 计算原始序列、一阶差分序列、二阶差分序列的单位根检验结果
ChinaBank['diff_1'] = ChinaBank['diff_1'].fillna(0)
ChinaBank['diff_2'] = ChinaBank['diff_2'].fillna(0)

timeseries_adf = ADF(ChinaBank['Close'].tolist())
timeseries_diff1_adf = ADF(ChinaBank['diff_1'].tolist())
timeseries_diff2_adf = ADF(ChinaBank['diff_2'].tolist())


# 打印单位根检验结果
print('timeseries_adf : ', timeseries_adf)
print('timeseries_diff1_adf : ', timeseries_diff1_adf)
print('timeseries_diff2_adf : ', timeseries_diff2_adf)

确定参数

import statsmodels.api as sm

#绘制
fig = plt.figure(figsize=(12,7))

ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(train, lags=20,ax=ax1)
ax1.xaxis.set_ticks_position('bottom') # 设置坐标轴上的数字显示的位置，top:显示在顶部  bottom:显示在底部
#fig.tight_layout()

ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(train, lags=20, ax=ax2)
ax2.xaxis.set_ticks_position('bottom')
#fig.tight_layout()
plt.show()

建立模型

#遍历，寻找适宜的参数
import itertools
import numpy as np
import seaborn as sns

#确定pq的取值范围
p_min = 0
d_min = 0
q_min = 0
p_max = 5
d_max = 0
q_max = 5

#Initialize a DataFrame to store the results,，以BIC准则
results_bic = pd.DataFrame(index=['AR{}'.format(i) for i in range(p_min,p_max+1)],
                           columns=['MA{}'.format(i) for i in range(q_min,q_max+1)])

for p,d,q in itertools.product(range(p_min,p_max+1),
                               range(d_min,d_max+1),
                               range(q_min,q_max+1)):
    if p==0 and d==0 and q==0:
        results_bic.loc['AR{}'.format(p), 'MA{}'.format(q)] = np.nan
        continue
    try:
        model = sm.tsa.ARIMA(train, order=(p, d, q),
                               #enforce_stationarity=False,
                               #enforce_invertibility=False,
                              )
        results = model.fit()
        results_bic.loc['AR{}'.format(p), 'MA{}'.format(q)] = results.bic
    except:
        continue

results_bic

#得到结果后进行浮点型转换
results_bic = results_bic[results_bic.columns].astype(float)

#绘制热力图
fig, ax = plt.subplots(figsize=(10, 8))
ax = sns.heatmap(results_bic,
                 mask=results_bic.isnull(),
                 ax=ax,
                 annot=True,
                 fmt='.2f',
                 cmap="Purples"
                 )

ax.set_title('BIC')
plt.show()

results_bic.stack().idxmin()

ADF检验取p、q的最优值

train_results = sm.tsa.arma_order_select_ic(train, ic=['aic', 'bic'], trend='n', max_ar=8, max_ma=8)

print('AIC', train_results.aic_min_order)
print('BIC', train_results.bic_min_order)

模型检验

#根据以上求得
p = 1
d = 0
q = 0

model = sm.tsa.ARIMA(train, order=(p,d,q))
results = model.fit()
resid = results.resid #获取残差

#绘制
#查看测试集的时间序列与数据(只包含测试集)
fig, ax = plt.subplots(figsize=(12, 5))

ax = sm.graphics.tsa.plot_acf(resid, lags=40,ax=ax)

plt.show()

模型预测

predict_sunspots = results.predict(dynamic=False)
print(predict_sunspots)

#查看测试集的时间序列与数据(只包含测试集)
plt.figure(figsize=(12,6))
plt.plot(train)
plt.xticks(rotation=45) #旋转45度
plt.plot(predict_sunspots)
plt.show()

#绘图
fig, ax = plt.subplots(figsize=(12, 6))
ax = sub.plot(ax=ax)
#预测数据
predict_sunspots.plot(ax=ax)
plt.show()