13 自回归滑动平均模型
13.1 ARMA(\(p,q\))模型及其平稳解
定义13.1 (ARMA模型) 设\(\{\varepsilon_t\}\) 是 \(\text{WN}(0,\sigma^2)\), 实系数多项式\(A(z)\) 和 \(B(z)\)没有公共根, 满足 \(b_0=1,\ a_p b_q\neq 0\) 和 \[\begin{align} A(z)=& 1-\sum^{p}_{j=1} a_j z^j\neq 0 , |z| \leq 1, \\ B(z)=& \sum^{q}_{j=0} b_j z^j\neq 0 ,\; |z|< 1, \tag{13.1} \end{align}\] 就称差分方程： \[\begin{align} X_t = \sum^{p}_{j=1} a_j X_{t-j} + \sum^{q}_{j=0} b_j \varepsilon_{t-j}, \quad t \in \mathbb Z, \tag{13.2} \end{align}\] 是一个自回归滑动平均模型, 简称为ARMA(\(p,q\))模型. 称满足(13.2)的平稳序列\(\{X_t\}\)为平稳解或ARMA(\(p,q\))序列.
模型写成 \[\begin{align} A(\mathscr B) X_t = B(\mathscr B) \varepsilon_t, \quad t \in \mathbb Z \tag{13.3} \end{align}\] \(A^{-1}(z) B(z)\)在\(|z| < \rho\) 解析(\(1 < \rho < \min\{|z_j|\}\), \(\{z_j\}\)为\(A(z)\)的所有根)， 可以Taylor展开 \[\begin{align} \Psi(z) \stackrel{\triangle}{=}& A^{-1}(z) B(z) = \sum_{j=0}^\infty \psi_j z^j, \quad |z| \leq \rho \tag{13.4} \end{align}\]
易见\(\psi_j = o(\rho^{-j})\), \[ A^{-1}(\mathscr B) B(\mathscr B)\varepsilon_t = \Psi(\mathscr B)\varepsilon_t = \sum_{j=0}^\infty \psi_j \varepsilon_{t-j} \] 是线性平稳列。两边用\(A(\mathscr B)\)作用，根据7中线性滤波的逆的定理7.2知 \[\begin{aligned} A(\mathscr B) \Psi(\mathscr B)\varepsilon_t = A(\mathscr B) A^{-1}(\mathscr B) B(\mathscr B) \varepsilon_t = B(\mathscr B) \varepsilon_t \end{aligned}\] 即\(\Psi(\mathscr B)\varepsilon_t\)是ARMA(\(p,q\))模型(13.2)的解。
反之，若\(\{Y_t\}\)是(13.2)的一个平稳解， 在(13.2)两边作用\(A^{-1}(\mathscr B)\)即得 \[\begin{aligned} A^{-1}(\mathscr B) A(\mathscr B) Y_t = Y_t = A^{-1}(\mathscr B) B(\mathscr B) \varepsilon_t = \Psi(\mathscr B)\varepsilon_t \end{aligned}\] 即 \[\begin{align} X_t = A^{-1}(\mathscr B) B(\mathscr B) \varepsilon_t = \Psi(\mathscr B)\varepsilon_t \tag{13.5} \end{align}\] 是ARMA(\(p,q\))模型(13.2)的唯一平稳解。
称(13.5)中的\(\{\psi_j\}\)为\(\{X_t\}\)的Wold系数。
模型(13.2)的任意解可以写成 \[\begin{align} Y_t = X_t + \sum_{j=1}^k \sum_{l=0}^{r(j)-1} V_{l,j} t^l \rho_j^{-t} \cos(\lambda_j t - \theta_{l,j}), t \in \mathbb Z, \tag{13.6} \end{align}\] 其中\(\{X_t\}\)为平稳解(13.5)， \(z_1,z_2,\dots,z_k\)为\(A(z)\)的全体互不相同的零点, \(z_j = \rho_j e^{i\lambda_j}\) 有重数\(r(j)\). 随机变量\(V_{j,l}, \theta_{l,j}\) 由\(Y_0-X_0,Y_1-X_1,\dots,Y_{p-1}-X_{p-1}\)惟一决定.
13.2 ARMA模型的模拟生成
(13.6)中的\(\{ Y_t \}\)与\(\{X_t \}\)当\(t \to \infty\)时无限接近： \[\begin{align} |Y_t - X_t| \leq \sum_{j=1}^k \sum_{l=0}^{r(j)-1} |V_{l,j}| t^l \rho_j^{-t} , t \to \infty \tag{13.7} \end{align}\] 可以据此模拟ARMA模型: 取初值\(Y_{-(p-1)} = \dots = Y_{-1} = Y_0 = 0\)，递推得 \[\begin{aligned} Y_t = \sum_{j=1}^p a_j Y_{t-j} + \sum_{j=0}^q b_j \varepsilon_{t-j}, t=1,2,\dots, m+n \end{aligned}\] 当\(m\)较大时取后一段\(Y_t, t=m+1,m+2,\dots,m+n\)作为ARMA(\(p,q\))模型的模拟数据。 当\(A(z)\)有靠近单位圆的根时\(m\)要取得较大。
13.3 ARMA(\(p,q\))序列的自协方差函数
13.3.1 用Wold系数表示
\(\{\gamma_k\}\)可由Wold系数表示: \[\begin{align} \gamma_k = \sigma^2 \sum_{j=0}^\infty \psi_j \psi_{j+k}, \quad k=0,1,2,\dots \tag{13.8} \end{align}\] 由于\(\psi_j = o(\rho^{-j}), j\to\infty\)， 由(13.8)可得\(\gamma_k = o(\rho^{-j}), j\to\infty\)。
13.3.2 Wold系数递推公式
记\(b_j=0, j<0\)或\(j>q\), \(b_0=1\)； \(\psi_j=0, j<0\)。 由参数\(\boldsymbol{a}_p = (a_1, \dots, a_p)^T\), \(\boldsymbol{b}_p = (b_1, \dots, b_q)^T\)计算\(\{\psi_j\}\)时可以递推 \[\begin{align} \psi_j = \begin{cases} 1, \quad & j=0, \\ b_j + \sum_{k=1}^p a_k \psi_{j-k}, & j=1,2,\dots \end{cases} \tag{13.9} \end{align}\]
证明:
记\(A(z) = 1 - \sum_{j=1}^p a_j z^j = \sum_{j=0}^p \phi_j z^j\)。 注意 \[\begin{aligned} A(z)\Psi(z) =& \sum_{k=0}^p \phi_k z^k \sum_{j=0}^\infty \psi_j z^j \\ =& \sum_{j=0}^\infty \sum_{k=0}^p \phi_k \psi_{j-k} z^j = B(z) \end{aligned}\] 比较系数得 \[\begin{aligned} \sum_{k=0}^p \phi_k \psi_{j-k} =& b_j, \quad j \geq 1 \\ \psi_j =& \sum_{k=1}^p a_k \psi_{j-k} + b_j, \quad j \geq 1 \end{aligned}\] 即(13.9)成立。
○○○○○○
13.4 ARMA(\(p,q\))模型的可识别性
我们将证明: 由ARMA(\(p,q\))模型的自协方差函数\(\{\gamma_k\}\)可以决定ARMA(\(p,q\)) 模型的参数 \[\begin{aligned} (\boldsymbol{a}^T, \boldsymbol{b}^T, \sigma^2) = (a_1, \dots, a_p, b_1, \dots, b_q, \sigma^2) \end{aligned}\]
引理13.1 设\(\{X_t\}\)是(13.2)的平稳解. 如果又有白噪声\(\{\eta_t\}\)和实系数多项式\(C(\mathscr B)\),\(D(\mathscr B)\)使得 \[C(\mathscr B)X_t = D(\mathscr B)\eta_t, \ \ t\in \mathbb Z,\] 成立, 则\(C(z)\)的阶数\(\geq p\), \(D(z)\)的阶数\(\geq q\).
这主要因为我们要求多项式\(A(z)\)和\(B(z)\)互素。 证明略。
13.4.1 ARMA序列的Y-W方程
ARMA模型的平稳解为 \[\begin{aligned} X_t = \sum_{j=0}^\infty \psi_j \varepsilon_{t-j} \end{aligned}\] 所以 \[\begin{aligned} E (\varepsilon_{t+k} X_t) = 0, \quad k>0 \end{aligned}\]
类似AR模型可推导ARMA模型的Y-W方程: 在模型方程两边同乘以\(X_{t-k}\)求期望得 \[\begin{aligned} E(X_t X_{t-k}) =& \sum_{j=1}^p a_j E(X_{t-j} X_{t-k}) + \sum_{j=0}^q b_j E(\varepsilon_{t-j} X_{t-k})\\ \text{即} & \\ \gamma_k =& \sum_{j=1}^p a_j \gamma_{k-j} + \sum_{j=0}^q b_j E(\varepsilon_{t-j} \sum_{l=0}^\infty \psi_l \varepsilon_{t-k-l}) \\ =& \sum_{j=1}^p a_j \gamma_{k-j} + \sum_{j=0}^q b_j \psi_{j-k} \sigma^2, \quad k \in \mathbb Z \end{aligned}\]
当\(k>q\)时\(\psi_{j-k}=0, j=0,1,\dots,q\)，上式为 \[\begin{aligned} \gamma_k =& \sum_{j=1}^p a_j \gamma_{k-j}, \quad k \geq q+1 \end{aligned}\]
总之 \[\begin{align} \gamma_k - \sum_{j=1}^p a_j \gamma_{k-j} = \begin{cases} \sigma^2 \sum_{j=\max(0,k)}^q b_j \psi_{j-k}, \ & k<q \\ \sigma^2 b_q, & k=q \\ 0 & k>q \end{cases} \tag{13.10} \end{align}\]
对\(k>q\)的Y-W方程可以写成矩阵形式： \[\begin{align} \left[ \begin{array}{c} \gamma_{q+1} \\ \gamma_{q+2} \\ \vdots \\ \gamma_{q+p} \end{array} \right] = \left[ \begin{array}{cccc} \gamma_{q} & \gamma_{q-1} & \cdots & \gamma_{q-p+1} \\ \gamma_{q+1} & \gamma_{q} & \cdots & \gamma_{q-p+2} \\ \vdots & \vdots & & \vdots \\ \gamma_{q+p-1} & \gamma_{q+p-2} & \cdots & \gamma_{q} \\ \end{array} \right] \left[ \begin{array}{c} a_1 \\ a_2 \\ \vdots \\ a_p \end{array} \right] \tag{13.11} \end{align}\]
对比AR(\(p\))的Y-W方程，相当于\(\Gamma_p\)的\((i,j)\)元素写成\(\gamma_{i-j}\)后 给所有\(\gamma_{\cdot}\)的下标都加上\(q\)。
把系数矩阵记为\(\Gamma_{p,q}\): \[\begin{aligned} \Gamma_{p,q} =&(\gamma_{|q+i-j|})_{i,j=1,2,\dots,p} \\ =& \left[ \begin{array}{cccc} \gamma_{q} & \gamma_{q-1} & \cdots & \gamma_{q-p+1} \\ \gamma_{q+1} & \gamma_{q} & \cdots & \gamma_{q-p+2} \\ \vdots & \vdots & & \vdots \\ \gamma_{q+p-1} & \gamma_{q+p-2} & \cdots & \gamma_{q} \\ \end{array} \right] \end{aligned}\] 只要\(\Gamma_{p,q}\)可逆则可解出\(a_1,\dots, a_p\)。
解出\(a_1,\dots,a_p\)后令 \[\begin{aligned} Y_t = A(\mathscr B) X_t = B(\mathscr B)\varepsilon_t, \quad t \in \mathbb Z \end{aligned}\] 则\(\{Y_t\}\)是一个MA(\(q\))序列，其自协方差函数为\(q\)步截尾，且 \[\begin{aligned} \gamma_y(k) =& E(Y_t Y_{t-k}) \\ =& \sum_{j=0}^p \sum_{l=0}^p \phi_j \phi_l E(X_{t-j} X_{t-l}) \\ =& \sum_{j=0}^p \sum_{l=0}^p \phi_j \phi_l \gamma_{k+l-j}, \quad 0\leq k \leq q \end{aligned}\] 可以用§12的方法唯一解出\(b_1,\dots, b_q, \sigma^2\)。
于是，只要\(\Gamma_{p,q}\)可逆, 则ARMA(\(p,q\))序列的自协方差函数和 ARMA(\(p,q\))模型的参数\((\boldsymbol{a}_p^T, \boldsymbol{b}_q^T, \sigma^2)\) 相互惟一决定。
13.4.2 ARMA模型中AR部分的参数求解
如果\(\Gamma_{p,q}\)可逆则由(13.11)可以解出\(a_1, \dots, a_p\)。
定理13.2 设\(\{\gamma_k\}\)为ARMA(\(p,q\))序列\(\{X_t\}\) 的自协方差函数列,则\(m\geq p\)时\(\Gamma_{m,q}\)可逆。
证明:
用反证法然后由引理13.1导出矛盾。
设\(\Gamma_{m,q}\)(\(m\times m\)矩阵)不满秩， 则存在\(\boldsymbol{\beta}=(\beta_0,\beta_1, \dots, \beta_{m-1})^T \neq 0\)使得 \(\Gamma_{m,q} \boldsymbol{\beta}=0\)，即 \[\begin{align} \sum_{l=0}^{m-1} \beta_l \gamma_{q+k-l} = 0, \quad k=0,1,\dots, m-1 \tag{13.12} \end{align}\] 注意当\(k\geq m\)时\(q+k-l > q\)，所以这时 \(\gamma_{q+k-l} = \sum_{j=1}^p a_j \gamma_{q+k-l-j}\)， 所以取\(k=m\)有 \[\begin{aligned} \sum_{l=0}^{m-1} \beta_l \; \gamma_{q+k-l} =& \sum_{l=0}^{m-1} \beta_l \sum_{j=1}^p a_j \gamma_{q+k-l-j} \\ =& \sum_{j=1}^p a_j \sum_{l=0}^{m-1} \beta_l \; \gamma_{q+(k-j)-l} \\ =& 0 \quad (\text{由(13.12)及}0 \leq k-j = m-j \leq m-1) \end{aligned}\]
递推得上式当\(k>m\)时也成立。因此 \[\begin{aligned} \sum_{l=0}^{m-1} \beta_l \; \gamma_{q+k-l} = 0, \quad k \geq 0. \end{aligned}\]
令\(Y_t = \sum_{l=0}^{m-1} \beta_l X_{t-l}\)则\(\{Y_t\}\)是零均值平稳列， 利用 \[\begin{aligned} E(Y_t X_{t-q-k}) = \sum_{l=0}^{m-1} \beta_l \gamma_{q+k-l} = 0, \quad k \geq 0 \end{aligned}\] 可知\(\{Y_t\}\)的自协方差\(q-1\)步截尾: \[\begin{aligned} E(Y_t Y_{t-q-k}) = 0, \quad k \geq 0 \end{aligned}\] 所以\(\{Y_t\}\)是MA(\(q'\))(\(q' \leq q-1\))序列， 存在\(\{\eta_t\} \sim \text{WN}(0,s^2)\)使得 \[\begin{aligned} \sum_{l=0}^{m-1} \beta_l X_{t-l} = \sum_{j=0}^{q'} d_j \eta_{t-j} \end{aligned}\] 与引理13.1矛盾。
○○○○○○
13.4.3 ARMA模型的一个充分条件
定理13.3 设零均值平稳序列\(\{X_t\}\)有自协方差函数\(\{\gamma_k\}\). 又设实数\(a_1,a_2,\cdots,a_p\) \((a_p\neq 0)\) 使得\(A(z)=1-\sum_{j=1}^p a_j z^j\) 满足最小相位条件, 另外 \[\begin{align} \gamma_k - \sum^{p}_{j=1} a_j\gamma_{k-j}=\begin{cases} c \neq 0 , \ & k=q, \\ 0, & k>q, \end{cases} \tag{13.13} \end{align}\] 则\(\{X_t\}\)是一个ARMA\((p',q')\)序列, 其中\(p'\leq p , \ q' \leq q\). \
证明: 设\(Y_t=A(\mathscr B) X_t = X_t-\sum^{p}_{j=1}a_j X_{t-j}\). 则\(\{Y_t\}\)是零均值平稳序列, 满足 \[ E(Y_t X_{t-k})=\gamma_k -\sum^{p}_{j=1} a_j \gamma_{k-j} =\begin{cases} c \neq 0, & k=q ,\\ 0, & k>q. \end{cases} \] 所以有 \[\begin{aligned} \gamma_y(k) =& E(Y_t Y_{t-k}) =E\big [ Y_t (X_{t-k} - \sum_{j=1}^p a_j X_{t-k-j}) \big]\\ =& \begin{cases} c \neq 0, \ & k=q ,\\ 0, & k>q. \end{cases} \end{aligned}\] 说明\(\{Y_t\}\)的自协方差函数是\(q\)后截尾的.
由定理12.2知道, \(\{Y_t\}\)为一个MA\((q)\)序列, 即存在单位圆内没有根的\(q\)阶实系数多项式\(B(z)\)使得\(B(0)=b_0=1\) 和 \[\begin{align} A(\mathscr B) X_{t}=Y_t = B(\mathscr B)\varepsilon_t, \quad t\in \mathbb Z, \tag{13.14} \end{align}\] 其中\(\{\varepsilon_t\}\) 是 \(\text{WN}(0,\sigma^2)\).
如果\(A(z)\)和\(B(z)\)没有公因子, 上述模型就是所需要的ARMA(\(p,q\))模型. 否则设公因子是\(C(z)\), 则有 \(A(z)=C(z)A'(z)\), \(B(z)=C(z)B'(z).\) 这时(13.14)变成 \[ C(\mathscr B)A'(\mathscr B)X_t =C(\mathscr B)B'(\mathscr B)\varepsilon_t. \] 两边乘\(C^{-1}(\mathscr B)\)(显然\(C(z)\)也满足最小相位条件) 后得到所需ARMA(\(p',q'\))模型: \[ A'(\mathscr B)X_t =B'(\mathscr B)\varepsilon_t. \]
○○○○○○
13.5 ARMA序列的谱密度
由于ARMA序列的\(\{\gamma_k\}\)绝对可和，以及平稳解的线性序列表达式， 可得ARMA(\(p,q\))序列(2.6)有谱密度 \[\begin{align} f(\lambda) =& \frac{1}{2\pi} \sum_{k=-\infty}^\infty \gamma_k e^{-ik\lambda} \\ =& \frac{\sigma^2}{2\pi} \left| \sum_{j=0}^\infty \psi_j e^{ij\lambda} \right|^2 = \frac{\sigma^2}{2\pi} \left| \frac{B(e^{i\lambda})}{A(e^{i\lambda})} \right|^2 \tag{13.15} \end{align}\] 形如(13.15)的谱密度被称为有理谱密度.
13.6 可逆ARMA序列
定义13.2 在ARMA\((p,q)\)模型的定义 13.1 中, 如果进一步要求\(B(z)\)在单位圆上无根: \[\begin{align} B(z)=1 + \sum^{q}_{j=1} b_j z^j \neq 0, |z|\leq 1 \tag{13.16} \end{align}\] 则称ARMA(\(p,q\))模型(13.16)为可逆的ARMA模型, 称相应的平稳解为可逆的ARMA(\(p,q\))序列.
从定理12.3(最小序列的谱密度条件）知道可逆的ARMA(\(p,q\))序列是最小序列.
对于可逆的ARMA(\(p,q\))模型(13.16)， 由于\(B^{-1}(z)A(z)\)在\(\{z: |z|\leq\rho\}\) \((\rho>1)\) 内解析, 所以有Taylor展式： \[\begin{align} B^{-1}(z)A(z)=\sum^{\infty}_{j=0} \varphi_j z^j,\ \ |z|\leq\rho, \tag{13.17} \end{align}\] 其中\(|\varphi_j|=o(\rho^{-j})\),当\(j\to\infty\),从而可以定义 \(B^{-1}(\mathscr B)A(\mathscr B)=\sum^{\infty}_{j=0}\varphi_j \mathscr B^j\). 在(13.17)两边乘以\(B^{-1}(\mathscr B)\), 得到： \[\begin{align} \varepsilon_t =B^{-1}(\mathscr B)A(\mathscr B)X_t =\sum^{\infty}_{j=0} \varphi_j X_{t-j}, \ \ t\in \mathbb Z. \tag{13.18} \end{align}\] (13.18)是(13.5)的逆转形式, 表明可逆ARMA(\(p,q\))序列和它的噪声序列可以相互线性表示.
○○○○○○
13.7 ARMA模型例子
先给出一些有关的自定义R函数。
## ARMA theoretical spectrum given ARMA coefficients
arma.true.spectrum <- function(a, b, ngrid=256, sigma=1,
tit="True ARMA Spectral Density",
plot.it=TRUE){
p <- length(a)
q <- length(b)
freqs <- seq(from=0, to=pi, length=ngrid)
spec1 <- numeric(ngrid)
spec2 <- numeric(ngrid)
for(ii in seq(ngrid)){
spec1[ii] <- 1 + sum(complex(mod=b, arg=freqs[ii]*seq(q)))
spec2[ii] <- 1 - sum(complex(mod=a, arg=freqs[ii]*seq(p)))
}
spec = sigma^2 / (2*pi) * abs(spec1)^2 / abs(spec2)^2
if(plot.it){
plot(freqs, spec, type='l',
main=tit,
xlab="frequency", ylab="spectrum",
axes=FALSE)
axis(2)
axis(1, at=(0:6)/6*pi,
labels=c(0, expression(pi/6),
expression(pi/3), expression(pi/2),
expression(2*pi/3), expression(5*pi/6), expression(pi)))
}
box()
invisible(list(frequencies=freqs, spectrum=spec,
ar.coefficients=a, ma.coefficients=b,
sigma=sigma))
}
## simulate ARMA(p, q) model.
## a is vector a_1, \dots, a_p
## b is vector b_1, \dots, b_q
## Model is
## X_t = a_1 X_{t-1} + \dots + a_p X_{t-p}
## + \epsilon_t + b_1 \epsilon_{t-1} + \dots + b_q \epsilon_{t-q}
## \Var(\epsilon_t) = sigma^2
arma.gen <- function(n, a, b, sigma=1.0,
by.roots.ar=FALSE, by.roots.ma=FALSE,
plot.it=FALSE, n0=1000,
x0=numeric(length(a))){
n2 <- n0 + n
p <- length(a)
## first generate n0+n MA(q) series
eps <- ma.gen(n2, b, sigma=sigma,
by.roots=by.roots.ma,
plot.it=FALSE)
b <- attr(eps, "b")
if(by.roots.ar){
require(polynom)
cf <- Re(c(poly.calc(a)))
cf <- cf / cf[1]
a <- -cf[-1]
}
##set.seed(1)
x2 <- filter(eps, a, method="recursive", side=1, init=x0)
x <- x2[(n0+1):n2]
x <- ts(x)
attr(x, "model") <- "ARMA"
attr(x, "a") <- a
attr(x, "b") <- b
attr(x, "sigma") <- sigma
if(plot.it) plot(x)
x
}
## Wold coefficients for the ARMA model
arma.Wold <- function(n, a, b=numeric(0)){
p <- length(a)
q <- length(b)
arev <- rev(a)
psi <- numeric(n)
psi[1] <- 1
for(j in seq(n-1)){
if(j <= q) bj=b[j]
else bj=0
psis <- psi[max(1, j+1-p):j]
np <- length(psis)
if(np < p) psis <- c(rep(0,p-np), psis)
psi[j+1] <- bj + sum(arev * psis)
}
psi
}
## Calculate theoretical autocovariance function
## of ARMA model using Wold expansion
arma.gamma.by.Wold <- function(n, a, b=numeric(0), sigma=1){
nn <- n + 100
psi <- arma.Wold(nn, a, b)
gam <- numeric(n)
for(ii in seq(0, n-1)){
gam[ii+1] <- sum(psi[1:(nn-ii)] * psi[(ii+1):nn])
}
gam <- (sigma^2) * gam
gam
}
arma.gamma <- arma.gamma.by.Wold
## Solve ARMA parameters given autocovariance functions
arma.solve <- function(gms, p, q){
Gpq <- matrix(0, nrow=p, ncol=p)
for(ii in seq(p)) for(jj in seq(p)){
Gpq[ii,jj] <- gms[abs(q + ii - jj)+1]
}
gs <- gms[(q+1+1):(q+p+1)]
a <- solve(Gpq, gs)
aa <- c(-1, a)
gys <- numeric(q+1)
for(k in seq(0, q)){
gys[k+1] <- sum(c(outer(0:p,0:p,
function(ii,jj) aa[ii+1] * aa[jj+1]
* gms[abs(k+jj-ii)+1])))
}
res <- ma.solve(gys)
b <- res$b
sigma <- sqrt(res$s2)
list(a=a, b=b, sigma=sigma)
}
## simulate MA(q) model.
## a is vector a_1, \dots, a_q
## Model is X_t = eps_t + a_1 eps_{t-1} + \dots + a_q eps_{t-q}
## \Var(\epsilon_t) = sigma^2
## This version uses the filter function.
ma.gen <- function(n, a, sigma=1.0, by.roots=FALSE,
plot.it=FALSE){
if(by.roots){
require(polynom)
cf <- Re(c(poly.calc(a)))
cf <- cf / cf[1]
a <- cf
}
q <- length(a)
n2 <- n+q
eps <- rnorm(n2, 0, sigma)
x2 <- filter(eps, c(1,a), method="convolution", side=1)
x <- x2[(q+1):n2]
x <- ts(x)
attr(x, "model") <- "MA"
attr(x, "b") <- a
attr(x, "sigma") <- sigma
if(plot.it) plot(x)
x
}
## 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))
}13.7.1 ARMA(4,2)例子
ARMA(4,2): \[\begin{align} \begin{aligned} a_1 = -0.9, \quad & a_2 = -1.4, \\ a_3=-0.7, \quad & a_4=-0.6; \\ b_1=0.5, \quad & b_2=-0.4. \end{aligned} \tag{13.19} \end{align}\]
\(A(z)\)的根为\(1.1380e^{\pm 2.2062i}\), \(1.1344e^{\pm 1.4896i}\), \(B(z)\)的两个实根为\(2.3252, -1.0752\)。 此时间序列有两个频率成分，对应于AR部分的特征根的辐角。
模拟生成长度\(n=80\)的样本（结果见图13.1）：
set.seed(1)
n <- 80
a <- c(-0.9, -1.4, -0.7, -0.6)
b <- c(0.5, -0.4)
x <- arma.gen(n, a, b, plot.it=F)
ts.plot(x, main="ARMA(4,2) Series")
图13.1: ARMA(4,2)模拟数据
对模拟数据估计ACF并作图（结果见图13.2）：
图13.2: ARMA(4,2)数据的ACF
模型的理论谱密度（结果见图13.3）：
图13.3: ARMA(4,2)数据的理论谱密度
模型的理论自协方差函数（结果见图13.4）：
ng <- 21
gams <- arma.gamma(ng, a, b, sigma=1)
plot(0:(ng-1), gams, type="h",
main="Theoretical gamma by Wold",
xlab="k", ylab=expression(gamma[k]))
abline(h=0)
图13.4: ARMA(4,2)数据的理论自协方差函数
##
## ====== Autocovariances ========## 0 1 2 3 4 5 6 7 8 9
## 6.6708 -1.5078 -4.5792 2.4672 1.2433 -0.4630 -0.3035 -1.4293 1.2894 1.3309
## 10 11 12 13 14 15 16 17 18 19
## -1.8203 -0.2699 1.0861 -0.1239 -0.1279 -0.3097 -0.1071 0.6939 -0.1810 -0.5477
## 20
## 0.3249下面从理论自协方差函数反解模型参数。
##
## ===== Solve ARMA from ACV ======## $a
## [1] -0.9 -1.4 -0.7 -0.6
##
## $b
## [1] 0.4999999 -0.4000000
##
## $sigma
## [1] 1用样本中估计的自协方差函数反解模型参数并与用理论值反解的结果对比：
## $a
## [1] -1.0102537 -1.5340719 -0.8871690 -0.7297601
##
## $b
## [1] 0.63530816 0.03771917
##
## $sigma
## [1] 1.358952○○○○○○
13.7.2 ARMA(2,2)例子
设\(\{X_t\} \sim\) ARMA(2,2), 已知 \[\begin{aligned} (\gamma_0,\gamma_1,\cdots,\gamma_4) =(4.61,\ -1.06,\ 0.29,\ 0.69,\ -0.12) \end{aligned}\] 要反解ARMA参数。
## $a
## [1] 0.08939301 -0.62648682
##
## $b
## [1] -0.3334025 0.8157936
##
## $sigma
## [1] 2.002966z1 <- polyroot(c(1, -res3[["a"]]))
cat("AR roots: ", Mod(z1[1]), "exp(+- i ", Arg(z1[1]), ")", "\n")## AR roots: 1.263409 exp(+- i 1.514296 )## MA roots: 1.107159 exp(+- i 1.385167 )解出的模型为 \[ X_t = 0.0894 X_{t-1} - 0.6265 X_{t-2} + \varepsilon_t - 0.3334 \varepsilon_{t-1} + 0.8158 \varepsilon_{t-2}, \ \{ \varepsilon_t \} \sim \text{WN}(0, 2.0030^2) \]
AR部分的特征根为\(1.26e^{\pm i 1.51}\), MA部分的特征根为\(1.11 e^{\pm i 1.39}\)。
○○○○○○
13.8 附录：ARAM(1,1)例子
对ARMA(1,1) \[ X_t = a X_{t-1} + \varepsilon_t + b \varepsilon_{t-1} \] 由Wold系数递推公式得 \(\psi_0=1\), \[\begin{aligned} \psi_1 =& b + a \psi_{1-1} = a + b \\ \psi_j =& a \psi_{j-1} = \cdots = a^{j-1} \psi_1 = a^{j-1}(a + b), \ j=2,3,\dots \end{aligned}\] 于是ARMA(1,1)的平稳解可以写成 \[ X_t = \varepsilon_t + (a + b)\sum_{j=1}^\infty a^{j-1} \varepsilon_{t-j} \]
从Wold系数计算协方差函数得 \[\begin{aligned} \gamma_0 =& \sigma^2 \sum_{j=0}^\infty \psi_j^2 = \sigma^2\left\{ 1 + \sum_{j=1}^\infty (a + b)^2 a^{2(j-1)} \right\} \\ =& \sigma^2\left\{ 1 + \frac{(a + b)^2}{1 - a^2} \right\} = \sigma^2 \frac{1 + 2 a b + b^2}{1 - a^2} \\ \gamma_1 =& \sigma^2 \sum_{j=0}^\infty \psi_j \psi_{j+1} = \sigma^2 \left\{ \psi_1 + \sum_{j=1}^\infty a (a + b)^2 a^{2(j-1)} \right\} \\ =& \sigma^2 \frac{(a+b)(1+ab)}{1 - a^2} = \sigma^2 \frac{a + b + a^2 b + a b^2}{1 - a^2} \\ \gamma_k =& \sigma^2 \sum_{j=0}^\infty \psi_j \psi_{j+k} = \sigma^2 a \sum_{j=0}^\infty \psi_j \psi_{j+k-1} = a \gamma_{k-1} = \cdots \\ =& a^{k-1} \gamma_1 = a^{k-1} \sigma^2 \frac{(a+b)(1+ab)}{1 - a^2} \end{aligned}\]
也可以利用YW方程计算协方差函数。 对\(k=0,1,2,\dots\)， 模型方程两边同乘以\(X_{t-k}\)取期望， 因为\(E \varepsilon_t X_{t-k}=0\)对\(k\geq 1\)， 有 \[\begin{aligned} \gamma_0 =& a \gamma_1 + E(X_t \varepsilon_t) + b E(X_t \varepsilon_{t-1}) \\ \gamma_1 =& a \gamma_0 + 0 + b E(X_{t-1} \varepsilon_{t-1}) \\ \gamma_k =& a \gamma_{k-1} + 0 + b \cdot 0, k=2,3,\dots \end{aligned}\] 由Wold表示可知\(E(X_{t} \varepsilon_{t-j}) = \sigma^2 \psi_j\)， 所以 \[\begin{aligned} \gamma_0 =& a \gamma_1 + \sigma^2[1 + b(a+b)] \\ \gamma_1 =& a \gamma_0 + \sigma^2 b \end{aligned}\] 将第二式代入第一式就可求解\(\gamma_0, \gamma_1\)为 \[\begin{aligned} \gamma_0 =& \sigma^2 \frac{1 + 2ab + b^2}{1 - a^2}, \\ \gamma_1 =& \sigma^2 \frac{a+b + a^2 b + a b^2}{1 - a^2} \\ \gamma_k =& a \gamma_{k-1} = a^{k-1} \gamma_1, k=2,3,\dots \end{aligned}\]
自相关函数为 \[\begin{aligned} \rho_1 =& \frac{(a+b)(1 + ab)}{1 + 2ab + b^2} \\ \rho_k =& a \rho_{k-1} = a^{k-1} \rho_1, k=2,3,\dots \end{aligned}\]
下面求出ARMA(1,1)的偏相关函数前几项。 \[ a_{11}=\rho_1 = \frac{(a+b)(1 + ab)}{1 + 2ab + b^2} \] \(a_{22}\)满足方程 \[ \left(\begin{matrix} \gamma_0 & \gamma_1 \\ \gamma_1 & \gamma_0 \end{matrix}\right) \left(\begin{matrix} a_{21} \\ a_{22} \end{matrix}\right) = \left(\begin{matrix} \gamma_1 \\ \gamma_2 \end{matrix}\right) \] 两边除以\(\gamma_0\)得 \[ \left(\begin{matrix} \rho_0 & \rho_1 \\ \rho_1 & \rho_0 \end{matrix}\right) \left(\begin{matrix} a_{21} \\ a_{22} \end{matrix}\right) = \left(\begin{matrix} \rho_1 \\ \rho_2 \end{matrix}\right) \] 其中\(\rho_0=1\), \(\rho_2 = a \rho_1\)，即 \[ \left(\begin{matrix} 1 & \rho_1 \\ \rho_1 & 1 \end{matrix}\right) \left(\begin{matrix} a_{21} \\ a_{22} \end{matrix}\right) = \left(\begin{matrix} \rho_1 \\ \rho_2 \end{matrix}\right) \] 由Cramer法则， \[ a_{22} = \frac{\left| \begin{matrix} 1 & \rho_1 \\ \rho_1 & a \rho_1 \end{matrix}\right|}{\left| \begin{matrix} 1 & \rho_1 \\ \rho_1 & 1 \end{matrix}\right|} = \frac{\rho_1(a - \rho_1)}{1 - \rho_1^2} \] 类似地， \[\begin{aligned} a_{33} =& \frac{\left| \begin{matrix} 1 & \rho_1 & \rho_1 \\ \rho_1 & 1 & \rho_2 \\ \rho_2 & \rho_1 & \rho_3 \end{matrix}\right|}{\left| \begin{matrix} 1 & \rho_1 & \rho_2 \\ \rho_1 & 1 & \rho_1 \\ \rho_2 & \rho_1 & 1 \end{matrix}\right|} = \frac{\rho_1 (a - \rho_1)^2}{1 + 2 a \rho_1^3 - \rho_1^2(2 + a^2)} \end{aligned}\]
○○○○○○
13.9 附录：ARMA模型的谱条件
ARMA模型的谱表示: \[\begin{aligned} X_t = \int_{-\pi}^\pi e^{it\lambda} \frac{B(e^{-i\lambda})}{A(e^{-i\lambda})} d Z_\varepsilon(\lambda) \end{aligned}\] 其中\(\frac{B(z)}{A(z)}\)叫做ARMA模型的极大解析函数。
AMRA模型的谱密度:
设\(A(\cdot)\), \(B(\cdot)\)根都在单位圆外， \(\{\varepsilon_t\}\)是WN(\(0,\sigma^2\)), 则 平稳列\(\{X_t\}\)是可逆ARMA模型 \[\begin{aligned} A(\mathscr B) X_t = B(\mathscr B) \varepsilon_t \end{aligned}\] 的充分必要条件是它有谱密度 \[\begin{aligned} f(\lambda) = \frac{\sigma^2}{2\pi} \left| \frac{B(e^{-i\lambda})}{A(e^{-i\lambda})} \right|^2 \end{aligned}\] 见(谢衷洁 1990) P.76定理2.4。
证明
必要性教材中已证明。设充分性条件成立， 这时设 \[\begin{aligned} \frac{A(z)}{B(z)} = \sum_{j=0}^\infty d_j z^j, \ |z| \leq \rho_2 \ (\rho_2>1) \end{aligned}\] 令 \[\begin{aligned} \eta_t = \sum_{j=0}^\infty d_j X_{t-j} \end{aligned}\] 则 \[\begin{aligned} A(\mathscr B) X_t = B(\mathscr B)\eta_t \end{aligned}\] 且\(\{\eta_t \}\)的谱密度为 \[\begin{aligned} f_\eta(\lambda) =& \left|\sum_{j=0}^\infty d_j e^{-ij\lambda} \right|^2 \cdot \frac{\sigma^2}{2\pi} \left| \frac{B(e^{-i\lambda})}{A(e^{-i\lambda})} \right|^2 \\ =& \left| \frac{A(e^{-i\lambda})}{B(e^{-i\lambda})} \right|^2 \cdot \frac{\sigma^2}{2\pi} \left| \frac{B(e^{-i\lambda})}{A(e^{-i\lambda})} \right|^2 \\ =& \frac{\sigma^2}{2\pi} \end{aligned}\] 即\(\{\eta_t \}\)为白噪声， 所以\(\{X_t \}\)为可逆ARMA序列。
○○○○○○
13.10 附录：ARMA模型与Hilbert空间
参考(谢衷洁 1990) P.78。
考虑可逆ARMA模型 \[\begin{aligned} A(\mathscr B) X_t = B(\mathscr B) \varepsilon_t \end{aligned}\] 设 \[\begin{aligned} H_X =& {\mathcal L}\{X_t: t \in \mathbb Z \} \\ H_\varepsilon =& {\mathcal L}\{\varepsilon_t: t \in \mathbb Z \} \\ H_X(t) =& {\mathcal L}\{X_s: s \leq t, \ s \in \mathbb Z \} \\ H_\varepsilon(t) =& {\mathcal L}\{\varepsilon_s: s \leq t, \ s \in \mathbb Z \} \end{aligned}\] 其中\(\mathcal L\)表示线性闭包为Hilbert空间。 则\(\varepsilon_t \in H_X(t)\), \(\varepsilon_t \in H_X(t) \ominus H_X(t-1)\), \(\{\varepsilon_t/\sigma, t \in \mathbb Z \}\)是\(H_X\)的一组完备标准正交基。 把\(X_t\)用标准正交基展开成Wold表示 \[\begin{aligned} X_t = \sigma \sum_{j=0}^\infty \psi_j \frac{\varepsilon_{t-j}}{\sigma} \end{aligned}\] 展开的系数\(\sigma\psi_j\)是\(H_X\)中\(X_t\)对正交基\(\{\varepsilon_t/\sigma, t \in \mathbb Z \}\) 的广义傅立叶系数 \[\begin{aligned} \sigma\psi_j = < X_t, \varepsilon_{t-j}/\sigma >, \ j=0,1,2,\dots \end{aligned}\]
可逆ARMA模型中的白噪声\(\varepsilon_t\)是新息(Wold序列)： \[\begin{aligned} \varepsilon_t = X_t - L(X_t | H_X(t-1)) \end{aligned}\] (参考(谢衷洁 1990) P.81定理2.6).
设可逆ARMA模型Wold表示为 \[\begin{aligned} X_t = \sum_{j=0}^\infty \psi_j \varepsilon_{t-j} \end{aligned}\] 逆表示为 \[\begin{aligned} \varepsilon_{t} = \sum_{j=0}^\infty d_j X_{t-j} \end{aligned}\] 则\(\{d_j \}\)可解为 \[\begin{aligned} d_0 =& 1 \\ d_j =& - \sum_{k=1}^j \psi_k d_{j-k}, \ j=1,2,\dots \end{aligned}\] 这是因为 \[\begin{aligned} \left(\sum_{k=0}^\infty d_k z^k \right) \left(\sum_{l=0}^\infty \psi_l z^l \right) = 1 \end{aligned}\] 即 \[\begin{aligned} \sum_{j=0}^\infty \left(\sum_{k=0}^j \psi_k d_{j-k} \right) z^j = 1 \end{aligned}\]