何気なく読んでいて、途中で「?」と思った記事がありました。
何をやっているのかなー、と思って読み進めていったら一つ引っ掛かるところが。まず、この特集で扱っているのは「気温」と「電力消費(の日次最大値)」という時系列データなんですよね。なのに、4ページ目で普通に線形単回帰してます。
時系列をプロットしたのを眺めれば、どう見たって互いに相関しているのは丸分かりなのでどう計算しても構わないなんていうのは一目瞭然なんですが、それでも手法の説明のところで「時系列分析(ARMA / ARIMA)」とか言っているので、もうちょっと色々その辺を踏まえた何かがあっても良いのかなと思ったのでした。ということで、いつもながらRで見てみようと思います。
必要なRパッケージ
{tseries}, {vars}をインストールして展開して下さい。記事中で使われているデータですが、GitHubに多少直したものをテキストで上げてあります。適当に読み込んで、ここではdataとでも名付けておきます。
記事中でやっていることを再現してみる
まずデータをプロットしてみます。ちなみにこのやり方は@teramonagi氏のブログ記事を参考にしました。大仏様、ありがとう!
> plot(data$temp,lwd=2,type="l",col="red") > par(new=T) > plot(data$power,lwd=2,type="l",col="blue",axes=F) > axis(4)
次に、記事中にもあるように普通に電力消費データを気温データで単回帰してみましょう。
> summary(lm(power~temp,data)) Call: lm(formula = power ~ temp, data = data) Residuals: Min 1Q Median 3Q Max -537.52 -133.84 -15.36 176.80 481.02 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 952.508 198.962 4.787 1.14e-05 *** temp 104.072 7.285 14.286 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 227.4 on 60 degrees of freedom Multiple R-squared: 0.7728, Adjusted R-squared: 0.769 F-statistic: 204.1 on 1 and 60 DF, p-value: < 2.2e-16
ということで、気温が電力消費に影響を与えているらしいということが分かります。偏回帰係数のp値も十分に低いし、モデルのp値も十分に低いし、R^2値も調整済みで0.769と十分にでかい。言うことない感じです。
でもこれって時系列データなのだから...
ただ、やっぱりちょっと気になりました。最初に時系列分析の話をしているわけですから、その後にこんな時系列データが出てきたら時系列分析しないわけにはいかないでしょー、と。そこでadf.test(){tseries}を使って単位根検定してみました。
> adf.test(temp) Augmented Dickey-Fuller Test data: temp Dickey-Fuller = -2.9079, Lag order = 3, p-value = 0.2074 alternative hypothesis: stationary > adf.test(power) Augmented Dickey-Fuller Test data: power Dickey-Fuller = -3.2881, Lag order = 3, p-value = 0.08182 alternative hypothesis: stationary
p値の大きさ的には微妙なラインですが*1、このADF検定の結果から言えば気温データも電力消費データも単位根過程ということになります。
ということは、お互いに単位根過程なので、普通に計算すると見せかけの回帰を踏んでしまうと思うんですが、どうなんでしょうか? 少なくとも、互いの回帰(相関)関係を過大評価している可能性は否定できないような気がします。
一つ目の解決策:一階階差を取ってから回帰する
これはもはや見せかけの回帰の回避策で言うところの鉄板ですね。株価の世界では見せかけの回帰が頻発するので、「一階階差を取ってから回帰する」というのがルーチンになっているところもあるそうです。ということでやってみましょう。
> d_data<-as.data.frame(diff(ts(data),lag=1)) # diff()を使いたかったので一旦時系列に直すというゴミっぷり > summary(lm(power~temp,d_data)) Call: lm(formula = power ~ temp, data = d_data) Residuals: Min 1Q Median 3Q Max -783.92 -90.70 -2.53 63.74 674.05 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6.146 25.274 0.243 0.809 temp 38.875 8.435 4.609 2.22e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 197.1 on 59 degrees of freedom Multiple R-squared: 0.2647, Adjusted R-squared: 0.2522 F-statistic: 21.24 on 1 and 59 DF, p-value: 2.218e-05
R^2が大幅に下がりますが、これでも十分に有意でR^2が比較的大きい回帰関係が見出されます。ちなみに、一階階差した気温と電力消費の時系列データをプロットするとこんな感じになります。
後の方(50日目ぐらいのところ)に極端に相関が強くなっているところがまだ残ってますが、残りはまぁ穏当な感じがしますね。散布図は面倒なのでここでは割愛します。
もう一つの解決策:VARモデルで回帰関係を見る
見せかけの回帰が懸念される状況で、変数間の関係性を調べるためのもう一つの方法がVARモデリングです。これも以前の記事で取り上げた方法ですね。ただし、共和分があるとまずいので先にca.jo(){urca}で調べておきます。
> dts<-ts(data) > dts.vecm<-ca.jo(dts,type="eigen",ecdet="trend",K=2,spec="transitory") > summary(dts.vecm) ###################### # Johansen-Procedure # ###################### Test type: maximal eigenvalue statistic (lambda max) , with linear trend in cointegration Eigenvalues (lambda): [1] 3.625518e-01 1.373237e-01 -6.261268e-18 Values of teststatistic and critical values of test: test 10pct 5pct 1pct r <= 1 | 8.86 10.49 12.25 16.26 r = 0 | 27.02 16.85 18.96 23.65 # p < 0.05でr = 0ということ Eigenvectors, normalised to first column: (These are the cointegration relations) temp.l1 power.l1 trend.l1 temp.l1 1.00000000 1.00000000 1.000000000 power.l1 -0.00777772 0.07035145 -0.003611764 trend.l1 0.01950542 -1.53983860 -0.610249967 Weights W: (This is the loading matrix) temp.l1 power.l1 trend.l1 temp.d -0.9869833 -0.02703177 2.065799e-16 power.d -2.4793758 -4.01793803 1.445091e-14
ということで、共和分関係は存在しないことが分かりました。なら安心して普通に生VARでも一階階差VARでも推定してOK*2ということで、まず生VARを推定してみましょう。
> VARselect(dts,lag.max=5,type="trend") $selection AIC(n) HQ(n) SC(n) FPE(n) 1 1 1 1 $criteria 1 2 3 4 AIC(n) 12.09845 12.21901 12.29833 12.42976 HQ(n) 12.18203 12.35831 12.49335 12.68050 SC(n) 12.31351 12.57744 12.80013 13.07494 FPE(n) 179627.79862 202787.92444 219876.80075 251472.72611 5 AIC(n) 12.53215 HQ(n) 12.83861 SC(n) 13.32070 FPE(n) 279834.99586 > dts.var<-VAR(dts,p=1,type="trend") > summary(dts.var) VAR Estimation Results: ========================= Endogenous variables: temp, power Deterministic variables: trend Sample size: 61 Log Likelihood: -539.972 Roots of the characteristic polynomial: 0.9973 0.04856 Call: VAR(y = dts, p = 1, type = "trend") Estimation results for equation temp: ===================================== temp = temp.l1 + power.l1 + trend Estimate Std. Error t value Pr(>|t|) temp.l1 0.039771 0.163263 0.244 0.808 power.l1 0.006916 0.001179 5.868 2.25e-07 *** # ココに注目 trend 0.001048 0.021341 0.049 0.961 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.429 on 58 degrees of freedom Multiple R-Squared: 0.9925, Adjusted R-squared: 0.9922 F-statistic: 2574 on 3 and 58 DF, p-value: < 2.2e-16 Estimation results for equation power: ====================================== power = temp.l1 + power.l1 + trend Estimate Std. Error t value Pr(>|t|) temp.l1 -1.2166 15.5934 -0.078 0.938 power.l1 1.0061 0.1126 8.937 1.68e-12 *** # ココに注目 trend 0.4652 2.0383 0.228 0.820 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 232 on 58 degrees of freedom Multiple R-Squared: 0.9965, Adjusted R-squared: 0.9963 F-statistic: 5443 on 3 and 58 DF, p-value: < 2.2e-16 Covariance matrix of residuals: temp power temp 5.897 363.3 power 363.341 53750.7 Correlation matrix of residuals: temp power temp 1.0000 0.6454 power 0.6454 1.0000
生VARモデルでは、気温・電力消費ともに相互にラグ次数1で有意に回帰しているということが分かりました。ちなみに、ついでなので一階階差VARモデルを推定した上でGranger因果も見てみましょう。
> d_dts<-diff(dts,lag=1) > plot(d_dts) > VARselect(d_dts,lag.max=5) $selection AIC(n) HQ(n) SC(n) FPE(n) 1 1 1 1 $criteria 1 2 3 4 AIC(n) 12.51636 12.57999 12.62192 12.73335 HQ(n) 12.60049 12.72021 12.81823 12.98574 SC(n) 12.73336 12.94166 13.12826 13.38435 FPE(n) 272818.42663 290960.02670 303928.70799 340771.50760 5 AIC(n) 12.69615 HQ(n) 13.00463 SC(n) 13.49183 FPE(n) 329884.44597 > d_dts.var<-VAR(d_dts,p=1) > summary(d_dts.var) VAR Estimation Results: ========================= Endogenous variables: temp, power Deterministic variables: const Sample size: 60 Log Likelihood: -543.158 Roots of the characteristic polynomial: 0.519 0.09515 Call: VAR(y = d_dts, p = 1) Estimation results for equation temp: ===================================== temp = temp.l1 + power.l1 + const Estimate Std. Error t value Pr(>|t|) temp.l1 -0.507708 0.139542 -3.638 0.000592 *** power.l1 0.003031 0.001840 1.647 0.105007 const 0.193946 0.360079 0.539 0.592245 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.784 on 57 degrees of freedom Multiple R-Squared: 0.1893, Adjusted R-squared: 0.1608 F-statistic: 6.653 on 2 and 57 DF, p-value: 0.002529 Estimation results for equation power: ====================================== power = temp.l1 + power.l1 + const Estimate Std. Error t value Pr(>|t|) temp.l1 1.5306 11.6497 0.131 0.896 power.l1 -0.1064 0.1536 -0.693 0.491 const 15.3058 30.0611 0.509 0.613 Residual standard error: 232.5 on 57 degrees of freedom Multiple R-Squared: 0.009533, Adjusted R-squared: -0.02522 F-statistic: 0.2743 on 2 and 57 DF, p-value: 0.7611 Covariance matrix of residuals: temp power temp 7.753 376.7 power 376.664 54037.6 Correlation matrix of residuals: temp power temp 1.0000 0.5819 power 0.5819 1.0000 > causality(d_dts.var,cause="temp") $Granger Granger causality H0: temp do not Granger-cause power data: VAR object d_dts.var F-Test = 0.0173, df1 = 1, df2 = 114, p-value = 0.8957 # Granger因果は認められない $Instant H0: No instantaneous causality between: temp and power data: VAR object d_dts.var Chi-squared = 15.1781, df = 1, p-value = 9.783e-05 # 同時影響は強い
という感じで、日次データでは気温⇒電力消費というようなGranger因果は存在せず、同時影響が強い(≒相関が強い)という割と当たり前の結果になりました。特に最初の単回帰の結果から大きくは乖離してませんね、ということで。
結論:いきなり回帰しても今回は実害はなかったけど、時系列データなんだからそれなりに気を付けるべき
そういうわけで、今回のデータに関して言えばいきなり単回帰してもそんなに問題はなかったらしい、ということになりました。まぁ、夏場の気温と電力消費なんていう「どう見ても相関するに決まっている」データを使っているので、当たり前っちゃ当たり前ですが。
なのですが、記事の冒頭で時系列分析云々ARMA/ARIMA云々と書いているからには、やっぱり時系列データの扱いにはもう少し神経を使っても良かったのでは?とも思ってしまいます。せめて時系列データではない、もっとありふれたクロスセクションデータでやった方が無難だったのでは?という気も。。。
という、超絶お節介を呟いたところで今回はおしまい。お粗末さまでした~。
おまけ
実はPP検定だと今回の単位根検定の結果は変わるんですよね。試しにpp.test(){tseries}でやってみるとこうなります。
> pp.test(temp) Phillips-Perron Unit Root Test data: temp Dickey-Fuller Z(alpha) = -34.5444, Truncation lag parameter = 3, p-value = 0.01 # 対立仮説が「単位根過程ではない」である点に注意 alternative hypothesis: stationary Warning message: In pp.test(temp) : p-value smaller than printed p-value > pp.test(power) Phillips-Perron Unit Root Test data: power Dickey-Fuller Z(alpha) = -20.2483, Truncation lag parameter = 3, p-value = 0.04765 # 上に同じ alternative hypothesis: stationary
なお、この問題については群馬大青木先生のサイト内の「統計学関連なんでもあり」の過去ログ内に詳しく解説したやり取りを見つけましたので、興味のある方はこちらもご覧ください。
ちなみに、今回用いた時系列は途中でズーンと跳ね上がるフェーズがある*3ので、そこを避けた方が良いのでは?という気もします。実際にその分をカットして改めて単位根検定をかけてみると、
> adf.test(temp[1:45]) Augmented Dickey-Fuller Test data: temp[1:45] Dickey-Fuller = -3.5968, Lag order = 3, p-value = 0.04417 alternative hypothesis: stationary > adf.test(power[1:45]) Augmented Dickey-Fuller Test data: power[1:45] Dickey-Fuller = -2.8804, Lag order = 3, p-value = 0.2246 alternative hypothesis: stationary
ということになって少なくとも気温データは単位根ではないらしいということになるので、そのまま単回帰しても良いだろうという結論に達します。こういうこともあるので、時系列データを扱う時は色々神経を使う必要があって面倒だなぁと思ってます*4。
何かこの記事のタイトルが裏RjpWikiっぽくなっていて微妙だと思ったんですが(笑)、このままにしておきます。すいません。