About connecting the dots.

data science related trivial things

みせかけの回帰 (3)

前回までで,みせかけの回帰に関する基礎概念の説明を終えました.今回は,実際にRを使ってみせかけの回帰の再現と,単位根検定の実施をしていきたいと思います.そのために,あらかじめ{tseries}パッケージを入れておきましょう.

MSCIデータ

データの読み込みとプロット

まず,定番の経済データとして,ネタ本の沖本先生のWebサイトでダウンロードできる,msci_day.xlsを読み込んでみましょう.これはG7先進国の2003年から2008年のMSCI指数です.エクセルそのままだと読み込めないので,いったんテキストファイルに持ってきてあげます。その上で,日本とアメリカとイギリスのデータをグラフにプロットしてみましょう.

# MSCIインデックス(日本および米国)
msci <- read.table("~/20130502_msci.txt", sep="\t", header=T)

# MSCIグラフプロット
span <- as.Date("2003-01-01") + 1:1391
par(new=F)
par(xaxt="n")
plot(span, msci$jp, type = "l", col="red", xlab="Time", ylab="Score", ylim=c(500, 4000))
par(new=T)
plot(span, msci$us, type = "l", col="blue", xlab="Time", ylab="Score", ylim=c(500, 4000))
par(new=T)
plot(span, msci$uk, type = "l", col="green", xlab="Time", ylab="Score", ylim=c(500, 4000))
par(xaxt="s")
axis.Date(1, at=seq(min(span), max(span), "years"))
legend("topleft", legend=c("jp", "usa", "uk"), col=c("red", "blue", "green"), lty=c(1, 1, 1))

f:id:SAM:20130504104649j:plain

一般的にこの手の経済指標は単位根過程に従うことが多いわけですが,たしかにグラフを見ると原系列は非定常ですね.3国とも時間経過とともにスコアが上昇トレンドといえます.では,1次差分をとってみると,どうなるでしょう.

# MSCIグラフプロット(1次差分)
span <- as.Date("2003-01-01") + 1:1390
par(new=F)
par(xaxt="n")
plot(span, diff(msci$jp), type = "l", col="red", xlab="Time", ylab="Score", ylim=c(-200, 200))
par(new=T)
plot(span, diff(msci$us), type = "l", col="blue", xlab="Time", ylab="Score", ylim=c(-200, 200))
par(new=T)
plot(span, diff(msci$uk), type = "l", col="green", xlab="Time", ylab="Score", ylim=c(-200, 200))
par(xaxt="s")
axis.Date(1, at=seq(min(span), max(span), "years"))
legend("topleft", legend=c("jp", "usa", "uk"), col=c("red", "blue", "green"), lty=c(1, 1, 1))

f:id:SAM:20130504104704j:plain

ぱっとみてわかるレベルで,定常ですね.ということで,単位根過程とみて問題なさそうです.

単位根検定

では実際に単位根過程なのかどうかを,前回説明したADF検定で示してみましょう.ADF検定の実施法は以下の通りです.

> adf.test(msci$jp)

	Augmented Dickey-Fuller Test

data:  msci$jp
Dickey-Fuller = -1.5432, Lag order = 11, p-value = 0.7717
alternative hypothesis: stationary

> adf.test(msci$us)

	Augmented Dickey-Fuller Test

data:  msci$us
Dickey-Fuller = -2.2458, Lag order = 11, p-value = 0.4742
alternative hypothesis: stationary

> adf.test(msci$uk)

	Augmented Dickey-Fuller Test

data:  msci$uk
Dickey-Fuller = -2.1418, Lag order = 11, p-value = 0.5183
alternative hypothesis: stationary

前回述べた通り,ADF検定は帰無仮説を,ターゲットの時系列データが単位根過程である,としているので,これが棄却されない(p値が一定基準を下回らない)場合に,ターゲットが単位根過程であるという帰無仮説を棄却できないということになります.

今回の場合,すべての指標でp値はとても大きな値をとっており,仮に基準のp値を0.05としたところで,これを大きく上回っており,どれも単位根過程と考えるのが妥当といえるでしょう.

みせかけの回帰

では,実際に上記の変数群を使って回帰を行った場合,みせかけの回帰が生じるのでしょうか.

> summary(lm(msci$jp~msci$us))

Call:
lm(formula = msci$jp ~ msci$us)

Residuals:
    Min      1Q  Median      3Q     Max 
-455.37 -141.26  -35.19  109.25  636.45 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -975.13069   37.17272  -26.23   <2e-16 ***
msci$us        3.08110    0.03182   96.83   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 201.2 on 1389 degrees of freedom
Multiple R-squared:  0.871,	Adjusted R-squared:  0.8709 
F-statistic:  9375 on 1 and 1389 DF,  p-value: < 2.2e-16

ということで,アメリカのMSCIから日本のMSCIを予測した場合,きれいに高いr-squaredと,有意な回帰係数が得られました.さて,ではこのような場合に実際どうすればみせかけの回帰を回避できるでしょうか.簡単なのは,変数の差分を取ってあげることです.念のため差分のADF検定も行っておきましょう.

> adf.test(diff(msci$jp))

	Augmented Dickey-Fuller Test

data:  diff(msci$jp)
Dickey-Fuller = -11.8327, Lag order = 11, p-value = 0.01
alternative hypothesis: stationary

Warning message:
In adf.test(diff(msci$jp)) : p-value smaller than printed p-value
> adf.test(diff(msci$us))

	Augmented Dickey-Fuller Test

data:  diff(msci$us)
Dickey-Fuller = -11.2011, Lag order = 11, p-value = 0.01
alternative hypothesis: stationary

Warning message:
In adf.test(diff(msci$us)) : p-value smaller than printed p-value
> adf.test(diff(msci$uk))

	Augmented Dickey-Fuller Test

data:  diff(msci$uk)
Dickey-Fuller = -11.8295, Lag order = 11, p-value = 0.01
alternative hypothesis: stationary

Warning message:
In adf.test(diff(msci$uk)) : p-value smaller than printed p-value

ということで,今度は検定結果をみると,すべてp-valueはすべて0.01未満ということですので,一般的な観点から帰無仮説は棄却され,各差分時系列は単位根過程ではないとみなせるでしょう.となると,差分同士で回帰分析を行った場合どうなるでしょうか.

> summary(lm(diff(msci$jp)~diff(msci$us)))

Call:
lm(formula = diff(msci$jp) ~ diff(msci$us))

Residuals:
     Min       1Q   Median       3Q      Max 
-187.721  -18.613    0.032   18.443  127.053 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept)     1.0004     0.9014   1.110    0.267
diff(msci$us)   0.0652     0.0895   0.729    0.466

Residual standard error: 33.58 on 1388 degrees of freedom
Multiple R-squared:  0.0003823,	Adjusted R-squared:  -0.0003379 
F-statistic: 0.5308 on 1 and 1388 DF,  p-value: 0.4664

ということで,先ほどみられていた,極めて高い相関関係はばっちり消え去って,まったくの無相関ということになりました.まさしくみせかけの回帰,ということになりますね.

いろいろなテストデータ

ということで,これに加えて様々な形状の時系列データで,実際にどのようなスコアが得られるのかをみていきたいと思います.

データの作成

arima.simでARIMA過程が作り出せますが,ここではわかりやすくするためAR(1)を3個ほど作成してみたいと思います.同様にMA(1)も作っておきましょう.続いてrnorm()を用いて正規分布に従う乱数を発生させて,これを時系列データとします.同様にして,1次関数,2次関数,sin波から生成した時系列データ,以上を用いて同様のプロット〜単位根検定を行ってみたいと思います.ちなみに各関数を生成する際の微妙な係数は,グラフプロットで見やすくするために適当に調整したもので,深い意味はありません.

# AR(1)
ar1_1 <- arima.sim(list(order=c(1, 1, 0), ar=0.8), n=99)
ar1_2 <- arima.sim(list(order=c(1, 1, 0), ar=0.8), n=99)
ar1_3 <- arima.sim(list(order=c(1, 1, 0), ar=0.8), n=99)

# MA(1)
ma1_1 <- arima.sim(list(order=c(0, 1, 1), ma=0.8), n=99)

# 正規分布
gd <- as.ts(rnorm(100)*10)

# 1次関数
f <- function(x) 0.5*x+4
linear <- as.ts(f(c(1:100)))

# 2次関数
f <- function(x) -0.02*(x+10)*(x-80)
quadric <- as.ts(f(c(1:100)))

# sin波
wave <- as.ts(sin(c(1:100)/3)*50)

グラフのプロット

作成したデータと,その1次差分をプロットしてみました.明らかに単位根過程っぽいものと,全然違いそうなもの,さらにはどっちなのかよくわからないものが混在していますね.

# 作成した変数のグラフプロット
par(new=F)
plot(ar1_1, type = "l", col="red", xlab="Time", ylab="Score", ylim=c(-100, 100))
par(new=T)
plot(ar1_2, type = "l", col="orange", xlab="Time", ylab="Score", ylim=c(-100, 100))
par(new=T)
plot(ar1_3, type = "l", col="brown", xlab="Time", ylab="Score", ylim=c(-100, 100))
par(new=T)
plot(ma1_1, type = "l", col="blue", xlab="Time", ylab="Score", ylim=c(-100, 100))
par(new=T)
plot(gd, type = "l", col="green", xlab="Time", ylab="Score", ylim=c(-100, 100))
par(new=T)
plot(linear, type = "l", col="purple", xlab="Time", ylab="Score", ylim=c(-100, 100))
par(new=T)
plot(quadric, type = "l", col="darkblue", xlab="Time", ylab="Score", ylim=c(-100, 100))
par(new=T)
plot(wave, type = "l", col="pink", xlab="Time", ylab="Score", ylim=c(-100, 100))
legend("topleft",
       legend=c("AR(1) 1", "AR(1) 2", "AR(1) 3", "MA(1)", "Gaussian", "Linear", "Quadric", "Sin wave"),
       col=c("red", "orange", "brown", "blue", "green", "purple", "darkblue", "pink"),
       lty=c(1, 1, 1, 1, 1, 1, 1, 1))

# 作成した変数のグラフプロット(差分)
par(new=F)
plot(diff(ar1_1), type = "l", col="red", xlab="Time", ylab="Score", ylim=c(-50, 50))
par(new=T)
plot(diff(ar1_2), type = "l", col="orange", xlab="Time", ylab="Score", ylim=c(-50, 50))
par(new=T)
plot(diff(ar1_3), type = "l", col="brown", xlab="Time", ylab="Score", ylim=c(-50, 50))
par(new=T)
plot(diff(ma1_1), type = "l", col="blue", xlab="Time", ylab="Score", ylim=c(-50, 50))
par(new=T)
plot(diff(gd), type = "l", col="green", xlab="Time", ylab="Score", ylim=c(-50, 50))
par(new=T)
plot(diff(linear), type = "l", col="purple", xlab="Time", ylab="Score", ylim=c(-50, 50))
par(new=T)
plot(diff(quadric), type = "l", col="darkblue", xlab="Time", ylab="Score", ylim=c(-50, 50))
par(new=T)
plot(diff(wave), type = "l", col="pink", xlab="Time", ylab="Score", ylim=c(-50, 50))
legend("topleft",
       legend=c("AR(1) 1", "AR(1) 2", "AR(1) 3", "MA(1)", "Gaussian", "Linear", "Quadric", "Sin wave"),
       col=c("red", "orange", "brown", "blue", "green", "purple", "darkblue", "pink"),
       lty=c(1, 1, 1, 1, 1, 1, 1, 1))

f:id:SAM:20130504104724j:plain
f:id:SAM:20130508002133j:plain

単位根検定

ということで,先ほどと同様に単位根検定を行ってみましょう.

> adf.test(ar1_1)

	Augmented Dickey-Fuller Test

data:  ar1_1
Dickey-Fuller = -1.1402, Lag order = 4, p-value = 0.9124
alternative hypothesis: stationary

> adf.test(ar1_2)

	Augmented Dickey-Fuller Test

data:  ar1_2
Dickey-Fuller = -3.0059, Lag order = 4, p-value = 0.1599
alternative hypothesis: stationary

> adf.test(ar1_3)

	Augmented Dickey-Fuller Test

data:  ar1_3
Dickey-Fuller = -1.9909, Lag order = 4, p-value = 0.5804
alternative hypothesis: stationary

> adf.test(ma1_1)

	Augmented Dickey-Fuller Test

data:  ma1_1
Dickey-Fuller = -2.9419, Lag order = 4, p-value = 0.1865
alternative hypothesis: stationary

> adf.test(gd)

	Augmented Dickey-Fuller Test

data:  gd
Dickey-Fuller = -4.3002, Lag order = 4, p-value = 0.01
alternative hypothesis: stationary

Warning message:
In adf.test(gd) : p-value smaller than printed p-value
> adf.test(linear)

	Augmented Dickey-Fuller Test

data:  linear
Dickey-Fuller = 1.7321, Lag order = 4, p-value = 0.99
alternative hypothesis: stationary

Warning message:
In adf.test(linear) : p-value greater than printed p-value
> adf.test(quadric)

	Augmented Dickey-Fuller Test

data:  quadric
Dickey-Fuller = 0.3704, Lag order = 4, p-value = 0.99
alternative hypothesis: stationary

Warning message:
In adf.test(quadric) : p-value greater than printed p-value
> adf.test(wave)

	Augmented Dickey-Fuller Test

data:  wave
Dickey-Fuller = -4.552979e+14, Lag order = 4, p-value = 0.01
alternative hypothesis: stationary

Warning message:
In adf.test(wave) : p-value smaller than printed p-value

結果をみてみると,正規分布とsin波のみが帰無仮説を棄却されましたが,それ以外はすべて単位根検定という判定がなされました.正規分布とsin波は,そもそも原系列からして明らかに定常過程なので,当然といえば当然かもしれませんね.

みせかけの回帰

ということで,最後にみせかけの回帰が生じるかどうか,みてみましょう.従属変数には,単位根過程と示されたar1_1を使って,他の変数を独立変数にして回帰を行ってみたいと思います.ちょっと結果が冗長になってしまうのはご勘弁ください.

> summary(lm(ar1_1~ar1_2))

Call:
lm(formula = ar1_1 ~ ar1_2)

Residuals:
    Min      1Q  Median      3Q     Max 
-27.760 -18.095  -3.478  19.462  30.272 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  15.9312     4.5492   3.502 0.000697 ***
ar1_2         0.7528     0.2727   2.761 0.006885 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 18.83 on 98 degrees of freedom
Multiple R-squared:  0.07216,	Adjusted R-squared:  0.06269 
F-statistic: 7.622 on 1 and 98 DF,  p-value: 0.006885

> summary(lm(ar1_1~ar1_3))

Call:
lm(formula = ar1_1 ~ ar1_3)

Residuals:
   Min     1Q Median     3Q    Max 
-37.40 -14.70  10.06  14.74  21.98 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  33.2140     2.3168  14.336  < 2e-16 ***
ar1_3         0.4624     0.1144   4.044 0.000105 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 18.1 on 98 degrees of freedom
Multiple R-squared:  0.143,	Adjusted R-squared:  0.1342 
F-statistic: 16.35 on 1 and 98 DF,  p-value: 0.0001049

> summary(lm(ar1_1~ma1_1))

Call:
lm(formula = ar1_1 ~ ma1_1)

Residuals:
    Min      1Q  Median      3Q     Max 
-33.084 -13.422   2.425  10.831  38.403 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  20.6680     2.0037  10.315  < 2e-16 ***
ma1_1         2.9642     0.4909   6.038  2.8e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 16.69 on 98 degrees of freedom
Multiple R-squared:  0.2712,	Adjusted R-squared:  0.2637 
F-statistic: 36.46 on 1 and 98 DF,  p-value: 2.797e-08

> summary(lm(ar1_1~gd))

Call:
lm(formula = ar1_1 ~ gd)

Residuals:
    Min      1Q  Median      3Q     Max 
-32.239 -19.073   1.822  16.069  28.424 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  27.2544     1.9594  13.910   <2e-16 ***
gd           -0.1272     0.2091  -0.608    0.544    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 19.51 on 98 degrees of freedom
Multiple R-squared:  0.003764,	Adjusted R-squared:  -0.006402 
F-statistic: 0.3703 on 1 and 98 DF,  p-value: 0.5443

> summary(lm(ar1_1~linear))

Call:
lm(formula = ar1_1 ~ linear)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.1295  -2.7227  -0.1812   2.5837   8.0285 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -11.03271    0.89877  -12.28   <2e-16 ***
linear        1.31271    0.02756   47.64   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.977 on 98 degrees of freedom
Multiple R-squared:  0.9586,	Adjusted R-squared:  0.9582 
F-statistic:  2270 on 1 and 98 DF,  p-value: < 2.2e-16

> summary(lm(ar1_1~quadric))

Call:
lm(formula = ar1_1 ~ quadric)

Residuals:
    Min      1Q  Median      3Q     Max 
-30.378 -10.685   2.599  11.477  19.625 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 38.54710    1.78456  21.600   <2e-16 ***
quadric     -0.58766    0.05934  -9.904   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 13.82 on 98 degrees of freedom
Multiple R-squared:  0.5002,	Adjusted R-squared:  0.4951 
F-statistic: 98.09 on 1 and 98 DF,  p-value: < 2.2e-16

> summary(lm(ar1_1~wave))

Call:
lm(formula = ar1_1 ~ wave)

Residuals:
    Min      1Q  Median      3Q     Max 
-29.849 -18.509   1.819  16.528  25.972 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 27.45283    1.95329  14.055   <2e-16 ***
wave        -0.03991    0.05475  -0.729    0.468    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 19.49 on 98 degrees of freedom
Multiple R-squared:  0.005391,	Adjusted R-squared:  -0.004758 
F-statistic: 0.5312 on 1 and 98 DF,  p-value: 0.4678

ということで,奇麗に単位根過程と示された変数はすべて有意な回帰係数とr-squared,それに対して単位根過程ではないと示された正規分布とsin波は有意ではない値が得られました.

まとめ

ということで,いろんなデータを使ってみせかけの回帰を行って,実際のところをみてみました.またみせかけの回帰を回避するための,差分を取った場合の回帰も行って,きちんと回避できていることが示せました.本当はこの上に共和分の検討をする必要があるのですが,今回の話の範囲からやや外れる(のと,きちんとやろうとするとかなり応用的な部分まで踏み込まないといけない)ので,今回は省略.ということでおしまいです.