【时序预测】一文梳理时间序列预测——ARMA模型

【目录】

• 1. 时间序列背景
  • 1.1. 数据分类
  • 1.2. 统计学分析方法
  • 1.3. 整体框架
• 2. 数学基础
  • 2.1. 几个基本概念
  • 2.2. 自相关
  • 2.3. 偏自相关
  • 2.4. 随机性
  • 2.5. 平稳性
• 3. 时间序列分析建模
• 4. 平稳时间序列之ARMA建模
  • 4.1. 模型介绍
  • 4.2. 模型识别
  • 4.3. 参数估计
  • 4.4. 模型检验

1. 时间序列背景

1.1. 数据分类

在统计研究中，常用按时间顺序排列的一组随机变量X1​,X2​,⋯,Xt​,⋯来表示一个随机事件的时间序列，简记为{Xt,t∈T}。在时间的角度上来说，数据类型可分为两类：横截面数据和时间序列[1]。横截面数据是值在某一时间点搜集来自不同对象的数据，时间序列是一组按照时间排序的数据；横截面数据与时间序列的组合在计量经济学上构成了面板数据集。

时间序列按照时间排序的事实意味着他们存在特殊的性质以及存在一些特定的方式去分析时间序列。

1.2. 统计学分析方法

早期人在认识自然和改造自然的过程中，对于很多自然现象发现只要观察时间长，就能通过直观的图示和简单的数据汇总就能发现自然规律。早期的时序分析通常都是通过直观的数据比较或绘图观测,寻找序列中蕴含的波动特征和发展规律,这种分析方法就称为描述性时序分析。随着时间和社会的发展，想通过对序列（比如金融、人口、互联网等相关行业）简单的观察和描述,总结出随机变量发展变化的规律,并准确预测出它们将来的走势通常是非常困难的。因此，时间序列分析将研究的重心从总结表面现象转移到分析序列值内在的相关关系上。

时间序列分析的方法在统计学上可以分为频域分析法和时域分析法[2]。

• 频域分析法

主要思想是任何时间序列都可以分解成若干个不同周期（频率）的波动，各分量的周期变化可通过傅里叶变换利用正弦函数表示，从而无限逼近原始时间序列。借助Fourier分析或小波变换从频率的角度解释时间序列，需要较强的数学基础，且结果抽象难以解释。

• 时域分析法

主要思想是事件的发展通常有一定的惯性，用统计学语言来说就是序列值之间存在一定的相关关系，而这种相关关系具有某种统计规律性，可以利用拟合法构建数学模型来描述时间序列。该方法理论基础扎实，操作规范简单，分析结果容易解释。因此，时域分析法衍生出的经典时间序列分析方法是本文所讲的重点。目前主流的时间序列分析软件有SAS、Eviews、S-plus、Matlab、R、Python、TSP等。

1.3. 整体框架

2. 数学基础

2.1. 几个基本概念

1. 均值：表示时间序列在各个时刻取值的平均值，其定义如下：

1. 方差：表示时间序列在各个时刻围绕其均值波动的平均程度，其定义如下：

1. 协方差：表示时间序列任意两个时刻直接的相关性，任取t,s∈Tt , s \in Tt,s∈T，则其定义如下：

2.2. 自相关

理论上，基于真实数据的生成过程来计算自相关函数。在公式上对其定义如下：

然而，在估计时我们不知道真正的数据生成过程。事实上，我们有数据，也就是时间序列，我们相用ARMA(p, q)模型估计数据的生成过程。如下，用自相关函数的样本来计算样本自相关函数（样本ACF）

2.3. 偏自相关

自相关函数虽然反映了时间序列在两个不同时刻Xt和Xs的相依程度，但是这种相关板含量Xs通过Xt和Xs之间的其他变量Xt+1, Xt+2, ...Xt-s传递到Xt(s<t)的影响，也就是说自相关函数实际上掺杂了其他变量的影响。为了剔除中间变量的影响，可引入偏自相关函数PACF，换句话说，

PACF-消除AR(p)过程的间接相关性。其公式定义如下：

2.4. 随机性

• 随机性的意义

纯随机性的时间序列也称作白噪声序列。纯随机性的时间序列是指序列值之间没有相关关系，具有高度随机性即无记忆序列，历史数据或过去的行为对将来的发展没有任何影响。纯随机性序列不具有研究的价值。

纯随机性时间序列是独立同分布、均值为0、协方差为0，满足以下三个性质：

• 随机性的检验

1、图检验法

根据纯随机性序列的均值为0、协方差为0的性质，可以观察时序图是否满足均值为零、具有高度纯随机性的条件；观察ACF自相关图，如果序列不与历史时刻序列存在自相关关系，则可认为是随机性序列。这种方法具有较为粗糙，不太准确，常用作预检验。

比如下图，随机波动，不具有任何明显或不明显的历史相关性，故为纯随机性时间序列。

QQ图：若残差满足正态性假设，残差序列为正态白噪声过程或高斯白噪声过程。残差序列的Q-Q图近似为过原点的一条直线，则残差服从正态分布且均值为零。

Python语言
方法一：参考链接：https://zhuanlan.zhihu.com/p/77786430
import statsmodels.api as sm
sm.qqplot(data, line='s', ax=qq_ax)
方法二：# 参考链接：https://www.jianshu.com/p/6250e60fa28a
model = sm.tsa.statespace.SARIMAX()
mini_aic_model = model.fit()
mini_aic_model.plot_diagnostics(figsize=(16, 8))  # 可画出Q-Q图

2、统计检验法

随机性的检验方法是由数学家Box 和Ljung提出的Box-Ljung检验，观察QLB统计量是否服从卡方分布。当P值大于显著性水平α，可认为序列是白噪声序列。Box和Pierce提出的QBP统计量检验适用于大样本场合（n>=50），在小样本场合检验效果不太精确。基于QLB检验统计量的Box-Ljung检验适用性更好。

Q统计量和LB统计量，其计算公式分别如下：

其中，nnn为序列的观察期数，mmm为指定延迟期数，kkk为延迟阶数。

Python语言
方法一：参考链接：https://mlln.cn/2017/10/26/python%E6%97%B6%E9%97%B4%E5%BA%8F%E5%88%97%E5%88%86%E6%9E%90-%E7%BA%AF%E9%9A%8F%E6%9C%BA%E6%80%A7%E6%A3%80%E9%AA%8C/
from statsmodels.sandbox.stats.diagnostic import acorr_ljungbox
qljungbox, pval, qboxpierce, pvalbp=acorr_ljungbox(data, boxpierce=True)
 
方法二：版本最新
from statsmodels.stats.diagnostic import acorr_ljungbox

经验表明，纯随机性检验时滞后阶数不必取太大，因为平稳序列通常具有短期相关性。如果序列值之间存在显著的相关关系，通常指存在于延迟时期较短的序列值之间；若短期延迟都不存在显著的相关关系，那么长期延迟之间就更不会存在显著的相关关系了。如果存在短期相关性，那么就一定不是纯随机序列。对于平稳时间序列，自相关函数随着延迟阶数的增长而逐渐趋于零，如果考虑的延迟阶数太长，反而会淹没序列的短期相关性[2]。

2.5. 平稳性

• 严平稳与宽平稳

严平稳：随机过程的所有统计性质不受时间起点的影响，不随时间的平移而变化，相同时间步长/间隔内的序列的联合概率分布函数完全相同。

严平稳在数学上的要求较高，在实际场景的运用中往往难以达到，因而研究最多的是宽平稳/弱平稳时间序列分析。需要满足以下三个条件即可认为时间序列是宽平稳（平稳）的：

1. 均值为常数
2. 方差为常数
3. 协方差为常数，且仅依赖于时间跨度s而与时间起点t无关

对于平稳时间序列有一个均值回归特性：如果一个时间序列的平稳的，那么发生在时间t的任何冲击，随着时间的推移都会有一个递减效应，最后消失在t+s，s趋近于无穷。非平稳时间序列中，冲击的影响要么在未来所有的时间持续呈现相同的规模，要么出现增大的趋势。任何白噪声序列都是平稳的。

• 平稳性的意义

时间序列认为每个时间点都是随机变量，而每个变量只有一个样本观察值，我们要运用数理统计方法——样本估计总体，从一个总体也就是一个随机变量中抽取若干个样本观测——就毫无办法。而时间序列又可列多个变量，每个变量在观测的历史进程中都只出现一次，如果时间序列不平稳（随时间的平移而变化），每个变量的均值、方差都只能用一个观察值去估计，效果肯定不佳。但如果时间序列是平稳的，就可以将多个随机变量看做来自相同的总体，他们的观察值就可以看做从同一个总体中抽出的样本，这样就可以利用经典的数理统计方法进行处理[2]。

引用：一般来说，对于每个变量希望样本容量越大越好，因为维度越小分析过程越简单，样本容量越大，分析结果越可靠[3]。

• 平稳性的检验

时序图检验和ACF图检验等图检验法常用于平稳性的初步判定，主观性和经验性较大；ADF单位根检验和PP/KPSS检验等统计检验法应用较多，是较为客观的检验方法，较为准确可靠。

1、图检验法

• 时序图检验

根据平稳时间序列均值、方差均为常数的性质，平稳时间序列的时序图（折线图）应该始终在一个常数值附近随机波动无明显的趋势性或周期性，而且波动的范围有界。

如左下图所示，如果时间序列是平稳的，那么序列应该是围绕某一个均值上下随机波动，而下图中的序列明显具有一定的增长趋势，因此，可以断定该序列肯定不是平稳时间序列。如右下图所示，数据基本上是围绕着0轴上下波动，在0轴上方和下方的波动幅度基本一致。我们可以初步认定，该时间序列是平稳的。

• PP图检验

Q-Q图的结果与P-P图非常相似，只是P-P图是用分布的累计比，而Q-Q图用的是分布的分位数来做检验。QQ图和P-P图一样，如果数据为正态分布，则在Q-Q正态分布图中，数据点应基本在图中对角线上

方法一：参考链接：https://zhuanlan.zhihu.com/p/77786430.
import scipy.stats as scs
scs.probplot(data, sparams=(data.mean(), data.std()), plot=pp_ax)

• ACF图检验

自相关系数和之后阶数之间的关系图。ACF长期位于零轴的一边——单调趋势序列；ACF呈现正弦波动规律——周期变化规律；一直较小且在两倍的标准差范围内——平稳性。

如左下图所示，该时间序列随着阶数的递增，自相关系数先减后增，因此，可以判断该时间序列不是平稳时间序列。如右下图所示，，自相关系数数迅速由1减小到0附近，之后随着阶数的上升在0轴上下小幅波动，这基本上是符合平稳性的要求的。

方法一：参考链接：https://zhuanlan.zhihu.com/p/77786430.
import scipy.stats as scs
scs.probplot(data, sparams=(data.mean(), data.std()), plot=pp_ax)

2、统计检验法

1979年，Dickey和Fuller提出的增广DF(augmented Dickey-Fuller)检验，即ADF检验。

1988年，Phillips和Perron提出的单位根检验方法，即PP检验法。

1992年，Kwiatkowski、Phillips、Schmidt和Shin提出的KPSS检验

ADF检验
方法一：参考链接：https://zhuanlan.zhihu.com/p/35069699
from statsmodels.tsa.stattools import adfuller
tpa_adf = adfuller(tpa_train, autolag='AIC')
方法二：参考链接：https://zhuanlan.zhihu.com/p/77786430
#使用stats子模块中diagnostic（模型诊断）单位根检验unitroot_adf
from statsmodels.stats.diagnostic import unitroot_adf
unitroot_adf(data.close)
 
KPSS检验：
from statsmodels.tsa.stattools import kpss

3. 时间序列分析建模

流程框图如下：

• 步骤1 对时间序列做白噪声检验。通常来说，只有那些序列值之间具有密切的相关关系，历史数据对未来的发展具有一定的影响，才可以用来建模以挖掘历史数据中的有效信息，并预测未来的发展。纯随机性序列不具有任何分析价值。
• 步骤2     对序列做平稳性检验，只有平稳的时间序列才可以用ARMA建模；若序列是非平稳的，则对其平稳化处理。
• 步骤3     ARMA模型建模，包含模型识别、参数估计、模型检验等三个过程。
• 步骤4 利用构建的拟合模型，可以进行动态结构分析以考察研究对象的运行规律，调整输入变量以控制研究对象的发展、对未来变化进行序列预测。

4. 平稳时间序列之ARMA建模

4.1. 模型介绍

ARMA模型是时间序列分析的基本模型，于1970年由美国统计学家Box和英国统计学家Jenkins共同提出，是用来估计平稳的不规则波动或时间序列季节性变动的最常见模型。后续在ARMA模型的基础上引入差分得到ARIMA模型，为了纪念两人在时间序列分析上的突出贡献，ARIMA模型也被称为Box-Jenkins模型。

（1）AR 模型研究第t期的序列值受t-1, t-2, …期序列值以及当期随机干扰值的影响；AR模型反应记忆特征和记忆衰减情况；

（2）MA模型研究序列在t时刻的值与t-1, t-2, … 时刻随机干扰值的相关关系；MA模型考察外部冲击对变量的影响情况和外部冲击的记忆期限。

（3）ARMA模型是将自回归模型与移动平均模型相结合。在这个公式中，p与q分别为自回归模型与移动平均模型的阶数，是需要人为定义的。γi与θi分别是两个模型的相关系数，是需要求解的。

对于平稳非随机的时间序列，ARMA建模主要包含三个步骤：模型识别、参数估计和模型检验。

# 参考链接
from statsmodels.tsa.arima_model import AR, ARMA, ARIMA
 
# MA模型
ma_model = ARMA(ts, order=(0, w))
 
# MA模型
y_hat_avg['moving_avg_forecast'] = train['Count'].rolling(60).mean().iloc[-1]
plt.figure(figsize=(16,8))
plt.plot(train['Count'], label='Train')
plt.plot(test['Count'], label='Test')
plt.plot(y_hat_avg['moving_avg_forecast'], label='Moving Average Forecast')
plt.legend(loc='best')
plt.show()
 
from sklearn.metrics import mean_squared_error
from math import sqrt
rms = sqrt(mean_squared_error(test['Count'], y_hat_avg['moving_avg_forecast'])) #计算均方根误差
print(rms)

4.2. 模型识别

模型识别问题是指选择合适的ARMA(p, q)模型，也称作模型定阶过程，估计自相关阶数p和移动平均阶数q。模型识别或模型定阶的主要方法有两种：一种是利用自相关系数和偏自相关系数，一种是最优信息准则

• 利用自相关系数和偏自相关系数

截尾：在最初d阶时明显大于2倍的标准差范围，而后面几乎95.5%以上的都落在2倍标准差范围之内，并且由非零值衰减为小值波动的过程非常突然，可认为d阶截尾。

拖尾：有超过5%的值落入2倍标准差范围之外，或者显著衰减为小值波动的过程比较缓慢，则拖尾。另一方面，如果值帅见到2倍标准差范围内呈现“伪正弦波动”或其他有明显规律的波动是，也认为是拖尾现象。

如图所示，自相关图里自相关系数没有快速的减为0（一般认为自相关系数低于2倍标准差即图中蓝色虚线一下时即为0），而是呈现出拖尾的特征，故判断序列为非平稳序列。

经验上，寻找模型的滞后阶数时，从1开始到样本大小的1/4比较合理。

# 参考链接：https://zhuanlan.zhihu.com/p/87609082
方法一
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
 
方法二
from statsmodels.tsa.stattools import acf, pacf
 
方法三(可能是老版本)
import statsmodels.tsa.api as smt
smt.graphics.plot_acf(data, lags=lags, ax=acf_ax, alpha=0.5)
smt.graphics.plot_pacf(data, lags=lags, ax=pacf_ax, alpha=0.5)

• 最优信息准则

根据最简约原则，常见的做法是实现选择多组(p, q)组合，通过参数估计建立多个ARMA模型，同时经过模型检验确定机组都能通过检验的(p, q)组合，然后根据AIC、BIC、HQIC等信息准则确定一组最好的(p, q)值。信息准则值越小，模型越好。

AIC=-2ln(模型的极大似然函数值)+2(模型中未知参数的个数)，前者刻画了一个随机事件的信息量以及模型的拟合效果，后者兼顾了未知参数的个数。
BIC= 2ln(模型的极大似然函数值)+ln(模型中未知参数的个数。

1. 补充1：信息准则：统计上来说，如果一个事件的发生给人们带来了信息，那么应该认为该事件是一个随机事件；显然，一个为人们所完全预料的事件，不会给人们带来信息；小概率的事件能带给人们更多信息。）
2. 补充2：最简约的模型是指数据拟合的非常好，与此同时使用最少的参数。使用较少参数的合理的拟合，是估计的一般原则；一个极其过参数化的模型并没有告诉我们很多关于生成数据的过程和时间的性质。
3. 补充3：AIC不能给出模型阶的相合估计，即当样本容量趋于无穷大时，AIC准则确定的模型阶数不能收敛到真实阶，会比真实阶更高。BIC准则弥补了这一缺陷，且可由贝叶斯理论推导出，也叫BSC准则。在样本数目非常大时，AIC偏向于过参数化模型，BIC加重了对过参数化模型的惩罚，而HQIC结合了两者，效果更准确。
4. 补充4：拟合优度中常见的衡量标准有R2、均方误差等。只有在比较有相同数目的解释变量的现行模型时才可以使用R2，由于比较模型的滞后阶数不同，不能用R2。而比较拟合模型与真实数据之间的均方误差也可以选择最优的模型。此外还可以采用统计学上的F检验法：通过比较ARMA(p, q)模型与ARMA(p-1, q-1)模型的残差平方和，判定模型阶数降低后与原来模型之间是否存在显著性差异的方法。

方法一：参考链接：https://zhuanlan.zhihu.com/p/54985638
import statsmodels.api as sm
sm.tsa.arma_order_select_ic(data,max_ar=6,max_ma=4,ic='aic')['aic_min_order']  # AIC
 
方法二：https://zhuanlan.zhihu.com/p/52084533
from statsmodels.tsa.stattools import adfuller
adfuller(timeseries, regression='c', autolag='AIC')

4.3. 参数估计

模型识别后，利用时间序列观察值确定模型的未知参数。模型的参数决定了不同时刻随机变量之间的相依关系，反映了随机变量随时间变化的记忆性大小和记忆期的长短。常用的有矩估计、最小二乘估计、极大似然估计。参数估计的过程在不同时间序列分析软件中有许多可用的标准构件，只需要调用接口即可，不需要读者掌握太多的细节和数学求解过程。

补充[2]：

1. 矩估计利用Yule-Walker方程求解。一般来说，低阶ARMA模型的矩估计具有计算量小、估计思想简单直观，且不需要假设总体分布的优点。但是矩估计只用到了p+q个样本自相关系数，即样本二阶矩的信息，时间序列的其他信息都被忽略了，导致了它是一种比较粗糙的估计方法，常用作极大似然估计和最小二乘估计的迭代计算初始值
2. 极大似然估计的思想是样本来自使该样本出现概率最大的总体，因此未知参数的极大似然估计就是使得似然函数（即联合密度函数）达到最大的参数值。极大似然估计充分应用了每一个观察值，估计精度高，同时还具有估计一致性、渐近正态性和渐进有效性等许多优良的统计性质，是一种非常优良的参数估计方法。
3. 最小二乘估计的思想是使残差平方和达到最小的那组参数值即为最小二乘估计值。最小二乘估计也充分应用了每一个观察值，估计精度高。

4.4. 模型检验

一个好的拟合模型应该确保提取了时间序列中几乎所有的样本相关信息，换言之，拟合残差项中不在蕴含任何相关信息，残差序列应为白噪声序列。换而言之，如果残差序列为非白噪声序列，则残差序列中还残留在相关信息未被提取，说明模型不够有效。因此，模型的显著性检验即为残差序列的随机性检验（白噪声检验），可采用观察残差序列时序图、残差序列Q-Q图、残差序列自相关函数图、以及Box-Ljung检验等方法检验。

模型识别时，为了简化问题提出来很多假设，这些假设往往因人而异带有主管因素，因此必须对模型本身进行检验。如果拟合模型未能通过检验，说明模型不够有效，通常需要选择其他模型重新拟合。

此外，还可以对模型的参数进行显著性检验。检验模型的每一个参数是否显著异于零，对于贡献不大的参数可以剔除以减小模型冗余度使模型精简。模型参数的显著性检验常用t检验法。

参考链接：

[1]时间序列分析方法与应用，elements of time series econometrics an applied approach.

[2]时间序列分析。潘雄锋等

[3]应用时间序列分析。白晓东

[4]https://blog.csdn.net/linchuhai/article/details/87920764   

[5]https://zhuanlan.zhihu.com/p/56093524 

[6]http://km.oa.com/group/34294/articles/show/388706?kmref=search&from_page=1&no=1

[7]https://otexts.com/fppcn/stationarity.html

[8]https://blog.csdn.net/Yuting_Sunshine/article/details/95317735?utm_medium=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-5.edu_weight&depth_1-utm_source=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-5.edu_weight  如何理解自相关和偏自相关图(最全面的讲解)