渋谷駅前で働くデータサイエンティストのブログ

元祖「六本木で働くデータサイエンティスト」です / 道玄坂→銀座→東京→六本木→渋谷駅前

何も考えずに線形回帰すると怖いので、計量時系列分析でダメ押ししてみる

何気なく読んでいて、途中で「?」と思った記事がありました。


何をやっているのかなー、と思って読み進めていったら一つ引っ掛かるところが。まず、この特集で扱っているのは「気温」と「電力消費(の日次最大値)」という時系列データなんですよね。なのに、4ページ目で普通に線形単回帰してます。


時系列をプロットしたのを眺めれば、どう見たって互いに相関しているのは丸分かりなのでどう計算しても構わないなんていうのは一目瞭然なんですが、それでも手法の説明のところで「時系列分析(ARMA / ARIMA)」とか言っているので、もうちょっと色々その辺を踏まえた何かがあっても良いのかなと思ったのでした。ということで、いつもながらRで見てみようと思います。


必要なRパッケージ


{tseries}, {vars}をインストールして展開して下さい。記事中で使われているデータですが、GitHubに多少直したものをテキストで上げてあります。適当に読み込んで、ここではdataとでも名付けておきます。


記事中でやっていることを再現してみる


まずデータをプロットしてみます。ちなみにこのやり方は@氏のブログ記事を参考にしました。大仏様、ありがとう!

> 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)


f:id:TJO:20130826121517p:plain


次に、記事中にもあるように普通に電力消費データを気温データで単回帰してみましょう。

> 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が比較的大きい回帰関係が見出されます。ちなみに、一階階差した気温と電力消費の時系列データをプロットするとこんな感じになります。


f:id:TJO:20130826144738p:plain


後の方(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っぽくなっていて微妙だと思ったんですが(笑)、このままにしておきます。すいません。

*1:最後の「おまけ」のセクションも参照のこと

*2:共和分関係が存在する場合には一階階差VAR表現は存在しないことになっている

*3:さらについでに{MSwM}とかでレジーム転換のチェックを行うと、結構良い感じの合致度のモデルで尚且つその「ズーンと跳ね上がる」ところでレジームが変わっているという結果になります

*4:あれだけ計量時系列分析推しでシリーズ記事書いてる奴が何言ってんだって感じですが