ARMA模型的概念

AR模型有偏自相关函数截尾性质；
MA模型有相关函数截尾性质。
有些因果线性时间序列有与AR和MA类似的表现，
但是不能在低阶实现偏自相关函数截尾或者相关函数截尾。

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

ARMA(1,1)模型为

Xt=ϕ0+ϕ1Xt−1+εt+θ1εt−1

或

Xt−ϕ1Xt−1=ϕ0+εt+θ1εt−1

其中|ϕ1|<1, |θ1|<1,
−ϕ1≠θ1。
{εt}是独立同分布零均值白噪声列,
εt与Xt−1,Xt−2,…独立。

一般的ARMA(p,q)类似。

AR(p)可以看成ARMA(p, 0),
MA(q)可以看成是ARMA(0, q)。

在ARMA(1,1)的系数条件中，
|ϕ1|<1是平稳解条件，
|θ1|<1是可逆性条件，
−ϕ1≠θ1是为了模型不至于退化：
模型也可以写成

(1−ϕ1B)Xt=ϕ0+(1+θ1B)εt

如果不加这个条件，两边的滞后算子的多项式就可以消去。

ARMA模型的性质

形式地，

11−ϕ1z=11−ϕ1B=1+θ1z1−ϕ1z=1+θ1B1−ϕ1B=∑j=0∞ϕj1zj,∑j=0∞ϕj1Bj,1+(ϕ1+θ1)∑j=1∞ϕj−11zj,1+(ϕ1+θ1)∑j=1∞ϕj−11Bj,

于是

Xt===11−ϕ1B{ϕ0+(1+θ1B)εt}ϕ01−ϕ1+1+θ1B1−ϕ1Bεtϕ01−ϕ1+εt+(ϕ1+θ1)∑j=1∞ϕj−11εt−j(6.1)

这是因果型线性时间序列，是弱平稳的，满足上述ARMA(1,1)模型方程。

EXt=Var(Xt)=ϕ01−ϕ1,σ2(1+(ϕ1+θ1)2∑j=1∞(ϕ21)j−1)=σ21+θ21+2ϕ1θ11−ϕ21

由方差的表达式可知|ϕ1|<1是平稳性的必要条件。

求自协方差函数和自相关函数。
因为协方差与均值无关，
而ϕ0只影响到均值，
所以求协方差与自相关函数时可以假定ϕ0=0。
对

Xt−ϕ1Xt−1=εt+θ1εt−1(6.2)

在(6.2)两边乘以Xt−1并取期望得

γ1−ϕ1γ0=θ1E(εt−1Xt−1)

由(*)可知E(εtXt)=σ2,
所以E(εt−1Xt−1)=σ2

γ1−ϕ1γ0=σ2θ1,

γ1=σ2(ϕ1+θ1)(1+θ1ϕ1)1−ϕ21

在(6.2)两边乘以Xt−k(k≥2)并取期望，得

γk−ϕ1γk−1=0,γk=ϕ1γk−1=ϕk−11γ1

所以ARMA(1,1)的自相关函数为

ρk=⎧⎩⎨⎪⎪(ϕ1+θ1)(1+ϕ1θ1)1+2ϕ1θ1+θ21,ϕ1ρk−1=ϕk−11ρ1,k=1k≥2

所以ARMA(1,1)的ACF与AR(1)的ACF很相似，
但是从k=2处才开始负指数衰减。
与AR类似，
自相关函数不能有限步截尾。

ARMA(1,1)的偏自相关函数与MA(1)的偏自相关函数类似，
但负指数衰减从k=2开始，
也不能在有限步截尾。

总之，
ARMA(1,1)的平稳性条件与AR(1)相同，
自相关函数与偏自相关函数均不能有限步截尾
(设ϕ1≠0, θ1≠0)。

一般ARMA模型

一般ARMA模型为

Xt=ϕ0+ϕ1Xt−1+⋯+ϕpXt−p+εt+θ1εt−1+⋯+θqεt−q

其中{ε}为独立同分布零均值白噪声列，
εt与Xt−1,Xt−2,…独立。
1−ϕ1z−⋯−ϕpzp称为特征多项式，
特征多项式的根都在单位园外，这是平稳性条件。
一般要求1+θ1z+⋯+θqzq的根也都在单位圆外，
这个条件称为可逆性条件。
两个多项式没有公共根，
否则同一模型可能会有不同的表示。

平稳解的均值为

EXt=ϕ01−ϕ1−⋯−ϕp

ARMA模型辨识

可以逐个从低阶模型尝试，
p+q越小越好，
找到AIC最小的选择，
用精确最大似然或者条件最大似然方法估计参数。
对残差进行白噪声检验以验证模型是否充分。

R的forecast包提供了一个auto.arima()函数，
可以自动进行模型选择。
TSA包提供了一个armasubsets()函数用于模型选择。

Tsay和Tiao(1984)提出了一个对ARMA定阶的辅助工具EACF，
其结果可以用与(p,q)有关的二维表格表示，
结果包含由字母“O”组成的三角形，
锐角顶点在(p,q)位置。如

AR012345MA0XX∗∗∗∗1XOX∗∗∗2XOOX∗∗3XOOOX∗4XOOOOX5XOOOOO6XOOOOO7XOOOOO

例6.1 考虑3M公司股票从1946年2月到2008年12月的月对数收益率，
共有755个观测。

d <- read_table2(
  "m-3m4608.txt",
  col_types=cols(.default=col_double(),
                 date=col_date(format="%Y%m%d")))
## Warning: `read_table2()` was deprecated in readr 2.0.0.
## ℹ Please use `read_table()` instead.
mmm <- xts(log(1 + d[["rtn"]]), d$date)
rm(d)
tclass(mmm) <- "yearmon"
ts.3m <- ts(coredata(mmm), start=c(1946,2), frequency=12)
head(ts.3m)
##          Series 1
## [1,] -0.081125460
## [2,]  0.018421282
## [3,] -0.105360516
## [4,]  0.190518702
## [5,]  0.005114897
## [6,]  0.073743834
plot(ts.3m, main="3M Monthly Log Return")
abline(h=0, col="gray")


3M公司月对数收益率

图6.1: 3M公司月对数收益率

ACF图形：

forecast::Acf(ts.3m, main="")


3M公司月对数收益率的ACF

图6.2: 3M公司月对数收益率的ACF

ACF很接近于白噪声。

PACF图形：


3M公司月对数收益率的PACF

图6.3: 3M公司月对数收益率的PACF

PACF也比较接近于白噪声但是有比较多的超出界限的值，
尽管超出量不大。

用TSA包提供的eacf()函数辨识模型：

## Registered S3 methods overwritten by 'TSA':
##   method       from    
##   fitted.Arima forecast
##   plot.Arima   forecast
## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12
## 0 o o x o o x o o o x o  x  o 
## 1 x o x o o x o o o o o  x  o 
## 2 x x x o o x o o o o o  o  o 
## 3 x x x o o o o o o o o  o  o 
## 4 x o x o o o o o o o o  o  o 
## 5 x x x o x o o o o o o  o  o 
## 6 x x x x x o o o o o o  o  o

从上面的结果不易选择ARMA阶。
函数的默认检验系数非零的检验水平是0.05水平。
如果用0.01水平，
则选择p=q=0。

TSA包还提供了armasubsets()函数用来选择ARMA模型阶，
办法是用长阶自回归获得新息εt的估计，
然后用普通最小二乘估计ARMA系数，
用回归的自变量选择方法进行模型选择。
如：

resr <- TSA::armasubsets(ts.3m, nar = 6, nma = 12)
plot(resr)

北京大学金融时间序列分析讲义第6章： ARMA模型

选项nar和nma分别是AR部分和MA部分最多考虑的阶数，
图形中方格的每一行是一种回归子集选择，
默认采用BIC比较不同模型，
可以用scale选项选择其它比较准则，
结果第一行对应于BIC最优的模型，
从结果看，
是仅有截距项和θ12的模型。

用forecasts包的auto.arima()函数定阶：

forecast::auto.arima(ts.3m, max.p = 6, max.q = 6, max.P = 1, max.Q = 1)
## Series: ts.3m 
## ARIMA(3,0,1)(1,0,1)[12] with non-zero mean 
## 
## Coefficients:
##          ar1      ar2      ar3      ma1    sar1     sma1    mean
##       0.0453  -0.0285  -0.0837  -0.1124  0.5319  -0.4435  0.0103
## s.e.  0.3146   0.0417   0.0387   0.3147  0.2885   0.3049  0.0023
## 
## sigma^2 = 0.003998:  log likelihood = 1016.63
## AIC=-2017.25   AICc=-2017.06   BIC=-1980.24

auto.arima()函数选择了一个季节ARMA模型。
参数max.p, max.q, max.P, max.Q分别指定要考虑的模型各个最大阶数。

ARMA模型预测

ARMA(p,q)的预测与AR和MA的预测都有关系。
对AR部分仍是递推地向前预测，
对MA部分，
需要估计新息εt的值。

超前一步预测：

x̂ h(1)=ϕ0+∑j=1pϕjxh+1−j+∑j=1qθjε̂ h+1−j

其中ε̂ i=E(εi|X1,…,Xh)。

对超前k步预测有

x̂ h(k)=ϕ0+∑j=1pϕjx̂ h+k−j+∑j=1qθjε̂ h+k−j

其中h+k−j≤h时x̂ h+k−j=xh+k−j,
h+k−j>h时ε̂ h+k−j=0。
预测误差为eh(k)=xh+k−x̂ h+k。

ARMA模型的三种表示

设ARMA(p,q)模型的系数满足平稳性条件与可逆性条件。

第一种表示：

(1−ϕ1B−⋯−ϕpBp)Xt=ϕ0+(1+θ1B+⋯+θqBq)εt

ARMA模型的MA表示

用滞后算子的性质，
令P(z)=1−∑pj=1ϕpzj,
Q(z)=1+∑qj=1θqzj，
则

Ψ(z)=Q(z)P(z)=∑j=0∞ψjzj,

其中{ψj}绝对可和（实际上，j→∞时cj以负指数速度衰减）。
所以

Xt=μ+Ψ(B)εt=μ+∑j=0∞ψjεt−j, t∈ℤ

这称为ARMA模型的MA表示或者Wold表示。
μ=EXt=ϕ0/P(1)=ϕ0/(1−ϕ1−⋯−ϕp)。
ψ0=1。
系数{ψj,j=0,1,…}称为ARMA模型的脉冲响应函数或者Wold系数。
ψj是脉冲响应函数的意义是，
如果ϵt=1，
它将给Xt+j施加一个ψj的增量影响。

ψj随j→+∞以负指数速度衰减，
这样的性质是合理的，
即一个时刻发生的事情的影响应该随历史向前推进而迅速降低。

MA表示对估计多步预测的均方误差有影响。
理论上

E(Xh+k|ℱh)=μ+∑i=0∞ψi+kεh−i(6.3)

其中ℱh是包含{Xs,s≤t}的最小σ代数，
表示到时刻h为止的观测信息，
E(εh+k|ℱh)=0, k≥1。
于是理论的预测误差为

eh(k)=Xh+k−E(Xh+k|ℱh)=∑j=0k−1ψjεh+k−j

均方误差为

Var(eh(k))=σ2∑j=0k−1ψ2j

因为{ψj}绝对可和，
所以(6.3)中的级数当k→∞时趋于0，
于是长期预测值趋于μ，
这称为均值反转。
当k→∞是多步预测的均方误差趋于
σ2∑∞j=0ψ2j=Var(Xt)，
即用均值作为预测的均方误差。

ARMA模型的AR表示

当平稳性与可逆性条件都成立时，
令

π(z)=P(z)Q(z)

则

π(z)=∑j=0∞πjzj,

{πj}绝对可和，
当j→∞时πj以负指数速度收敛到0。
π0=1。
于是ARMA模型有如下的AR表示

π(B)Xt=ϕ0Q(1)+εt

或

Xt=ϕ̃ 0−π1Xt−1−π2Xt−2−⋯+εt

其中ϕ̃ 0=ϕ0/Q(1)=ϕ0/(1+θ1+⋯+θq)。
这称为ARMA的AR表示或者长阶自回归形式。
这说明平稳可逆ARMA序列可以用自回归模型近似表示。
πj称为ARMA模型的π权重。

附录

文献和历史：

  • EVGENIJ EVGENIEVICH SLUTZKY和GEORGE UDNY YULE在19世纪初提出了MA模型和AR模型。
  • HERMAN WOLD(1938),
    A Study in the Analysis of Stationary Time Series,
    Almquist and Wicksell, Stockholm. 博士论文。
    完善了线性时间序列模型。
  • GEORGE E.P. BOX and GWILYM M. JENKINS(1970),
    Time Series Analysis: Forecasting and Control,
    Holden Day, San Francisco et al.; 2nd enlarged edition 1976.
    本书给出了应用线性时间序列模型的理论与方法。
  • 下面的文献提出，
    一元时间序列模型短期预测常常优于大量经济变量的联立方程的预测：
    CLIVE W.J. GRANGER and PAUL NEWBOLD(1975),
    Economic Forecasting: The Atheist’s Viewpoint, in:
    G.A. RENTON (ed.), Modelling the Economy, Heinemann, London
    pp. 131 – 148.