时间序列 ARMA 模型实战！

ARMA全称Autoregressive moving average model（自回归滑动平均模型），由美国统计学家博克斯（G.E.P.Box）和英国统计学家詹金斯（G.M.Jenkins）在二十世纪七十年代提出，也称B-J方法。

ARMA模型有三种基本形式：

1. 自回归模型（AR，Auto-regressive）
   自回归模型根据历史观测值进行预测。
   
   其中，是预测误差，为回归系数。

2. 移动平均模型（MA，Moving Average）
   移动平均模型根据历史预测误差进行预测。
   
   其中，为前q期的随机扰动项，误差项是当前期的随机干扰（零均值白噪声序列），为平滑系数。

3. 混合模型（ARMA，Auto-regressive Moving Average）
   混合模型为同时包含AR模型和MA模型。

自回归不是多元线性回归

小标题中的多元线性回归这里特指平时常说的带有多个自变量的非时间序列多元线性回归模型。

回归模型：
其中，为因变量，为自变量，为待估计的参数，为误差项。

自回归模型：
其中，为序列本身前个观测值。

回归是一种分析自变量和因变量关系的模型，当因变量和自变量为线性关系时叫线性回归模型，当有多个变量时叫多元线性回归模型。自回归（AR）模型是一种特殊的多元线性回归模型，和平时用的多元线性回归模型略有不同。

多元线性回归模型有多个自变量，但时间序列因其特殊性，一个时刻只有一个观测值，在无其他维度数据补充的前提下，没有可辅助预测的自变量，因而只能取其自身延迟p个时刻的观测值作为自变量，所以称自回归。自回归使用自身的数据进行预测，序列前后应当有相关性。所以多元线性回归模型预测时并不要求序列平稳，但AR模型则要求序列必须是平稳的。

AR模型对平稳性要求体现在如下几个方面：

下面是AR(2)模型对不同取值画出ACF的典型图形，供参考：

偏自相关函数p阶截尾表明，变量只与前p个变量有关，因而也可以借助PACF确定AR模型的阶数p。

比如可以画出PACF图的上下界限，以此判断PACF在哪里截尾。
上图中PACF虽然在k=3,9,14,16等位置超出界限，但是超出不多，可考虑用 AR(3) 建模。

移动平均不是移动平滑

在之前的文章中我们已经讲过了简单移动平均（SMA）、加权移动平均（WMA）、指数移动平均（EMA），他们和指数平滑一样都可以认为是一种平滑方法，利用修匀技术消弱短期随机波动对序列的影响使序列平滑化。

而ARMA中的MA模型，关注的是自回归模型中误差项的累加，基于移动平均的特性，能有效地消除预测中的随机波动。MA模型是用于预测未来值的方法，而移动平均平滑法是用来估计历史值的循环趋势。

MA(q)模型结构如下：

其中，为前q期的随机扰动项，误差项是当前期的随机干扰（零均值白噪声序列），为序列均值，为平滑系数且为待优化的参数。

残差序列怎么来的？如下公式递归计算出来的。

其中，为序列均值，需要估算。

贴个参考链接：https://stats.stackexchange.com/questions/26024/moving-average-model-error-terms

MA模型认为主要是受过去q期的误差项的影响。对于高阶的AR模型，有些可以用低阶的MA模型更好地描述。一般的AR模型也可以用高阶MA模型近似。MA模型同样要求序列宽平稳。

MA模型的自相关函数具有截尾性，即自相关函数在k>q后为零。

示例：
上图ACF在k=1很大，在k=3和k=9也比较明显，可以考虑拟合 MA(3) 或 MA(9)。

回归配平均效果翻倍

AR模型有偏自相关函数截尾性质，MA模型有相关函数截尾性质。ARMA模型结合了AR和MA模型，在对数据拟合优度相近的情况下往往可以得到更好的模型，而且不要求偏自相关函数截尾也不要求相关函数截尾。

ARMA模型的公式也很直接，AR模型+MA模型：

以上可以看出的取值是前p期和前q期的多元线性函数。

q=0时为AR(p)模型可以看成ARMA(p, 0)，p=0时为MA(q)模型可以看成是ARMA(0, q)。

众里寻参千百度

实际用ARMA建模时，阶数p、q都是未知的，确定p、q的问题称为定阶，一般画自相关偏自相关图或者根据最优信息准则定阶。

上文中已经给出了第一种判断方式，即根据ACF和PACF图确定。但是有时肉眼判断会产生困扰，看不清晰多少阶截尾且需要主观判断。此时就可以根据最优信息准则选择p、q，筛选策略为尽可能选择更简单的模型。常用最优信息准则有AIC、BIC。

AIC（Akaike information criterion，赤池信息量准则），由日本统计学家Akaike于1973年提出，综合考虑拟合精度和模型参数个数，值越小越好。



其中，k表示模型参数的个数，L表示似然函数。

k好说为参数个数为p+q+1；L似然函数怎么理解，在机器学习中经常用到这个概念，希望预测结果和实际值越接近越好（长得越像越好），当前场景下可以简记为负的残差平方和，L越大拟合越精确。

网上看到的一段公式：

BIC（Bayesian information criterion，Bayesian information criterion）是由Schwarz在1978年提出，公式和AIC基本一致，仅在惩罚项中增加了样本数量的考量。

为何？因为AIC虽好却有不足，就是样本量大的场景下，AIC准则通常愿意选择更复杂的模型，即便实际模型复杂度很低。原因体现在AIC的惩罚项为参数个数和样本容量没关系，但是似然函数却会受到样本容量的放大。



BIC考虑了样本个数，惩罚项比AIC更大，样本数量过多时，可以防止模型精度过高造成的模型复杂度过高。

系数内自省也

ARMA中的系数，模型训练时会自动拟合最优参数，平时用的时候不需要关心，但还是应该知道大体是怎么来的。

AR(p)模型的参数估计：

• Yule-Walker估计：用Levinson递推公式解Y-W方程，能够保证最小相位性，计算简单。
• 最小二乘估计：最小化残差平方和，根据线性模型理论解正规方程，计算简单。
• 最大似然估计：联合密度函数最大化，一般精度较高。

MA(q)模型的参数估计：

• 矩估计：利用参数与自协方差函数关系，用解非线性方程组的方法求解，误差较大。
• 逆相关函数法：把MA变成一个AR再用Yule-Warker方法估计参数。
• 新息方法：计算样本自协方差函数，根据新息预报递推公式计算参数。

ARMA(p,q)模型的参数估计

• 矩估计方法：用Yule-Walker方程估计AR模型参数，然后可以用逆相关函数法估计MA模型参数。在中小样本， 矩估计不是一种实用的方法， 即使作为最大似然估计的初值也不太合适。
• 自回归逼近法：可以看成是回归系数，先拟合一个近似自回归模型。然后算出残差序列，对残差平方和极小化最终得到估计参数。
• 最大似然估计：也可以说是最小平方和估计，平方和指的是误差项。用最大似然估计可以大大改善MA部分的参数估计精度。

以上内容参考自北大李东风的应用时间序列分析备课笔记：
https://www.math.pku.edu.cn/teachers/lidf/course/atsa/atsanotes/html/_atsanotes/atsa-estarma.html

短期预言家

ARMA模型适合短期预测，对于长期预测，预测结果偏向于均值，预测均方误差趋于序列的方差。对于AR模型来说，预测得越远，均方误差越大；对于MA模型来说，到q步之后就没了预测能力。

具体ARMA模型的预测公式是怎样的，我们不妨一起看看AR模型和MA模型的预测。

AR模型的一步预测：
其中，为h+1时刻的误差项未知。
所以超前一步的预测公式为：
预测误差为：

AR模型的多步预测：

其中，

所以计算需要递推计算。
超前k步的预测误差为：
当超前步数k趋近无穷大时，趋近于序列均值，称为均值回归，和前一步的越强（越大），多步预报的效果也越差。

MA模型的一步预测：
根据，

MA模型的多步预测：

其中，，即MA(2)中

因为MA(q)序列在间隔超过q步以后就独立，所以超前多步预测，只能预测到步，从q+1步开始就只能用均值预测了。

ARMA模型的一步预测：

ARMA模型的多步预测： 

其中，
；

输入必要平稳

序列平稳是ARMA分析的前提，所以输入序列必须是平稳的。时间序列的平稳性检验主要有三种方法：

1. 图形分析方法：绘制时间序列的折线图，看曲线是否稳定；看自相关图中自相关系数是否能很快退化到零。缺点是需要主观判断，没有量化依据。
2. 简单统计方法：根据宽平稳均值、方差不随时间变化的特性，将序列前后拆分成两段，分别计算均值、方差，对比看是否差异明显。
3. 假设检验方法：当前主流单位根检验，检验序列中是否存在单位根，若存在，则为非平稳序列，不存在则为平稳序列。检验方法有ADF检验、PP检验、KPSS检验等。

具体实现可参考“时间序列的平稳性检验方法汇总篇”。

然事常与愿违，现实中的时间序列数据多是非平稳的，检验后发现非平稳怎么办，难道就不能使用ARMA了吗？并不，事总在人为，可以先对序列进行平稳化处理。

常用平稳化的方法有：

1. 差分：差分可以去除序列中的趋势和季节性。一阶差分可以去除线性趋势，如果还有二次趋势，还可以继续二阶差分。二阶差分后还未平稳的话就要注意了，继续差分即便最终平稳了，但是多次差分后解释力下降，且有可能造成过度差分，最差为差分后的序列为白噪声，后面也没法分析了。对于周期型序列也可以用季节差分的方式去除时间序列季节性。
2. 平滑：对当前序列值减去平滑值得到一个残差序列，当平滑结果能比较好的描述原始序列趋势特征的时候，残差序列一般是平稳的，后续可对残差序列进行建模预测。计算平滑值的方法可以用简单移动平均、加权移动平均、一次指数平滑、二次指数平滑等。同类思想，还可以拟合一个回归方程，用回归方程描述原始序列的趋势特征。
3. 变换：如对数变换，能够去除方差随时间增长的趋势。对数据进行取log处理，变换前的序列必须满足大于0（小于0的话可以对序列用+x的方式变为正数）。取对数后，原数据越大，缩小的幅度越大，可以使得方差随时间波动大的时间序列的方差变得更稳定，从而一定程度上使得序列平稳。但也不一定变换后即平稳，比如呈指数趋势的序列，变换后只能将指数趋势转化为线性趋势，此时再使用一阶差分即可将序列变得平稳，同时变换后的数据可以看成增长率的对数，解释性强。其它还有开根号、Box-Cox变换、Yeo-Johonson变换等。这些变换试图将数据转换为正态分布，虽然对于平稳性来说并不总是必要的，但通常能够纠正序列的非线性问题。
4. 分解：可以将时间序列分解成3部分：长期趋势、季节变动、不规则波动，3种成分相加叫加法模型，3种成分相乘叫乘法模型，1加1乘叫混合模型。分解目的为去除季节性的影响，分解后可对分解出的趋势项、季节项和余项分别进行预测。常用时间序列分解方法有朴素分解、X11分解、SEATS分解、STL分解等，其中STL分解用的较多。

平稳化之后也莫急着训练模型，还需检验序列是否白噪声。因为白噪声是平稳的，但是却是完全随机的，没法预测。所以如果序列是白噪声的话，也就可以回家吃饭，莫要费功夫钻研了。

输出还可再检查

如果模型拟合充分，则残差序列服从零均值正态分布且完全随机无相关性。所以还需要分析残差的正态性和无关性，以此检验模型拟合是否良好。如果检验未通过，说明模型不合适或者参数没调好。

检验残差的正态性，可以查看残差的正态概率图、QQ图，也可以使用normaltest正态性检验等方法。

• 分布直方图：以比较直观的方式看数据分布是否服从正态分布。
• QQ图：Quantile-Quantile Plot，即分位数-分位数图，QQ图中若残差基本完全落在45°线上即为符合正态性假设。
• normaltest：scipy中normaltest根据数据的峰度偏度检验是否符合正态分布，p值大于显著性水平0.05认为样本数据符合正态分布，检验较严格。

检验残差的无关性（白噪声检验），可以看残差的自相关图，也可以使用假设检验方法D-W检验、Box-Ljung检验、Ljung-Box检验。

• 自相关图：白噪声完全无自相关性，除0阶自相关系数为1外，延迟k阶的样本自相关系数均为0或在0值附近。
• D-W检验：杜宾-瓦特森检验，检验序列是否存在一阶自相关。DW值接近2时不存在一阶自相关性。
• Box-Ljung检验：原假设为滞后m阶序列值之间相互独立，即序列为独立同分布的白噪声，如果p值大于显著性水平如0.05，不能拒绝原假设，序列为白噪声序列。
• Ljung-Box检验：Box-Ljung检验在小样本量下不太精确，Ljung-Box检验弥补了这一缺陷。同Box-Ljung检验，p值大于显著性水平如0.05，不能拒绝原假设，序列为白噪声序列。

详细白噪声检验方法可以看下“白噪声检验”。

放码过来

读取数据

使用statsmodels中自带的美国夏威夷莫纳罗亚天文台大气二氧化碳数据，并转换为月度数据。

import numpy as np
import pandas as pd
from statsmodels.datasets import co2
from statsmodels.tsa.stattools import adfuller
from matplotlib import pyplot as plt

# 加载数据
data = co2.load(as_pandas=True).data
# 月度数据
data = data['co2']
data = data.resample('M').mean().ffill()

# 可视化
data.plot(figsize=(12,4))
plt.title('co2 data')
plt.show()

平稳性检验

能够很直观的看出来，存在趋势和周期，必然非平稳，ADF检验p值应该很大。

res = adfuller(data)
print('p value:', res[1])

p value: 0.9989453312516823
p值显著大于0.05，序列非平稳。

平稳化

先来使用差分的方式进行平稳化处理，因为有周期直接使用季节差分。

# # 一阶差分
# data_diff1 = data.diff()
# # 二阶差分
# data_diff2 = data_diff1.diff()

# 季节差分
data_diff = data.diff(12).dropna()

data_diff.plot(figsize=(12,4))
plt.title('co2 - seasonal difference')
plt.show()

# ADF检验
res = adfuller(data_diff)
print('p value:', res[1])

季节差分后，ADF检验中p值为0.0007小于0.05，故而差分后序列平稳。

白噪声检验

检验差分后的数据是否为白噪声。

from statsmodels.stats.diagnostic import acorr_ljungbox
res = acorr_ljungbox(data, lags=[6,12,24], return_df=True)
print(res)

各滞后期数下p值都为0均小于0.05，差分后的序列非白噪声，可以进行预测。

模型定阶

根据最小化BIC准则确定p、q。计算比较耗时，为了控制计算量，这里限制AR最大阶不超过6，MA最大阶不超过4。

from statsmodels.tsa.stattools import arma_order_select_ic
bic_min_order = arma_order_select_ic(data_diff, max_ar=6, max_ma=4, ic='bic')['bic_min_order']
print(bic_min_order)

求得ARMA的阶数为(6,3)。

模型训练

from statsmodels.tsa.arima_model import ARMA
model = ARMA(data_diff, order=bic_min_order).fit(disp=-1)
print(model.summary())

系数p值均小于0.05，说明系数均显著不为0。如果存在某个系数的p值较大，说明t检验中系数不显著，贡献不大，可以剔除，减小模型冗余度。

模型检验

可视化残差直方图、QQ图，看是否正态分布，可视化ACF看是否仍存有自相关性。

import statsmodels.api as sm
from statsmodels.graphics.tsaplots import plot_acf

fig, axs = plt.subplots(2, 2)
fig.subplots_adjust(hspace=0.3)

model.resid.plot(ax=axs[0][0])
axs[0][0].set_title('residual')

model.resid.plot(kind='hist', ax=axs[0][1])
axs[0][1].set_title('histogram')

sm.qqplot(model.resid, line='45', fit=True, ax=axs[1][0])
axs[1][0].set_title('Normal Q-Q')

plot_acf(model.resid, ax=axs[1][1])

plt.show()

上图可以看出残差基本服从正态分布，但ACF图显示滞后12期时仍存在自相关性，所以模型还可优化。猜想是否有可能模型定阶时限制了AR和MA的最大阶数，所以直接根据差分后样本数据的ACF图和PACF图定阶。

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

fig, axs = plt.subplots(2, 1)
plot_acf(data_diff, ax=axs[0])
plot_pacf(data_diff, ax=axs[1])
plt.show()

根据以上ACF、PACF图，取p=2,q=11训练模型。

model2 = ARMA(data_diff, order=(2,11)).fit(disp=-1)

fig, axs = plt.subplots(2, 2)
fig.subplots_adjust(hspace=0.3)

model.resid.plot(ax=axs[0][0])
axs[0][0].set_title('residual')

model.resid.plot(kind='hist', ax=axs[0][1])
axs[0][1].set_title('histogram')

sm.qqplot(model.resid, line='45', fit=True, ax=axs[1][0])
axs[1][0].set_title('Normal Q-Q')

plot_acf(model.resid, ax=axs[1][1])

plt.show()

效果有所改善，符合预期。

模型预测

预测未来6期数据（注意ARMA模型预测时，预测期数超过MA模型最大阶数后，MA部分的预测值实际为序列均值不再变化，AR部分倒还会有持续变化但效果下降）。

preds = model.predict(0, len(data_diff)+6)

# 也可只取未来预测值
# fcast = model2.forecast(6)

plt.figure(figsize=(12, 4))
data_diff.plot(color='g', label='data_diff')
preds.plot(color='r', label='predict')
plt.legend()
plt.show()

差分还原

把曾经减去的项再加回来。

df1 = pd.DataFrame(data)
df2 = pd.DataFrame(preds, columns=['predict'])
df = pd.concat([df1, df2], axis=1)
df['result'] = df['predict'] + df['co2'].shift(12)

plt.figure(figsize=(12, 4))
df['result'].plot(color='r', label='arma result')
df['co2'].plot(color='g', label='co2 data')
plt.legend()
plt.show()