ARIMA模型的定阶原理与建模分析

前言

ARIMA模型是很经典的自回归模型，这篇文章将全面的讲述ARIMA的建模步骤。从定阶原理解释到实际数据代码编写模型来进行回归预测。基于理论推导和代码编写一气呵成！

岁月如云，匪我思存，写作不易，望路过的朋友们点赞收藏加关注哈，在此表示感谢！

一：AR$(p)$模型的定阶原理

AR模型是一个线性模型，p阶自回归模型的一般表达式为：

$x_t={\varphi }_0+{\varphi }_1{x}_{t-1}+{\varphi }_2{x}_{t-2}+...+{\varphi }_p{x}_{t-p}+{\epsilon }_t(\#)$ ，

其中 { ε t } \{\varepsilon_t\} {εt​} 是一个白噪声序列，既然AR模型被建立，此AR模型是满足弱平稳条件的，则存在 $\left|{\varphi }_p\right|<1$和自相关系数，以及$E({\epsilon }_t)=0;Var({\epsilon }_t)={\sigma }^2;E({\epsilon }_s{\epsilon }_t)=0,\forall s\ne t$。

• 首先我们先建立AR(2)模型

$x_t={\varphi }_0+{\varphi }_1{x}_{t-1}+{\varphi }_2{x}_{t-2}+{\epsilon }_t,\left|{\varphi }_1\right|<1,\left|{\varphi }_2\right|<1(\ast )$

那么我们对上式 (*) 左右两边各减去 $u$ 得：
$x_t-\mu ={\varphi }_0+{\varphi }_1(x_{t-1}-u)+{\varphi }_2(x_{t-2}-u)+({\varphi }_1+{\varphi }_2-1)\mu +{\epsilon }_t(1)$ ，

又由弱平稳性质， $(\ast )$ 两边取均值可得：

$\mu ={\varphi }_0+{\varphi }_{1}E(x_{t-1})+{\varphi }_{2}E(x_{t-2})+E({\epsilon }_t)$ ,

即 $\mu ={\varphi }_0+{\varphi }_1\mu +{\varphi }_2\mu +0\Rightarrow \mu =\frac{{\varphi }_0}{1-{\varphi }_1-{\varphi }_2}$ ，

那么把此结果带入 $(1)$ 式可得：

$x_t-\mu ={\varphi }_1(x_{t-1}-\mu )+{\varphi }_2(x_{t-2}-u)+{\epsilon }_t(2)$ ，

$(2)$ 式两边乘以 $(x_{t-1}-u)$:

$(x_t-\mu )(x_{t-1}-u)={\varphi }_1(x_{t-1}-\mu )(x_{t-1}-u)+{\varphi }_2(x_{t-2}-u)(x_{t-1}-u)+{\epsilon }_t(x_{t-1}-u)(3)$

$(3)$ 式两边同时取期望后得：

$E\left[(x_t-\mu )(x_{t-1}-u)\right]={\varphi }_{1}E\left[(x_{t-1}-\mu )(x_{t-1}-u)\right]+{\varphi }_{2}E\left[(x_{t-2}-u)(x_{t-1}-u)\right](4)$

对 $(4)$ 式两边再除以方差 ${\sigma }_0$ 之后得 ${\rho }_1={\varphi }_1+{\varphi }_2{\rho }_1$ ，这里的 ${\rho }_1$ 为自相关系数。
则可得 ${\rho }_1={\varphi }_1/(1-{\varphi }_2)$ ，同理 (2) 式两边同时乘以 $(x_{t-2}-\mu )$

可得 ${\rho }_2={\varphi }_1{\rho }_1+{\varphi }_2$。

同理我们推广 $(2)$ 式两边乘以 $(x_{t-k}-\mu )(k\ge 3)$ ，

可得 ${\rho }_k={\varphi }_1{\rho }_{k-1}+{\varphi }_2{\rho }_{k-2}$ 。

从 ${\rho }_k$ 的表达式我们容易发现尽管 $\left|{\varphi }_1\right|<1,\left|{\varphi }_2\right|<1$ ，但 ${\rho }_k$ 永远不会为0，所以会出现拖尾现象。

• 接着我们拓展至AR$(p)$模型

在平稳的前提下，我们容易得 $\mu =\frac{{\varphi }_0}{1-{\varphi }_1-{\varphi }_2-...-{\varphi }_p}$ ，

将 $(\#)$ 两边减去均值 $\mu$ 可得： $x_t-\mu ={\varphi }_1(x_{t-1}-\mu )+{\varphi }_2(x_{t-2}-u)+...+{\varphi }_p(x_{t-p}-\mu )+{\epsilon }_t(5)$ ，

那么对 $(5)$ 左右两边同乘以 $(x_t-\mu )、(x_{t-1}-\mu )、...、$ 并除以方差 ${\sigma }_0$ 可得：

$1={\varphi }_1{\rho }_1+{\varphi }_2{\rho }_2+...+{\varphi }_p{\rho }_p,{\rho }_1={\varphi }_1+{\varphi }_2{\rho }_2+{\varphi }_3{\rho }_3+...+{\varphi }_p{\rho }_{p-1},{\rho }_2={\varphi }_1{\rho }_1+{\varphi }_2+{\varphi }_3{\rho }_3+...+{\varphi }_p{\rho }_{p-2},...,{\rho }_p={\varphi }_1{\rho }_{p-1}+{\varphi }_2{\rho }_{p-2}+{\varphi }_3{\rho }_{p-3}+...+{\varphi }_p$ ，

根据这些关系式，模仿AR(2)的递推关系式可得：
${\rho }_k={\varphi }_1{\rho }_{k-1}+{\varphi }_2{\rho }_{k-2}+...+{\varphi }_k{\rho }_{k-p}(k\ge p)$ ，

因此符合$AR(p)$的平稳序列模型，其自相关系数在$p$ 阶之后一直不会为0，存在所谓拖尾现象。

二：MA$(q)$模型的定阶原理

$MA(q)$模型被称为移动平均模型，一个 $q$ 阶的移动平均模型可以用数学式表达为：
$x_t=\mu +{\epsilon }_t+{\theta }_1{\epsilon }_{t-1}+{\theta }_2{\epsilon }_{t-2}+...+{\theta }_q{\epsilon }_{t-q}$ ，

那么满足的性质有

$E({\epsilon }_t)=0;Var({\epsilon }_t)={\sigma }^2;E({\epsilon }_s{\epsilon }_t)=0,\forall s\ne t$；

$E(x_t)=\mu ;Var(x_t)=(1+{\theta }_1^2+{\theta }_2^2+...+{\theta }_q^2){\sigma }^2$

• 首先我们还是建立$MA(2)$模型

如$x_t=\mu +{\epsilon }_t+{\theta }_1{\epsilon }_{t-1}+{\theta }_2{\epsilon }_{t-2}$ ，

则 $E(x_t)=\mu ,Var(x_t)=(1+{\theta }_1^2+{\theta }_2^2){\sigma }_{ϵ}^2$

对于 $Var(x_t)$ ，两边同时被 ${\sigma }_1=Cov(x_t,x_{t-1})$ 相除有

${\rho }_1=\frac{Cov(x_t,x_{t-1})}{(1+{\theta }_1^2+{\theta }_2^2){\sigma }^2}$ ，

又$Cov(x_t,x_{t-1})=E\left[(x_t-\mu )(x_{t-1}-\mu )\right]=E\left[({\epsilon }_t+{\theta }_1{\epsilon }_{t-1}+{\theta }_2{\epsilon }_{t-2})({\epsilon }_{t-1}+{\theta }_1{\epsilon }_{t-2}+{\theta }_2{\epsilon }_{t-3})\right]=$

$E\left({\epsilon }_{t-1}{\epsilon }_t+{\theta }_1{\epsilon }_{t-2}{\epsilon }_t+{\theta }_2{\epsilon }_{t-3}{\epsilon }_t+{\theta }_1{\epsilon }_{t-2}{\epsilon }_{t-1}+{\theta }_1{\theta }_2{\epsilon }_{t-3}{\epsilon }_{t-1}+{\theta }_2{\epsilon }_{t-2}{\epsilon }_{t-1}+{\theta }_2{\epsilon }_{t-3}{\epsilon }_{t-2}+{\theta }_1{\epsilon }_{t-1}^2+{\theta }_1{\theta }_2{\epsilon }_{t-2}^2\right)$
$={\theta }_{1}E({\epsilon }_{t-1}^2)+{\theta }_1{\theta }_{2}E({\epsilon }_{t-2}^2)$

那么最终 ${\rho }_1=\frac{{\theta }_1+{\theta }_1{\theta }_2}{1+{\theta }_1^2+{\theta }_2^2}$。

如果我们相同的方法求解 ${\rho }_2=\frac{{\theta }_2}{1+{\theta }_1^2+{\theta }_2^2}$ ，那么 ${\rho }_3=0$ 这是显然的。

• 接着我们建立$MA(q)$模型

同理对于$MA(q)$模型，我们经过相同的运算可得最终表达式

${\rho }_l=\frac{{\theta }_l+{\theta }_1{\theta }_{l+1}+{\theta }_2{\theta }_{l+2}+...+{\theta }_{q-l}{\theta }_q}{1+{\theta }_1^2+{\theta }_2^2+...+{\theta }_q^2}$ ，

那么当 $l>q$ 时同理可得 ${\rho }_l=0$ 。

所以，通过上述推导我们有理由相信：$MA(q)$模型的自相关系数 $q$ 阶截尾。所谓$q$ 阶截尾意思是在 $q$ 阶以后$MA(q)$模型的自相关系数立马截止， $q+1$ 阶起就为0。

以上就是通过理论解释了AR§模型和MA(q)模型的拖尾和截尾的底层逻辑。

三：ARMA模型

• 参数估计过程

当把$AR(p)$模型和$MA(q)$模型相结合时，我们得到$ARMA(p,q)$模型如下：

$x_t={\varphi }_0+{\varphi }_1{x}_{t-1}+...+{\varphi }_p{x}_{t-p}+{\epsilon }_t+{\theta }_1{\epsilon }_{t-1}+...+{\theta }_q{\epsilon }_{t-q}$

相较于前两个模型，此模型是更具有普遍性。首先我们通过一些定阶模型确定 $p,q$，当阶数确定后，可以根据最小二乘最大似然估计或者梯度下降法更新所有方程系数。根据模型的表达式一直迭代下去即可完成“无穷的”预测。但是作为长期预测，理论上是可行的，实际确实长期预测所受的干扰因素太多了，除非你的预测数据是周期性、趋势性或者季节性的，那长期还是有点实际意义，否则任何回归模型，还是作为短期预测才有更大的实际意义。

• 建模过程

1：序列判断

(a)：判断我们需要建立的模型数据是否为平稳序列，若非平稳序列我们要对其进行变换处理（一般用差分方法即可）至平稳序列，

(b)：接着再判断平稳序列时候为白噪声序列，若为白噪声序列则建模结束（白噪声序列无法构成ARMA模型），否则进行下一步。

2：模型估计与建立

(a)：判断 $p$ 和 $q$ 的值。当我们建立好自回归模型时，为了得到最优的模型结构，我们需要定下 $p$ 和 $q$ 值。这里的定阶一是可以通过自相关系数$ACF$和偏自相关系数$PACF$大致决定。由上面的理论分析，我们知道$AR(p$)将出现 $p$ 阶拖尾，$MA(q)$将出现 $q$ 阶截尾，

(b)：如果序列的$ACF$和$PACF$不是很明确的话，我们可以用其他模型来定阶。其中就包括AIC和BIC信息准备判别。AIC是一种用于模型选择的指标，同时考虑模型的拟合程度以及简单性，BIC是对AIC的改进，一般来说较小的AIC或者BIC表示在保持模型简单的同时，能够更好的对时间序列进行拟合。

3：模型诊断

即对模型残差进行验证，确保其为服从正态分布的白噪声序列，当模型的残差为白噪声时，说明我们已经将序列的信息充分提取到模型中，建模彻底结束。

在上一篇文章我们对于ARMA模型 $x_t={\sum }_{i=1}^q{\theta }_i{\epsilon }_{t-i}+{\varphi }_0+{\sum }_{j=1}^p{\varphi }_j{L}^j{x}_t$ 分析发现，ARMA其实和AR模型在平稳性上的判断是一样的，都有这相同的特征方程，同样可以通过单位根方法判断是否平稳性成立。

上篇文章地址。。。

四：实际建模运用

我们接下来基于实际销量数据开始建立时序模型，首先观察下销量数据可视化结果，由曲线图发现销量的变化明显具有上涨的趋势性，符合自回归移动平均模型的建模直观要求。

图1

from statsmodels.graphics.tsaplots import plot_acf,plot_pacf
from statsmodels.tsa.stattools import adfuller as ADF
from statsmodels.stats.diagnostic import acorr_ljungbox
from statsmodels.tsa.arima_model import ARIMA
import statsmodels.api as sm
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

yv = np.array([2800,2811,2832,2850,2880,2910,2960,3023,3039,3056,3138,3150,3198,3100,3029,2950,2989,3012,3050,3142,3252,3342,3365,3385,3340,3410,3443,3428,3554,3615,3646,3614,3574,3635,3738,3764,3788,3820,3840,3875,3900,3942,4000,4021,4055])
yv_serie = pd.Series(yv[:-10])##样本外数据

def testwhitenoise(data):
    m = 10# 检验10个自相关系数
    acf,q,p = sm.tsa.acf(data,nlags=m,qstat=True)
    out = np.c_[range(1,m+1),acf[1:],q,p]
    output = pd.DataFrame(out,columns=['lag','自相关系数','统计量Q值','p_values'])
    output = output.set_index('lag')# 设置第一列索引名称,可省略重复索引列1
    print(output)

def teststeady(data,count=0):
    res_ADF = ADF(data)
    print('ADF检验结果为：', res_ADF)
    Pv = res_ADF[1]
    if Pv > 0.05:
        print('\033[1;31mP值：%s，原始序列不平稳，要进行差分！\033[0m' % round(Pv,5))
        count = count + 1
        print('\033[1;32m进行了%s阶差分后的结果如下\033[0m' % count)
        data = data.diff(1).dropna()
        teststeady(data,count)
    else:
        print('\033[1;34mP值：%s，原始序列平稳，继续建模\033[0m'% round(Pv,5))
testwhitenoise(yv_serie)
teststeady(yv_serie)

图2

图2就是平稳性和自相关性（白噪声）检验的结果，我们发现当进行一阶差分后序列平稳，按照建模步骤我们接下来开始定阶。

def confirm_p_q(data):
    fig = plt.figure(figsize=(8,6))
    testwhitenoise(data)
    train = teststeady(data)
    ax1 = fig.add_subplot(211)
    fig = sm.graphics.tsa.plot_pacf(train, lags=10, ax=ax1)
    ax2 = fig.add_subplot(212)
    fig = sm.graphics.tsa.plot_acf(train, lags=10, ax=ax2)
    plt.show()  ###可视化定阶

    pmax = int(len(data) / 10)
    qmax = int(len(data) / 10)
    AIC = sm.tsa.arma_order_select_ic(train,max_ar=pmax,max_ma=qmax,ic='aic')['aic_min_order']
    BIC = sm.tsa.arma_order_select_ic(train,max_ar=pmax,max_ma=qmax,ic='bic')['bic_min_order']
    HQIC = sm.tsa.arma_order_select_ic(train,max_ar=pmax,max_ma=qmax,ic='hqic')['hqic_min_order']
    print('AIC：',AIC)
    print('BIC：',BIC)
    print('HQIC：',HQIC)
    return AIC
pq = confirm_p_q(yv_serie)##返回p,q值

图3

图4

由上面图4自相关函数图可知，定阶在 $p,q=(1,1)$ 阶比较合理，再由相应的信息准则，我们最终定阶 $p,q=(2,2)$ 也是合理的。

这里的定阶结果都是理论给的结果，实际中的定阶还是要根据模型表现不断调整，一般阶数越高越复杂，拟合效果越强，但过拟合概率也越高，所以要不断尝试不断调整。

接着我们正式开始预测

def prediction(data):
    tempmodel = ARMA(teststeady(data),pq).fit(disp=-1)
    print(tempmodel.summary())
    #num = 10
    #predictoutside1 = tempmodel.forecast(num)[0]#预测样本外的
    predictoutside2 = tempmodel.predict(len(tempmodel.predict()),len(tempmodel.predict()) + 9,dynamic=True)##也是样本外预测，预测结果一致
    predictinside = tempmodel.predict()##样本内预测
    init_value = yv[0]

    fig = plt.figure(figsize=(8, 6))
    predictinside = predictinside.cumsum()##差分还原
    pretrueinside = init_value + predictinside
    startprevalue = list(pretrueinside)[-1]
    predictoutside2 = predictoutside2.cumsum()##差分还原
    pretrueoutside = startprevalue + predictoutside2
    
    ##作图
    plt.plot(yv,label='原始值')
    plt.plot([init_value] + list(pretrueinside),label='样本内预测值')
    X = [i for i in range(len(yv)-11,len(yv))]
    plt.plot(X,[startprevalue] + list(pretrueoutside), label='样本外预测值')
    allpredata = [init_value] + list(pretrueinside) + list(pretrueoutside)
    plt.legend()
    plt.show()
    return tempmodel,allpredata
preres = prediction(yv_serie)

最后我们对模型进行评价

def evaluate_model(model,apd):
    delta = model.fittedvalues - tsres
    score = 1 - delta.var() / tsres.var()
    print('R^2：', score)
    allmse = mean_squared_error(apd,yv)##所有预测值跟所有原始值的MSE
    print('ALLMSE:',allmse)

    ###残差白噪声检验
    testwhitenoise(delta)
evaluate_model(preres[0],preres[1])

五：建模结果比较分析

• 当我们选择阶数 $p,q=(1,1)$ 时看下建模效果：

图5：p,q=(1,1)

图6

注：这里涉及两个评价指标，一个是拟合优度 $R^2$ 值，公式如下： $R^2=1-\frac{Va{r}_{残差}}{Va{r}_{样本内x_t}}$ ， $Var$ 是方差意思， $R^2$ 越接近1，说明拟合越好。 另一个是均方误差，公式如下： $MSE=\frac{1}{n}{\sum }_{i=1}^n(x_i-\overline{x_i}{)}^2$ ， $\overline{x_i}$ 是样本估计量（预测值），此实验中，预测值指的是样本内预测值+样本外预测值，样本值是全体原数据值。

从残差检验是白噪声序列后，我们完整的建模算正式结束！

当我们选择阶数 $p,q=(2,2)$ 时看下建模效果：

图7：p,q=(2,2)

图8

由图5和图6比较，直观上感觉图6总体拟合效果更好，再观察理论评价指标，也是 $p,q=(2,2)$ 表现的更好，所以具体定阶时，我们不妨多个指标一起观察。

六：总结

• 此篇文章涉及的内容很多，有详细的理论推导解释AR§拖尾和MA(q)截尾的缘故，并最终一步一步建立ARMA模型来解决实际问题，
• 在上一篇文章我们也谈到ARMA对趋势性，周期性和季节性数据做短期预测是非常有效的，这篇文章主要是对趋势性数据做预测，周期性和季节性当然也是同理而得，
• 对于阶数 $p,q$ 的取定，一直是个非常重要的步骤，所以在实际中，我们一定要结合实际结合多种方法综合定阶。