北大金融时间序列分析讲义第4章： 自回归模型

自回归模型的概念

如果ρ1≠0，
则Xt与Xt−1相关，
可以用Xt−1预测Xt。
最简单的预测为线性组合，
如下模型：

Xt=ϕ0+ϕ1Xt−1+εt(4.1)

称为一阶自回归模型(Autoregression model)，记作AR(1)模型。
其中{εt}是零均值独立同分布白噪声序列，
方差为σ2，
并设εt与Xt−1,Xt−2,…独立。
系数|ϕ1|<1。
更一般的定义中仅要求{εt}是零均值白噪声，
不要求独立同分布。

这个模型与一元线性回归模型Yi=ϕ0+ϕ1xi+εi在某些方面类似，
比如，
εt都是起到误差或者扰动的作用。
但是，
在自回归模型中，自变量Xt−1在时刻t−1时作为因变量，
所以自回归模型中因变量和自变量不是两个变量而是同一个变量的不同时刻。
回归模型的理论结果不能直接应用到自回归模型当中。

AR(1)模型也是马尔可夫(Markov)过程：
Xt在Xt−1,Xt−2,…条件下的条件分布，
只与Xt−1有关。
已知Xt−1后，
用Xt−1,Xt−2,…去预测Xt，
与仅用Xt−1去预测的效果相同。
这种性质称为马氏性。

条件期望和条件方差：

E(Xt|Xt−1)=ϕ0+ϕ1Xt−1,Var(Xt|Xt−1)=σ2

即在Xt−1=xt−1已知后，
Xt条件的条件分布是期望为ϕ0+ϕ1xt−1,
方差为σ2的分布。
可以证明

Var(Xt)=σ21−ϕ21

因为|ϕ1|<1，
所以Xt在Xt−1=xt−1已知条件下的条件方差小于其无条件方差，
也就是说用Xt−1的信息去预测Xt，
可以使得Xt的波动减小，
能够达到预测的效果。

AR(1)模型的推广是AR(p)模型：

Xt=ϕ0+ϕ1Xt−1+⋯+ϕpXt−p+εt(4.2)

其中{εt}是零均值独立同分布白噪声序列，方差为σ2，
且εt与Xt−1,Xt−2,…独立。
系数ϕ1,…,ϕp需要满足如下条件：方程

1−ϕ1z−⋯−ϕpzp=0

的所有复根z∗都满足|z∗|>1，
上述方程的左边的多项式称为AR(p)模型的特征多项式，
特征多项式的所有复根称为特征根，
对系数的条件称为“特征根都在单位圆外”。
在有的教材中将特征多项式的根的倒数定义为特征根，
在(Tsay 2013)中就是以倒数为特征根。
对AR(1)就是要求|ϕ1|<1。

更一般的定义中仅要求{εt}是零均值白噪声，
不要求独立同分布。

滞后算子

设{ξt,t∈ℤ}为弱平稳时间序列或者常数列，
定义如下的滞后算子B:

Bξt=ξt−1,Bjξt=ξt−j, j∈ℤ

考虑复变函数P(z)=∑∞j=0ajzj，
其中{aj}为实数列，
设∑∞j=0|aj|<∞（称为系数绝对可和），
有限阶多项式为P(z)的特例。
定义

P(B)ξt=∑j=0∞ajξj−t.

P(B)也称为一个（常系数时齐）线性滤波器，
P(B)ξt是{ξt}序列的一个滑动平均变换。

设Q(z)=∑∞j=0bjzj，{bj}为实数列，绝对可和，
则C(z)=P(z)Q(z)=∑∞j=0cjzj，
其中

cj=∑i=0jaibj−i, j=0,1,…

且{cj}绝对可和，
称数列{cj}是数列{aj}和数列{bj}的离散卷积。
对任意弱平稳列或者常数列{ξt}均有

P(B)Q(B)ξt=Q(B)P(B)ξt=C(B)ξt, t∈ℤ

若ξt≡ξ与t无关，
则Bjξ=ξ。
比如，
Bj1=1，
P(B)1=P(1)。

令D(z)=P(z)Q(z)，
若Q(z)≠0对任意满足|z|≤1的复数z均成立，
则D(z)=∑∞j=0djzj,
{dj}绝对可和，且

P(B)ξt=Q(B)D(B)ξt=D(B)Q(B)ξt, t∈ℤ

AR(1)模型的性质

研究AR模型的性质，
可以在对实际数据选择适当模型时，
知道那些情况下使用AR模型，
使用何种AR模型是合适的。

对理论证明感兴趣的同学请参见(何书元 2003)的论述。

AR(1)模型(4.1)要求|ϕ1|<1，
这是(4.1)中的{Xt}有弱平稳解的充分必要条件。

充分性：
当|ϕ1|<1且ϕ0=0时，
因为

11−ϕ1z=∑j=0∞ϕj1zj,11−ϕ1B=∑j=0∞ϕj1Bj

取

X∗t=∑j=0∞ϕj1εt−j

右侧的随机变量级数均方收敛且a.s.收敛。
将上式代入Xt=ϕ1Xt−1+εt，
直接验证或者利用滞后算子性质可以证明X∗t满足方程。
所以ϕ0=0, |ϕ1|<1时AR(1)模型有X∗t这样的因果线性时间序列形式的弱平稳解。
对于一般的ϕ0，因为平稳要求EXt=EXt−1=μ所以要求
μ=ϕ0+ϕ1μ即μ=ϕ01−ϕ1，
令

X∗t=μ+∑j=0∞ϕj1εt−j

可以验证其满足AR(1)模型，
是因果线性时间序列形式的弱平稳解。

必要性证明：
如果模型有弱平稳解Xt，
因为要求εt与Xt−1,Xt−2,…独立，
在模型两边求方差得

Var(Xt)=Var(ϕ0+ϕ1Xt−1+εt)=ϕ21Var(Xt−1)+σ2,σ2X=ϕ21σ2X+σ2,σ2X=σ21−ϕ21

因为σ2X>0所以1−ϕ21>0, |ϕ1|<1。

上述证明过程也证明了

μ=EXt=ϕ01−ϕ1,σX=Var(Xt)=σ21−ϕ21.

注：可以证明，
ϕ1=1时，
一定没有弱平稳列满足方程Xt=ϕ0+ϕ1Xt−1+εt。
当|ϕ1|>1时，
取

Xt=μ−∑j=1∞ϕ−j1εt+j, t∈ℤ

则{Xt}是满足Xt=ϕ0+ϕ1Xt−1+εt的弱平稳时间序列，
但εt不再与Xt−1,Xt−2,…相互独立。

由μ=ϕ0/(1−ϕ1)得ϕ0=(1−ϕ1)μ，
AR(1)模型可以写成

Xt=(1−ϕ1)μ+ϕ1Xt−1+εt

即t时刻的值，
是历史均值与t−1时刻值的加权平均，
加上一个与历史值相互独立的扰动。

当ϕ1=0时，AR(1)模型就退化成白噪声，
所以AR(1)模型通常都认为ϕ1≠0。
|ϕ1|越接近于1，
说明前后的相关性越强。

AR(1)模型的自相关函数

AR(1)模型的平稳解满足ϵt与Xt−1,Xt−2,…独立。
将模型写成

Xt−μ=ϕ1(Xt−1−μ)+εt(4.3)

在(4.3)两边乘以εt后取期望，有

E[εt(Xt−μ)]=ϕ1E[εt(Xt−1−μ)]+σ2=σ2

在(4.3)两边乘以Xt−μ后取期望，有

γ0=ϕ1γ1+σ2

在(4.3)两边乘以Xt−j−μ后取期望，有

γj=ϕ1γj−1, j=1,2,…

即有

γ0=γj=σ2X=σ21−ϕ21=ϕ1γ1+σ2ϕ1γj−1, j=1,2,…

于是ACF为

ρj=γjγ0=ϕj1, j=1,2,…

当0<ϕ1<1时，
ACF为正的单调下降序列，
以负指数速度（几何级数）下降。
当−1<ϕ1<0时，
ACF为负正交替的序列，
绝对值以负指数速度下降。

例：ϕ1=0.8的ACF图形：

tmp.x <- 0:15
tmp.y <- 0.8^tmp.x
plot(tmp.x, tmp.y, type="h", xlab="lag", ylab="ACF",
     ylim=c(-1, 1))
abline(h=0)

例：ϕ1=−0.8的ACF图形：

tmp.x <- 0:15
tmp.y <- (-0.8)^tmp.x
plot(tmp.x, tmp.y, type="h", xlab="lag", ylab="ACF",
     ylim=c(-1, 1))
abline(h=0)

AR(2)模型的性质

AR(2)模型为

Xt=ϕ0+ϕ1Xt−1+ϕ2Xt−2+εt(4.4)

其中{εt}为独立同分布的零均值白噪声,
εt与Xt−1,Xt−2,…独立。
系数满足特征方程的根都在单位圆外要求。
特征方程为

1−ϕ1z−ϕ2z2=0(4.5)

这个模型有因果线性时间序列形式的弱平稳解。
对(4.4)两边取期望得μ=EXt的方程

μ=ϕ0+ϕ1μ+ϕ2μ

所以

μ=ϕ01−ϕ1−ϕ2

将(4.4)中的Xt减去μ，
方程变成

Xt−μ=ϕ1(Xt−1−μ)+ϕ2(Xt−2−μ)+εt(4.6)

在(4.6)两边乘以Xt−j−μ，
并利用εt与Xt−1,Xt−2,…独立的性质，可得

γj=ϕ1γj−1+ϕ2γj−2, j=1,2,…

两边除以γ0得

ρj=ϕ1ρj−1+ϕ2ρj−2, j=1,2,…(4.7)

这种形式的递推称为齐次线性差分方程。
(4.7)可以用滞后算子写成

(1−ϕ1B−ϕ2B2)ρj=0, j=1,2,…(4.8)

可以看出，
递推中用到的1−ϕ1B−ϕ2B2
就是将特征多项式中的z替换成了滞后算子B，
这样递推得到的ρj的性质自然就是由特征多项式决定的。

(4.7)的递推需要初值ρ0和ρ1。
ρ0=1，而

ρ1=ϕ1ρ0+ϕ2ρ1,ρ1=ϕ11−ϕ2

这样就可以给定ϕ1,ϕ2后递推计算出ρj,j=1,2,3,…。
特征根都在单位圆外的条件可以保证limj→∞ρj=0
（参见(何书元 2003)§2.1），且衰减速度为负指数速度。
limj→∞ρj=0也是弱平稳列的必要条件。

特征多项式1−ϕ1z−ϕ2z2=0的根的判别式为
Δ=ϕ21+4ϕ2，
当ϕ2≥−4ϕ21时为实根，
当ϕ2<−4ϕ21时为复根。
两个特征根都在单位圆外的充分必要条件是

|ϕ1|<2, |ϕ2|<1, ϕ2<1±ϕ1

图4.1的三角形阴影部分代表了(ϕ1,ϕ2)使得模型有平稳解的区域，
不包含三角形边界；
蓝色区域以及蓝色和红色边界线为两个实特征根的区域，
红色区域为两个共轭复根的区域。

有两个实根时，两个根为

z1,z2=ϕ1±ϕ21+4ϕ2‾‾‾‾‾‾‾‾‾√−2ϕ2

这两个根的绝对值都大于1。
付过两个根都是正的，则ρj为正，呈现出负指数速度单调衰减。
如果两个根正负不同，则绝对值仍以负指数速度衰减但是不单调，
当负根绝对值更小时会正负交替衰减。

有共轭复根时，两个根为

z1,z2=ϕ1−2ϕ2±i|ϕ21+4ϕ2‾‾‾‾‾‾‾‾‾√|−2ϕ2

复根的模为

|zi|=1−ϕ2‾‾‾‾√

辐角为

ω=cos−1ϕ12−ϕ2‾‾‾‾√

反过来有

ϕ2=−1|z1|2, ϕ2=2|z1|cosω

ω这个辐角对应的周期为

P=2πω=2πcos−1ϕ12−ϕ2√

这样的AR(2)模型会体现出随机地平均以周期P的波动；
ρj序列呈现出以周期P震荡的幅度负指数衰减的变化。

下面对不同的ϕ1,ϕ2取值画出ACF的典型图形。

上图中的模型AR(2)分别为：

Xt=Xt=Xt=Xt=Xt=Xt=Xt=Xt=0.1Xt−1+0.5Xt−2+εt−0.1Xt−1+0.5Xt−2+εt0.8Xt−1+εtXt−1−0.1Xt−2+εt−Xt−1−0.1Xt−2+εt(21.1cosπ6)Xt−1−11.12Xt−2+εt(21.1cosπ2)Xt−1−11.12Xt−2+εt(21.1cos2π3)Xt−1−11.12Xt−2+εt

例4.1 考虑美国的国民生产总值(GNP)经过季节调整后的季度增长率。
时间是1947年第二季度到2010年第一季度，总计252个观测值。

读入数据并计算对数增长率：

da <- read_table("q-gnp4710.txt", col_types=cols(.default = col_double()))
gnp <- ts(da[["VALUE"]], start=c(1947, 1), frequency=4)

这个序列从1947年第一季度到2010年第一季度共253个观测。

GNP的对数值的图形：

plot(log(gnp), xlab="year", ylab="log(GNP)")

以上的GNP是经过季节调整的，所以看不出季节性。

计算GNP的对数值的差分作为对数增长率，并作图：

rate <- diff(log(gnp))
plot(rate, xlab="Year", ylab="Growth")
abline(h=0, col="gray")

增长率的ACF估计：

forecast::Acf(rate, xlab="Years", main="")

这个图形呈现出一定周期的震荡衰减，
可尝试AR建模。

拟合的模型为

(1−0.438B−0.206B2+0.156B3)(Xt−0.016)=εt, σ2=Var(ϵt)=9.55×10−5

R函数polyroot(x)可以求出多项式的所有复根，
x是多项式按升幂排列的系数。

roots <- polyroot(c(1, -0.438, -0.206, 0.156)); roots

## [1]  1.614790+0.866168i -1.909068+0.000000i  1.614790-0.866168i

## [1] 1.832429 1.909068 1.832429

有一个实根−1.9091和一对共轭复根1.6148±0.8662i,
复根的模为1.8324，
复根对应的周期为

## [1] 12.7619

12.8个季度，约为3年。这是美国GNP的随机周期，
也是ACF震荡的周期。

时间序列的周期研究是很常用的工具。

○○○○○

AR(p)模型的性质

对于一般的AR(p)模型，
其ACF的性质以及序列的随机周期，
也由其特征根决定。
ACF可以是单调衰减、震荡衰减、正负交替衰减、呈周期震荡衰减。
在有复特征根根或者有接近−1的特征根时时间序列呈现出一定的随机周期变化。

由平稳性得

μ=ϕ01−ϕ1−⋯−ϕp

自相关函数(ACF)满足如下的递推（差分方程）

(1−ϕ1B−⋯−ϕpBp)ρj=0, j=1,2,…

AR(p)模型的平稳解是线性时间序列。

偏自相关函数

实际数据用AR模型建模时，阶数p是未知的，
确定p的问题称为定阶。
一般常用偏自相关函数和AIC准则。

设X1,…,Xn,Y为随机变量，

L(Y|X1,…,Xn)=argminŶ =b0+b1X1+⋯+bnXnE(Y−Ŷ )2

称为用X1,…,Xn对Y的最优线性预测。
Y−L(Y|X1,…,Xn)与Z−L(Z|X1,…,Xn)
的相关系数称为Y和Z在扣除X1,…,Xn影响后的偏相关系数。

对平稳线性时间序列，
对n=1,2,…，
有

L(Xt|Xt−1,…,Xt−n)=ϕn0+ϕn1Xt−1+⋯+ϕnnXt−n

其中ϕnj,j=0,1,…,n与t无关。
称ϕnn为时间序列{Xt}的偏自相关系数，
{ϕnn}序列称为时间序列{Xt}的偏自相关函数(PACF)。

ϕnn实际是Xt与Xt−n在扣除Xt−2,…,Xt−n+1的影响后的偏相关系数。
ϕ11就是ρ1。

ϕnn用样本进行估计，
得到的估计值ϕ̂ nn,n=1,2,…称为样本偏自相关函数。
在R软件中用pacf(x)估计并作图。

如果{Xt}服从如下AR(p)模型：

Xt=ϕ0+ϕ1Xt−1+⋯+ϕpXt−p+εt, ϕp≠0

这意味着用Xt−1,Xt−2,…的线性组合预测Xt时，
只需要用到Xt−1,…,Xt−p，
增加Xt−p−1,Xt−p−2,…不能改进预测。
这意味着ϕkk=0, k>p。
这种性质叫做AR模型的偏自相关函数截尾性。

AR(p)序列的样本偏自相关函数ϕ̂ kk满足如下性质：

• T→∞时ϕ̂ pp→ϕp≠0。
• 对k>p, ϕ̂ kk→0(T→∞)。
• 对k>p，ϕ̂ kk渐近方差为1T。

这样，可以用类似对ACF的白噪声检验那样给PACF图画出±2T√的上下界限，
以此判断PACF在哪里截尾。

例4.2 考虑CRSP价值加权指数月度收益率，
时间从1926-1到2008-12。

d <- read_table(
  "m-ibm3dx2608.txt",
  col_types=cols(.default=col_double(),
                 date=col_date(format="%Y%m%d")))
ibmind <- xts(as.matrix(d[,-1]), d$date)
tclass(ibmind) <- "yearmon"
vw <- ts(coredata(ibmind)[,"vwrtn"], start=c(1926,1), frequency=12)
head(ibmind)

##             ibmrtn     vwrtn     ewrtn     sprtn
## 1月 1926 -0.010381  0.000724  0.023174  0.022472
## 2月 1926 -0.024476 -0.033374 -0.053510 -0.043956
## 3月 1926 -0.115591 -0.064341 -0.096824 -0.059113
## 4月 1926  0.089783  0.038358  0.032946  0.022688
## 5月 1926  0.036932  0.012172  0.001035  0.007679
## 6月 1926  0.068493  0.056888  0.050487  0.043184

plot(vw, main="CRSP Value Weighted Index Monthly Return")
abline(h=0, col="gray")

forecast::Acf(vw, main="")

ACF到k=21还没有截尾。
PACF图:

两倍的渐近标准差为:

## [1] 0.06337243

PACF值的列表：

tmp <- pacf(vw, main="", plot=FALSE, lag.max=12)
cbind(k=round(tmp$lag*12), pacf=round(tmp$acf, 4))

##        k    pacf
##  [1,]  1  0.1154
##  [2,]  2 -0.0304
##  [3,]  3 -0.1025
##  [4,]  4  0.0326
##  [5,]  5  0.0618
##  [6,]  6 -0.0502
##  [7,]  7  0.0312
##  [8,]  8  0.0517
##  [9,]  9  0.0635
## [10,] 10  0.0053
## [11,] 11 -0.0052
## [12,] 12  0.0109

如果仅看前12个PACF，
取p=3应该可以，
但是实际上PACF在k=17仍未截尾。

例4.3 考虑例4.1中的美国经过季节调整后的GNP季度对数增长率数据。

forecast::Acf(rate, xlab="Years", main="")

ACF明显不截尾。

pacf(rate, xlab="Years", main="")

PACF虽然在k=3,9,14,16等位置超出界限，
但是超出不多，
可考虑用AR(3)建模。

信息准则

信息准则是统计建模中常用的模型比较工具，
其基本思想是模型拟合数据的拟合优度与模型简单化的折衷。

AIC准则（Akaike’s Information Criterion）：

AIC=−2Tln(似然函数值)+2T(参数个数)

其中似然函数值是在参数最大似然估计处的似然函数值。
当模型为高斯AR(p)，即{ϵt}是独立同N(0,σ2)序列时的AR(p)模型时，
AIC公式为

AIC(k)=lnσ̃ 2k+2kT

其中k是模型的阶，
σ̃ 2k是阶为k的条件下εt的方差的最大似然估计。
lnσ̃ 2k代表了模型对数据的拟合优劣，
此值越大拟合越差；
2kT是对模型复杂程度的惩罚，
此值越大，
模型越复杂，
稳定性越差，
对未来的情况的适应性也越差。
在某个范围内取k使得AIC(k)最小，
就达成了拟合优度与模型简单程度的折衷。

另一个常用的信息准则是BIC准则(Bayesian Information Criterion)，
高斯AR模型为：

BIC(k)=lnσ̃ 2k+klnTT

BIC倾向于取比AIC更低阶的模型。

可以取k=0,1,…,P0计算AIC或BIC，
取最小值点的k。
P0可取为10log10T。

参考：

• HIROTUGU AKAIKE,
  Fitting Autoregressive Models for Prediction,
  Annals of the Institute of Statistical Mathematics AC-19 (1974), pp. 364 – 385
• GIDEON SCHWARZ,
  Estimating the Dimensions of a Model,
  Annals of Statistics 6(1978), pp. 461 – 464

例4.4 考虑例4.1中的美国经过季节调整后的GNP季度对数增长率数据的AR建模定阶。

基本R软件stats包中的ar()函数可以对时间序列样本进行AR建模，
默认采用AIC准则定阶。
用选项aic=FALSE, order.max=p可以指定p阶模型。

gnprate <- diff(log(gnp))
resm <- ar(gnprate, method="mle"); resm

## 
## Call:
## ar(x = gnprate, method = "mle")
## 
## Coefficients:
##       1        2        3        4        5        6        7        8  
##  0.4318   0.1985  -0.1180   0.0189  -0.1607   0.0900   0.0615  -0.0814  
##       9  
##  0.1940  
## 
## Order selected 9  sigma^2 estimated as  8.918e-05

产生了一个9阶的模型。
作其AIC的图形：

plot(as.numeric(names(resm$aic)), resm$aic, type="h",
     xlab="k", ylab="AIC")

从AIC图形可以看出，如果取较低的阶，
3阶也是可以的。

虽然ar()函数没有提供BIC的值，
但是因为BIC(k)−AIC(k)=k(lnT−2)T,
可以自己计算：

tmp.T <- length(gnprate)
tmp.bic <- resm$aic + as.numeric(names(resm$aic)) * (log(tmp.T) - 2)/tmp.T
plot(as.numeric(names(resm$aic)), tmp.bic, type="h",
     xlab="k", ylab="BIC")

BIC也建议p=9。

AR模型参数估计方法

AR模型有多种估计方法，
比如，
用普通线性回归的最小二乘法估计，
假设正态分布用最大似然估计，
Yule-Walker递推计算，
Burg递推计算，等等。

在stats::ar()函数中，
method="ols"为最小二乘估计，
设ϕi的估计为ϕ̂ i，
则拟合值为

x̂ t=ϕ̂ 0+ϕ̂ 1xt−1+⋯+ϕ̂ pxt−p, t=p+1,…,T

残差为

et=xt−x̂ t, t=p+1,…,T

相应的新息方差σ2=Var(εt)的估计为

σ̂ 2=1T−2p−1∑t=p+1Te2t

如果使用高斯条件最大似然估计（认为x1,…,xp固定），
则ϕ̂ i估计不变，
但是新息方差得估计变成了

σ̃ 2=1T−p∑t=p+1Te2t=T−2p−1T−p∑t=p+1Te2t

例4.5 对例4.5的CRSP价值加权指数月度收益率用AR(3)建模。

先用AIC定阶：

resm <- ar(vw, method="mle"); resm

## 
## Call:
## ar(x = vw, method = "mle")
## 
## Coefficients:
##       1        2        3        4        5        6        7        8  
##  0.1167  -0.0112  -0.1126   0.0217   0.0735  -0.0452   0.0254   0.0462  
##       9  
##  0.0660  
## 
## Order selected 9  sigma^2 estimated as  0.002831

用AIC定阶为9阶。AIC数值：

##     0     1     2     3     4     5     6     7     8     9    10    11    12 
## 22.33 10.99 12.07  3.35  4.37  2.46  1.96  3.04  2.24  0.00  1.97  3.94  5.81

从AIC数值看出，选较低阶时p=3是一个可以接受的选择。
使用arima()函数估计，这个函数的输出比较容易解释。

mean.vw <- mean(vw); mean.vw

## [1] 0.00890439

resm2 <- arima(vw, order=c(3,0,0))
resm2

## 
## Call:
## arima(x = vw, order = c(3, 0, 0))
## 
## Coefficients:
##          ar1      ar2      ar3  intercept
##       0.1158  -0.0187  -0.1042     0.0089
## s.e.  0.0315   0.0317   0.0317     0.0017
## 
## sigma^2 estimated as 0.002875:  log likelihood = 1500.86,  aic = -2991.73

结果中的intercept实际是均值μ的估计，
μ̂ =0.0089。
σ̂ =0.054。
模型为:

(Xt−0.0089)=0.1158(Xt−1−0.0089)−0.0187(Xt−2−0.0089)−0.1042(Xt−3−0.0089)

结果还给出了系数的标准误差。
以θ̂ ±2SE(θ̂ )作为简单的近似95%置信区间，
结果标准误差说明μ,ρ1,ρ3是显著不等于0的，
而ρ2不显著。

μ̂ =0.008904，
这是CRSP价值加权指数的平均月度简单收益率。
折合成年收益率，再换算到1926-2008年共83年的总简单收益率：

## [1] 0.1122442

((1 + mean.vw)^12)^83 - 1

## [1] 6832.002

按平均月度简单增长率0.8904%计算，
得到年度简单增长率为11.22%；
83年可增长6832倍。

按实际数据83年的实际增长计算83年的总简单增长率，
然后折合到一年的复利：

## [1] 1591.953

## [1] 0.09290084

实际的复利年均增长率为9.29%，而不是11.22%;
实际的83年的简单增长率是1592倍，
而不是6832倍。

用μ̂ 计算的结果与用真实数据计算的结果的83年增长率有这么大差距，
是因为μ̂ >0，
按μ̂ 计算总是在按几何级数增长；
但是按∏(1+xt)计算则有时增长有时下降，
所以实际的总增长率要低得多。

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

tseries包的arma()函数可以用条件极小二乘方法估计ARMA模型参数，
使用了optim()函数来求极值。结果可以用summary()显示。
tseries::arma()的结果还支持print(), plot(), coef(), vcov(),
residuals(), fitted()等信息提取函数。

○○○○○

AR模型检验

拟合AR模型后，
如果模型是能够准确表示数据的规律的，
则拟合的残差应该近似为白噪声，
应该能通过白噪声检验。
因为残差是从模型估计计算得到的，
自由度有损失，
在Box.test()中用fitdf=p指定自由度减少个数。

例4.6 例如，对例4.5
的CRSP价值加权指数月简单收益率的AR(3)模型的残差进行白噪声检验：

resm2 <- arima(vw, order=c(3,0,0))
Box.test(resm2$residuals, lag=12, type="Ljung", fitdf=3)

## 
##  Box-Ljung test
## 
## data:  resm2$residuals
## X-squared = 16.352, df = 9, p-value = 0.05988

结果在0.05水平下不显著，说明模型拟合是充分的。

从例4.5的估计结果看出ϕ2不限制。
在arima()函数中可以用fixed=指定某些系数固定为给定值，
NA表示不规定。这样做的一个潜在问题是得到的模型可能不符合平稳性条件。

resm3 <- arima(vw, order=c(3,0,0), fixed=c(NA, 0, NA, NA)); resm3

## Warning in arima(vw, order = c(3, 0, 0), fixed = c(NA, 0, NA, NA)): some AR
## parameters were fixed: setting transform.pars = FALSE

## 
## Call:
## arima(x = vw, order = c(3, 0, 0), fixed = c(NA, 0, NA, NA))
## 
## Coefficients:
##          ar1  ar2      ar3  intercept
##       0.1136    0  -0.1063     0.0089
## s.e.  0.0313    0   0.0315     0.0017
## 
## sigma^2 estimated as 0.002876:  log likelihood = 1500.69,  aic = -2993.38

检验稳定性：

abs(polyroot(c(1, -coef(resm3)[1:3])))

## [1] 2.031585 2.279433 2.031585

三个根都在单位圆外，符合平稳性条件。

对固定ϕ2=0后的模型进行残差的白噪声检验：

Box.test(resm3$residuals, lag=12, type="Ljung", fitdf=2)

## 
##  Box-Ljung test
## 
## data:  resm3$residuals
## X-squared = 16.828, df = 10, p-value = 0.07827

不否认残差为白噪声。

也可以用forecast::checkresiduals()函数进行模型诊断：

forecast::checkresiduals(resm3)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(3,0,0) with non-zero mean
## Q* = 46.896, df = 21, p-value = 0.0009699
## 
## Model df: 3.   Total lags used: 24

结果图形中包括残差，
残差的ACF图，残差的直方图和正态密度拟合。
此函数自动进行了残差的Ljung-Box白噪声检验，
但它减去的自由度包括均值参数占用的一个自由度。
可以用lag选项指定Ljung-Box检验使用的平方和项数。

stats包的tsdiag()函数也可以对模型残差进行诊断，
包括标准化残差图、残差的ACF图、残差的Ljung-Box白噪声检验p值，
包括使用不同平方和项数的检验p值，
但是可能没有进行自由度调整。
程序如：

○○○○○

AR模型拟合优度指标

类似于线性回归模型的拟合优度判断，
在线性时间序列建模中，也可以定义如下的拟合优度R2统计量

R2=1−残差平方和序列的总离差平方和

比如，拟合了AR(p)模型后，可以计算R2为

R2=1−∑Tt=p+1e2t∑Tt=p+1(xt−x¯)2

其中x¯=1T−p∑Tt=p+1xt。
0≤R2≤1。
R2越大，
说明模型对数据拟合越好，
残差越小。
这样的指标仅对线性时间序列有意义。
如果是非平稳的单位根过程，
拟合AR模型的R2会趋于1。

只要数据不变，R2随阶数p的增加而增大。
根据统计学中的精简性原则（称为Ocam’s razor原则），
应尽可能选择简单模型。
所以提出调整的R2指标(adjusted R2)，定义为

R2Adj=1−σ̂ 2σ̂ 2x

其中σ̂ 2是新息方差σ2的估计，
σ̂ 2x是Xt的样本方差。
仍有0≤R2Adj≤1，
但R2Adj对阶数p的增加有惩罚。

用估计的AR模型进行预测

前面定义的L(Y|X1,…,Xn)是最佳线性预测。
还可以定义Ŷ 为X1,…,Xn的任意函数，
使得E(Y−Ŷ )2最小，
称为最佳预测（最优预测）。
最佳预测等于条件期望E(Y|X1,…,Xn)。

对AR(p)模型（新息为独立同分布白噪声），
最佳预测与最佳线性预测相同。

在时刻h，用截止到h位置的信息预测xh+ℓ，
得到的最佳预测记为x̂ h(ℓ)。

resm4 <- arima(vw[1:984], order=c(3,0,0))
pred4 <- predict(resm4, n.ahead=12, se.fit=TRUE)
cbind(Observed=round(c(vw[985:996]), 4), 
      Predict=round(c(pred4$pred), 4), 
      SE=round(c(pred4$se), 3))

##       Observed Predict    SE
##  [1,]  -0.0623  0.0075 0.053
##  [2,]  -0.0220  0.0160 0.054
##  [3,]  -0.0105  0.0117 0.054
##  [4,]   0.0511  0.0098 0.054
##  [5,]   0.0238  0.0088 0.054
##  [6,]  -0.0786  0.0092 0.054
##  [7,]  -0.0132  0.0094 0.054
##  [8,]   0.0110  0.0096 0.054
##  [9,]  -0.0981  0.0095 0.054
## [10,]  -0.1847  0.0095 0.054
## [11,]  -0.0852  0.0095 0.054
## [12,]   0.0215  0.0095 0.054

plot(window(vw, start=c(2007, 1), end=c(2008,12)),
     ylab="", type="b", ylim=c(-0.2, 0.2), lwd=2)
lines(ts(c(pred4$pred), start=c(2008,1), frequency = 12), 
      col="red", lwd=1, lty=2, type="b", pch=2)
lines(ts(c(pred4$pred) - 2*c(pred4$se), start=c(2008,1), frequency = 12), 
      col="green", lwd=1, lty=3, type="l")
lines(ts(c(pred4$pred) + 2*c(pred4$se), start=c(2008,1), frequency = 12), 
      col="green", lwd=1, lty=3, type="l")

#foreplot(pred=pred4, rt=vw, orig=984, start=(83-2)*12+1, p=0.95)