【时间序列分析】10.偏相关系数与Levinson递推公式

十、偏相关系数与Levinson递推公式

1.对Yule-Walker方程的推广——偏相关系数

在提出Yule-Walker方程时，我们提出过一个问题——既然$\mathrm{A}\mathrm{R}(p)$序列只需要$p$个自相关系数就能够建立Yule-Walker方程，为什么要将这个方程延拓到$n$阶呢？其实，在实际生活中，我们往往并不知道一个$\mathrm{A}\mathrm{R}(p)$序列的阶数，即我们不知道自回归系数一共有几个，在这种情况下，根本无法确定$p$，自然不能直接得到${\Gamma }_p,{\mathbf{\gamma }}_p$具体的形态。本篇文章会解答这个问题。

首先我们将Yule-Walker方程的适用情况，从$\mathrm{A}\mathrm{R}(p)$序列推广到一切平稳序列 { X t } \{X_t\} {Xt​}，假设平稳序列有自协方差函数 { γ k } \{\gamma_k\} {γk​}和自协方差矩阵${\Gamma }_n$，则关于${\mathbf{a}}_n$的方程
$${\mathbf{\gamma }}_n={\Gamma }_n{\mathbf{a}}_n$$
被称为 { γ k } \{\gamma_k\} {γk​}的$n$阶Yule-Walker方程，这里的${\mathbf{a}}_n=(a_{n,1},a_{n,2},\cdots ,a_{n,n}{)}^{\prime }$被称为$n$阶Yule-Walker系数。如果${\Gamma }_n$正定则${\mathbf{a}}_n={\Gamma }_n^{-1}{\mathbf{\gamma }}_n$，此时的Yule-Walker系数由自协方差函数完全决定了。

需要注意这里的Yule-Walker系数与$\mathrm{A}\mathrm{R}(p)$序列的不一样。在$\mathrm{A}\mathrm{R}(p)$序列中，如果$n\ge p$，$n$阶Yule-Walker系数${\mathbf{a}}_n$的前$p$个值就是其自回归系数，后面的值都是0，而在$n<p$的情况则不一定是自回归系数的前$n$个值。这里定义的$n$阶Yule-Walker系数，是对任何平稳序列 { X t } \{X_t\} {Xt​}或 { γ k } \{\gamma_k\} {γk​}都存在的，是方程的解向量，自然不存在自回归系数这一概念。

并且，$n$阶Yule-Walker系数是依赖于$n$的，在没有给出$n$时是无法对其进行求解的。而且，不同阶数的${\mathbf{a}}_n$的同位置分量往往数值也不同。比如$\mathrm{A}\mathrm{R}(p)$序列的$p$阶Yule-Walker系数是自回归系数，但小于$p$阶时，如${\mathbf{a}}_1={\gamma }_1/{\gamma }_0$，就不等于第一个自回归系数$a_1$。不过作为特例，$\mathrm{A}\mathrm{R}(p)$序列在$p$阶以上的Yule-Walker系数的同位置分量是一样的。

最后，由于${\Gamma }_n$正定时，${\mathbf{a}}_n$由 { γ k } \{\gamma_k\} {γk​}完全确定，所以我们给${\Gamma }_n$正定时的$n$阶Yule-Walker系数起一个名字——偏相关系数。我们接下来将探讨Yule-Walker系数的应用。

2.最优线性预测与Yule-Walker系数

要探究Yule-Walker系数的作用，我们先对上一篇最后提到的完全可预测性作出延伸。虽然对于极大多数平稳序列，$X_{n+1}$不能由$X_1,\cdots ,X_n$线性表示，但我们可以用这$n$个向量的线性组合对$X_{n+1}$作出估计，也就是用$X_1,\cdots ,X_n$的线性组合作为$X_{n+1}$的估计量。

不同的线性组合会产生不同的估计量，不同的估计量之间必然存在优劣，一个常用来评判估计量好坏的指标是均方误差，即$\mathrm{E}(X_{n+1}-\sum _{j=1}^n{b}_j{X}_i{)}^2$。在这么多$X_{n+1}$的估计量中，存在一个均方误差最小的线性组合，它被称为最优线性预测。

巧合的是，最优线性预测对应的系数，与Yule-Walker系数${\mathbf{a}}_n$完全一致，也就是满足${\mathbf{\gamma }}_n={\Gamma }_n{\mathbf{a}}_n$的这组系数，下面我们给出定理及其证明：

最优线性预测：如果${\mathbf{a}}_n$是 { X t } \{X_t\} {Xt​}的$n$阶Yule-Walker系数，则$X_{n+1}$的最优线性预测是
$${\hat{X}}_{n+1}=\sum _{j=1}^n{a}_{n,j}{X}_{t+1-j}={\mathbf{a}}_n^{\prime }{\mathbf{X}}_n.{\mathbf{X}}_n=(X_n,\cdots ,X_1{)}^{\prime }.$$
将预测的均方误差设置为
$${\sigma }_n^2\stackrel{\Delta }{=}\mathrm{E}(X_{n+1}-{\mathbf{a}}_n^{\prime }{\mathbf{X}}_n{)}^2.$$
现在对此定理作出证明。任取一组常数${\mathbf{b}}_n$，其线性预测为${\mathbf{b}}_n^{\prime }{\mathbf{X}}_n$，则有
$$\begin{array}{rl}& \mathrm{E}(X_{n+1}-{\mathbf{b}}_n^{\prime }{\mathbf{X}}_n{)}^2\\ =& \mathrm{E}(X_{n+1}-{\mathbf{a}}_n^{\prime }{\mathbf{X}}_n+({\mathbf{a}}_n^{\prime }-{\mathbf{b}}_n{)}^{\prime }{\mathbf{X}}_n{)}^2\\ =& \mathrm{E}(X_{n+1}-{\mathbf{a}}_n^{\prime }{\mathbf{X}}_n{)}^2+\mathrm{E}[({\mathbf{a}}_n-{\mathbf{b}}_n{)}^{\prime }{\mathbf{X}}_n{]}^2+2\mathrm{E}(X_{n+1}-{\mathbf{a}}_n^{\prime }{\mathbf{X}}_n)[({\mathbf{a}}_n-{\mathbf{b}}_n{)}^{\prime }{\mathbf{X}}_n]\\ =& \mathrm{E}(X_{n+1}-{\mathbf{a}}_n^{\prime }{\mathbf{X}}_n{)}^2+\mathrm{E}[({\mathbf{a}}_n-{\mathbf{b}}_n{)}^{\prime }{\mathbf{X}}_n{]}^2+2\mathrm{E}(X_{n+1}-{\mathbf{a}}_n^{\prime }{\mathbf{X}}_n){\mathbf{X}}_n^{\prime }({\mathbf{a}}_n-{\mathbf{b}}_n)\\ =& \mathrm{E}(X_{n+1}-{\mathbf{a}}_n^{\prime }{\mathbf{X}}_n{)}^2+\mathrm{E}[({\mathbf{a}}_n-{\mathbf{b}}_n{)}^{\prime }{\mathbf{X}}_n{]}^2+2({\mathbf{\gamma }}_n-{\mathbf{a}}_n^{\prime }{\Gamma }_n)({\mathbf{a}}_n-{\mathbf{b}}_n)\\ =& \mathrm{E}(X_{n+1}-{\mathbf{a}}_n^{\prime }{\mathbf{X}}_n{)}^2+\mathrm{E}[({\mathbf{a}}_n-{\mathbf{b}}_n{)}^{\prime }{\mathbf{X}}_n{]}^2\\ \ge & \mathrm{E}(X_{n+1}-{\mathbf{a}}_n^{\prime }{\mathbf{X}}_n{)}^2.\end{array}$$
这里第三行到第四行，是将$({\mathbf{a}}_n-{\mathbf{b}}_n{)}^{\prime }{\mathbf{X}}_n$进行转置，因为1阶矩阵转置等于自身；第四行到第五行是自然的期望向自协方差函数转化；第五行到第六行是因为${\mathbf{\gamma }}_n={\mathbf{a}}_n^{\prime }{\Gamma }_n$。

这个定理说明，偏相关系数的一个重要作用，是用历史的信息对未来做出预测，这种预测是最佳线性预测。

3.Yule-Walker系数的计算——Levinson递推

要求解Yule-Walker系数很简单，只要找到${\Gamma }_n$的逆矩阵，计算${\mathbf{a}}_n={\Gamma }_n^{-1}{\mathbf{\gamma }}_n$即可，但这种简单仅仅是理论上的，因为随着阶数$n$的增长，${\Gamma }_n$的求逆计算复杂度会很快增加。为了简化计算，我们通常采用Levinson递推公式，它用来递推地计算Yule-Walker系数与最优线性预测的均方误差。

Levinson递推公式：如果${\Gamma }_{n+1}$正定，则对$1\le k\le n$有
$$\{\begin{array}{l}{a}_{1,1}=\frac{{\gamma }_1}{{\gamma }_0},\\ {\sigma }_0^2={\gamma }_0,\\ ----\\ {\sigma }_k^2={\sigma }_{k-1}^2(1-a_{k,k}^2),\\ a_{k+1,k+1}=\frac{{\gamma }_{k+1}-{\sum }_{j=1}^k{\gamma }_j{a}_{k,k+1-j}}{{\gamma }_0-{\sum }_{j=1}^k{\gamma }_j{a}_{k,j}}=\frac{{\gamma }_{k+1}-{\gamma }_k{a}_{k,1}-\cdots -{\gamma }_1{a}_{k,k}}{{\gamma }_0-{\gamma }_1{a}_{k,1}-\cdots -{\gamma }_k{a}_{k,k}},\\ a_{k+1,j}=a_{k,j}-a_{k+1,k+1}{a}_{k,k+1-j},1\le j\le k.\end{array}$$
这里${\sigma }_k^2$是用${\mathbf{X}}_k$预测${\hat{X}}_{k+1}$时的均方误差。

看起来Levinson递推公式非常硕大，但具有非常好记的规律。首先，前两个式子是递推基础，可以直接通过计算得到；后三个式子是递推关系，即在获得上一行Yule-Walker系数和均方误差后，可以通过这三个式子计算下一行的Yule-Walker系数和均方误差。

然后看计算顺序，在获得上一行的所有Yule-Walker系数后，先计算的是下一行中多出来的位置的系数值，再计算下一行的其他Yule-Walker系数值，最后计算下一行的预测均方误差。在整个计算过程中的Yule-Walker系数计算顺序可以如下呈现（均方误差在获得某行所有数值后计算）：
$$a_{1,1}\to a_{2,2}\to a_{2,\cdot }\to a_{3,3}\to a_{3,\cdot }\to \cdots \to a_{k,k}\to a_{k,\cdot }\to a_{k+1,k+1}\to a_{k+1,\cdot }\to \cdots$$
最后观察计算形式，在$a_{k+1,k+1}$的计算中，分子是协方差倒着排列从${\gamma }_{k+1}$开始与顺序排列的Y-W相乘的，而分母是协方差正着排列从${\gamma }_0$开始与顺序排列的Y-W系数相乘的；在$a_{k+1,j}$的计算中，被减数是上一行对应位置的Y-W系数，减数则要用刚算出来的$a_{k+1,k+1}$与上一行对称位置的Y-W系数相乘；最后是预测均方误差，是在上一行的均方误差基础上乘上$(1-a_{k+1,k+1}^2)$，可以看出，随着阶数的增长，均方误差不会增长。

4.偏相关系数的性质

偏相关系数是${\Gamma }_n$正定时的Y-W系数，它满足最小相位条件，即如果${\Gamma }_{n+1}$正定，则
$$1-\sum _{j=1}^n{a}_{n,j}{z}^j\ne 0,|z|\le 1.$$
现在我们可以回到最初提出的问题：在$\mathrm{A}\mathrm{R}(p)$自回归系数不定阶时，如何确定其阶数并计算其自回归系数？归根到底，就是阶数的确定。为了解决这个问题，我们需要知道$\mathrm{A}\mathrm{R}(p)$序列偏相关系数的一个重要性质——$p$后截尾。

$p$后截尾性：对于平稳序列，如果在$n\ge p$时偏相关系数满足
$$a_{n,n}=\{\begin{array}{lc}{a}_p,& n=p,\\ 0,& n>p.\end{array}$$
就称它的偏相关系数是$p$后截尾的。

$\mathrm{A}\mathrm{R}(p)$的判定：零均值平稳序列 { X t } \{X_t\} {Xt​}是$\mathrm{A}\mathrm{R}(p)$序列的充要条件是，它的偏相关系数$a_{n,n}$在$p$后截尾。

这个定理给出了$\mathrm{A}\mathrm{R}(p)$序列的重要判定方式，即计算其$n$阶偏相关系数，如果观察到在某一位后都是0，则这个序列是$\mathrm{A}\mathrm{R}(p)$序列，并且自然地得到它的阶数。

证明其充分性，记${\mathbf{a}}_p=(a_{p,1},\cdots ,a_{p,p}{)}^{\prime }=(a_1,\cdots ,a_p{)}^{\prime }$，既然它是$p$截尾的，那么$a_{p+k,p+k}=0$，由Levinson递推公式，得到
$$\begin{gathered} a_{p+1,j}=a_{p,j}-a_{p+1,p+1}{a}_{p,p+1-j}=a_j,1\le j\le p, \\ {a}_{p+k,j}=a_{p+k-1,j}=\cdots =a_{p,j}=a_j,1\le j\le p,k\ge 2, \\ {a}_{p+k,j}=a_{j,j}=0,p<j\le p+k. \end{gathered}$$
这说明对$n\ge p$，总有
$${\mathbf{a}}_n^{\prime }=(a_{n,1},a_{n,2},\cdots ,a_{n,n})=(a_1,\cdots ,a_p,0,\cdots ,0).$$
那么由Yule-Walker方程，$k\ge 1$时${\gamma }_k=\sum _{j=1}^p{a}_j{\gamma }_{k-j}$。现定义${\epsilon }_t=X_t-\sum _{j=1}^p{a}_j{X}_{t-j}$，则如果 { ε t } \{\varepsilon_t\} {εt​}是白噪声，就说明$X_t$是$\mathrm{A}\mathrm{R}(p)$序列。
$$\begin{gathered} \mathrm{E}{\epsilon }_t=\mathrm{E}(X_t-\sum _{j=1}^p{a}_j{X}_{t-j})=0; \\ \mathrm{D}{\epsilon }_t=\mathrm{E}{\epsilon }_t^2={\sigma }_p^2>0. \end{gathered}$$
先证明${\epsilon }_t$与$t$时刻之前时刻$s$时候的$X_s$无关，也就是$\forall t>s$，
$$\mathrm{E}({\epsilon }_t{X}_s)=\mathrm{E}\left[\left(X_t-\sum _{j=1}^p{a}_j{X}_{t-j}\right)X_s\right]={\gamma }_{t-s}-\sum _{j=1}^p{a}_j{\gamma }_{t-s-j}=0,$$
于是有$\forall t>s$，
$$\mathrm{E}({\epsilon }_t{\epsilon }_s)=\mathrm{E}\left[{\epsilon }_t\left(X_s-\sum _{j=1}^p{a}_j{X}_{s-j}\right)\right]=0.$$
所以 { ε t } \{\varepsilon_t\} {εt​}是$\mathrm{W}\mathrm{N}(0,{\sigma }_p^2)$，并且我们之前证明了偏相关系数满足最小相位条件（稳定性条件），所以 { X t } \{X_t\} {Xt​}是$\mathrm{A}\mathrm{R}(p)$序列。

反过来，如果零均值平稳序列 { X t } \{X_t\} {Xt​}是$\mathrm{A}\mathrm{R}(p)$序列，则由Yule-Walker方程，当$n\ge p$时自然有
$${\mathbf{\gamma }}_n={\Gamma }_n{\mathbf{a}}_n,{\mathbf{a}}_n=(a_1,\cdots ,a_p,0,\cdots ,0).$$
这就说明偏相关系数是$p$后截尾的。

有了这个定理，我们验证某序列是平稳的就有了固定的模式：先根据其历史信息计算自协方差函数的估计值 { γ ^ k } \{\hat \gamma_k\} {γ^​k​}，如果它是收敛的，则根据${\hat{\gamma }}_k$计算偏相关系数的估计值 { a ^ n , n } \{\hat a_{n,n}\} {a^n,n​}，如果它是$p$后截尾的就说明原序列是平稳序列，且阶数就是截尾位置$p$。

回顾总结

1. Yule-Walker方程是对任何平稳序列的方程：${\mathbf{\gamma }}_n={\Gamma }_n{\mathbf{a}}_n$，方程的解${\mathbf{a}}_n$称为Yule-Walker系数。

2. 如果${\Gamma }_n$是正定矩阵，则${\mathbf{a}}_n={\Gamma }^{-n}{\mathbf{\gamma }}_n$，即Yule-Walker系数由自协方差函数唯一确定，此时的Yule-Walker系数被称为偏相关系数。

3. Yule-Walker系数可以用来预测平稳序列。如果已知前$n$个历史信息，则$X_{n+1}$的最优线性预测为
   $${\mathbf{a}}_n^{\prime }{\mathbf{X}}_n=\sum _{j=1}^{\infty }{a}_{n,j}{X}_{n+1-j}.$$
   均方误差为
   $${\sigma }_n^2=\mathrm{E}{\left(X_{n+1}-{\mathbf{a}}_n^{\prime }{\mathbf{X}}_n\right)}^2.$$

4. 在已知自协方差函数时，用来计算偏相关系数的常用方法不是求正定矩阵${\Gamma }_{n+1}$逆，而是用Levinson递推公式，其形式是
   $$\{\begin{array}{l}{a}_{1,1}={\gamma }_1/{\gamma }_0,\\ {\sigma }_1^2={\gamma }_0,\\ {\sigma }_k^2={\sigma }_{k-1}^2(1-a_{k,k}^2),\\ a_{k+1,k+1}=\frac{{\gamma }_{k+1}-a_{k,1}{\gamma }_k-\cdots -a_{k,k}{\gamma }_1}{{\gamma }_0-a_{k,1}{\gamma }_1-\cdots -a_{k,k}{\gamma }_k},\\ a_{k+1,j}=a_{k,j}-a_{k+1,k+1}{a}_{k,k+1-j}.\end{array}$$

5. $p$后截尾指的是$a_{n,n}$在$n>p$时恒为0。对于零均值平稳序列 { X t } \{X_t\} {Xt​}，它是$\mathrm{A}\mathrm{R}(p)$序列的充要条件是它的偏相关系数$p$后截尾。

6. 验证$\mathrm{A}\mathrm{R}(p)$序列的一般步骤：根据样本求自协方差函数的估计值、观察自协方差函数的收敛性、用Levinson递推公式计算偏相关系数的估计值、观察偏相关系数的后截尾性。