时间序列分析:Yule-Walker 估计示例
Contents [show]
1. AR(2)
考虑 AR(2) 过程
Xt=ϕ1Xt−1+ϕ2Xt−2+et
其中 et∼N(0,σ2e),et 独立于 Xt−1,Xt−2,⋯,并假设序列具有零均值。
1.1. 估计 ϕ
首先计算样本的自相关 r1 和 r2,由方程
[r1r2]=[1r1r11][ˆϕ1ˆϕ2]
可以解得参数 ˆϕ1 和 ˆϕ2。
1.2. 估计方差
对式 (1) 等号两边计算方差,得
Var(Xt)=ϕ21Var(Xt−1)+ϕ22Var(Xt−2)+2ϕ1ϕ2Cov(Xt−1,Xt−2)+σ2e
可以解得
σ2e=γ0–ϕ21γ1+ϕ22γ2–2ϕ1ϕ2γ1=γ0(1–ϕ21–ϕ22–2ϕ1ϕ2ρ1)
由 Yule-Walker 方程组
[ρ1ρ2]=[1ρ1ρ11][ϕ1ϕ2]
即
ρ1=ϕ1+ρ1ϕ2ρ2=ρ1ϕ1+ϕ2
于是
1–ϕ21–ϕ22–2ϕ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(1–p∑i=1ˆϕiri)
1.3. 代码示例
首先模拟 AR(2) 过程
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
的值为 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
的值为
1.0000000 0.6701659 0.6701659 1.0000000
b
的值为
0.6701659 0.7319296
解得 phi.hat
的值为
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)
解得 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)
然后进行估计
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
sd.hat
为 4.02255420924941
,可见估计结果与实际值都很接近。