Loading [MathJax]/jax/output/HTML-CSS/jax.js

时间序列分析:Yule-Walker 估计示例

1. AR(2)

  考虑 AR(2) 过程

Xt=ϕ1Xt1+ϕ2Xt2+et

其中 etN(0,σ2e)et 独立于 Xt1,Xt2,,并假设序列具有零均值。

1.1. 估计 ϕ

  首先计算样本的自相关 r1r2,由方程

[r1r2]=[1r1r11][ˆϕ1ˆϕ2]

可以解得参数 ˆϕ1ˆϕ2

1.2. 估计方差

  对式 (1) 等号两边计算方差,得

Var(Xt)=ϕ21Var(Xt1)+ϕ22Var(Xt2)+2ϕ1ϕ2Cov(Xt1,Xt2)+σ2e

可以解得

σ2e=γ0ϕ21γ1+ϕ22γ22ϕ1ϕ2γ1=γ0(1ϕ21ϕ222ϕ1ϕ2ρ1)

  由 Yule-Walker 方程组

[ρ1ρ2]=[1ρ1ρ11][ϕ1ϕ2]

ρ1=ϕ1+ρ1ϕ2ρ2=ρ1ϕ1+ϕ2

于是

1ϕ21ϕ222ϕ1ϕ2ρ1=1ϕ21ϕ1ϕ2ρ1ϕ22ϕ1ϕ2ρ1=1ϕ1(ϕ1+ρ1ϕ2)ϕ2(ϕ2+ρ1ϕ1)=1ϕ1ρ1ϕ2ρ2

将上式带入式 (3),得

σ2e=γ0(1ϕ1ρ1ϕ2ρ2)

σ2e 的 Yule-Walker 估计是

ˆσ2e=c0(1ˆϕ1r1ˆϕ2r2)

其中 c0 是样本协方差系数,r1,r2 是样本自相关系数,ˆϕ1,ˆϕ2 是估计的参数。

  更一般地,对 AR(2) 过程 σ2e 的 Yule-Walker 估计是

ˆσ2e=c0(1pi=1ˆϕiri)

1.3. 代码示例

  首先模拟 AR(2) 过程

set.seed(42)
ar.process=arima.sim(10000, model=list(ar=c(1/3, 1/2)), sd=4)
set.seed(42) ar.process=arima.sim(10000, model=list(ar=c(1/3, 1/2)), sd=4)
set.seed(42)
ar.process=arima.sim(10000, model=list(ar=c(1/3, 1/2)), sd=4)

使用的参数为 ϕ1=1/3ϕ1=1/2σ2e=42=16

  计算样本自相关系数

r <- acf(ar.process, plot=F)$acf[2:3]
r <- acf(ar.process, plot=F)$acf[2:3]
r <- acf(ar.process, plot=F)$acf[2:3]

这里 r 的值为 0.6701659498254 0.731929559010564

  构造 Yule-Walker 方程组并求解

R <- matrix(1, 2, 2)
R[1, 2] <- r[1]
R[2, 1] <- r[1]
b <- matrix(r, nrow=2, ncol=1)
phi.hat <- solve(R, b)
R <- matrix(1, 2, 2) R[1, 2] <- r[1] R[2, 1] <- r[1] b <- matrix(r, nrow=2, ncol=1) phi.hat <- solve(R, b)
R <- matrix(1, 2, 2)
R[1, 2] <- r[1]
R[2, 1] <- r[1]

b <- matrix(r, nrow=2, ncol=1)

phi.hat <- solve(R, b)

这里 R 的值为

1.0000000 0.6701659
0.6701659 1.0000000
1.0000000 0.6701659 0.6701659 1.0000000
1.0000000 0.6701659
0.6701659   1.0000000

b 的值为

0.6701659
0.7319296
0.6701659 0.7319296
0.6701659
0.7319296

解得 phi.hat 的值为

0.3261191
0.5133757
0.3261191 0.5133757
0.3261191
0.5133757

可见 phi.hat 的结果与实际参数(1/3,1/2)很接近。

  接着估计 ˆσ2e

c0 <- acf(ar.process, type='covariance', plot=F)$acf[1]
var.hat <- c0 * (1 - sum(phi.hat * r))
sd.hat <- sqrt(var.hat)
c0 <- acf(ar.process, type='covariance', plot=F)$acf[1] var.hat <- c0 * (1 - sum(phi.hat * r)) sd.hat <- sqrt(var.hat)
c0 <- acf(ar.process, type='covariance', plot=F)$acf[1]
var.hat <- c0 * (1 - sum(phi.hat * r))
sd.hat <- sqrt(var.hat)

解得 sd.hat 的值为 4.02197762014053,与实际值 σ2e=4 也很接近。

2. AR(3)

  AR(3) 过程的 Yule-Walker 估计方法与前面类似。首先模拟 AR(3) 过程如下

set.seed(42)
ar.process=arima.sim(10000, model=list(ar=c(1/3, 1/2, 7/100)), sd=4)
set.seed(42) ar.process=arima.sim(10000, model=list(ar=c(1/3, 1/2, 7/100)), sd=4)
set.seed(42)
ar.process=arima.sim(10000, model=list(ar=c(1/3, 1/2, 7/100)), sd=4)

然后进行估计

r <- acf(ar.process, plot=F)$acf[2:4]
R <- matrix(1, 3, 3)
R[1, 2] <- r[1]
R[2, 1] <- r[1]
R[1, 3] <- r[2]
R[3, 1] <- r[2]
R[2, 3] <- r[1]
R[3, 2] <- r[1]
b <- matrix(r, nrow=3, ncol=1)
phi.hat <- solve(R, b)
c0 <- acf(ar.process, type='covariance', plot=F)$acf[1]
var.hat <- c0 * (1 - sum(phi.hat * r))
sd.hat <- sqrt(var.hat)
r <- acf(ar.process, plot=F)acf[2:4]R<matrix(1,3,3)R[1,2]<r[1]R[2,1]<r[1]R[1,3]<r[2]R[3,1]<r[2]R[2,3]<r[1]R[3,2]<r[1]b<matrix(r,nrow=3,ncol=1)phi.hat<solve(R,b)c0<acf(ar.process,type=covariance,plot=F)acf[1] var.hat <- c0 * (1 - sum(phi.hat * r)) sd.hat <- sqrt(var.hat)
r <- acf(ar.process, plot=F)$acf[2:4]

R <- matrix(1, 3, 3)
R[1, 2] <- r[1]
R[2, 1] <- r[1]
R[1, 3] <- r[2]
R[3, 1] <- r[2]
R[2, 3] <- r[1]
R[3, 2] <- r[1]

b <- matrix(r, nrow=3, ncol=1)

phi.hat <- solve(R, b)

c0 <- acf(ar.process, type='covariance', plot=F)$acf[1]
var.hat <- c0 * (1 - sum(phi.hat * r))
sd.hat <- sqrt(var.hat)

得到的 phi.hat

0.3259706
0.5124173
0.0689914
0.3259706 0.5124173 0.0689914
0.3259706
0.5124173
0.0689914

sd.hat4.02255420924941,可见估计结果与实际值都很接近。