【时间序列分析】19.平稳序列递推预测实例

十九、平稳序列递推预测实例

1.$\mathrm{A}\mathrm{R}(1)$与$\mathrm{A}\mathrm{R}(2)$序列

在上一篇文章中，我们说到$\mathrm{A}\mathrm{R}(p)$序列只需要用前$p$个历史信息即可，这个结论与递推定理是否一致？现在假设$\mathrm{A}\mathrm{R}(1)$序列为
$$X_t=a{X}_{t-1}+{\epsilon }_t,{\epsilon }_t\sim \mathrm{W}\mathrm{N}(0,{\sigma }^2).$$
给定三个观测值$x_1,x_2,x_3$，基于这三个观测值预测$X_4$，应该有
$${\hat{X}}_4=a{x}_3.$$
是否如此？现在先计算自协方差函数之间的关系，显然有
$$\begin{gathered} {\gamma }_1=\mathrm{E}[(a{X}_{t-1}+{\epsilon }_t)X_{t-1}]=a{\gamma }_0, \\ {\gamma }_k=\mathrm{E}[a(X_{t-1}+{\epsilon }_t)X_{t-k}]=a{\gamma }_{k-1}. \end{gathered}$$
也就是${\gamma }_k=a^k{\gamma }_0$，现在开始递推。

第一轮递推：
$$\begin{gathered} {\nu }_0={\gamma }_0,Z_1=x_1 \\ {\theta }_{1,1}=\frac{{\gamma }_1}{{\gamma }_0}=a,{\hat{X}}_2=a{x}_1, \\ {Z}_2=X_2-{\hat{X}}_2=x_2-a{x}_1,{\nu }_1={\gamma }_0-{\theta }_{1,1}^2{\nu }_0=(1-a^2){\gamma }_0. \end{gathered}$$
第二轮递推：
$$\begin{gathered} {\theta }_{2,2}=\frac{{\gamma }_2}{{\gamma }_0}=a^2,{\theta }_{2,1}=\frac{{\gamma }_1-{\theta }_{1,1}{\theta }_{2,2}{\nu }_0}{{\nu }_1}=a, \\ {\hat{X}}_3={\theta }_{2,2}{Z}_1+{\theta }_{2,1}{Z}_2=a^2{x}_1+a(x_2-a{x}_1)=a{x}_2, \\ {Z}_3=x_3-a{x}_2,{\nu }_2={\gamma }_0-{\theta }_{2,2}^2{\nu }_0-{\theta }_{2,1}^2{\nu }_1=(1-a^2){\gamma }_0. \end{gathered}$$
第三轮递推：
$$\begin{gathered} {\theta }_{3,3}=\frac{{\gamma }_3}{{\gamma }_0}=a^3,{\theta }_{3,2}=\frac{{\gamma }_2-{\theta }_{1,1}{\theta }_{3,3}{\nu }_0}{{\nu }_1}=a^2, \\ {\theta }_{3,1}=\frac{{\gamma }_1-{\theta }_{2,2}{\theta }_{3,3}{\nu }_0-{\theta }_{2,1}{\theta }_{3,2}{\nu }_1}{{\nu }_2}=a, \\ {\hat{X}}_4={\theta }_{3,3}{Z}_1+{\theta }_{3,2}{Z}_2+{\theta }_{3,1}{Z}_3=a{x}_3, \\ {\nu }_3={\gamma }_0-{\theta }_{3,3}^2{\nu }_0-{\theta }_{3,2}^2{\nu }_1-{\theta }_{3,1}^2{\nu }_2=(1-a^2){\gamma }_0. \end{gathered}$$
可以看到使用递推得到的${\hat{X}}_4$与直接代入$\mathrm{A}\mathrm{R}(1)$模型得到的结果一样，都是$a{x}_3$。事实上，由Y-W方程与最佳线性预测的性质，可以直接得到
$$L(X_{n+1}|{\mathbf{X}}_n)={\Gamma }_n^{-1}{\gamma }_n{\mathbf{X}}_n,$$
这就说明${\mathbf{X}}_n$前面的系数是偏相关系数。从刚才的递推，我们还得到了其递推阵为
$$\begin{array}{cccc}{\nu }_0& & & \\ {\theta }_{1,1}& {\nu }_1& & \\ {\theta }_{2,2}& {\theta }_{2,1}& {\nu }_2& \\ {\theta }_{3,3}& {\theta }_{3,2}& {\theta }_{3,1}& {\nu }_3\\ \cdots & \cdots & \cdots & \cdots \end{array}=\begin{array}{cccc}{\gamma }_0& & & \\ a& (1-a^2){\gamma }_0& & \\ a^2& a& (1-a^2){\gamma }_0& \\ a^3& a^2& a& (1-a^2){\gamma }_0\\ \cdots & \cdots & \cdots & \cdots \end{array}$$
可以看到$\mathrm{A}\mathrm{R}(1)$序列的递推系数具有很好的规律，对于$\mathrm{A}\mathrm{R}(2)$序列$X_t=a_1{X}_{t-1}+a_2{X}_{t-2}+\epsilon$，则没有这么好的规律，但是$\mathrm{A}\mathrm{R}(p)$序列这种只依赖于前$p$个历史信息的预测问题，在计算的时候是十分方便的。

最后说说$\mathrm{A}\mathrm{R}(p)$模型的递推预测：$L(X_{n+k}|{\mathbf{X}}_n)$，由最佳线性预测的线性性，有
$$L(X_{n+k}|{\mathbf{X}}_n)=L\left(\sum _{j=1}^p{a}_j{X}_{n+k-j}|{\mathbf{X}}_n\right).$$
由此递推计算其$k$步预测。

2.$\mathrm{M}\mathrm{A}(1)$与$\mathrm{M}\mathrm{A}(2)$序列

$\mathrm{M}\mathrm{A}(q)$序列最重要的是它的$q$后截尾性质，这使得它的递推系数也具有类似“截尾”的效果。我们先考虑$\mathrm{M}\mathrm{A}(1)$序列：$X_t={\epsilon }_t+b{\epsilon }_{t-1}$，显然
$${\gamma }_0=(1+b^2){\sigma }^2,{\gamma }_1=b{\sigma }^2,$$
对于$k\ge 2,{\gamma }_k=0$。可惜的是，预测误差项与预测值并没有呈现什么关系，也不能由前$q$个历史信息得到未来信息（而是用前$q$个预测误差）。代入递推公式，有

第一轮递推：
$$\begin{gathered} {\nu }_0={\gamma }_0=(1+b^2){\sigma }^2, \\ {\theta }_{1,1}=\frac{{\gamma }_1}{{\nu }_0}, \\ {\nu }_1={\gamma }_0-{\theta }_{1,1}^2{\nu }_0. \end{gathered}$$
第二轮递推：
$$\begin{gathered} {\theta }_{2,2}=\frac{{\gamma }_2}{{\gamma }_0}=0,{\theta }_{2,1}=\frac{{\gamma }_1-{\theta }_{1,1}{\theta }_{2,2}{\nu }_0}{{\nu }_1}=\frac{{\gamma }_1}{{\nu }_1}, \\ {\nu }_2={\gamma }_0-{\theta }_{2,2}^2{\nu }_0-{\theta }_{2,1}^2{\nu }_1={\gamma }_0-{\theta }_{2,1}^2{\nu }_1. \end{gathered}$$

第三轮递推：
$$\begin{gathered} {\theta }_{3,3}=\frac{{\gamma }_3}{{\gamma }_0},{\theta }_{3,2}=\frac{{\gamma }_2-{\theta }_{1,1}{\theta }_{3,3}{\nu }_0}{{\nu }_1}=0,{\theta }_{3,1}=\frac{{\gamma }_1-{\theta }_{2,2}{\theta }_{3,3}{\nu }_0-{\theta }_{2,1}{\theta }_{3,2}{\nu }_1}{{\nu }_2}=\frac{{\gamma }_1}{{\nu }_2}, \\ {\nu }_3={\gamma }_0-{\theta }_{3,3}^2{\nu }_0-{\theta }_{3,2}^2{\nu }_1-{\theta }_{3,3}^2{\nu }_2. \end{gathered}$$
以此类推，我们可以得到以下的递推系数阵：
$$\begin{array}{cccc}{\nu }_0& & & \\ {\theta }_{1,1}& {\nu }_1& & \\ {\theta }_{2,2}& {\theta }_{2,1}& {\nu }_2& \\ {\theta }_{3,3}& {\theta }_{3,2}& {\theta }_{3,1}& {\nu }_3\\ \cdots & \cdots & \cdots & \cdots \end{array}=\begin{array}{cccc}{\gamma }_0& & & \\ \frac{{\gamma }_1}{{\nu }_0}& {\gamma }_0-{\theta }_{1,1}^2{\nu }_0& & \\ 0& \frac{{\gamma }_1}{{\nu }_1}& {\gamma }_0-{\theta }_{2,1}^2{\nu }_1& \\ 0& 0& \frac{{\gamma }_1}{{\nu }_2}& {\gamma }_0-{\theta }_{3,1}^2{\nu }_2\\ \cdots & \cdots & \cdots & \cdots \end{array}$$
可惜的是，我们只能得到系数的比较漂亮的递推式，但如果将其表示为$b$的函数则没有很好的形式；而且，代入$x_1,x_2$等实际观测值后，在新的预测中也没有好的形式，因为对未来的预测总会用到所有历史信息（尽管只用到一个预测误差，但这个预测误差包含所有历史信息）。总结起来，就是
$$\begin{gathered} {\nu }_0={\gamma }_0,{\theta }_{n,1}=\frac{{\gamma }_1}{{\nu }_{n-1}}; \\ {\nu }_n={\gamma }_0-{\theta }_{n,1}^2{\nu }_{n-1}, \\ \forall j\ge 2,{\theta }_{n,j}=0. \end{gathered}$$
以上结论也可以用两步数学归纳法证明。设以上结论对$n<k,k\ge 2$都成立，则
$${\theta }_{k,k}=\frac{{\gamma }_k}{{\nu }_0}=0\because k\ge 2.$$
假设${\theta }_{k,k-j}=0$对所有$j\le m<k-1$成立，则
$$\begin{gathered} {\theta }_{k,k-m}=\frac{{\gamma }_{k-m}-{\sum }_{j=0}^{m-1}{\theta }_{k,k-j}{\theta }_{m,m-j}{\nu }_j}{{\nu }_m}=0, \\ \because m<k-1,\therefore k-m>1,{\gamma }_{k-m}=0; \\ \because {\theta }_{k,k-j}=0,\forall j\le m,\therefore {\theta }_{k,k-j}=0,\forall j=0,1,\cdots ,m-1. \end{gathered}$$
这就证明对所有$m\ge 2$有${\theta }_{k,m}=0$，而
$${\theta }_{k,1}=\frac{{\gamma }_1-{\sum }_{j=0}^{m-1}{\theta }_{k,k-j}{\theta }_{m,m-j}{\nu }_j}{{\nu }_{k-1}}=\frac{{\gamma }_1}{{\nu }_{k-1}}.$$
这就完成了结论的证明。

对于$\mathrm{M}\mathrm{A}(2)$序列，我们猜想也有类似的递推式，假设$\mathrm{M}\mathrm{A}(2)$模型是
$$X_t={\epsilon }_t+b_1{\epsilon }_{t-1}+b_2{\epsilon }_{t-2}.$$
所以对于$k\ge 3,{\gamma }_3=0$，容易递推得到：
$$\begin{gathered} {\nu }_0={\gamma }_0,{\theta }_{1,1}=\frac{{\gamma }_1}{{\gamma }_0},{\nu }_1={\gamma }_0-{\theta }_{1,1}^2{\nu }_0, \\ \forall n\ge 2:{\theta }_{n,2}=\frac{{\gamma }_2}{{\nu }_{n-2}},{\theta }_{n,1}=\frac{{\gamma }_1-{\theta }_{n,2}{\theta }_{n-1,1}{\nu }_{n-2}}{{\nu }_{n-1}}. \\ \forall n\ge 2:{\nu }_n={\gamma }_0-{\theta }_{n,2}^2{\nu }_{n-2}-{\theta }_{n,1}^2{\nu }_{n-1}. \\ \forall n\ge 3,j\ge 3,j\le n:{\theta }_{n,j}=0. \end{gathered}$$
事实上，$\mathrm{M}\mathrm{A}(q)$序列的递推都以前$q$行作为递推基础，后面的行都可以用递推公式表达。

3.$\mathrm{A}\mathrm{R}\mathrm{M}\mathrm{A}(1,1)$序列

$\mathrm{A}\mathrm{R}\mathrm{M}\mathrm{A}(p,q)$序列是最繁琐的，但$\mathrm{A}\mathrm{R}\mathrm{M}\mathrm{A}(1,1)$序列由于阶数较低，还是可以尝试计算的，此时
$$m=\max (p,q)=1.$$
设$\mathrm{A}\mathrm{R}\mathrm{M}\mathrm{A}(1,1)$序列为
$$X_t=a{X}_{t-1}+{\epsilon }_t+b{\epsilon }_{t-1},$$
构造的辅助序列为：
$$Y_t=\{\begin{array}{lc}\frac{X_t}{\sigma },& t=1,\\ \frac{X_t-a{X}_{t-1}}{\sigma },& t>2.\end{array}$$
根据我们之前讨论的顺序，先从其协方差求起，以下将时间指标设为$s,t$，由于$Y_t$不平稳，需要分别讨论。由于$Y_t$依赖于$X_t$，而$X_t$是平稳序列，所以先求$\mathrm{A}\mathrm{R}\mathrm{M}\mathrm{A}(1,1)$序列的自协方差函数。有
$$\begin{gathered} {\gamma }_0=\mathrm{E}(X_t^2)=\mathrm{E}(a{X}_{t-1}+{\epsilon }_t+b{\epsilon }_{t-1}{)}^2=a^2{\gamma }_0+(1+b^2+2ab){\sigma }^2, \\ \Rightarrow {\gamma }_0=\frac{1+b^2+2ab}{1-a^2}{\sigma }^2, \\ {\gamma }_1=\mathrm{E}(X_t{X}_{t-1})=\mathrm{E}[(a{X}_{t-1}+{\epsilon }_t+b{\epsilon }_{t-1})X_{t-1}]=a{\gamma }_0+b{\sigma }^2 \\ \Rightarrow {\gamma }_1=a{\gamma }_0+b{\sigma }^2, \\ \forall k\ge 2:{\gamma }_k=\mathrm{E}(X_t{X}_{t-k})=\mathrm{E}[(a{X}_{t-1}+{\epsilon }_t+b{\epsilon }_{t-1})X_{t-k}]=a{\gamma }_{k-1}. \end{gathered}$$
当$s=t=1$时情况比较简单，就是
$$\mathrm{E}(Y_s{Y}_t)=\frac{1}{{\sigma }^2}\mathrm{E}(X_1^2)=\frac{1+2ab+b^2}{1-a^2}.$$
当$s\ge t>1$时，$Y_t$可以视为$\mathrm{M}\mathrm{A}(1)$序列$\sigma Y_t={\epsilon }_t+b{\epsilon }_{t-1}$，所以其自协方差函数显然是
$$\begin{gathered} \mathrm{E}(Y_t^2)=1+b^2,t\ge 2, \\ \mathrm{E}(Y_{t+1}{Y}_t)=b,t\ge 1. \\ \mathrm{E}(Y_{t+k}{Y}_t)=0,t\ge 1,k\ge 2. \end{gathered}$$
这就求得了$Y_t$的协方差。由于$t>1$时，$Y_t$是$\mathrm{M}\mathrm{A}(1)$序列，所以可以直接代入$\mathrm{M}\mathrm{A}(1)$的递推结论：
$$\begin{gathered} {\theta }_{n,1}=\frac{\mathrm{E}(Y_n{Y}_{n+1})}{{\nu }_{n-1}}=\frac{b}{{\nu }_{n-1}}, \\ \forall n\ge j\ge 2:{\theta }_{n,j}=0, \\ {\nu }_n=\mathrm{E}(Y_n^2)-{\theta }_{n,1}^2{\nu }_{n-1}=1+b^2-\frac{b^2}{{\nu }_{n-1}}. \end{gathered}$$
最后加上递推初始值即可：
$${\nu }_0={\gamma }_0=\frac{1+2ab+b^2}{1-a^2}.$$
这里要注意，在使用$\mathrm{M}\mathrm{A}(1)$序列递推时，应当用$\mathrm{M}\mathrm{A}(1)$序列的自协方差函数（而不是$\mathrm{A}\mathrm{R}\mathrm{M}\mathrm{A}(1,1)$序列的），如果搞不明白为什么，可以看看非平稳序列的递推公式（《十八、时间序列的递推预测》）。而初始递推值又应该使用$\mathrm{A}\mathrm{R}\mathrm{M}\mathrm{A}(1,1)$序列的自协方差函数，因为递推开始时有${\nu }_0=\mathrm{E}(Y_1^2)$。

此时我们得到了$Y_t$的线性预测：
$${\hat{Y}}_{n+1}={\theta }_{n,1}{W}_n,W_n=Y_n-L(Y_n|{\mathbf{Y}}_{n-1}).$$
于是由于$X_{n+1}=\sigma Y_{n+1}+a{X}_n$，就有
$${\hat{X}}_{n+1}={\theta }_{n,1}{Z}_n+a{X}_n,Z_n=X_n-L(X_n|{\mathbf{X}}_n).$$
从上面的讨论我们也可以发现，讨论$\mathrm{A}\mathrm{R}\mathrm{M}\mathrm{A}(p,q)$序列时，我们往往不直接计算$\mathrm{A}\mathrm{R}\mathrm{M}\mathrm{A}(p,q)$序列的递推系数，而是选了一条看起来更艰难的路，构造了一个非平稳序列进行递推。但是，这个非平稳序列的后半部分是一个$\mathrm{M}\mathrm{A}(q)$序列，因此递推系数也有“截尾”性，还可以套$\mathrm{M}\mathrm{A}(q)$序列的递推式，因此只需要计算前$m$项的递推基础即可。