时间序列分析：偏自相关函数

Author: nex3z 2019-07-21

1. 问题

　　在 $\mathrm{AR}(p)$ 过程中，每个时刻的值都与历史时刻相关，导致其自相关函数呈现逐渐衰减，而不会像 $\mathrm{MA}(q)$ 那样出现截断。我们希望能够单独分析两个时刻随机变量之间的相关性，而不受其他时刻的影响，这样就可以方便地确定 $\mathrm{AR}(p)$ 过程的阶数。

2. 偏自相关的一般例子

　　看一个更一般的例子。R 中 isdals 包的 bodyfat 数据包含了 20 个人的体脂（Fat）、三头肌皮脂厚度（Triceps）、大腿围（Thigh）和中臂围（Midarm）的数据，期望通过 Triceps、Thigh 和 Midarm 预测 Fat 的数值。

library(isdals)
data(bodyfat)
attach(bodyfat)
pairs(cbind(Fat, Triceps, Thigh, Midarm))

由图 1 可以看出，Fat、Triceps 和 Thigh 之间似乎都有着线性关系，计算几个特征之间的相关可以印证这一点。

cor(cbind(Fat, Triceps, Thigh, Midarm))

输出为

              Fat   Triceps     Thigh    Midarm
Fat     1.0000000 0.8432654 0.8780896 0.1424440
Triceps 0.8432654 1.0000000 0.9238425 0.4577772
Thigh   0.8780896 0.9238425 1.0000000 0.0846675
Midarm  0.1424440 0.4577772 0.0846675 1.0000000

这在直观上也容易理解：对于一个人来说，他的 Fat、Triceps 和 Thigh 是与他的体型紧密相关的。

　　我们希望能消除 Thigh 的影响，单独衡量 Fat 和 Triceps 的关系，即偏自相关。此时可以训练两个线性模型，使用 Thigh 分别对 Fat 和 Triceps 进行预测，然后衡量预测值之间的相关：

Fat.hat = predict(lm(Fat ~ Thigh))
Triceps.hat = predict(lm(Triceps ~ Thigh))
cor((Fat - Fat.hat), (Triceps - Triceps.hat))

上面代码的输出为 0.1749822，可见消除 Thigh 的影响后，Fat 和 Triceps 的相关降低了不少。

　　也可以使用类似地方法消除多个变量的影响，例如消除 Thigh 和 Midarm 的影响后，衡量 Fat 和 Triceps 的相关：

Fat.hat = predict(lm(Fat ~ Thigh+Midarm))
Triceps.hat = predict(lm(Triceps ~ Thigh+Midarm))
cor((Fat - Fat.hat), (Triceps - Triceps.hat))

输出为 0.33815。

　　也可以用 ppcor 包的 pcor 函数计算偏自相关

library(ppcor)
pcor(cbind(Fat, Triceps, Thigh, Midarm))

输出为

$estimate
               Fat   Triceps      Thigh     Midarm
Fat      1.0000000 0.3381500 -0.2665991 -0.3240520
Triceps  0.3381500 1.0000000  0.9963725  0.9955918
Thigh   -0.2665991 0.9963725  1.0000000 -0.9926612
Midarm  -0.3240520 0.9955918 -0.9926612  1.0000000

可见 Fat 与 Triceps 的偏自相关为 0.3381500，与前面的计算相符。

3. 时间序列的偏自相关函数

　　定义 $k$ 阶滞后的偏自相关函数为，消除中间介入变量 $Y_{t-1}, Y_{t-2}, \cdots, Y_{t-k+1}$ 的影响后 $Y_t$ 和 $Y_{t-k}$ 的相关系数，记为 $\phi_{kk}$。

　　如果 $\{Y_t\}$ 式正态分布的时间序列，则令

$$\begin{equation} \phi_{kk} = \mathrm{Corr}(Y_k, Y_{t-k}|Y_{t-1}, Y_{t-2}, \cdots, Y_{t-k+1}) \tag{1} \end{equation}$$

即 $\phi_{kk}$ 是二元分布 $Y_k, Y_{t-k}$ 的以 $Y_{t-1}, Y_{t-2}, \cdots, Y_{t-k+1}$ 为条件的相关系数。

　　更一般地，构建一个线性模型，使用中间介入变量 $Y_{t-1}, Y_{t-2}, \cdots, Y_{t-k+1}$ 对 $Y_t$ 进行预测，例如预测值为

$$\begin{equation} \widehat Y_t = \beta_1 Y_{t-1} + \beta_2 Y_{t-2} + \cdots + \beta_{k-1} Y_{t-k+1} \end{equation}$$

选择参数 $\beta$ 使得预测的均方误差最小。选定 $\beta$ 后，由平稳性可知，给定相同的数据 $Y_{t-1}, Y_{t-2}, \cdots, Y_{t-k+1}$，对 $Y_{t-k}$ 的最优预测是

$$\begin{equation} \widehat Y_{t-k} = \beta_1 Y_{t-k+1} + \beta_2 Y_{t-k+2} + \cdots + \beta_{t-1} Y_{t-k+1} \end{equation}$$

此时 $k$ 阶滞后的偏自相关函数定义为

$$\begin{align} \phi_{kk} = \mathrm{Corr}&(Y_t – \widehat Y_t, Y_{t-k} – \widehat Y_{t-k}) \\ = \mathrm{Corr}&[Y_t – (\beta_1 Y_{t-1} + \beta_2 Y_{t-2} + \cdots + \beta_{k-1} Y_{t-k+1}), \\ &\;Y_{t-k} – (\beta_1 Y_{t-k+1} + \beta_2 Y_{t-k+2} + \cdots + \beta_{t-1} Y_{t-k+1})] \end{align}$$

　　$Y_t$ 只基于 $Y_{t-1}$ 的最优线性预测为 $\rho_1 Y_{t-1}$，因此

$$\begin{align} &\mathrm{Cov}(Y_t – \rho_1 Y_{t-1}, Y_{t-2} – \rho_1 Y_{t-1}) \\ &= \mathrm{Cov}(Y_t, Y_{t-2}) – \rho_1\mathrm{Cov}(Y_t, Y_{t-1}) – \rho_1\mathrm{Cov}(Y_{t-1}, Y_{t-2}) + \rho_1^2\mathrm{Cov}(Y_{t-1}, Y_{t-1}) \\ &= \gamma_2 – \rho_1\gamma_1 – \rho_1\gamma_1 + \rho_1^2\gamma_0 \\ &= \gamma_0(\rho_2 – \rho_1^2 – \rho_1^2 + \rho_1^2) = \gamma_0(\rho_2 – \rho_1^2) \end{align}$$

又由

$$\begin{align} \mathrm{Var}(Y_t – \rho_1 Y_{t-1}) &= \mathrm{Var}(Y_{t-2} – \rho_1 Y_{t-1}) \\ &= \mathrm{Var}(Y_{t-2}) + \rho_1^2 \mathrm{Var}(Y_{t-1}) – 2\rho_1\mathrm{Cov}(Y_{t-2}, Y_{t-1}) \\ &= \gamma_0 + \rho_1^2\gamma_0 – 2\rho_1\gamma_1 = \gamma_0(1 + \rho_1^2 – 2\rho_1^2) = \gamma_0(1 – \rho_1^2) \end{align}$$

于是

$$\begin{equation} \phi_{22} = \frac{\gamma_0(\rho_2 – \rho_1^2)}{\gamma_0(1 – \rho_1^2)} = \frac{\rho_2 – \rho_1^2}{1 – \rho_1^2} \end{equation}$$

对于 $\mathrm{AR}(1)$ 过程，有 $\rho_k = \phi^k$，于是

$$\begin{equation} \phi_{22} = \frac{\phi^2 – \phi^2}{1 – \phi^2} = 0 \end{equation}$$

即 $\mathrm{AR}(1)$ 的滞后为 $2$ 的偏自相关为 $0$。对于 $\mathrm{AR}(1)$ 过程，当 $k = 1$ 时，有 $\phi_{11} \neq 0$；当 $k > 1$ 时，有 $\phi_{kk} = 0$。

　　一般地，$\mathrm{AR}(p)$ 过程的偏自相关函数图像在滞后超过 $p$ 后截尾。