R Time Series Analysis 時系列解析(12) 単位根検定1

単位根

$y_{t} = \phi y_{t-1} + a_{t},a_t \sim N(0, \sigma)$

において $\phi = 1$ となる過程を「ランダム・ウォーク(RW, 酔歩)」と呼ぶ。この場合、この過程は、特性方程式 $z - \phi = 0$ が $\phi = 1$ のとき1(単位根)を解にもつため、単位根過程と言う事が出来る。ちなみに、 $|\phi|$ <1ならこの過程は定常過程となり、 $|\phi|$ >1なら確率過程とはならない。

単位根検定が必要とされる理由

確定的トレンドと確率的トレンド

例えば、 $\triangle y_{t} = \mu$ の場合、 $y_{t}$ は毎期一定値ずつ増減する。一方、 $\triangle y_{t} = a_{t}$ とすると、 $y_{t}$ は毎期確率変数で増減が決まる。前者を確定的トレンドと言い、後者を確率的トレンドという。 (ただし、 $\triangle y_{t} = y_{t} - y_{t-1}$ )

見分けがつかない例

以下の2つの例はよくある例。

1つ目：

$y_{t}=\mu+\beta t +\psi (B) a_{t},\psi (B)=\psi_{0}+\psi_{1} B_{1}+\psi_{2} B^2 +\cdots,|\psi_{j}|$ <1, $for \forall j$

ただし、 Bはラグ演算子 $By_{t} \equiv y_{t-1}$ を表す。

2つ目：

$y_{t} = y_{t-1} + a_{t} + \mu$

1つ目は確定的トレンド $\mu + \beta t$ と定常過程MA(∞) $\psi (B) a_{t}$ が合わさった形、2つ目は確定的なトレンド $\mu$ とRW $y_{t-1} + a_{t}$ が合わさった形をしている。これらを図示する例を以下に示す。

n <- 100
d <- 0
sig2 <- 1.2
mu <- 0
beta <- 0.07

phi <- c(0.5,0.3);p <- length(phi)
theta <-  c(0.3,0.2);q <- length(theta)
time <- 1:n; trend <- mu+beta*time
pdq <- list(order=c(p,d,q),ar=phi,ma=theta)
wt <- arima.sim(n,model=pdq,sd=sqrt(sig2))
y1 <- trend + wt
y2 <- cumsum(beta+rnorm(n,sd=sqrt(sig2)))
yl <- c(min(y1,y2),max(y1,y2))
ttl <- "確率的トレンドとドリフト付きRW"
plot(y1,type="l",main=ttl,ylim=yl)
lines(y2,lty=2)

hanrei <- c("確定的トレンド+定常過程","ドリフト付きRW")
xx <- 0
yy <- max(y1,y2)
legend(xx,yy,legend=hanrei,lty=c(1,2),box.lty=0,col=0)

f:id:nakhirot:20130707143310p:plain

実際の時系列データには様々な"仕組"が内在し、見た目だけではどの過程かは見分けがつかない。時系列データにトレンドがあるのか、トレンドを除いたらそれは定常過程なのか、単位根を含むのかについて調べる手段が必要と思われる。

見せかけの回帰

全く独立にRWで発生させた系列 $y_{t} = y_{t-1} + a_{t},x_{t} = x_{t-1} + b_{t}$ を用意し、これを単回帰 $y_{t} = \alpha + \beta x_{t} + a_{t},a_t \sim N(0, \sigma)$ により推定する場合、一般に回帰係数が有意とである場合がある。

#2つの独立なRWを、単回帰分析で推定し、係数が有意となる割合を調べる
misekake <- function(n) {
  p <- numeric(n); cnta <- 0; r <- 0.05
  for (j in 1:n) {
   y <- cumsum(rnorm(100)); x <- cumsum(rnorm(100))
   kekka <- summary(lm(y~x))
   q <- kekka$coefficients[[2,4]]
   p[j] <- q; if(q<r){cnta <- cnta + 1}
  }
  rate <- cnta/n*100
  print(paste("n =",n,"rate = ",rate))
}

for(j in 1:10){misekake(100)}

結果は、

[1] "n = 100 rate =  78"
[1] "n = 100 rate =  74"
[1] "n = 100 rate =  77"
[1] "n = 100 rate =  76"
[1] "n = 100 rate =  70"
[1] "n = 100 rate =  76"
[1] "n = 100 rate =  73"
[1] "n = 100 rate =  77"
[1] "n = 100 rate =  75"
[1] "n = 100 rate =  75"

となり、7割強の場合で単回帰係数が有意となることが確認される。このように、単位根を含む系列で回帰を行うと、見せかけの回帰が発生する場合がある。

単位根検定

単位根を持つか否かを判定する方法として有名なのがADF検定（Augmented DF test)。まず、基本的なアイデアを整理する。検定式で最もシンプルなのはAR(1)の場合。

$y_{t} = \phi y_{t-1} + a_{t}$

において、これが単位根を持つには $\hat{\phi} = 1$ の帰無仮説をt検定し、棄却されれば単位根を持たないと判断すればよい。ただし、 $\phi = 1$ の場合、求める推定量 $\hat{\phi}$ のt値がt分布をしないことが証明されている(Fuller 1976)ため、以下のように変形する:

$\triangle y_{t} = (\phi -1) y_{t-1}+a_{t}$

と変形し、 $\pi \equiv \phi -1 = 0$ を帰無仮説としてt検定を行う。帰無仮説が棄却されない場合は、単位根を持つ可能性を否定出来ない。

上記をp次に拡張する。AR(p)モデル

$y_{t} = \phi_{1} y_{t-1} + \phi_{2} y_{t-2} + \cdots+\phi_{p} y_{t-p} +a_{t}$

が単位根を1つ持つとすると、これは下記のように変形できる。

$(1-\phi_{1} B - \phi_{2} B^2 - \cdots-\phi_{p} B^p)y_{t}=a_{t}$

$(1-\psi_{1}^* B - \psi_{2}^* B^2 - \cdots - \psi_{p-1}^* B^{p-1})(1 - B)y_{t}=a_{t}$

$(1-\psi_{1}^* B - \psi_{2}^* B^2 - \cdots - \psi_{p-1}^* B^{p-1})\triangle y_{t}=a_{t}$

よって、

$\triangle y_{t} = \psi_{1}^* \triangle y_{t-1} + \psi_{2}^* \triangle y_{t-2} + \cdots+\psi_{p-1}^* \triangle y_{t-p+1}+a_{t}$

が成り立つことを踏まえ、検定式を

$\triangle y_{t} = (\phi -1) y_{t-1} + \psi_{1}^* \triangle y_{t-1} + \psi_{2}^* \triangle y_{t-2} +\cdots+\psi_{p-1}^* \triangle y_{t-p+1}+a_{t}$

と置く。実際はトレンドを含む可能性もあるので、

$\triangle y_{t} = \beta_{1} + \beta_{2} t + \pi y_{t-1} + \psi_{1}^* \triangle y_{t-1} + \psi_{2}^* \triangle y_{t-2} + \cdots+\psi_{p-1}^* \triangle y_{t-p+1}+a_{t}$

などとトレンドの項(上記では線形トレンド $\beta_{1} + \beta_{2}$ を含めて置くことが多い。ここによると、ADF検定の論点は3つ。

1. 次数pをいくつにすべきか?

2. トレンド項を加えるべきか?

3. $\pi = 0$ が棄却出来ない場合、どのように考えるべきか?

Rではurcaパッケージのur.df()関数が上記を扱う。株価のデータを使って早速やってみる。

#終値に対して、トレンド項(線形項)を含む4次(p=4)のACF検定
stockACF <- summary(ur.df(stocksub[,2],type="trend",lag=4))
stockACF

以下の結果が表示される。

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 
Test regression trend
Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)

Residuals:
    Min      1Q  Median      3Q     Max 
-80.421 -11.975  -0.543  11.078  67.279 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.554065   2.719037   0.204    0.839    
z.lag.1     -0.975098   0.143974  -6.773 9.88e-11 ***
tt           0.006029   0.018799   0.321    0.749    
z.diff.lag1 -0.025066   0.128846  -0.195    0.846    
z.diff.lag2  0.039974   0.110123   0.363    0.717    
z.diff.lag3 -0.057675   0.092208  -0.625    0.532    
z.diff.lag4 -0.009274   0.065302  -0.142    0.887    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 20.66 on 237 degrees of freedom
Multiple R-squared:  0.5141,	Adjusted R-squared:  0.5018 
F-statistic:  41.8 on 6 and 237 DF,  p-value: < 2.2e-16

Value of test-statistic is: -6.7727 15.2971 22.9349 

Critical values for test statistics: 
      1pct  5pct 10pct
tau3 -3.99 -3.43 -3.13
phi2  6.22  4.75  4.07
phi3  8.43  6.49  5.47

少しずつ見ていく。

Test regression trend
Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)

これはz.diffは $\triangle y_{t}$ 、z.lag.1は $y_{t-1}$ の項、1は定数項、ttは線形項、z.diff.lagは $\psi_{1}^* \triangle y_{t-1} + \psi_{2}^* \triangle y_{t-2} + \cdots+\psi_{p-1}^* \triangle y_{t-p+1}+a_{t}$ を表しているのであろう。

Residuals:
    Min      1Q  Median      3Q     Max 
-80.421 -11.975  -0.543  11.078  67.279

これは残差の正規分布性のチェックのためのものか。通常の線形回帰と同様、残差の正規分布性が前提となっているのであろう。

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.554065   2.719037   0.204    0.839    
z.lag.1     -0.975098   0.143974  -6.773 9.88e-11 ***
tt           0.006029   0.018799   0.321    0.749    
z.diff.lag1 -0.025066   0.128846  -0.195    0.846    
z.diff.lag2  0.039974   0.110123   0.363    0.717    
z.diff.lag3 -0.057675   0.092208  -0.625    0.532    
z.diff.lag4 -0.009274   0.065302  -0.142    0.887    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

z.lag.1は $y_{t-1}$ の項以外有意水準5%で係数が有意ではない。つまり、

$\triangle y_{t} = -0.975098 y_{t-1} + a_{t}$

$\Leftrightarrow \triangle y_{t} = 0.024092 y_{t-1} + a_{t}$

と推定出来るということか？実質はAR(1)ということか。

Residual standard error: 20.66 on 237 degrees of freedom
Multiple R-squared:  0.5141,	Adjusted R-squared:  0.5018 
F-statistic:  41.8 on 6 and 237 DF,  p-value: < 2.2e-16

順に残差の標準誤差、決定係数、自由度調整済み決定係数、F値、p.value。回帰と読み方は同様のようだ。

Value of test-statistic is: -6.7727 15.2971 22.9349

一番左は $\hat{\pi}$ のt値らしい。一番右の値は帰無仮説 $\pi=0,\beta_{2}=0$ のF検定のF値らしい。真ん中の値は帰無仮説 $\pi=0,\beta_{1}=0,\beta_{2}=0$ のF検定のF値とのこと。

Critical values for test statistics: 
      1pct  5pct 10pct
tau3 -3.99 -3.43 -3.13
phi2  6.22  4.75  4.07
phi3  8.43  6.49  5.47

tau3の5pctは、帰無仮説 $\pi = 0$ に対応しているとのこと。t値は-6.7727で5pctを下回っているので、帰無仮説 $\pi = 0$ は棄却される。つまり、単位根はない。同様にphi2, phi3はそれぞれ $\pi=0,\beta_{1}=0,\beta_{2}=0$ 、 $\pi=0,\beta_{2}=0$ に対応し、F値はそれぞれ15.2971, 22.9349だから棄却される。

ur.df関数のオプション

なお、ur.df関数のtype引数と各モデルは下記のように対応する。

stockACF <- summary(ur.df(stocksub[,2],type="drift"))

だと、

$\triangle y_{t} = \beta_{1} + \pi y_{t-1} + \psi_{1}^* \triangle y_{t-1} + \psi_{2}^* \triangle y_{t-2} + \cdots+\psi_{p-1}^* \triangle y_{t-p+1}+a_{t}$

に対応し( $\beta_{1}$ をdriftと呼ぶことに由来するようだ)、

stockACF <- summary(ur.df(stocksub[,2],type="none"))

だと、

$\triangle y_{t} = \pi y_{t-1} + \psi_{1}^* \triangle y_{t-1} + \psi_{2}^* \triangle y_{t-2} + \cdots+\psi_{p-1}^* \triangle y_{t-p+1}+a_{t}$

に対応する。

それぞれやってみる。上側はtype="drift"の場合、下側はtype="none"の場合。(一部抜粋)

Call:
lm(formula = z.diff ~ z.lag.1 - 1)
Residuals:
   Min     1Q Median     3Q    Max 
-79.09 -10.96   0.00  12.92  65.00 
Coefficients:
        Estimate Std. Error t value Pr(>|t|)    
z.lag.1 -1.00673    0.06373   -15.8   <2e-16 ***

Call:
lm(formula = z.diff ~ z.lag.1 + 1)

Residuals:
    Min      1Q  Median      3Q     Max 
-80.379 -12.177  -1.246  11.633  63.754 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.24599    1.30384   0.956     0.34    
z.lag.1     -1.01026    0.06385 -15.822   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

最終的な結果はtype="trend"の場合と $y_{t-1}$ の係数が異なる。

まとめ

ur.df関数により、単位根を含むか否かの判定が出来る。単位根を含まないと判断が出来れば、その結果の次数応じてARIMAモデルを当てはめることが出来る。単位根が存在する場合は、階差をとって再びur.df関数を用いて単位根を含むかどうかの判定を行う。

残課題

ADF検定で推定しているモデルは単位根の項、トレンドの項、系列相関の項(ARモデル)を含んでいるが、MAモデルの項は含んでいないのは何故か。
type引数は、"trend", "drift", "none"の順番で検定を行い、いずれでも単位根の存在が棄却出来ない場合、階差系列を取り、再び同じプロセスで単位根の検定を行う。単位根の存在が棄却出来た段階で、その時点の階差次数dに対して、ARIMA(p,d,q)を適用する、という理解だが、確信が無い。