【統計検定準1級】時系列解析（4）系列相関

はじめに

この記事では時系列モデルにおける系列相関（serial correlation）の検定方法と、系列相関の疑いがあるときの対処法について扱います。ダービー・ワトソン（Durbin-Watoson）比（検定）、コクラン・オーカット法に関する、統計検定準1級レベルの内容です。

はじめに
ダービン・ワトソン比の定義
ダービン・ワトソン比の近似と解釈
系列相関のもとでの回帰モデルの推定（コクラン・オーカット法）
Rを使って可視化する
参考文献

ダービン・ワトソン比の定義

回帰分析では、誤差項について以下の仮定が置かれています。

誤差項の共分散がゼロ
誤差項はそれぞれ独立に平均 $0$ 、分散 $σ^2$ の正規分布に従う

時系列データを扱うとき、誤差項の系列 $u_1,u_2, \cdots , u_i$ の符号がプラスもしくはマイナスが連続で続く場合は正の系列相関がある、プラスとマイナスが交互で続く場合は負の系列相関がある、と言います。系列相関の特徴は、 $i$ 期の誤差項 $u_i$ の符号が $u_1,u_2 \cdots , u_{i-1}$ から予想されるということです。これは「誤差項はそれぞれ独立に正規分布に従う」という仮定に反します。この仮定が満たされない場合、不均一分散の場合と同じく、最小二乗推定量（OLSE）は最良線形不偏推定量ではなくなります。

ダービン・ワトソン比とは、回帰モデルと誤差項を仮定し、

$\begin{eqnarray} Y_i &=& α + β X_{i} + u_i \tag{1} \\ \\ u_i &=& ρ u_{i-1} + ε_i \tag{2} \end{eqnarray}$

以下の検定*1を行うための統計量です。

$\begin{eqnarray} H_o : ρ=0 \ vs \ H_1 : |ρ| < 1 \tag{3} \end{eqnarray}$

(2)の誤差項 $u_i$ は互いに独立に正規分布 $(0,σ^2)$ に従います。(2)のような誤差項の関数形をAR(1)（1st order autoregressive model）と呼びます。
(3)は、 $ρ=0$ ならば系列相関なし、 $ρ\neq 0$ ならば系列相関が存在します。 $ρ>0$ の場合は正の系列相関、 $ρ<0$ の場合は負の系列相関が存在することになります。ダービン・ワトソン比の定義は、

$\begin{eqnarray} DW= \frac{\sum_{i=2}^{n} ( \hat{u}_i - \hat{u}_{i-1} )^2}{\sum_{i=1}^{n} \hat{u}^2_i} \tag{4} \end{eqnarray}$

です。 $\hat{u_i}$ は(1)の推定結果を使って導出した残差です。

ダービン・ワトソン比の近似と解釈

ダービー・ワトソン比は、誤差項の一次の自己相関係数 $ρ$ を使って近似することができます。(4)の右辺を展開すると、

$\begin{eqnarray} DW &=& \frac{\sum_{i=2}^{n} \hat{u}^2_i + \sum_{i=2}^{n} \hat{u}^2_{i-1} - 2\sum_{i=2}^{n} \hat{u}^2_i \hat{u}^2_{i-1} }{\sum_{i=1}^{n} \hat{u}^2_i} \\ & \simeq & \frac{2 \sum_{i=2}^{n} \hat{u}^2_i - 2\sum_{i=2}^{n} \hat{u}^2_i \hat{u}^2_{i-1} }{\sum_{i=1}^{n} \hat{u}^2_i} \\ &=& 2(1-\hat{ρ})　\tag{5} \end{eqnarray}$

(5)ではサンプルサイズが十分に大きいとき、誤差項の二乗和は近しい値になると考えます。

$\begin{eqnarray} \sum_{i=1}^{n} \hat{u}^2_i \simeq \sum_{i=2}^{n} \hat{u}^2_i \simeq \sum_{i=2}^{n} \hat{u}^2_{i-1} \tag{6} \end{eqnarray}$

$\begin{eqnarray} \hat{ρ} = \frac{\sum_{i=2}^{n} \hat{u}^2_i \hat{u}^2_{i-1}}{\sum_{i=1}^{n} \hat{u}^2_i}　\tag{7} \end{eqnarray}$

$\hat{ρ}$ を求めるには、まず(1)より誤差項 $u_i$ の推定値を求め、それらを(2)に代入します。ダービン・ワトソン比の解釈は、

2前後のときは系列相関なし（ $\hat{ρ}=0$ のとき、 $DW \simeq 2$ ）
2より小さいときは正の系列相関
2より大きいときは負の系列相関

ただし、ダービン・ワトソン比は(2)の仮定が限定的であるという批判もあります。そのため、最近は参考指標として2に近いかどうかをざっくりチェックする、という程度で使われているようです。

系列相関のもとでの回帰モデルの推定（コクラン・オーカット法）

誤差項に系列相関が存在すると判断された場合、コクランオーカット法を用いて回帰モデルを推定することができます。
(1)の添え字を $i$ から $i-1$ に変えると、

$\begin{eqnarray} Y_{i-1} = α + β X_{i-1} + u_{i-1} \tag{8} \end{eqnarray}$

(1)から $ρ$ ×(8)を引くと、

$\begin{eqnarray} Y_{i} - ρY_{i-1} &=& (1-ρ)α + β(X_i - ρX_{i-1}) + (u_{i} -ρu_{i-1}) \\ &=& γ + β(X_i - ρX_{i-1}) + ε_i \tag{9} \end{eqnarray}$

(9)の変形は(2)の関係を用いています。 $γ=(1-ρ)α$ です。もし $ρ$ の値が既知であれば、(9)を最小二乗法で推定すれば効率的な推定量が得られます。しかし、実際には未知なので、 $ρ$ の推計値を用いて推定を行います。

具体的な手順としてはまず、(1)の回帰モデルを推定し、残差を $\hat{u_i}$ とします。次に、 $\hat{u_i}$ を $\hat{u_{i-1}}$ に回帰して、 $ρ$ の推計量 $\hat{ρ}$ を得ます。
そして、(9)に従い、

$\begin{eqnarray} Y^{*}_i &=& Y_i-\hat{ρ}Y_{i-1} \tag{10}　\\ X^{*}_i &=& X_i - \hat{ρ}X_{i-1} \tag{11} \end{eqnarray}$

を新たな変数として

$\begin{eqnarray} Y^{*}_{i} &=& γ + β X^{*}_i + ε_{i} \tag{12} \end{eqnarray}$

(11)を推定します。(11)の誤差項は、標準的仮定を満たすため、最良線形不偏推定量である $\hat{β}$ を求めることができます。

Rを使って可視化する

それではRの標準データセットを使って、実際に(1)を推定します。以下、2種類の150期の時系列データを使います。

売上高 $Y_i, i = 1, \cdots 150$
売上高の先行指標 $X_i, i = 1, \cdots 150$

(1)を推定します。線形な関数だと当てはまりが良いようです。

残差をプロットします。

上図を見る限り、残差は一期前の値に影響を受けている、すなわち系列相関が考えられます。ちなみに(1)のパラメータは $\hat{α}=30.88, \hat{β}=16.81$ です。DW比を計算すると、0.68828となり、やはり正の系列相関があることが分かります。

自己相関係数 $ρ$ のコレログラムを描くと、

青線は95%信頼区間です。

回帰モデル(1)には系列相関が発生しているため、コクラン・オーカット法を用いて(1)の $β$ を推定します。以下のresには(1)の回帰モデルの推定結果が格納されています。

#(1)の推定結果を使って残差を求める
resid1 <- resid(res)

#残差の一階ラグをとる
resid1_lag1 <- c(0,resid(res)[1:length(resid1)-1])

Sales_lag1 <- c(0,df[,2][1:length(df[,1])-1])
Sales.lead_lag1 <- c(0,df[,3][1:length(df[,2])-1])
df_dw <- cbind(df,resid1,resid1_lag1,Sales_lag1,Sales.lead_lag1)

#(2)を推定してρを求める
rho <- coef(summary(lm(resid1 ~ resid1_lag1)))[1,1]

#(10)(11)を求める
Sales_star<-cbind(c(df[,2][2:length(df[,1])]-rho*Sales_lag1[2:length(df[,1])]))
Sales.lead_star<-cbind(c(df[,3][2:length(df[,2])]-rho*Sales.lead_lag1[2:length(df[,2])]))

#(12)を推定してβを求める                  
df_dw2 <- cbind(Sales_star,Sales.lead_star)
cochrane_orcutt<-lm(Sales_star ~ Sales.lead_star)

#推定結果とDW検定の結果
summary(cochrane_orcutt)

コクラン・オーカット法を用いることで、(1)のパラメータは $\hat{α}=30.27, \hat{β}=16.78$ だと求まりました。実はこれらの値、系列相関を仮定せずに(1)を推定した場合の推定値とほとんど変わりません。主な理由は(1)の推定が当てはまりがよく、決定係数が $R = 0.905$ と十分に高いためです。そもそもモデルの当てはまりがよいので、誤差項にAR(1)を仮定したところで、推定結果はほとんど変わらない、ということです。