日本橋濱町Weblog(日々酔亭)

Quality Economic Analyses Produces Winning Markets

R言語:系列相関の除去・・・FGLSでの推定

さて、回帰分析(lm)を実行して、ダービンワトソン検定を行ったら見事に系列相関ありという結果となったところまでが前回まで。
まず系列相関があると何がまずいのかだが・・・一致性は担保されるが、有効性(効率性)が担保されないということだった。有効性とは、最小分散であるということだが、最小分散が満たされないとどういう問題が起こるのであろうか。
Cov \(u_t, u_s | \bf X) \neq 0
この式、共分散がゼロという仮定が成り立たないのが系列相関がある場合・・・この時、t統計量、F統計量、LM統計量を用いての仮設検定を行うことが妥当ではなくなるということだ。なぜ、共分散行列がゼロではないとこれらの検定が妥当でなくなるのかは別途。
今回はこの系列相関を取り除く作業をRで行う。
教科書には2通りの系列相関の除去の方法が出ているが、今回はFGLSでの推定を採用*1
モデルは以下の通り・・・って、前に示したもの・・・同じだ。
net.expenditure = f(income, income^2, trend, m02, m04, m06, m07, m08, m09, m10, d1403)
関数型はこんな感じ。あえて書く必要もないかもしれないが、これもTeXで書いてみた(長すぎると全体が表示されない。これも今後の課題)。
前回と異なっているのは、対数線形モデルではなく、線形モデルになっているという点。
net.expenditure_t = \alpha + \beta_1_t \cdot income + \beta_2_t \cdot income^2 + \beta_3_t \cdot trend + \beta_4_t \cdot m02 + \beta_5_t \cdot m04 + \beta_6_t \cdot m06 + \beta_7_t \cdot m07 + \beta_8_t \cdot m08 + \beta_9_t \cdot m09 + \beta_10_t \cdot m10 + \beta_11_t \cdot d1403)
これでいいのか?久しぶりのこういう世界・・・表現方法が違ってたりして。
コマンドは以下の通り。
> results_gls <- gls(net.expenditure ~ income + incomesq+ trend + m02 + m02 + m04 + m06 + m07 + m08 + m09 + m10 + d1403, corr = corARMA(p=1), data = nedatats)
前半部分と最後のオブジェクトの指定はは回帰分析(lm)の時と同じだが、一階の自己相関に従うという条件を入れないといけない。それが、corr = corARMA(p=1)の部分。
これで推定すると以下の結果を得る。

Generalized least squares fit by REML
Model: net.expenditure ~ income + incomesq + trend + m02 + m04 + m06 + m07 + m08 + m09 + m10 + d1403
Data: nedatats
AIC BIC logLik
1988.81 2029.69 -980.4051

Correlation Structure: AR(1)
Formula: ~1
Parameter estimate(s):
Phi
0.3812407

Coefficients:
Value Std.Error t-value p-value
(Intercept) -2202.6800 676.9292 -3.25393 0.0014
income 0.0114 0.0027 4.19858 0.0000
incomesq 0.0000 0.0000 -3.60734 0.0004
trend 35.2387 0.6352 55.47597 0.0000
m02 -398.8778 62.2992 -6.40261 0.0000
m04 -277.7479 66.4095 -4.18235 0.0001
m06 -918.0920 139.7650 -6.56882 0.0000
m07 -503.3101 117.5292 -4.28243 0.0000
m08 -237.0749 80.8551 -2.93210 0.0039
m09 -305.2581 69.9174 -4.36598 0.0000
m10 -329.8766 73.9142 -4.46296 0.0000
d1403 1647.6444 195.7975 8.41505 0.0000

Correlation:
(Intr) income incmsq trend m02 m04 m06 m07 m08 m09 m10
income -0.996
incomesq 0.993 -0.999
trend -0.217 0.148 -0.146
m02 0.462 -0.473 0.472 -0.073
m04 0.552 -0.563 0.562 -0.092 0.306
m06 0.872 -0.879 0.870 -0.127 0.445 0.523
m07 0.789 -0.804 0.800 -0.112 0.427 0.496 0.832
m08 0.457 -0.479 0.480 -0.063 0.298 0.336 0.517 0.636
m09 -0.015 -0.008 0.012 0.003 0.081 0.077 0.074 0.171 0.432
m10 0.507 -0.523 0.523 -0.079 0.301 0.345 0.506 0.510 0.443 0.380
d1403 0.042 -0.042 0.043 -0.083 0.118 0.116 0.049 0.057 0.057 0.042 0.050

Standardized residuals:
Min Q1 Med Q3 Max

  • 2.64971947 -0.61196371 -0.03390092 0.51201045 3.08647732

Residual standard error: 222.6958
Degrees of freedom: 149 total; 137 residual

推定できたので、推定結果を再現してみたいのだが、incomesqが下4桁しか表示してくれない・・・0.0000・・・これ。
実際は0ではなく、所得の二乗ということで非常に大きな数値のため、小数点以下8桁か9桁ぐらいなのだがこれを直接取り出す方法がわからない。
そこで止まってしまった。
一応、ここまで系列相関がある場合、FGLSで推定すればよいというところの再推定までできたことになる。実際に系列相関が除去できているかを推定結果から得られる残差をプロットして確認してみたいが、そこはincomesqの値がないとできないので、そこを取り出す方法を見つけないといけない・・・そうか、なんてことはない・・・R言語上でグラフ化すればよいのでした。これは次回。

*1:実はもう一方のやり方はエラーが出て結果を得られなかった。これについては別途原因を調べる予定。