20 MA模型的参数估计
20.1 MA(1)模型的参数估计
考虑MA(1) \[\begin{aligned} X_t = \varepsilon_t + b \varepsilon_{t-1}, \quad t=1,2,\dots \end{aligned}\] 其中\(|b|<1\)。 \(\gamma_0=E X_t^2 = \sigma^2(1 + b^2)\), \(\gamma_1 = \sigma^2 b\). \(\rho_1 = \frac{b}{1+b^2}\), 关于\(b\)的方程为 \[\begin{aligned} \rho_1 b^2 - b + \rho_1 = 0 \end{aligned}\] 当\(|\rho_1|\leq 0.5\)时\(b\)才有实解: \[\begin{aligned} b = \frac{1 - \sqrt{1 - 4\rho_1^2}}{2\rho_1} \end{aligned}\] (另一个解使\(|b|>1\)，抛弃) 用\(\hat\rho_1 = \hat\gamma_1 / \hat\gamma_0\)代替理论值 可得\(b\)的矩估计 \[\begin{aligned} \hat b = \frac{1 - \sqrt{1 - 4\hat\rho_1^2}}{2\hat\rho_1} \end{aligned}\]
如果\(\{\varepsilon_t\}\)是独立同分布的白噪声, 由§16.2定理16.1知道\(\hat\rho_1\)是\(\rho_1\)的强相合估计: \(\lim_{ N\to \infty} \hat \rho_1=\rho_1\) , a.s., 于是 \[ \lim_{N\to \infty} \hat b = \frac{1-\sqrt{1-4 \rho^2_1} }{2 \rho_1} =b, \ a.s., \] 所以\(\hat b\)是\(b\)的强相合估计. 实际上利用§16.2定理16.2还可以证明, 当\(N \to \infty\)时, \(\sqrt{N}( \hat b - b)\) 依分布收敛到正态分布([27]). \[ \text{N}\left(0, \frac{1+b^2+b^4+b^6+b^8}{(1- b^2)^2}\right). \]
20.2 MA模型的矩估计及其计算
考虑可逆MA(\(q\))模型: \[\begin{align} X_t = \varepsilon_t + \sum_{j=1}^q b_j \varepsilon_{t-j}, \quad t \in \mathbb Z, \tag{20.1} \end{align}\] \(\{\varepsilon_t\}\sim\text{WN}(0,\sigma^2)\), 系数满足可逆条件： \[\begin{align} B(z) = 1 + \sum_{j=1}^q b_j z^j \neq 0, \quad |z|\leq 1, \tag{20.2} \end{align}\] 假定\(q\)已知，估计\(\boldsymbol{b} = (b_1,\dots, b_q)^T\)和\(\sigma^2\)。
参数与自协方差函数关系为 \[\begin{align} \gamma_k = \sigma^2(b_0 b_k + b_1 b_{k+1} + \dots + b_{q-k}b_q), \quad 0 \leq k \leq q. \tag{20.3} \end{align}\] (其中\(b_0=1\))
从样本估计出\(\hat\gamma_0, \hat\gamma_1, \dots, \hat\gamma_q\)后， 可以用解非线性方程组的方法求解\(\boldsymbol{b}\)和\(\sigma^2\)， 但不能保证解唯一，也不能保证可逆性条件。 求解方法包括线性迭代方法和Newton-Raphson迭代方法。
还可以根据§12.8计算矩估计。 由MA(\(q\))的\(q\)步截尾性， 可定义 \[\begin{aligned} \tilde \gamma_k = \begin{cases} \hat\gamma_k, & 0 \leq k \leq q,\\ 0, & k > q, \end{cases} \end{aligned}\] 用\(\{\tilde\gamma_k\}\)作为\(\{\gamma_k\}\)的估计代入§12.8的计算公式。 记 \[\begin{align} \begin{array}{llr} A=\left( \begin{array}{llllll} 0 & 1 & 0 & \cdots & 0 & 0 \\ 0 & 0 & 1 & \cdots & 0 & 0 \\ \cdots & \cdots & & \cdots & \cdots & \\ 0 & 0 & 0 & \cdots & 0 & 1 \\ 0 & 0 & 0 & \cdots & 0 & 0 \end{array} \right )_{q\times q}, & C=\left( \begin{array}{c} 1 \\ 0 \\ \vdots \\0 \end{array} \right)_{q \times 1} \\ \\ \Omega_k=\left( \begin{array}{llll} \tilde\gamma_1 & \tilde\gamma_2 & \cdots & \tilde\gamma_k \\ \tilde\gamma_2 & \tilde\gamma_3 & \cdots & \tilde\gamma_{k+1} \\ \cdots & \cdots & \cdots & \cdots \\ \tilde\gamma_q & \tilde\gamma_{q+1} & \cdots & \tilde\gamma_{q+k-1} \end{array} \right), & \boldsymbol{\gamma}_q=\left( \begin{array}{c} \tilde\gamma_1 \\ \tilde\gamma_2 \\ \vdots \\ \tilde\gamma_q \end{array} \right) & \end{array} \tag{20.4} \end{align}\] 矩估计为 \[\begin{align} (\hat b_1, \dots, \hat b_q)^T = \frac{1}{\hat\sigma^2}(\boldsymbol\gamma_q - A \hat\Pi C), \quad \hat\sigma^2 = \hat\gamma_0 - C^T \hat\Pi C \tag{20.5} \end{align}\] 其中 \[ \hat\Pi = \lim_{k\to\infty} \hat\Omega_k \tilde\Gamma_k^{-1} \hat\Omega_k \]
定理20.1 如果模型(20.1)中的\(\{\varepsilon_t\}\)是独立同分布的\(\text{WN}(0,\sigma^2)\), 则几乎必然地当\(N\)充分大后由(20.5)计算的 \(\hat b_1, \hat b_2,\dots,\hat b_q\)满足可逆条件(20.2).
证明 由于当\(N\to \infty\),\(\hat \gamma_k\to \gamma_k\) , a.s., 所以 \[\begin{aligned} \hat f(\lambda) =& \frac{1}{2\pi}\sum_{k=-q}^{q} \hat \gamma_k e^{-ik\lambda} \\ \to& f(\lambda) = \frac{1}{2\pi}\sum_{k=-q}^{q} \gamma_k e^{-ik\lambda}, \ \text{a.s.}, \ \hbox{当}\ N\to \infty \end{aligned}\] 在\([-\pi, \pi]\)上一致成立. 于是当\(N\)充分大后\(\hat f(\lambda)\)恒正. 利用§12.6的引理12.1知道有惟一的\((b_1',b_2',\dots,b_q')\)满足可逆性条件(20.2), 并且使得 \[ \hat f(\lambda) =\frac{\sigma^2_0}{2\pi} |1+\sum_{j=1}^q b_j'e^{-ij\lambda}|^2 \] 是 \[ Y_t = e_t+\sum_{j=1}^q b_j'e_{t-j}, \ \ \{e_t\} \sim \text{WN}(0,\sigma^2_0) \] 的谱密度. 这时, \(\tilde \gamma_k = E(Y_t Y_{t+k})\). 再利用§12.8知道 \[ (\hat b_1, \hat b_2,\dots,\hat b_q) =(b_1',b_2',\cdots,b_q'). \]
○○○○○○
从定理20.1的证明知道这样得到的 \((\hat b_1, \hat b_2,\cdots,\hat b_q)\)满足 \[ 1+\sum_{j=1}^q\hat b_j z^j \neq 0, \ |z|\leq 1 \] 的充分条件是 \[ \sum_{k=-q}^{q} \hat \gamma_k e^{-ik\lambda}>0, \ \lambda \in [-\pi, \pi]. \]
20.2.1 矩估计法的R程序
# Given \gamma_0, \gamma_1, \dots, \gamma_q,
## Solve MA(q) coefficients b_1, \dots, b_q, \sigma^2
## Using Li Lei's algorithm.
## Input: gms -- \gamma_0, \gamma_1, \dots, \gamma_q
ma.solve <- function(gms, k=100){
q <- length(gms)-1
if(q==1){
rho1 <- gms[2] / gms[1]
b <- (1 - sqrt(1 - 4*rho1^2))/(2*rho1)
s2 <- gms[1] / (1 + b^2)
return(list(b=b, s2=s2))
}
A <- matrix(0, nrow=q, ncol=q)
for(j in seq(2,q)){
A[j-1,j] <- 1
}
cc <- numeric(q); cc[1] <- 1
gamma0 <- gms[1]
gammas <- numeric(q+k)
gammas[1:(q+1)] <- gms
gamq <- gms[-1]
Gammak <- matrix(0, nrow=k, ncol=k)
for(ii in seq(k)){
for(jj in seq(k)){
Gammak[ii,jj] <- gammas[abs(ii-jj)+1]
}
}
Omk <- matrix(0, nrow=q, ncol=k)
for(ii in seq(q)){
for(jj in seq(k)){
Omk[ii,jj] <- gammas[ii+jj-1+1]
}
}
PI <- Omk %*% solve(Gammak, t(Omk))
s2 <- gamma0 - c(t(cc) %*% PI %*% cc)
b <- 1/s2 * c(gamq - A %*% PI %*% cc)
return(list(b=b, s2=s2))
}例20.1 考虑如下MA(2)模型： \[ X_t = \varepsilon_t - 0.36 \varepsilon_{t-1} + 0.85 \varepsilon_{t-2}, \ \varepsilon_t \sim \text{WN}(0, 2^2) \] 对样本量\(N=100, 300\)重复产生样本\(M=400\)组， 进行参数估计， 计算估计参数的平均值、标准差、根均方误差。
demo.ma2.mom <- function(m=400){
b <- c(-0.36, 0.85)
sig <- 2.0
ests.100 <- matrix(0, nrow=m, ncol=3)
ests.300 <- matrix(0, nrow=m, ncol=3)
for(ii in seq(m)){
x <- arima.sim(
model=list(ma=b), n=300,
rand.gen=function(n, ...) rnorm(n, 0, sig))
## 样本量100的结果
gms1 <- c(acf(x[1:100], type="covariance", lag.max=2, plot=FALSE)$acf)
res1 <- ma.solve(gms1)
## 样本量300的结果
gms2 <- c(acf(x[1:300], type="covariance", lag.max=2, plot=FALSE)$acf)
res2 <- ma.solve(gms2)
ests.100[ii,] <- c(res1$b, res1$s2)
ests.300[ii,] <- c(res2$b, res2$s2)
}
res <- matrix(0, nrow=7, ncol=3)
rownames(res) <- c("真实参数", "样本量100估计均值", "样本量300估计均值",
"样本量100估计标准差", "样本量300估计标准差",
"样本量100估计根均方误差", "样本量300估计根均方误差"
)
colnames(res) <- c("b1", "b2", "s2")
res[1,] <- c(b, sig^2)
res[2,] <- apply(ests.100, 2, mean)
res[3,] <- apply(ests.300, 2, mean)
res[4,] <- apply(ests.100, 2, sd)
res[5,] <- apply(ests.300, 2, sd)
res[6,] <- sqrt(1/m*( (m-1)*res[4,]^2 + m*(res[2,]-res[1,])^2))
res[7,] <- sqrt(1/m*( (m-1)*res[5,]^2 + m*(res[3,]-res[1,])^2))
print(round(res, 4))
invisible(res)
}
demo.ma2.mom(m=400)## b1 b2 s2
## 真实参数 -0.3600 0.8500 4.0000
## 样本量100估计均值 -0.4583 1.0040 6.4688
## 样本量300估计均值 -0.2798 0.3969 4.3447
## 样本量100估计标准差 2.1513 6.0264 31.4712
## 样本量300估计标准差 1.6003 5.0242 10.3337
## 样本量100估计根均方误差 2.1509 6.0209 31.5286
## 样本量300估计根均方误差 1.6003 5.0384 10.3265可以看出矩估计法的误差太大， 在中小样本情形无法接受。 经检查， 如果输入的是理论自协方差函数， 求解的结果是完全正确的， 所以并不是程序问题， 而是矩估计法的问题。 极大似然估计可以获得较满意的估计精度。
20.3 MA模型参数估计的逆相关函数法
因为AR的Yule-Warker估计能保证最小相位性所以想到把MA变成一个AR再估计。
定义20.1 设平稳序列\(\{X_t\}\)有恒正的谱密度\(f(\lambda)\). 通常称 \[\begin{align} f_y(\lambda) =\frac{1}{4\pi^2 f(\lambda)}. \tag{20.6} \end{align}\] 为\(\{X_t\}\)的逆谱密度. 称 \[ \gamma_y(k) \stackrel{\triangle}{=} \int_{-\pi}^{\pi} e^{ik\lambda} f_y(\lambda) \; d \lambda \] 为\(\{X_t\}\)的逆相关函数或逆自相关函数. 严格来说应该叫做逆自协方差函数。
设\(\{X_t\}\)为可逆MA(\(q\))序列。其谱密度为 \[\begin{aligned} f(\lambda) = \frac{\sigma^2}{2\pi} |B(e^{i\lambda})|^2 \end{aligned}\] 逆谱密度为 \[\begin{align} f_y(\lambda) = \frac{\sigma^{-2}}{2\pi} |B(e^{i\lambda})|^{-2} \tag{20.7} \end{align}\] \(f_y\)为如下AR(\(q\))序列的谱密度 \[\begin{align} Y_t = -\sum_{j=1}^q b_j Y_{t-j} + e_t, \quad t\in \mathbb Z \tag{20.8} \end{align}\] 其中\(\{e_t\}\sim\text{WN}(0,\sigma^{-2})\)。
\(\{Y_t\}\)的自协方差函数 \[\begin{align} \gamma_y(k) = \int_{-\pi}^\pi e^{ik\lambda} f_y(\lambda)\;d\lambda, \quad k=0,1,2,\dots \tag{20.9} \end{align}\] 是\(\{X_t\}\)的逆相关函数。
只要能估计\(\{\gamma_y(k)\}\)，设估计为\(\{\hat\gamma_y(k) \}\)， 就可以用Y-W方法估计\(\{Y_t\}\)的参数\(b_1,\dots,b_q\)和\(\sigma^{-2}\)： \[\begin{align} & \left(\begin{array}{cccc} \hat\gamma_y(0) & \hat\gamma_y(1) & \cdots & \hat\gamma_y(q-1) \\ \hat\gamma_y(1) & \hat\gamma_y(0) & \cdots & \hat\gamma_y(q-2) \\ \vdots & \vdots & \ddots & \vdots \\ \hat\gamma_y(q-1) & \hat\gamma_y(q-2) & \cdots & \hat\gamma_y(0) \end{array}\right) \left(\begin{array}{c} -\hat b_1 \\ -\hat b_2 \\ \vdots \\ -\hat b_q \end{array}\right) = \left(\begin{array}{c} \hat\gamma_y(1) \\ \hat\gamma_y(2) \\ \vdots \\ \hat\gamma_y(q) \end{array}\right) \tag{20.10} \\ & \hat\sigma^{-2} = \hat\gamma_y(0) + \hat b_1 \hat\gamma_y(1) + \dots + \hat b_q \hat\gamma_y(q) \tag{20.11} \end{align}\]
从而得到MA(\(q\)) 序列\(\{X_t\}\)的参数估计， 且系数满足可逆性条件。
引理20.1 如果 \((a_1,a_2,\dots,a_p)\) 和 \(\sigma^2\)分别是AR\((p)\)模型 \[ X_t=\sum_{j=1}^p a_j X_{t-j} + \varepsilon_t, \ \ t\in \mathbb Z, \] 的自回归系数和白噪声\(\{\varepsilon_t\}\)的方差, 则 \(\{X_t\}\)的逆相关函数为 \[ \gamma_y(k) = \begin{cases} \frac{1}{\sigma^2} \sum_{j=0}^{p-k}a_j a_{j+k}, &0 \leq k \leq p, \ \ a_0\stackrel{\triangle}{=} -1, \\ 0, \ &k>p. \end{cases} \]
证明 \(\{X_t\}\)有谱密度 \[ f(\lambda) =\frac{\sigma^2}{2\pi}|A(e^{i\lambda})|^{-2}, \ \hbox{ 其中 } \ A(z)=1-\sum_{j=1}^p a_j z^j, \] 和逆谱密度 \[ f_y(\lambda)=\frac{1}{ 4\pi^2 f(\lambda)} = \frac{1}{2\pi\sigma^2}|A(e^{i\lambda})|^2, \] 于是有逆相关函数 \[\begin{aligned} \gamma_y(k) = \int_{-\pi}^{\pi} e^{ik\lambda} f_y(\lambda) d \lambda = \begin{cases} \frac{1}{\sigma^2} \sum_{j=0}^{p-k} a_j a_{j+k}, & 0 \leq k \leq p, \\ 0, &k>p. \end{cases} \end{aligned}\]
○○○○○○
满足(20.1)的可逆MA\((q)\)序列\(\{X_t\}\)可以写成无穷阶自回归的形式 \[\begin{align} X_t - \sum_{j=1}^{\infty} a_j X_{t-j} =\varepsilon_t, \quad t=1,2,.., \tag{20.12} \end{align}\] 这里的回归系数\(\{a_j\}\) 由 \(1/B(z)\) 在单位圆内的Taylor级数决定: \[ \frac1{B(z)} = 1 - \sum_{j=1}^\infty a_j z^j , \ |z|\leq 1. \] 由于当 \(j\to \infty\), \(a_j\)是以负指数阶收敛到\(0\)的, 所以对较大的正整数 \(p\) 可以将无穷阶自回归模型(20.12)写成近似的长阶(\(p\)阶)自回归的形式 \[\begin{align} X_t\approx \sum_{j=1}^p a_j X_{t-j} + \varepsilon_t, \ t=1,2,.., \tag{20.13} \end{align}\] 利用引理20.1知道\(\{X_t\}\)的逆相关函数\(\gamma_y(k)\)满足 \[\begin{align} \gamma_y(k) \approx \frac1{\sigma^2} \sum_{j=0}^{p-k} a_j a_{j+k}, \ 0\leq k \leq q, \ a_0\stackrel{\triangle}{=} -1. \tag{20.14} \end{align}\]
用样本\(x_1,x_2,\dots,x_N\)计算逆相关函数\(\hat\gamma_y(k)\)和\(\hat b_1,\dots,\hat b_q, \hat\sigma^2\)的方法如下。
(1) 首先利用\(\{x_t\}\)的样本自协方差函数\(\hat \gamma_k\)建立一个AR\((p_N)\)模型, 这里\(p_N\)可以是AR模型的AIC定阶, 也可以取作\(K\ln (N)\)的整数部分, \(K\)是一个正数.
(2) 对\(p\equiv p_N\)解样本Yule-Walker方程(19.3), (19.4), 得到样本Yule-Walker 系数 \[ \hat{\boldsymbol a} = (\hat a_{1}, \hat a_{2},\dots,\hat a_{p})^T\ \hbox{ 和 } \ \hat\sigma^2_p. \]
(3) 计算样本逆相关函数 \[ \hat \gamma_y(k) = \frac{1}{\hat \sigma^2_p} \sum _{j=0}^{p-k} \hat a_{j} \hat a_{j+k}, \ k=0,1,2,\dots,q, \ \hat a_{0} \stackrel{\triangle}{=} -1. \]
(4) 利用样本Yule-Walker方程(20.10)和(20.11)计算出 MA\((q)\)系数的估计 \(\boldsymbol{\hat b}_q=(\hat b_1, \hat b_2,..,\hat b_q)^T\) 和 \(\hat \sigma^2\).
这种办法是用一个长阶自回归来近似\(\{X_t\}\)， 自回归序列的逆自相关函数很容易计算，这样得到逆相关函数。
逆相关函数法估计MA参数的R程序：
## 用输入的$\gamma_0, \dots, \gamma_p$求解AR(p)模型参数
ar.yw <- function(gam){
p <- length(gam) - 1
G <- matrix(0, p, p)
for(i in 1:p){
for(j in 1:p){
G[i,j] = gam[1 + abs(i-j)]
}
}
gri <- gam[-1]
a <- solve(G, gri)
sig2 <- sum(gam * c(1, -a))
list(ar = a, var.pred = sig2)
}
## 从AR模型参数计算逆相关函数
ar_racv <- function(a, sig2){
p <- length(a)
a <- c(-1, a)
a2 <- c(a, numeric(p+1))
res <- numeric(p+1)
for(k in seq(0, p, by=1)){
res[k+1] <- sum(a * a2[(k+1):(k+p+1)])
}
res/sig2
}
## 用长阶自回归方法和逆相关函数方法估计MA模型
ma.solve.racv <- function(x, q=1){
## 拟合长阶自回归，用来估计逆自协方差函数
n <- length(x)
plar <- round(sqrt(n))
if(plar < q) plar <- q
mod1 <- ar(x, aic = FALSE, order.max = plar, method="yule-walker")
a <- mod1$ar
sigsq <- mod1$var.pred
gamr <- ar_racv(a, sigsq)[1:(q+1)]
res2 <- ar.yw(gamr)
b <- -res2$ar
sig2 <- 1/res2$var.pred
list(ma = b, var.pred = sig2)
}20.4 MA模型参数估计的新息方法
用新息预报公式可以计算\(\{X_t\}\)的样本新息。 \[\begin{aligned} \hat\varepsilon_{t+1} =& X_{t+1} - \hat X_{t+1} \stackrel{\triangle}{=} X_{t+1} - L(X_{t+1} | X_1,\dots,X_{t}) \\ =& X_{t+1} - L(X_{t+1} | \hat\varepsilon_1,\dots,\hat\varepsilon_{t}) \\ =& X_{t+1} - \sum_{j=1}^t \theta_{t,j} \hat\varepsilon_{t+1-j}, \quad t=1,2,\dots \end{aligned}\] 其中\(\hat X_1 \stackrel{\triangle}{=} 0\)，\(\{\theta_{t,j}\}\)可递推计算， 预报误差\(\nu_t = E \hat\varepsilon_{t+1}^2\)可递推计算。 对MA(\(q\))序列\(\{X_t\}\)上述新息预报公式中当\(t\geq q, j>q\)时\(\theta_{t,j}=0\), 新息预报公式变成 \[\begin{align} \hat X_{t+1} = \sum_{j=1}^q \theta_{t,j} \hat\varepsilon_{t+1-j}, \quad t\geq q \tag{20.15} \end{align}\]
对MA(\(q\))序列\(\{X_t\}\)，其白噪声项\(\{\varepsilon_t\}\)是新息, \(t\to\infty\)时: \[\begin{aligned} \varepsilon_{t+1} =& X_{t+1} - L(X_{t+1} | X_t, X_{t-1}, X_{t-2}, \dots) \\ \approx& X_{t+1} - L(X_{t+1} | X_t, X_{t-1}, X_{t-2}, \dots, X_1) \\ =& \hat\varepsilon_{t+1} \end{aligned}\] 这是因为§23.2.2定理23.3， 无穷长历史最佳线性预测是有限历史预测的极限。 于是\(t\)较大时 \[\begin{aligned} X_t =& \varepsilon_t + b_1 \varepsilon_{t-1} + \dots + b_q \varepsilon_{t-q} \\ \approx& \hat\varepsilon_t + b_1 \hat\varepsilon_{t-1} + \dots + b_q \hat\varepsilon_{t-q} \\ =& X_t - \hat X_t + b_1 \hat\varepsilon_{t-1} + \dots + b_q \hat\varepsilon_{t-q} \end{aligned}\] 说明 \[\begin{aligned} \hat X_t \approx b_1 \hat\varepsilon_{t-1} + \dots + b_q \hat\varepsilon_{t-q} \end{aligned}\] 上式与(20.15))的新息预报公式比较可知当\(t\)较大时 可以用\(\hat b_j = \theta_{t,j}\)估计\(b_j\)， 用\(\nu_t\)估计\(\sigma^2\)。 这种估计称为新息估计。
MA参数的新息估计算法如下:
- 给定观测数据\(x_1,x_2,\dots,x_N\), 取\(m=o(N^{1/3})\).
- 计算样本自协方差函数 \(\hat\gamma_0, \hat \gamma_1, \cdots, \hat \gamma_m\).
- \(\boldsymbol{b}\) 和 \(\sigma^2\)的新息估计 \[\begin{align} (\hat b_1, \hat b_2,\dots,\hat b_q) = (\hat \theta_{m,1}, \hat \theta_{m,2},\dots,\hat \theta_{m,q}), \ \ \hat \sigma^2 = \hat \nu_m. \tag{20.16} \end{align}\] 由下面的新息预报递推公式得到. \[\begin{align} \begin{cases} \hat \nu_0=\hat \gamma_0, &\\ \hat \theta_{n,n-k} = \hat \nu_k^{-1}[\hat \gamma_{n-k} - \sum_{j=0}^{k-1} \hat \theta_{k,k-j} \hat \theta_{n,n-j}\hat \nu_j], & 0\leq k \leq n-1,\\ \hat \nu_n=\hat \gamma_0-\sum_{j=0}^{n-1} \hat \theta^2_{n,n-j}\hat \nu_j, & 1\leq n \leq m, \end{cases} \tag{20.17} \end{align}\] 其中\(\sum_{j=0}^{-1}(\cdot)\stackrel{\triangle}{=} 0\).
- 递推次序是 \[ \hat \nu_0; \ \hat \theta_{1,1}, \hat \nu_1; \ \hat \theta_{2,2}, \hat \theta_{2,1}, \hat \nu_2; \hat \theta_{3,3},\hat \theta_{3,2}, \hat \theta_{3,1}, \hat \nu_3; \cdots. \]
下面的定理与推论给出了新息估计方法的相合性。
定理20.2 设\(\{X_t\}\)是可逆ARMA(\(p,q\))序列, 满足 \[ A(\mathscr B)X_t=B(\mathscr B)\varepsilon_{t}, \ \ t\in \mathbb Z. \] \(\{\psi_j\}\)是\(B(z)/A(z)\)的Taylor级数系数. 如果\(\{\varepsilon_t\}\)是4阶矩有限的独立同分布\(\text{WN}(0,\sigma^2)\), 正整数列\(m=m(N)<N\) 满足当\(N\to \infty\) 时, \(m\to \infty\) 和 \(m=o(N^{1/3})\). 则对任何正整数\(q\), 当\(N\to \infty\) \[ \sqrt{N}(\hat \theta_{m,1}- \psi_1, \hat \theta_{m,2}- \psi_2, \dots,\hat\theta_{m,q} - \psi_q) \] 依分布收敛到\(q\)维正态分布\(N(0,A)\), 其中\(q\times q\)矩阵 \[ A=(a_{i,j}), \ \ a_{i,j}= \sum_{k=1}^{\min(i,j)}\psi_{i-k}\psi_{j-k}. \] 并且\(\hat \nu_m\) 依概率收敛到\(\sigma^2\).
推论20.1 对于MA(\(q\))序列\(\{X_t\}\), 当模型(20.1)中的白噪声是4阶矩有限的独立同分布序列时, 新息估计(20.16)是相合估计, 当\(N\)趋于无穷时, \[ \sqrt{N}(\hat b_{1}- b_1, \hat b_2 - b_2,\cdots,\hat b_q - b_q) \] 依分布收敛到\(q\)维正态分布\(N(0,A)\), 其中\(q\times q\)矩阵 \[ A=(a_{i,j}), \ \ a_{i,j}= \sum_{k=1}^{\min(i,j)}b_{i-k}b_{j-k}, \ \ b_0\stackrel{\triangle}{=} 1. \] 并且\(\hat \nu_m\) 依概率收敛到\(\sigma^2\).
注意需要\(m\to\infty\)时新息方法得到的参数估计才是相合的， 不能用\(\theta_{q,1},\dots,\theta_{q,q}\)作为\(b_1,\dots,b_q\)的估计。 实际中，\(m\)也不能取得太大，因为新息预测的递推公式中用到\(\hat\gamma_{m-1}\)。
20.5 MA模型的定阶方法
由于MA(\(q\))序列的特征是自相关系数\(q\)后截尾, 所以当样本自相关系数 \(\hat\rho_k = \hat \gamma_k/ \hat \gamma_0\) 从某一点\(\hat q\)后变得很小时, 可以\(\hat q\)作为\(q\)的估计. 从§16.3的定理16.2及§16.3.1知道, 对于\(m>q\), \(\sqrt{N} \hat \rho_m\)依分布收敛到期望为0, 方差为 \[ 1 + 2\rho_1^2 + 2\rho_2^2+ \dots + 2\rho_q^2 \] 的正态分布. 这样由样本自相关系数\(\hat \rho_k, k=1,2,\cdots,m\)的图形可以大致得到\(q\)的估计, 同时也可以判断采用MA(\(q\))模型的合理与否.
还可以用AIC定阶方法: 如果根据问题 的背景或数据的特性能够判定MA(\(q\))模型阶数\(q\)的上界是\(Q_0\). 对于\(m=0,1,2,\dots,Q_0\)按前述的方法逐个拟合MA(\(m\))模型. 白噪声方差\(\sigma^2\)的估计量记做\(\hat \sigma^2_m\). 定义AIC函数 \[ \text{AIC}(m) = \ln (\hat \sigma^2_m) + 2m/N, \quad m=0,1,\dots,Q_0. \] 这里 \(N\)是样本个数. AIC(m)的最小值点\(\hat q\) (如不惟一, 应取小的) 称为MA(\(q\))模型的AIC定阶.
20.6 MA模型的拟合检验
从观测数据\(x_1,x_2,\cdots,x_N\) 得到模型的参数估计\(\hat q\), \(\hat b_1,\hat b_2,\dots,\hat b_{\hat q}\) 和\(\hat \sigma^2\)后, 取 \[\begin{aligned} \hat \varepsilon_{1-\hat q} & = \hat \varepsilon_{2-\hat q} =\dots=\hat \varepsilon_0 =0,\\ y_t =& x_t-\bar x_N,\\ \hat \varepsilon_t =& y_t - \sum_{j=1}^{\hat q} \hat b_j \hat \varepsilon_{t-j}, \ t=1,2,\dots,N. \end{aligned}\] 对\(L=O(N^{1/3})\), 如果\(\{\hat \varepsilon_t: t=L, L+1,\cdots,N\}\) 能够通过白噪声检验, 就认为模型的选择合适. 否则改变\(\hat q\)的取值, 拟合新的MA模型或改用其他的模型, 例如改用ARMA模型等.
20.7 MA谱密度估计
如果从数据得到了MA 模型的参数估计, 模型的检验也已经通过, 可以利用 \[ \hat f(\lambda) = \frac{\hat \sigma^2}{2\pi} | 1+ \sum_{j=1}^{\hat q} \hat b_j e^{ij\lambda}|^2 \] 作为所关心的平稳序列的谱密度的估计. 这是因为如果观测数据确实是MA\((q)\)序列(20.1)时, 它的谱密度是 \[ f(\lambda) = \frac{\sigma^2}{2\pi} | 1+ \sum_{j=1}^{q} b_j e^{ij\lambda}|^2. \] 不难看出, 如果\(\hat q\), \(\hat b_1, \hat b_2,\dots,\hat b_q\)和 \(\hat\sigma^2\)分别是\(q\), \(b_1,b_2,..,b_q\) 和 \(\sigma^2\)的相合估计, 则\(\hat f(\lambda)\) 是 \(f(\lambda)\)的相合估计.