みせかけの回帰 (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))
一般的にこの手の経済指標は単位根過程に従うことが多いわけですが，たしかにグラフを見ると原系列は非定常ですね．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))
ぱっとみてわかるレベルで，定常ですね．ということで，単位根過程とみて問題なさそうです．
単位根検定
では実際に単位根過程なのかどうかを，前回説明した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))
単位根検定
ということで，先ほどと同様に単位根検定を行ってみましょう．
> 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波は有意ではない値が得られました．
まとめ
ということで，いろんなデータを使ってみせかけの回帰を行って，実際のところをみてみました．またみせかけの回帰を回避するための，差分を取った場合の回帰も行って，きちんと回避できていることが示せました．本当はこの上に共和分の検討をする必要があるのですが，今回の話の範囲からやや外れる（のと，きちんとやろうとするとかなり応用的な部分まで踏み込まないといけない）ので，今回は省略．ということでおしまいです．
