在上一篇，我们介绍了时间序列统计建模的几个基本概念：平稳性、自相关和白噪声。时间序列数据的建模是在平稳性的基础上进行的。然而，即使平稳的时间序列依然有潜在无限多的参数（各阶自相关函数），我们需要的是更加简约的（也就是参数个数有限的）平稳时间序列模型，ARIMA模型就是这样一种简约的模型。白噪声是用来建立ARIMA的基本模块。

1、滞后算子

首先介绍一种符号：滞后算子B（也称后向算子或延迟算子，有时候也记做L），它在时间序列的模型表示中很有用。滞后算子的作用是将序列值向后移动（倒退）一个时期：

滞后算子也可以多次运用：

比如假设

是某个月数据，

就表示一年前这个月的数据。

2、AR(p)模型

对于具有显著自相关的时间序列，滞后项​很可能会对预测​有用。因此，我们可以将序列{​

}建模为过去观测（也就是滞后项）的加权平均再加上一个白噪声误差的形式，这就得到了自回归模型AR(p)，其中p称为模型的阶数，它对应于模型中包括进来的滞后项的个数。满足AR(p)模型的序列也可以叫做AR(p)过程。我们从最简单自回归模型AR(1)和AR(2)开始讨论AR(p)模型。

2.1 AR(1)模型

基本的AR(1)模型如下：

​

(1）

其中，系数​

称为权重，{​

}是均值为0，方差为

​的白噪声序列。

公式（1）也可以用滞后算子的形式来表示：

​

（2）

2.1.1 AR模型的性质

AR模型具有很多性质。作为以数据分析为目标的教程，我们更关心与模型建立和数据特征有关的性质：平稳性和自相关性。

(1) 平稳性条件

AR(1)模型的平稳条件取决于​

。

假定满足AR(1)模型的序列平稳，由公式（1）有，该序列的期望和方差

​

​

因为平稳序列的方差是非负的常数，可知：​

这就是AR(1)平稳的条件。

(2) 自相关函数

由公式（1），在满足平稳性的条件下，AR(1)的自相关函数满足：

​

​，l=2,3,...,t

所以，根据​

的取值，AR(1)的自相关函数ACF图像有两种模式：

（a）对正的

​,呈现为指数衰减的图像

（b）对负的

​ ,呈现为横轴上下两组以​衰减的图像。

AR(1)的自相关函数只依赖​

这一个参数，这是是非常节约的，代价就是它的ACF只具有上述两种模式，它只能描述具有这两种自相关模式的时间序列。如果我们从某个序列的数据中看到的样本ACF不满足以上两种模式，AR(1)就是不合适于这个序列。此时，对这个序列的建模需要添加更多的参数，或者切换到其他类型的模型。

例：对AR(1)过程的模拟

使用base R中stats包的函数arima.sim()模拟一个AR(1)过程：

​

观察ACF图像（图-1）中的震荡衰减现象（为便于观察，只取了前100个模拟值和前20阶的ACF进行展示）。

####AR模型的模拟——stats::arima.sim
library(tidyverse)
library(tsibble)
library(cowplot)
library(feasts)####ACF函数（自制）
gg_acf<-function(ACF,col,T){ACF<-bind_rows(as_tibble(ACF),as_tibble_row(c(lag=0,acf=1)))ACF%>%ggplot(aes(x=lag,y=acf))+geom_ribbon(aes(ymin=-2/sqrt(T),ymax=2/sqrt(T)),fill=col,alpha=0.4)+geom_segment(aes(x=lag,xend=lag,y=0,yend=acf),color=col)+geom_point(shape=21,color=col,fill=col)+theme_bw()}set.seed(1234)
ar1_sim<-arima.sim(n = 100, list(ar = -0.8))
ar1_sim_tbl<-ar1_sim%>%as_tsibble()p1<-ar1_sim_tbl%>%
ggplot(aes(x=index,y=value))+
geom_line()+theme_bw()p2<-ar1_sim_tbl%>%
ACF(lag_max=50)%>%gg_acf(col='red',T=100)cowplot::plot_grid(p1,p2,nrow=1)

由样本值（即观测值或模拟值）计算出的ACF称为样本ACF。由于观测或抽样的随机性，计算出的样本ACF与由公式推导的理论结果不一定十分吻合。

（3）非平稳的AR(1)

若​

，AR(1)过程是非平稳的，其均值、方差和自相关都不是常数。通常可分为两种情况：

• 随机行走

满足

​时的AR(1)过程称为随机行走。随机行走是一种著名的随机过程，在很多领域，特别是经济和金融中非常有用。

随机行走的均值依赖于（任意的）起始点，而方差随时间线性增加。正是渐增的方差使随机行走过程实现“漫步”，而不是呈现出均值回复的形态。

• ​
  的AR(1)序列

在形态上，会具有爆炸性的行为。

例：模拟随机行走

为了更好地展示随机行走的特点，下面我们采用交互的可视化包plotly模拟了20个从0出发的随机行走序列（图-2）。模拟是这样进行的，考虑到随机行走模型的迭代特点，序列​可以表示为一系列正态分布随机数的累加和。

library(plotly)set.seed(1234)p1<-plot_ly()for(i in 1:20){rw<-cumsum(c(0,rnorm(500,0,1)))p1<-p1%>%add_lines(x=0:(length(rw)-1),y=as.numeric(rw))}
p1

2.2 AR(2)

AR(2)模型为：

​

(3)

（1）特征方程与自相关函数

假设AR(2)过程平稳，由公式（3），它的自相关函数方程满足：

​

​

​

(​​

)

将上述第三个等式用滞后算子进行表示，得到：

​​

这是一个二阶差分方程，它决定了平稳AR(2)的自相关函数。

与上述二阶差分方程对应，一元二次方程：

​

称为AR(2)的特征方程，它的两个解可以由求根公式得出。

这两个解的倒数称为AR(2)模型的特征根（characteristic root）。

如果两个特征根都是实数值的，则模型的二阶差分方程能分解为​

，这时AR(2)可以看作一个AR(1)在另一个AR(1)上的运算。此时AR(2)的ACF图像是两个指数衰减的AR(1)的混合。如果两个解都是复根（共轭复根），AR(2)的ACF图像将呈现出幅度减少的正弦和余弦图像。

（2）平稳性

由差分方程的数学性质（参见[1]），AR(2)平稳的条件是两个特征根的模都小于1（若是两个实根，就是特征根的绝对值小于1）。具体来说，​

和​

满足：

从几何的观点，AR(2)平稳的条件是它的两个特征根都位于一个单位圆内。这个几何性质可以用来检验平稳性。

例：模拟AR(2)过程

使用arima.sim()，对如下AR(2)过程进行模拟（图-3）：

这个过程有两个复数特征根：

ar2_sim<-arima.sim(n = 1000, list(ar = c(1.3,-0.7)))
ar2_sim_tbl<-ar2_sim%>%as_tsibble()p1<-ar2_sim_tbl[1:100,]%>%
ggplot(aes(x=index,y=value))+
geom_line()+theme_bw()p2<-ar2_sim_tbl%>%
ACF(lag_max=20)%>%gg_acf(col='red',T=100)cowplot::plot_grid(p1,p2,nrow=1)

2.3 AR(p)

将上述两个模型进行推广，对于p>2，可以建立AR(p)模型：

​

它的特征方程是：​

AR(p)的平稳条件：特征根（特征方程解的倒数）的模都小于1，或者等价的说，这些特征根都位于一个单位圆内。

平稳AR(p)的ACF图像呈现减幅的正弦、余弦和指数衰减的混合形式，具体形式取决于特征根的性质。

2.3.1偏自相关函数

对于一个平稳的AR(p)过程，k阶滞后的自相关系数​

并不仅仅只反映了

​和​

之间的相关关系。比如对于AR(1)过程，​

和

​之间的自相关系数​

等于​​

与

​的自相关系数​

和​

与​

的自相关系数（依然是

）的乘积。因此，为了只衡量​对于​的影响，可以在控制了​

的取值之后（也就是不再把它们看作随机变量），再计算​

与​

之间的相关系数，这样得到的相关系数叫做

​与

的偏相关系数，这些偏自相关系数的集合叫做{​

}的偏自相关函数（PACF）

PACF可以这样来计算：在

​的AR(k-1)模型中引入新变量

​，可以得到一个新的模型AR(k)。这个新模型中​

（第k项）的系数

就是​

与

​的偏相关系数。它表示了添加​

之后，这个新变量对​

的影响或贡献。

对一个AR(p)模型，滞后为p的PACF应不为零，而对所有j>p的PACF应接近于0。这个性质叫做AR(p)模型的偏自相关函数截尾性。可以利用这一点，帮助确定AR(p)的阶数。

2.3.2 AR模型的建立和筛选

建立AR模型

可以使用极大似然估计或者条件最小二乘估计得到AR模型的参数。估计的数学细节不做讨论。R有多个函数可用于 估计AR模型，特别的，fable包提供了函数AR()，用于估计AR模型。

筛选AR模型

在建立AR(p)模型的时候，我们并不知道真实的阶数p是多少，那么该如何确定一个恰当的阶数p？这需要在模型的复杂度和简洁性之间进行选择。增加滞后项p的个数，可以使残差平方和减少，提升拟合的优度，然而这也会导致过拟合，降低预测的效果。时间序列模型筛选的标准和其他统计模型或机器学习模型的筛选思路是一致的：在保证一定复杂度的前提下，使模型尽可能的简洁。

根据AR(p)模型偏自相关函数的截尾性质，我们可以通过偏自相关图像来判断AR模型的阶数。然而在实际运用中，由于数据并不一定“完美”地符合理论模型，偏自相关图像的截尾性可能表现相当复杂的情况。因此，这一性质可以作为一种经验法则识别AR模型，要筛选出最简洁的模型，信息准则还是更可靠的方法。

信息准则基于似然函数，对模型的参数个数（复杂度）进行惩罚，根据惩罚函数的不同有多种具体的形式。我们常用的信息准则包括赤池信息量AIC和贝叶斯信息准则BIC。

（1）赤池信息量

AIC=-2/T*ln(似然函数)+2*（参数个数）/T

其中，T是样本容量（样本的大小）。AIC的第一项表示了模型的拟合优度（拟合的好坏程度），第二项是罚函数，对过多的参数加以惩罚。

（2）贝叶斯信息准则

将AIC公式第二项中的2替换为ln(T)，可以得到贝叶斯信息准则BIC

AIC对每个参数的惩罚为2，BIC对参数的惩罚为ln(T)。所以与AIC相比，在样本容量较大时，BIC倾向于选择一个低阶的AR模型。

例：采用TSstudio包中的数据集Coffee_Prices，这个数据集提供了两种咖啡豆的价格序列。我们选择其中的Robusta（罗布斯塔）的数据建立AR模型。首先，选取原数据集的一个子集并转化为tsibble格式，再对Robusta的价格序列计算差分（考虑到Robusta的价格序列可能并不平稳）。

接下来，通过观察Robusta的价格序列（value）与差分序列（value_diff）的自相关图（图-4），可见value的自相关图呈现缓慢衰减，可以判断是不平稳的，因此选择差分序列进行建模。观察差分序列value_diff的偏自相关图（图-5），一阶滞后的偏自相关是显著不为0的，作为练习，使用AR()函数对p=1到10的情况拟合AR模型。

函数AR的用法如下：

AR(formula, ic = c("aicc", "aic", "bic"), ...)

其中，formula是进行模型的指定，可以这样来完成：order(p = 0:15, fixed = list())。p参数指定模型的阶数，如果同时指定多个阶数，AR()会根据信息准则的最小值进行筛选（信息准则分别为AIC、AICc和BIC，默认是AIC）。fixed参数用来指定具体的模型系数，如fixed = list(ar1 = 0.3, ar3 = 0)。我们对拟合的模型分别用AIC和BIC准则进行了筛选。代码和运行结果如下：

library(fable)
library(TSstudio)coffee_prices<-as_tsibble(Coffee_Prices)%>%filter(key=="Robusta")%>%mutate(value_diff=difference(value))###价格序列的自相关
coffee_prices%>%
ACF(lag_max=100)%>%gg_acf(col='red',T=nrow(coffee_prices))###差分序列的偏自相关
coffee_prices%>%
select(-value)%>%
ACF(type="partial",lag_max=50)%>%gg_acf(col='red',T=nrow(coffee_prices)-1)+ylab("pacf")   ###拟合，分别根据aic和bic选择模型
fit_ar1<-series[2:nrow(coffee_prices),]%>%
model(AR(value_diff~order(1:10),ic="aic"))fit_ar2<-series[2:nrow(series),]%>%
model(AR(value_diff~order(1:10),ic="bic"))report(fit_ar1)
##Series: value_diff
##Model: AR(1) ##Coefficients:
##     ar1
##  0.2778##sigma^2 estimated as 0.007137
##AIC = -16.57 AICc = -16.55  BIC = -13.18report(fit_ar2)
##Series: value_diff
##Model: AR(1) ##Coefficients:
##     ar1
##  0.2778##sigma^2 estimated as 0.007137
##AIC = -16.57 AICc = -16.55  BIC = -13.18

由结果可知，AIC和BIC准则都选择了AR(1)模型。

2.3.3 模型检验

如果拟合的模型是充分的，其残差序列应该是白噪声（即不存在自相关），为了防止模型的不充分性，需要对拟合模型的残差序列进行白噪声检验。使用feasts包的函数gg_tsresiduals()可以对拟合模型的残差序列进行可视化，包括残差的时间序列图，残差序列的自相关图和残差的直方图（以便和正态分布进行比较）（图-6）。使用fable包函数residuals()提取拟合模型的残差序列，并使用函数ljung_box()进行白噪声检验：

###残差的图形
gg_tsresiduals(fit_ar1)###残差的白噪声检验
fit_res<-residuals(fit_ar1)
ljung_box(fit_res[,3])
##lb_stat    lb_pvalue
##3.930085e-05 9.949981e-01 

从检验结果可知，拟合模型fit_ar1的残差序列是一个白噪声。

如果残差序列的ACF中显示出显著的相关性，应该考虑对模型加以扩展。

说明：

（1）R中提供了很多时间序列拟合和检验的方法。我们这个系列文章主要采用tidyverst的tsibble包、feasts包和fable包提供的方法。需要注意的是，tidyverst的主要目标是时间序列预测，而且还在开发之中，功能还不够全面。很多其他R工具或许会给出更全面的结果，比如说作为tidyverst前身forecast包。然而，forecast并不支持tidyverse语法。

（2）关于AR模型的检验。对于拟合的AR(k)模型，除了检验残差序列符合白噪声之外，还应该检验模型系数的显著性，如果某个系数不显著，可以在拟合模型中去掉该系数对应的项，这样可以简化拟合模型。具体做法可以参考[2]

参考：

[1]、汉密尔顿，《时间序列分析》，中国人民大学出版社，2015.

[2]、蔡瑞胸，《金融数据分析导论》，机械工业出版社，2016.