Goodな生活

Goodな生活

それは、いきあたりばっちりな人生。being good and haphazard.

MENU

【統計検定準1級】系列相関の検定方法と対処法

はじめに

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

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

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

  • 誤差項の共分散がゼロ
  • 誤差項はそれぞれ独立に平均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

f:id:good_na_life:20200516214102p:plainf:id:good_na_life:20200516214142p:plain

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

f:id:good_na_life:20200517084524p:plain

残差をプロットします。

f:id:good_na_life:20200516214637p:plain

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

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

f:id:good_na_life:20200516220017p:plain

青線は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)を仮定したところで、推定結果はほとんど変わらない、ということです。

読んでいただき、ありがとうございました。

参考文献

日本統計学会公式認定 統計検定 1級・準1級 公式問題集[2016〜2017年]

日本統計学会公式認定 統計検定 1級・準1級 公式問題集[2016〜2017年]

  • 発売日: 2018/03/14
  • メディア: 単行本(ソフトカバー)

計量経済学 (サピエンティア)

計量経済学 (サピエンティア)

*1:対立仮説でφ > 1の場合を考えないのはデータ分析において発散的な確率過程を扱うことがほとんどないためです。