Rを使って計量経済分析
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
検索
|
最終更新
|
ヘルプ
|
ログイン
]
開始行:
|このページについて|
''概要'':このページは、統計言語Rで計量経済分析する方法に...
関連するサイトは、[[リンク:http://www.sugi-shun.com/econw...
''親ページ'':このページの親ページは[[R]]です。
|目次|
#contents
----
*古典的線形回帰モデル [#yed9bf5c]
**最小二乗法による推定 [#u30dda11]
とりあえず、RでOLSする方法を学びます。
例として、こんな問題を考えて見ましょう。
|投資関数は利子率の減少関数って習ったけど、これって本当？|
いや、あくまで例ですよ。
以下のデータを使います。
#ref(inv.csv)
このデータは、日本の企業設備投資と国内総支出と長期貸し出...
見出しの意味は、
|企業設備投資|国内総支出|Long-term Prime Lending Rate|
|INV | GDE |R|
です。
とりあえず、
> read.table("http://sugi-shun.com/econwiki/index.php?pl...
と入力して、Rにこのデータを読み込ませましょう。COLOR(RED)...
次に、
> attach(i)
と入力しましょう。こうすると、見出しの変数名を読み込んで...
> INV
などと打つと、企業設備投資のデータが出力されることを確認...
***単回帰 [#b8aa6a65]
とりあえず、単回帰してみましょう。
投資(Inv)と国内支出(GDE)の対前年同期比伸び率をまず求め、...
|INVの対前年同期比伸び率 | GDE の対前年同期比伸び率 |
|I|Y|
どうやったらいいかとというと、以下の関数を使います。
q.growth <- function(x) {
n <- length(x)
x1 <- x[-(1:4)]
x4 <- x[-((n-3):n)]
z1 <- (x1-x4)/x4*100
z1
}
この関数をつかって、
> q.growth(INV)->I
> q.growth(GDE)->Y
と入力すれば、四半期の対前年同期比伸び率データをつくれま...
さて次に単回帰分析してみましょう。その前に、それぞれの変...
> length(R)
[1] 109
> length(I)
[1] 105
という結果からもわかるとおり、このままだと変数ごとに長さ...
> R[5:length(R)]->r
と入力しておきましょう。これで長さがそろました。
COLOR(RED){（注意）：Rでは、大文字と小文字を区別します。...
次に、単回帰式&mimetex(I_t=a+bY_t);, &mimetex(I_t=a+br_t...
まず、&mimetex(I_t=a+bY_t);を推定してみましょう。
> lm(I ~ Y)
と入力してみましょう。すると、
Call:
lm(formula = I ~ Y)
Coefficients:
(Intercept) Y
-2.894 1.947
と出力されます。
このコマンドの意味は、
> lm(被説明変数 ~ 説明変数)
です。結果を見ればわかるとおり、定数項ありモデルを推定し...
> lm(I ~ 0 + Y)
と入力します。
この結果は、&mimetex(I_t=-2.894+1.947Y_t);と推定されたこ...
> summary(lm(I ~ Y))
と入力します。
すると、
Call:
lm(formula = I ~ Y)
Residuals:
Min 1Q Median 3Q Max
-15.130 -4.498 0.697 5.139 14.235
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.8943 1.1516 -2.513 0.0135 *
Y 1.9465 0.2659 7.321 5.62e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1...
Residual standard error: 6.195 on 103 degrees of freedom
Multiple R-Squared: 0.3422, Adjusted R-squared: 0.33...
F-statistic: 53.59 on 1 and 103 DF, p-value: 5.623e-11
という結果が返ってきます。
では次に、&mimetex(I_t=a+br_t);を推定してみましょう。
> summary(lm(I ~ r))
と入力すれば、
Call:
lm(formula = I ~ r)
Residuals:
Min 1Q Median 3Q Max
-16.5949 -4.6235 0.0347 6.0004 13.9655
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.0113 2.9102 2.409 0.0178 *
r -0.3836 0.3954 -0.970 0.3343
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1...
Residual standard error: 7.604 on 103 degrees of freedom
Multiple R-Squared: 0.009055, Adjusted R-squared: -0.0...
F-statistic: 0.9411 on 1 and 103 DF, p-value: 0.3343
と出力されます。
尚，どうして伸び率データを使うかというと，いま使っている...
***重回帰分析 [#f9aaed03]
単回帰のやり方を学んだので、次に重回帰のやり方を覚えまし...
たとえば、&mimetex(I_t=a+bY_t+cr_t);を推定することを考え...
以下のように入力してやりましょう。
> summary(lm(I ~ Y + r))
意味は、
> lm(被説明変数 ~ 説明変数その１ + 説明変数その２)
ということです。説明変数を複数にする場合、""+""でつなげれ...
出力結果は以下のようになります。
Call:
lm(formula = I ~ Y + r)
Residuals:
Min 1Q Median 3Q Max
-14.6601 -4.6180 0.4353 4.5796 16.2559
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.8428 2.4228 0.761 0.4486
Y 2.0229 0.2633 7.683 9.86e-12 ***
r -0.7051 0.3190 -2.211 0.0293 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1...
Residual standard error: 6.081 on 102 degrees of freedom
Multiple R-Squared: 0.3723, Adjusted R-squared: 0.36
F-statistic: 30.25 on 2 and 102 DF, p-value: 4.838e-11
なお、
X<-cbind(Y,r)
と説明変数行列を定義した上で（cbindは、列方向にべクトルを...
summary(lm(I ~ X))
とやっても同じ結果が得られます。これは便利です。
summary(lm(I ~0+ X))
で定数項なしモデルだって対応してくれます。
ちなみに、lm()とは、Linear Modelの略です。lm()はstatsパッ...
***予測 [#ae2821b1]
演習として、予測もやってみましょう。
> predict(lm(I ~ Y + r))
とすると、Iの予測値（理論値）を出力できます。
さらに
> predict(lm(I ~ Y + r),interval="confidence",level=0.95)
とすると、予測値の95%信頼区間も計算できます。
例えば、第1番目のIの理論値（すなわち、1971年第1四半期の理...
> predict(lm(I ~ Y + r),interval="confidence",level=0.95...
fit lwr upr
5.646562 4.104640 7.188485
として出てきます。
単回帰の場合の予測をしてみましょう。
predict(lm(I ~ Y ))->yosoku1
plot(Y,yosoku1)
とやれば、ちゃんと理論値はOLS推定された直線を表してくれる...
**ここまでの推定結果の解釈 [#f089e914]
+まず、投資関数とかいう名前で習うような&mimetex(I_t=a+br_...
+もう一つの単回帰式&mimetex(I_t=a+bY_t);については、&mime...
+重回帰式&mimetex(I_t=a+bY_t+cr_t);の推定結果を見ると、&m...
+さらに、&mimetex(r_t);の係数の大きさを見てみると、重回帰...
**より良い計量モデルを考えてみる [#m8485961]
[[ここまでの推定結果の解釈:http://www.sugi-shun.com/econw...
「企業が今年実行する設備投資計画は、去年、計画したものだ...
というわけで、&mimetex(I_t=a+br_t);ではなく、&mimetex(I_t...
まず、"&mimetex(r_{t-4});"に対応する変数を定義しないとい...
> R[1:(length(R)-4)]->r1
そして、回帰分析してみましょう。
> summary(lm(I ~ r1))
Call:
lm(formula = I ~ r1)
Residuals:
Min 1Q Median 3Q Max
-17.4018 -3.6172 0.1533 5.2735 14.1096
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.9412 3.1914 4.055 9.76e-05 ***
r1 -1.1820 0.4244 -2.785 0.00637 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1...
Residual standard error: 7.366 on 103 degrees of freedom
Multiple R-Squared: 0.07003, Adjusted R-squared: 0.061
F-statistic: 7.756 on 1 and 103 DF, p-value: 0.006371
ほら、やっぱり予想通りです。t値は有意ですし、係数の符号も...
さきほどは、利子率だけでなく、GDPも回帰式に入れてみたら、...
> summary(lm(I ~ Y + r1))
Call:
lm(formula = I ~ Y + r1)
Residuals:
Min 1Q Median 3Q Max
-13.8819 -3.4922 -0.1570 4.1544 15.5820
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.8212 2.5194 3.104 0.00247 **
Y 2.0903 0.2444 8.554 1.27e-13 ***
r1 -1.5349 0.3280 -4.679 8.88e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1...
Residual standard error: 5.649 on 102 degrees of freedom
Multiple R-Squared: 0.4585, Adjusted R-squared: 0.44...
F-statistic: 43.18 on 2 and 102 DF, p-value: 2.597e-14
ほら、t値はさらに有意になりました。符合条件も、理論通りで...
**仮説検定 [#i997fb8b]
RIGHT:参考：浅野・中村(2000)[pp.78-80]、"[[Econometrics i...
ここでは、LSE(Linear Squares Estimation)における仮説検定...
***lm()が自動出力するt検定量とF検定量を自分で計算してみよ...
lm()が自動的にやってくれる仮説検定は、係数が有意に0になる...
waldtest
パッケージ：lmtest
というコマンドがあります。Wald検定はこれを使うはずなんで...
COLOR(Fuchsia){（保留）：waldtest（パッケージ：lmtest）の...
代わりのコマンドを紹介します。
linear.hypothesis()
パッケージ：car
もしくは、
lht()
パッケージ：car
を使う方法を記述します。
> summary(lm(I ~ Y + r1))
を例にとりましょう。すでに、"H0:回帰係数はすべて0"が真の...
F-statistic: 43.18 on 2 and 102 DF, p-value: 2.597e-14
という結果を得ています。これと同じ結果を、自分で計算して...
ここでは、"[[Econometrics in R:http://cran.r-project.org/...
次のように入力します。
> unrestricted<-lm(I ~ Y + r1)
> rhs<-c(0,0)
> hm <- rbind(c(0,1,0),c(0,0,1))
> lht(unrestricted,hm,rhs)
一般に、線形制約式は、''Rβ''=''r''とかけます。ここで、''R...
rhsが''r''に対応します。hmが''R''に対応します。
結果は、
Hypothesis:
Y = 0
r1 = 0
Model 1: I ~ Y + r1
Model 2: restricted model
Res.Df RSS Df Sum of Sq F Pr(>F)
1 102 3254.5
2 104 6009.8 -2 -2755.4 43.178 2.597e-14 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1...
となりました。つまり、"H0:Yの係数が0かつr1の係数が0"は棄...
ちゃんとF統計量がlm()が返した値と一緒になっていることが確...
同じようにt検定を行うことが出来ます。例えば、Yのt値を自分...
> rhs2<-c(0)
> hm2 <- rbind(c(0,1,0))
> lht(unrestricted,hm2,rhs2,test=c("Chisq"))
と入力しましょう。こうすると、制約がたった一つで"H0:Yの係...
Linear hypothesis test
Hypothesis:
Y = 0
Model 1: I ~ Y + r1
Model 2: restricted model
Res.Df RSS Df Sum of Sq Chisq Pr(>Chisq)
1 102 3254.5
2 103 5589.0 -1 -2334.5 73.166 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1...
となります。ここでわれわれが興味があるのは、&mimetex(\chi...
> lht(unrestricted,hm2,rhs2,test=c("Chisq"))$Chi[2]
とやればよいです。これは自由度が1の&mimetex(\chi^2);分布...
RIGHT:参考：Greene(2003)[pp.851]
平方根を計算してやると、
> sqrt(lht(unrestricted,hm2,rhs2,test=c("Chisq"))$Chi[2])
[1] 8.553713
を得ますが、これはYの係数のt値にほかなりません。
ここまでやってきたことで、実は次のようなことがわかりまし...
***WALD検定 [#n7d120c5]
さて、すでに述べたとおり、Wald test, LM test, 制約付残差...
RIGHT:参考：Greene(2003)[pp.853-854]
結果のみ示すと、
&mimetex(\chi^2);分布に従う統計量 = 制約の数 × F分布に...
となります。
先ほどは、何も考えずにF検定量が出ました。これは、Rコマン...
> lht(unrestricted,hm,rhs,test=c("Chisq"))
とします。すると
Linear hypothesis test
Hypothesis:
Y = 0
r1 = 0
Model 1: I ~ Y + r1
Model 2: restricted model
Res.Df RSS Df Sum of Sq Chisq Pr(>Chisq)
1 102 3254.5
2 104 6009.8 -2 -2755.4 86.357 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1...
と出力されます。いま制約の数は2ですから、確かに&mimetex(\...
*回帰診断 [#j659cc8a]
OLS推定したら、「この回帰分析は成功なのか、失敗なのか」と...
以下の表の「判断方法」に記された検定などの結果、全ての仮...
ここで工夫とは、「古典的仮定が満たされていなかった場合、...
その上で、どう工夫しても改善できなかった場合、表の一番右...
|仮定 |内容 |判断方法|対処法|
|仮定1|説明変数は非確率的である(or,説明変数と撹乱項は直交...
|仮定2|撹乱項の期待値は0である。 |残差の期待値はいつも0...
|仮定3|撹乱項の分散は一定である。 |Breush Pagan test, Wh...
|仮定4|撹乱項は自己相関しない。 |Durbin-Watson, Breusc...
|仮定5|撹乱項は正規分布に従う。 |残差に対する各種の正...
表では、古典的仮定については、[[浅野・中村(2000):http://w...
なお、仮定１～仮定４のすべてが崩れた場合は、GMMという最終...
それでは、実際にこれらのチェックをしてみましょう。
とりあえず、&mimetex(I_t=a+bY_t+Cr_{t-4});というモデリン...
**説明変数の外生性のチェック [#fd0a5ff1]
hausman.systemfit {systemfit}
を使います．
COLOR(Fuchsia){（保留）：（説明変数と撹乱項の直交条件）。...
**分散均一性のチェック [#t8c26c13]
***図を描いて目視 [#x372ae75]
RIGHT:参考：浅野・中村(2000)[pp.132]
横軸に説明変数、縦軸に残差をプロットしてみて、目でまずは...
で、残差をプロットする方法を説明します。その前に、残差系...
> ?lm
などと入力してみましょう。すると、lmについてのヘルプがぽ...
Value:
という項目があります。これは、「lm()というコマンドを打っ...
residuals
という変数がありまます。たとえば、回帰分析したときの残差...
> lm(I ~ Y + r1)$residuals
のように、$(ドルマーク）で繋いで、入力します。すると、残...
Rは賢く、
> lm(I ~ Y + r1)$res
でも同じ結果を返します。Value一覧のなかで、アルファベット...
> lm(I ~ Y + r1)$r
とすれば、
NULL
となります。これは、Valueのうち、residualsではなく、rank...
> names(lm(I ~ Y + r1))
と入力しても、Value一覧を見ることができます。これを見ると...
> names(summary(lm(I ~ Y + r1)))
と打ってみましょう。はい、ちゃんとありましたね。
さらに余談を続けます。「r1のt値だけとりだす」にはどうした...
とりえあず、
> summary(lm(I ~ Y + r1))$coefficients
と入力してみましょう。すると、
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.821201 2.5193943 3.104397 2.468089e-03
Y 2.090334 0.2443774 8.553713 1.268153e-13
r1 -1.534889 0.3280416 -4.678946 8.883377e-06
と出力されます。これは、(3*4)の行列とRは認識しているよう...
> summary(lm(I ~ Y + r1))$coefficients[3,3]
とすればよいことがわかります。
ちなみに、「F検定のp値だけとりだす」にはどうしたらいいで...
> summary(lm(I ~ Y + r1))$fstatistic
と入力してみましょう。
value numdf dendf
43.17844 2.00000 102.00000
と出力されますが、ここにはなぜかp値がない！
苦肉の策として、
> BAG<-summary(lm(I ~ Y + r1))$fstatistic
> pf(BAG[1], BAG[2], BAG[3], lower.tail=FALSE)
value
2.597481e-14
が考えられますが、どうもいけてないですね、これは。RjpWiki...
さて、話をもとに戻します。横軸に説明変数、縦軸に残差をプ...
> plot(r1,lm(I ~ Y + r1)$res)
あるいは、
> plot(Y,lm(I ~ Y + r1)$res)
とすればよいです。この図の結果を見てみると、まぁ微妙。分...
***Breush Pagan test [#i533be5c]
RIGHT:参考：浅野・中村(2000)[pp.135]
コマンドは、
bptest()
(パッケージ：lmtest)
です。
具体的に、&mimetex(I_t=a+bY_t+Cr_{t-4});というモデルにお...
> bptest(I ~ Y + r1)
あるいは、
> bptest(lm(I ~ Y + r1))
> bptest(summary((lm(I ~ Y + r1))))
でも同じ結果が返されるようです。
COLOR(RED){（注意）：パッケージの読み込みを忘れずに。}
これを実行すると、
studentized Breusch-Pagan test
data: summary((lm(I ~ Y + r1)))
BP = 1.8808, df = 2, p-value = 0.3905
という結果が返ってきます。よって、"H0:分散は均一である"は...
COLOR(Fuchsia){（保留）：studentizedってなんですか？疑問...
***White Test [#o13a8fec]
RIGHT:参考：浅野・中村(2000)[pp.135]
コマンドは、
white.test()
(パッケージ：tseries)
です。
具体的に、&mimetex(I_t=a+bY_t+cr_{t-4});というモデルにお...
***Goldfeld-Quandt test [written by [[Akihiko Noda:http:/...
RIGHT:参考：浅野・中村(2000)[pp.134]
コマンドは、
gqtest()
(パッケージ：lmtest)
です。
基本的な考え方としては、存在するデータを2つに分割して、そ...
具体的に、&mimetex(I_t=a+bY_t+Cr_{t-4});というモデルにお...
> gqtest(I ~ Y + r1)
あるいは、
> gqtest(lm(I ~ Y + r1))
> gqtest(summary((lm(I ~ Y + r1))))
でも同じ結果が返されるようです。
COLOR(RED){（注意）：パッケージの読み込みを忘れずに。}
これを実行すると、
Goldfeld-Quandt test
data: I ~ Y + r1
GQ = 0.6811, df1 = 50, df2 = 49, p-value = 0.9102
という結果が返ってきます。よって、"H0:2つの分散は等しい"...
**系列相関のチェック [#q9258075]
***図を描いて目視 [#i2020d58]
RIGHT:参考：浅野・中村(2000)[pp.146-147]
図の描き方として、二つあります。一つ目としては、横軸に残...
このプロットの仕方は、[[分散均一性のチェック:http://www.s...
> plot(lm(I ~ Y + r1)$res[1:(length(lm(I ~ Y + r1)$res)-...
と入力してみてください。ちょっとコマンドが複雑で見通しが...
> lm(I ~ Y + r1)$res->zansa
> plot(zansa[1:(length(zansa)-1)],zansa[2:length(zansa)])
と入力しても同じです。
プロットされた図は、円状ではありません。何か、右上がりの...
二つ目の図の描き方としては、横軸の時間(t),縦軸の残差をと...
> plot(c(1:length(zansa)),zansa)
とかやれば目的は果たせます。
> plot(zansa)
でも出来るみたいです。
出てきた図を見ると、ブリッジ上に反り返っているので、やは...
横軸に時間(t)をとる方法も説明しておきます。
> plot(ts(zansa,start=c(1971,1),frequency=4),type="p")
もしくは
> plot(ts(zansa,start=c(1971,1),frequency=4),type="o")
などと入力すればよいです。いま、もともとのデータの観測期...
そこで、これを統計的にちゃんと判断する方法を以下で説明し...
***Durbin-Watson [#e59ae9ad]
RIGHT:参考：浅野・中村(2000)[pp.148]
コマンドは、
dwtest()
(パッケージ：lmtest)
もしくは、
durbin.watson
(パッケージ：car)
です。
具体的に、&mimetex(I_t=a+bY_t+cr_{t-4});というモデルにお...
> dwtest(I ~ Y + r1)
あるいは、
> dwtest(lm(I ~ Y + r1))
> dwtest(summary((lm(I ~ Y + r1))))
でも同じ結果が返されるようです。
COLOR(RED){（注意）パッケージの読み込みを忘れずに。}
これを実行すると、
Durbin-Watson test
data: summary((lm(I ~ Y + r1)))
DW = 0.2924, p-value < 2.2e-16
alternative hypothesis: true autocorrelation is greater ...
という結果が返されます。したがって、"H0:系列相関はない"が...
しかし、DWテストでは正規分布のAR(1)しか想定していません。...
***Breusch-Godfrey test [#c12150d4]
RIGHT:参考：浅野・中村(2000)[pp.152]
コマンドは、
bgtest()
(パッケージ：lmtest)
です。
具体的に、&mimetex(I_t=a+bY_t+cr_{t-4});というモデルにお...
> bgtest(I ~ Y + r1)
あるいは、
> bgtest(lm(I ~ Y + r1))
> bgtest(summary((lm(I ~ Y + r1))))
でも同じ結果が返されるようです。
COLOR(RED){（注意）パッケージの読み込みを忘れずに。}
これを実行すると、
Breusch-Godfrey test for serial correlation of o...
data: lm(I ~ Y + r1)
LM test = 77.1659, df = 1, p-value < 2.2e-16
という結果が返されます。したがって、"H0:系列相関はない"が...
***Q Test(Box-Pierce, Ljung-Box test) [#e6d7a41c]
コマンドは
> Box.test()
です。[[管理人(SUGIYAMA Shunsuke):http://sugi-shun.com]]...
**RESETテスト [#j1869d35]
resettest {lmtest}
を使います．
COLOR(Fuchsia){（保留）：そのうち加筆予定。関数形が正しい...
**ここまでの回帰診断のまとめ [#f1a7ac05]
どうやら、分散不均一性は発生していないようです。しかし、...
しかし系列相関が発生してもOLS推定量は最小分散性を失うだけ...
COLOR(Fuchsia){（保留）：Hausman test, RESET についてその...
*どう工夫しても仮定が崩れるときの対処法 [#d23942c0]
[[回帰診断:http://www.sugi-shun.com/econwiki/index.php?R%...
仮定1,仮定2は、OLS推定量が一致性を持つために必要です。仮...
仮定3,仮定4は、OLS推定量が最小分散性（有効性＝効率性）を...
仮定5はなくても漸近的に何も問題ありません。直感的には、ど...
考慮すべきは、仮定1,仮定3,仮定4の三つということになります...
というわけで、古典的仮定が崩れたときの対処法をコンパクト...
|崩れる仮定|対処法|
|仮定1|2SLS(IV)|
|仮定3|HCSE(White), GLS(WLSなど)|
|仮定4|HCSE(Newey-West), GLS(Cochrane-Orcutt法など)|
|仮定3and仮定4|GLS|
|仮定1and仮定3and仮定4|GMM|
これらの対処法をRでどうやってやるのかを以下で説明したいと...
**2SLS(IV):2Step Least Squares(Instrument Variables)[#qa5...
RIGHT:参考：Wooldridge(2002)[Ch5]
日本語では二段階最小二乗法(操作変数法)のことです（IV≠TSLS...
例として、ここではWooldridge(2002)[Problem 5.4]をとりあげ...
#ref(card.csv)
です。このデータについての説明は
#ref(card.txt)
を見てください。
このデータは、[[Wooldridgeのウェブサイト:http://www.msu.e...
とりあえず、
> read.table("card.csv",header=T,sep=",")->c
> attach(c)
などとしてデータを読み込みましょう。
まず、考えたいモデルは以下であるとします。
lwage = A + B*educ + C*exper + D*expersq + E*black + F*so...
とりあえず、OLS推定してみましょう。
> summary(lm(lwage ~ educ + exper + expersq + black + so...
+ + smsa + reg661 + reg662+reg663+reg664+reg665+reg666+r...
その結果として、
Call:
lm(formula = lwage ~ educ + exper + expersq + black + so...
smsa + reg661 + reg662 + reg663 + reg664 + reg665 + ...
reg667 + reg668 + smsa66)
Residuals:
Min 1Q Median 3Q Max
-1.62326 -0.22141 0.02001 0.23932 1.33340
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.7393766 0.0715282 66.259 < 2e-16 ***
educ 0.0746933 0.0034983 21.351 < 2e-16 ***
exper 0.0848320 0.0066242 12.806 < 2e-16 ***
expersq -0.0022870 0.0003166 -7.223 6.41e-13 ***
black -0.1990123 0.0182483 -10.906 < 2e-16 ***
south -0.1479550 0.0259799 -5.695 1.35e-08 ***
smsa 0.1363846 0.0201005 6.785 1.39e-11 ***
reg661 -0.1185697 0.0388301 -3.054 0.002281 **
reg662 -0.0222026 0.0282575 -0.786 0.432092
reg663 0.0259703 0.0273644 0.949 0.342670
reg664 -0.0634942 0.0356803 -1.780 0.075254 .
reg665 0.0094551 0.0361174 0.262 0.793504
reg666 0.0219476 0.0400984 0.547 0.584182
reg667 -0.0005887 0.0393793 -0.015 0.988073
reg668 -0.1750058 0.0463394 -3.777 0.000162 ***
smsa66 0.0262417 0.0194477 1.349 0.177327
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1...
Residual standard error: 0.3723 on 2994 degrees of freedom
Multiple R-Squared: 0.2998, Adjusted R-squared: 0.29...
F-statistic: 85.48 on 15 and 2994 DF, p-value: < 2.2e-16
を得ます。
しかしこのような賃金関数の推定において、educ(教育年数)は...
このときOLS推定量は一致性を失いますから、上記の推定結果は...
そこで操作変数の登場です。ここでは、nearc4（4年制大学から...
+説明変数educとnearc4は相関するか？
+εとnearc4は無相関であるか？
の二つです（これは、厳密な言い方ではありません）。
***操作変数(nearc4)は，説明変数(educ)と相関するか？ [#g15...
まず一つ目のチェックをします。これは、より厳密には、
educ = a + b*nearc4 + c*exper + d*expersq + e*black + f*s...
というモデルにおいて、b≠0である、というのが正確です。
これをチェックしてみましょう。
> summary(lm(educ ~ nearc4 + exper + expersq + black + s...
+ + reg661 + reg662+reg663+reg664+reg665+reg666+reg667+r...
と入力し、
Call:
lm(formula = educ ~ nearc4 + exper + expersq + black + s...
smsa + reg661 + reg662 + reg663 + reg664 + reg665 + ...
reg667 + reg668 + smsa66)
Residuals:
Min 1Q Median 3Q Max
-7.54513 -1.36996 -0.09103 1.27836 6.23847
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 16.8485239 0.2111222 79.805 < 2e-16 ***
nearc4 0.3198989 0.0878638 3.641 0.000276 ***
exper -0.4125334 0.0336996 -12.241 < 2e-16 ***
expersq 0.0008686 0.0016504 0.526 0.598728
black -0.9355287 0.0937348 -9.981 < 2e-16 ***
south -0.0516126 0.1354284 -0.381 0.703152
smsa 0.4021825 0.1048112 3.837 0.000127 ***
reg661 -0.2102710 0.2024568 -1.039 0.299076
reg662 -0.2889073 0.1473395 -1.961 0.049992 *
reg663 -0.2382099 0.1426357 -1.670 0.095012 .
reg664 -0.0930890 0.1859827 -0.501 0.616742
reg665 -0.4828875 0.1881872 -2.566 0.010336 *
reg666 -0.5130857 0.2096352 -2.448 0.014442 *
reg667 -0.4270887 0.2056208 -2.077 0.037880 *
reg668 0.3136204 0.2416739 1.298 0.194490
smsa66 0.0254805 0.1057692 0.241 0.809644
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1...
Residual standard error: 1.941 on 2994 degrees of freedom
Multiple R-Squared: 0.4771, Adjusted R-squared: 0.47...
F-statistic: 182.1 on 15 and 2994 DF, p-value: < 2.2e-16
を得ます。nearc4のt値は3.641です。ちゃんと有意です。これ...
***操作変数(nearc4)は，撹乱項(ε)と関係しないか？ [#vd6434...
次にnearc4が操作変数としての条件2をクリアするか検討しまし...
educ(教育年数)⇔能力(観測されない)⇔撹乱項
なので、操作変数nearc4を使います。nearc4が撹乱項と直交条...
nearc4(大学から近いところに住んでいるか？)⇔能力(観測され...
において、nearc4と能力（の代理変数であるIQ)が相関しないか...
(能力とnearc4が相関していなくても、撹乱項とnearc4に相関が...
すなわち、IQがnearc4と関係しないということを確認しないと...
・・・確認しようと思ったら、iqに欠損値があるようで、これ...
(read.tableにはオプションにna.stringsがあって、読み込みの...
とりあえず，欠損値をすべて「空白」に置き換えたcsvファイル...
#ref(card2.csv)
空白にしておけば，Rに読み込むとこれをNA，すなわち欠損値と...
> read.table("card2.csv",header=T,sep=",")->c
> attach(c)
で改めてデータを読み込みましょう．
> summary(lm(nearc4~ iq))
とやります．結果は，以下です．
Call:
lm(formula = nearc4 ~ iq)
Residuals:
Min 1Q Median 3Q Max
-0.8021 -0.6690 0.2701 0.3016 0.4031
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.4773200 0.0671001 7.114 1.55e-12 ***
iq 0.0022555 0.0006477 3.483 0.000507 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1...
Residual standard error: 0.4534 on 2059 degrees of freedom
(949 observations deleted due to missingness)
Multiple R-Squared: 0.005856, Adjusted R-squared: 0.00...
F-statistic: 12.13 on 1 and 2059 DF, p-value: 0.0005071
nearc4 に対してiqで単回帰すると、iqのt値は有意になりまし...
> summary(lm(nearc4 ~ iq+smsa66+reg661+reg662+reg669))
とやれば，結果は，以下です．
Call:
lm(formula = nearc4 ~ iq + smsa66 + reg661 + reg662 + re...
Residuals:
Min 1Q Median 3Q Max
-0.9659 -0.3916 0.1208 0.2293 0.6265
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.3346827 0.0613386 5.456 5.45e-08 ***
iq 0.0006253 0.0005919 1.056 0.290884
smsa66 0.3747183 0.0199138 18.817 < 2e-16 ***
reg661 0.1433780 0.0414904 3.456 0.000560 ***
reg662 0.1745655 0.0241401 7.231 6.72e-13 ***
reg669 0.1029124 0.0308465 3.336 0.000864 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1...
Residual standard error: 0.4082 on 2055 degrees of freedom
(949 observations deleted due to missingness)
Multiple R-Squared: 0.1959, Adjusted R-squared: 0.19...
F-statistic: 100.1 on 5 and 2055 DF, p-value: < 2.2e-16
もともと，単回帰分析を考えていたわけではにので，このよう...
条件2もクリアした，ということが確認されました．
***2SLSする [#raaf0f44]
そこでいよいよ2SLSを行います。
使うコマンドは、
>tsls()
パッケージ:sem
です。尚，
EquationsModelling {fMultivar}
や
nlsystemfit {systemfit}
や
systemfit {systemfit}
などでは，サブコマンドでOLS,2SLS,WLS,SUR,3SLSなどが網羅的...
さて，
>tsls()
パッケージ:sem
の場合，基本的に、X1のための操作変数をZとした場合、
summary(tsls(Y ~ X1 + X2 + X3, ~ Z + X2 +X3 ))
とやります。操作変数が複数の場合は、例えばX1のための操作...
summary(tsls(Y ~ X1 + X2 + X3, ~ Z1 + Z2 + X2 +X3 ))
とやればよいです。
ここで．
#ref(card.csv)
のデータにまた戻ります．Rを再起動して，
> read.table("card2.csv",header=T,sep=",")->c
> attach(c)
で再びデータを読み込んでおきましょう．
以下のように入力しましょう。
> summary(tsls(lwage ~ educ + exper + expersq + black + ...
+ + smsa + reg661 + reg662+reg663+reg664+reg665+reg666+r...
+ +reg668+smsa66,
+ ~ nearc4 + exper + expersq + black + south + smsa + re...
+ + reg662+reg663+reg664+reg665+reg666+reg667+reg668+sms...
結果は以下のようになります。
2SLS Estimates
Model Formula: lwage ~ educ + exper + expersq + black + ...
reg662 + reg663 + reg664 + reg665 + reg666 + reg667 ...
smsa66
Instruments: ~nearc4 + exper + expersq + black + south +...
reg663 + reg664 + reg665 + reg666 + reg667 + reg668 ...
Residuals:
Min. 1st Qu. Median Mean 3rd Qu. M...
-1.83e+00 -2.41e-01 2.43e-02 -5.84e-12 2.52e-01 1.43e...
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.773966 0.9349469 4.0366 5.560e-05
educ 0.131504 0.0549637 2.3926 1.679e-02
exper 0.108271 0.0236586 4.5764 4.923e-06
expersq -0.002335 0.0003335 -7.0014 3.116e-12
black -0.146776 0.0538999 -2.7231 6.504e-03
south -0.144671 0.0272846 -5.3023 1.227e-07
smsa 0.111808 0.0316620 3.5313 4.197e-04
reg661 -0.107814 0.0418137 -2.5784 9.972e-03
reg662 -0.007046 0.0329073 -0.2141 8.305e-01
reg663 0.040445 0.0317806 1.2726 2.033e-01
reg664 -0.057917 0.0376059 -1.5401 1.236e-01
reg665 0.038458 0.0469387 0.8193 4.127e-01
reg666 0.055089 0.0526597 1.0461 2.956e-01
reg667 0.026758 0.0488287 0.5480 5.837e-01
reg668 -0.190891 0.0507113 -3.7643 1.702e-04
smsa66 0.018531 0.0216086 0.8576 3.912e-01
Residual standard error: 0.3883 on 2994 degrees of freedom
ここで、2SLS推定のほうが、OLS推定よりも標準誤差が大きいこ...
ここまでは操作変数がただ一つでした。このとき、IV推定量と2...
> summary(tsls(lwage ~ educ + exper + expersq + black + ...
+ + smsa + reg661 + reg662+reg663+reg664+reg665+reg666+...
+ +reg668+smsa66,
+ ~ nearc2+nearc4 + exper + expersq + black + south + s...
+ + reg662+reg663+reg664+reg665+reg666+reg667+reg668+sm...
とするだけです。簡単ですね。
2SLS Estimates
Model Formula: lwage ~ educ + exper + expersq + black + ...
reg662 + reg663 + reg664 + reg665 + reg666 + reg667 ...
smsa66
Instruments: ~nearc2 + nearc4 + exper + expersq + black ...
reg662 + reg663 + reg664 + reg665 + reg666 + reg667 ...
smsa66
Residuals:
Min. 1st Qu. Median Mean 3rd Qu. M...
-1.94e+00 -2.51e-01 1.93e-02 -3.00e-12 2.65e-01 1.47e...
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.3396875 0.8945377 3.733423 1.924e-04
educ 0.1570593 0.0525782 2.987155 2.839e-03
exper 0.1188149 0.0228061 5.209792 2.019e-07
expersq -0.0023565 0.0003475 -6.780907 1.433e-11
black -0.1232778 0.0521500 -2.363907 1.815e-02
south -0.1431945 0.0284448 -5.034120 5.084e-07
smsa 0.1007530 0.0315193 3.196547 1.405e-03
reg661 -0.1029760 0.0434224 -2.371497 1.778e-02
reg662 -0.0002287 0.0337943 -0.006766 9.946e-01
reg663 0.0469556 0.0326490 1.438194 1.505e-01
reg664 -0.0554084 0.0391828 -1.414098 1.574e-01
reg665 0.0515041 0.0475678 1.082752 2.790e-01
reg666 0.0699968 0.0533049 1.313139 1.892e-01
reg667 0.0390596 0.0497499 0.785119 4.324e-01
reg668 -0.1980371 0.0525350 -3.769623 1.667e-04
smsa66 0.0150626 0.0223360 0.674364 5.001e-01
Residual standard error: 0.4053 on 2994 degrees of freedom
以上で2SLSによる推定がRで出来るようになりました。
尚、上記の推定結果は、Woolridge(2002)[Problem 5.4]の出題...
***データに関する注意 [#iac47a73]
COLOR(RED){（注意）：データに関する注意．}
expersqとlwageについて、Wooldridgeが用意したデータを使う...
> exper2<-exper^2
> logwage<-log(wage)
とかやってRで加工した系列を使うと、微妙に異なる推定結果を...
具体的には，
#ref(card_3.csv)
のように，expersq, lwageが入っていない状態のデータのcsvを...
read.table("card_3.csv",header=T,sep=",")->c
attach(c)
library(sem)
exper2<-exper^2
logwage<-log(wage)
summary(tsls(logwage ~ educ + exper + exper2 + black + s...
+ smsa + reg661 + reg662+reg663+reg664+reg665+reg666+reg...
+reg668+smsa66,
~ nearc4 + exper + exper2 + black + south + smsa + reg661
+ reg662+reg663+reg664+reg665+reg666+reg667+reg668+smsa6...
とかやれば，
2SLS Estimates
Model Formula: logwage ~ educ + exper + exper2 + black +...
reg662 + reg663 + reg664 + reg665 + reg666 + reg667 ...
smsa66
Instruments: ~nearc4 + exper + exper2 + black + south + ...
reg663 + reg664 + reg665 + reg666 + reg667 + reg668 ...
Residuals:
Min. 1st Qu. Median Mean 3rd Qu. M...
-1.84e+00 -2.45e-01 2.58e-02 2.17e-11 2.54e-01 1.43e...
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.742871 0.9520939 3.9312 8.644e-05
educ 0.134309 0.0563767 2.3824 1.726e-02
exper 0.106333 0.0230529 4.6126 4.144e-06
exper2 -0.002209 0.0003321 -6.6495 3.483e-11
black -0.138156 0.0572038 -2.4152 1.579e-02
south -0.143702 0.0275365 -5.2186 1.926e-07
smsa 0.108207 0.0329625 3.2827 1.040e-03
reg661 -0.106913 0.0422370 -2.5313 1.142e-02
reg662 -0.006620 0.0332152 -0.1993 8.420e-01
reg663 0.042308 0.0322907 1.3102 1.902e-01
reg664 -0.058365 0.0378677 -1.5413 1.234e-01
reg665 0.040148 0.0476794 0.8420 3.998e-01
reg666 0.052881 0.0526310 1.0048 3.151e-01
reg667 0.026308 0.0491505 0.5352 5.925e-01
reg668 -0.191616 0.0511800 -3.7440 1.845e-04
smsa66 0.018182 0.0217943 0.8343 4.042e-01
Residual standard error: 0.3913 on 2994 degrees of freedom
となり，微妙に異なる値が得られている．
WooldridgeのウェブサイトでDLできるexpersqとlwageのデータ...
・・・このように，計量分析では，論文の追試をやると，微妙...
***補足 [#nf96aed0]
systemfit
パッケージ：systemfit
でも2SLS推定が出来るようです。2SLSだけでなく、3SLSとかも...
?systemfit
で見てください。
***過剰識別制約テスト [#u19330c8]
COLOR(Fuchsia){（保留）：操作変数の過剰識別制約テストの使...
**HCSE(Heteroskedasticity Consistene Standard Error) [#x0...
RIGHT:参考："[[Econometrics in R:http://cran.r-project.or...
***White[#i69bbbc1]
分散不均一が発生してもOLS推定量は不偏性（一致性）は満たし...
コマンドは、
vcovHC
パッケージ:sandwich
もしくは
hccm
パッケージ：car
です。どちらを使っても同じ結果が得られることをちゃんと確...
使うデータは、[[最小二乗法による推定:http://www.sugi-shun...
[[lm()が自動出力するt検定量とF検定量を自分で計算してみよ...
COLOR(RED){（注意）パッケージ：car, sandwichを読み込みま...
> unrestricted<-lm(I ~ Y + r1)
> rhs2<-c(0)
> hm2 <- rbind(c(0,1,0))
> lht(unrestricted,hm2,rhs2,test=c("Chisq"))
と入力すれば、
Linear hypothesis test
Hypothesis:
Y = 0
Model 1: I ~ Y + r1
Model 2: restricted model
Res.Df RSS Df Sum of Sq Chisq Pr(>Chisq)
1 102 3254.5
2 103 5589.0 -1 -2334.5 73.166 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1...
と出力されます。欲しい情報はカイ二乗統計量なので、
> lht(unrestricted,hm2,rhs2,test=c("Chisq"))$Chi[2]
です。t値はこれの平方根なのでsqrtをとります。
> sqrt(lht(unrestricted,hm2,rhs2,test=c("Chisq"))$Chi[2])
[1] 8.553713
そしてこれはYのt値の他ならないことをすでに見ました。
このt値は分散均一の仮定のもとで計算されたものです。この仮...
まず、修正した共分散行列を使って仮説検定する方法は
> lht(unrestricted,hm2,rhs2,vcov=vcovHC(unrestricted) )
あるいは
> lht(unrestricted,hm2,rhs2,vcov=hccm(unrestricted) )
です。欲しいのはカイ二乗統計量の平方根なので、
> lht(unrestricted,hm2,rhs2,vcov=hccm(unrestricted),test...
> sqrt(whitese)
もしくは
> lht(unrestricted,hm2,rhs2,vcov=vcovHC(unrestricted),te...
> sqrt(whitese2)
と入力します。hccm()でもvcovHC()でも
[1] 6.390953
が得られました。これがWhiteの一致標準誤差を使った場合のt...
尚、White修正にはいろいろなバージョンがあります。
詳しくは、
> ?hccm
> ?vcovHC
で見てください。ここの"type"というサブコマンドをどのよう...
hccm()を使った場合は以下です。
> sqrt(lht(unrestricted,hm2,rhs2,vcov=hccm(unrestricted,...
[1] 6.741862
> sqrt(lht(unrestricted,hm2,rhs2,vcov=hccm(unrestricted,...
[1] 6.644852
> sqrt(lht(unrestricted,hm2,rhs2,vcov=hccm(unrestricted,...
[1] 6.564739
> sqrt(lht(unrestricted,hm2,rhs2,vcov=hccm(unrestricted,...
[1] 6.390953
> sqrt(lht(unrestricted,hm2,rhs2,vcov=hccm(unrestricted,...
[1] 6.374456
vcovHC()を使った場合は以下です。
> sqrt(lht(unrestricted,hm2,rhs2,vcov=vcovHC(unrestricte...
[1] 6.741862
> sqrt(lht(unrestricted,hm2,rhs2,vcov=vcovHC(unrestricte...
[1] 6.741862
> sqrt(lht(unrestricted,hm2,rhs2,vcov=vcovHC(unrestricte...
[1] 6.644852
> sqrt(lht(unrestricted,hm2,rhs2,vcov=vcovHC(unrestricte...
[1] 6.564739
> sqrt(lht(unrestricted,hm2,rhs2,vcov=vcovHC(unrestricte...
[1] 6.390953
> sqrt(lht(unrestricted,hm2,rhs2,vcov=vcovHC(unrestricte...
[1] 6.374456
ここで、"hc0"(hccm()を使った場合),"HC"もしくは"HC0"（vcov...
尚、vcovを使った場合、typeとして"const"を指定できます。こ...
> sqrt(lht(unrestricted,hm2,rhs2,vcov=vcovHC(unrestricte...
[1] 8.553713
ここの項目は、"[[Econometrics in R:http://cran.r-project....
を参考にしました。
***Newey-West [#s56bb41a]
Newey-Westの方法は、White修正の「系列相関があっても大丈夫...
コマンドは、
NeweyWest
パッケージ：sandwich
です。
**GLS(Generalized Least Squares) [#sc894fcf]
RIGHT:参考：浅野・中村(2000)[pp.130]
日本語では一般化最小二乗法です。[[管理人(SUGIYAMA Shunsuk...
そのうち、自分でプログラムを書くなり、組み込み関数を探す...
以下は覚書です。関連しそうなのは、
lm.gls
パッケージ:MASS
gls
パッケージ:nlme
glsD
パッケージ:Design
あたりか？
***WLS:Weighted Least Squares [#m1633509]
RIGHT:参考：浅野・中村(2000)[pp.127-130]
日本語では、加重最小二乗法です。
そのうち、自分でプログラムを書くなり、組み込み関数を探す...
***Cochrane Orcutt法 [#k6400586]
RIGHT:参考：浅野・中村(2000)[pp.153]
そのうち、自分でプログラムを書くなり、組み込み関数を探す...
**GMM [#u08f350f]
[[GMM:http://www.sugi-shun.com/econwiki/index.php?R%A4%C7...
*時系列分析 [#mf5cad5a]
[[Rで時系列分析]]を見てください。
*計量ミクロ [#ced0ba2a]
[[Rで計量ミクロ]]を見てください。
*OLS以外の推定方法 [#x0baa89f]
[[RでOLS以外の推定方法]]を見てください。
*Reference [#h7afec48]
**書籍 [#h8a36313]
-[[浅野・中村(2000):http://www.amazon.co.jp/gp/product/46...
-[[Greene(2003):http://www.amazon.co.jp/gp/product/013110...
-[[Hamilton(1994):http://www.amazon.co.jp/gp/product/0691...
-[[Wooldridge(2002):http://www.amazon.co.jp/Econometric-A...
**リンク [#f4df0207]
-[[Econometrics in R:http://cran.r-project.org/doc/contri...
-[[R for economists:http://www.mayin.org/ajayshah/KB/R/R_...
-[[A Brief Guide to R for Beginners in Econometrics:http:...
-[[CRAN Task View: Computational Econometrics:http://cran...
-[[パッケージ一覧:http://cran.md.tsukuba.ac.jp/src/contri...
終了行:
|このページについて|
''概要'':このページは、統計言語Rで計量経済分析する方法に...
関連するサイトは、[[リンク:http://www.sugi-shun.com/econw...
''親ページ'':このページの親ページは[[R]]です。
|目次|
#contents
----
*古典的線形回帰モデル [#yed9bf5c]
**最小二乗法による推定 [#u30dda11]
とりあえず、RでOLSする方法を学びます。
例として、こんな問題を考えて見ましょう。
|投資関数は利子率の減少関数って習ったけど、これって本当？|
いや、あくまで例ですよ。
以下のデータを使います。
#ref(inv.csv)
このデータは、日本の企業設備投資と国内総支出と長期貸し出...
見出しの意味は、
|企業設備投資|国内総支出|Long-term Prime Lending Rate|
|INV | GDE |R|
です。
とりあえず、
> read.table("http://sugi-shun.com/econwiki/index.php?pl...
と入力して、Rにこのデータを読み込ませましょう。COLOR(RED)...
次に、
> attach(i)
と入力しましょう。こうすると、見出しの変数名を読み込んで...
> INV
などと打つと、企業設備投資のデータが出力されることを確認...
***単回帰 [#b8aa6a65]
とりあえず、単回帰してみましょう。
投資(Inv)と国内支出(GDE)の対前年同期比伸び率をまず求め、...
|INVの対前年同期比伸び率 | GDE の対前年同期比伸び率 |
|I|Y|
どうやったらいいかとというと、以下の関数を使います。
q.growth <- function(x) {
n <- length(x)
x1 <- x[-(1:4)]
x4 <- x[-((n-3):n)]
z1 <- (x1-x4)/x4*100
z1
}
この関数をつかって、
> q.growth(INV)->I
> q.growth(GDE)->Y
と入力すれば、四半期の対前年同期比伸び率データをつくれま...
さて次に単回帰分析してみましょう。その前に、それぞれの変...
> length(R)
[1] 109
> length(I)
[1] 105
という結果からもわかるとおり、このままだと変数ごとに長さ...
> R[5:length(R)]->r
と入力しておきましょう。これで長さがそろました。
COLOR(RED){（注意）：Rでは、大文字と小文字を区別します。...
次に、単回帰式&mimetex(I_t=a+bY_t);, &mimetex(I_t=a+br_t...
まず、&mimetex(I_t=a+bY_t);を推定してみましょう。
> lm(I ~ Y)
と入力してみましょう。すると、
Call:
lm(formula = I ~ Y)
Coefficients:
(Intercept) Y
-2.894 1.947
と出力されます。
このコマンドの意味は、
> lm(被説明変数 ~ 説明変数)
です。結果を見ればわかるとおり、定数項ありモデルを推定し...
> lm(I ~ 0 + Y)
と入力します。
この結果は、&mimetex(I_t=-2.894+1.947Y_t);と推定されたこ...
> summary(lm(I ~ Y))
と入力します。
すると、
Call:
lm(formula = I ~ Y)
Residuals:
Min 1Q Median 3Q Max
-15.130 -4.498 0.697 5.139 14.235
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.8943 1.1516 -2.513 0.0135 *
Y 1.9465 0.2659 7.321 5.62e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1...
Residual standard error: 6.195 on 103 degrees of freedom
Multiple R-Squared: 0.3422, Adjusted R-squared: 0.33...
F-statistic: 53.59 on 1 and 103 DF, p-value: 5.623e-11
という結果が返ってきます。
では次に、&mimetex(I_t=a+br_t);を推定してみましょう。
> summary(lm(I ~ r))
と入力すれば、
Call:
lm(formula = I ~ r)
Residuals:
Min 1Q Median 3Q Max
-16.5949 -4.6235 0.0347 6.0004 13.9655
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.0113 2.9102 2.409 0.0178 *
r -0.3836 0.3954 -0.970 0.3343
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1...
Residual standard error: 7.604 on 103 degrees of freedom
Multiple R-Squared: 0.009055, Adjusted R-squared: -0.0...
F-statistic: 0.9411 on 1 and 103 DF, p-value: 0.3343
と出力されます。
尚，どうして伸び率データを使うかというと，いま使っている...
***重回帰分析 [#f9aaed03]
単回帰のやり方を学んだので、次に重回帰のやり方を覚えまし...
たとえば、&mimetex(I_t=a+bY_t+cr_t);を推定することを考え...
以下のように入力してやりましょう。
> summary(lm(I ~ Y + r))
意味は、
> lm(被説明変数 ~ 説明変数その１ + 説明変数その２)
ということです。説明変数を複数にする場合、""+""でつなげれ...
出力結果は以下のようになります。
Call:
lm(formula = I ~ Y + r)
Residuals:
Min 1Q Median 3Q Max
-14.6601 -4.6180 0.4353 4.5796 16.2559
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.8428 2.4228 0.761 0.4486
Y 2.0229 0.2633 7.683 9.86e-12 ***
r -0.7051 0.3190 -2.211 0.0293 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1...
Residual standard error: 6.081 on 102 degrees of freedom
Multiple R-Squared: 0.3723, Adjusted R-squared: 0.36
F-statistic: 30.25 on 2 and 102 DF, p-value: 4.838e-11
なお、
X<-cbind(Y,r)
と説明変数行列を定義した上で（cbindは、列方向にべクトルを...
summary(lm(I ~ X))
とやっても同じ結果が得られます。これは便利です。
summary(lm(I ~0+ X))
で定数項なしモデルだって対応してくれます。
ちなみに、lm()とは、Linear Modelの略です。lm()はstatsパッ...
***予測 [#ae2821b1]
演習として、予測もやってみましょう。
> predict(lm(I ~ Y + r))
とすると、Iの予測値（理論値）を出力できます。
さらに
> predict(lm(I ~ Y + r),interval="confidence",level=0.95)
とすると、予測値の95%信頼区間も計算できます。
例えば、第1番目のIの理論値（すなわち、1971年第1四半期の理...
> predict(lm(I ~ Y + r),interval="confidence",level=0.95...
fit lwr upr
5.646562 4.104640 7.188485
として出てきます。
単回帰の場合の予測をしてみましょう。
predict(lm(I ~ Y ))->yosoku1
plot(Y,yosoku1)
とやれば、ちゃんと理論値はOLS推定された直線を表してくれる...
**ここまでの推定結果の解釈 [#f089e914]
+まず、投資関数とかいう名前で習うような&mimetex(I_t=a+br_...
+もう一つの単回帰式&mimetex(I_t=a+bY_t);については、&mime...
+重回帰式&mimetex(I_t=a+bY_t+cr_t);の推定結果を見ると、&m...
+さらに、&mimetex(r_t);の係数の大きさを見てみると、重回帰...
**より良い計量モデルを考えてみる [#m8485961]
[[ここまでの推定結果の解釈:http://www.sugi-shun.com/econw...
「企業が今年実行する設備投資計画は、去年、計画したものだ...
というわけで、&mimetex(I_t=a+br_t);ではなく、&mimetex(I_t...
まず、"&mimetex(r_{t-4});"に対応する変数を定義しないとい...
> R[1:(length(R)-4)]->r1
そして、回帰分析してみましょう。
> summary(lm(I ~ r1))
Call:
lm(formula = I ~ r1)
Residuals:
Min 1Q Median 3Q Max
-17.4018 -3.6172 0.1533 5.2735 14.1096
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.9412 3.1914 4.055 9.76e-05 ***
r1 -1.1820 0.4244 -2.785 0.00637 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1...
Residual standard error: 7.366 on 103 degrees of freedom
Multiple R-Squared: 0.07003, Adjusted R-squared: 0.061
F-statistic: 7.756 on 1 and 103 DF, p-value: 0.006371
ほら、やっぱり予想通りです。t値は有意ですし、係数の符号も...
さきほどは、利子率だけでなく、GDPも回帰式に入れてみたら、...
> summary(lm(I ~ Y + r1))
Call:
lm(formula = I ~ Y + r1)
Residuals:
Min 1Q Median 3Q Max
-13.8819 -3.4922 -0.1570 4.1544 15.5820
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.8212 2.5194 3.104 0.00247 **
Y 2.0903 0.2444 8.554 1.27e-13 ***
r1 -1.5349 0.3280 -4.679 8.88e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1...
Residual standard error: 5.649 on 102 degrees of freedom
Multiple R-Squared: 0.4585, Adjusted R-squared: 0.44...
F-statistic: 43.18 on 2 and 102 DF, p-value: 2.597e-14
ほら、t値はさらに有意になりました。符合条件も、理論通りで...
**仮説検定 [#i997fb8b]
RIGHT:参考：浅野・中村(2000)[pp.78-80]、"[[Econometrics i...
ここでは、LSE(Linear Squares Estimation)における仮説検定...
***lm()が自動出力するt検定量とF検定量を自分で計算してみよ...
lm()が自動的にやってくれる仮説検定は、係数が有意に0になる...
waldtest
パッケージ：lmtest
というコマンドがあります。Wald検定はこれを使うはずなんで...
COLOR(Fuchsia){（保留）：waldtest（パッケージ：lmtest）の...
代わりのコマンドを紹介します。
linear.hypothesis()
パッケージ：car
もしくは、
lht()
パッケージ：car
を使う方法を記述します。
> summary(lm(I ~ Y + r1))
を例にとりましょう。すでに、"H0:回帰係数はすべて0"が真の...
F-statistic: 43.18 on 2 and 102 DF, p-value: 2.597e-14
という結果を得ています。これと同じ結果を、自分で計算して...
ここでは、"[[Econometrics in R:http://cran.r-project.org/...
次のように入力します。
> unrestricted<-lm(I ~ Y + r1)
> rhs<-c(0,0)
> hm <- rbind(c(0,1,0),c(0,0,1))
> lht(unrestricted,hm,rhs)
一般に、線形制約式は、''Rβ''=''r''とかけます。ここで、''R...
rhsが''r''に対応します。hmが''R''に対応します。
結果は、
Hypothesis:
Y = 0
r1 = 0
Model 1: I ~ Y + r1
Model 2: restricted model
Res.Df RSS Df Sum of Sq F Pr(>F)
1 102 3254.5
2 104 6009.8 -2 -2755.4 43.178 2.597e-14 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1...
となりました。つまり、"H0:Yの係数が0かつr1の係数が0"は棄...
ちゃんとF統計量がlm()が返した値と一緒になっていることが確...
同じようにt検定を行うことが出来ます。例えば、Yのt値を自分...
> rhs2<-c(0)
> hm2 <- rbind(c(0,1,0))
> lht(unrestricted,hm2,rhs2,test=c("Chisq"))
と入力しましょう。こうすると、制約がたった一つで"H0:Yの係...
Linear hypothesis test
Hypothesis:
Y = 0
Model 1: I ~ Y + r1
Model 2: restricted model
Res.Df RSS Df Sum of Sq Chisq Pr(>Chisq)
1 102 3254.5
2 103 5589.0 -1 -2334.5 73.166 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1...
となります。ここでわれわれが興味があるのは、&mimetex(\chi...
> lht(unrestricted,hm2,rhs2,test=c("Chisq"))$Chi[2]
とやればよいです。これは自由度が1の&mimetex(\chi^2);分布...
RIGHT:参考：Greene(2003)[pp.851]
平方根を計算してやると、
> sqrt(lht(unrestricted,hm2,rhs2,test=c("Chisq"))$Chi[2])
[1] 8.553713
を得ますが、これはYの係数のt値にほかなりません。
ここまでやってきたことで、実は次のようなことがわかりまし...
***WALD検定 [#n7d120c5]
さて、すでに述べたとおり、Wald test, LM test, 制約付残差...
RIGHT:参考：Greene(2003)[pp.853-854]
結果のみ示すと、
&mimetex(\chi^2);分布に従う統計量 = 制約の数 × F分布に...
となります。
先ほどは、何も考えずにF検定量が出ました。これは、Rコマン...
> lht(unrestricted,hm,rhs,test=c("Chisq"))
とします。すると
Linear hypothesis test
Hypothesis:
Y = 0
r1 = 0
Model 1: I ~ Y + r1
Model 2: restricted model
Res.Df RSS Df Sum of Sq Chisq Pr(>Chisq)
1 102 3254.5
2 104 6009.8 -2 -2755.4 86.357 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1...
と出力されます。いま制約の数は2ですから、確かに&mimetex(\...
*回帰診断 [#j659cc8a]
OLS推定したら、「この回帰分析は成功なのか、失敗なのか」と...
以下の表の「判断方法」に記された検定などの結果、全ての仮...
ここで工夫とは、「古典的仮定が満たされていなかった場合、...
その上で、どう工夫しても改善できなかった場合、表の一番右...
|仮定 |内容 |判断方法|対処法|
|仮定1|説明変数は非確率的である(or,説明変数と撹乱項は直交...
|仮定2|撹乱項の期待値は0である。 |残差の期待値はいつも0...
|仮定3|撹乱項の分散は一定である。 |Breush Pagan test, Wh...
|仮定4|撹乱項は自己相関しない。 |Durbin-Watson, Breusc...
|仮定5|撹乱項は正規分布に従う。 |残差に対する各種の正...
表では、古典的仮定については、[[浅野・中村(2000):http://w...
なお、仮定１～仮定４のすべてが崩れた場合は、GMMという最終...
それでは、実際にこれらのチェックをしてみましょう。
とりあえず、&mimetex(I_t=a+bY_t+Cr_{t-4});というモデリン...
**説明変数の外生性のチェック [#fd0a5ff1]
hausman.systemfit {systemfit}
を使います．
COLOR(Fuchsia){（保留）：（説明変数と撹乱項の直交条件）。...
**分散均一性のチェック [#t8c26c13]
***図を描いて目視 [#x372ae75]
RIGHT:参考：浅野・中村(2000)[pp.132]
横軸に説明変数、縦軸に残差をプロットしてみて、目でまずは...
で、残差をプロットする方法を説明します。その前に、残差系...
> ?lm
などと入力してみましょう。すると、lmについてのヘルプがぽ...
Value:
という項目があります。これは、「lm()というコマンドを打っ...
residuals
という変数がありまます。たとえば、回帰分析したときの残差...
> lm(I ~ Y + r1)$residuals
のように、$(ドルマーク）で繋いで、入力します。すると、残...
Rは賢く、
> lm(I ~ Y + r1)$res
でも同じ結果を返します。Value一覧のなかで、アルファベット...
> lm(I ~ Y + r1)$r
とすれば、
NULL
となります。これは、Valueのうち、residualsではなく、rank...
> names(lm(I ~ Y + r1))
と入力しても、Value一覧を見ることができます。これを見ると...
> names(summary(lm(I ~ Y + r1)))
と打ってみましょう。はい、ちゃんとありましたね。
さらに余談を続けます。「r1のt値だけとりだす」にはどうした...
とりえあず、
> summary(lm(I ~ Y + r1))$coefficients
と入力してみましょう。すると、
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.821201 2.5193943 3.104397 2.468089e-03
Y 2.090334 0.2443774 8.553713 1.268153e-13
r1 -1.534889 0.3280416 -4.678946 8.883377e-06
と出力されます。これは、(3*4)の行列とRは認識しているよう...
> summary(lm(I ~ Y + r1))$coefficients[3,3]
とすればよいことがわかります。
ちなみに、「F検定のp値だけとりだす」にはどうしたらいいで...
> summary(lm(I ~ Y + r1))$fstatistic
と入力してみましょう。
value numdf dendf
43.17844 2.00000 102.00000
と出力されますが、ここにはなぜかp値がない！
苦肉の策として、
> BAG<-summary(lm(I ~ Y + r1))$fstatistic
> pf(BAG[1], BAG[2], BAG[3], lower.tail=FALSE)
value
2.597481e-14
が考えられますが、どうもいけてないですね、これは。RjpWiki...
さて、話をもとに戻します。横軸に説明変数、縦軸に残差をプ...
> plot(r1,lm(I ~ Y + r1)$res)
あるいは、
> plot(Y,lm(I ~ Y + r1)$res)
とすればよいです。この図の結果を見てみると、まぁ微妙。分...
***Breush Pagan test [#i533be5c]
RIGHT:参考：浅野・中村(2000)[pp.135]
コマンドは、
bptest()
(パッケージ：lmtest)
です。
具体的に、&mimetex(I_t=a+bY_t+Cr_{t-4});というモデルにお...
> bptest(I ~ Y + r1)
あるいは、
> bptest(lm(I ~ Y + r1))
> bptest(summary((lm(I ~ Y + r1))))
でも同じ結果が返されるようです。
COLOR(RED){（注意）：パッケージの読み込みを忘れずに。}
これを実行すると、
studentized Breusch-Pagan test
data: summary((lm(I ~ Y + r1)))
BP = 1.8808, df = 2, p-value = 0.3905
という結果が返ってきます。よって、"H0:分散は均一である"は...
COLOR(Fuchsia){（保留）：studentizedってなんですか？疑問...
***White Test [#o13a8fec]
RIGHT:参考：浅野・中村(2000)[pp.135]
コマンドは、
white.test()
(パッケージ：tseries)
です。
具体的に、&mimetex(I_t=a+bY_t+cr_{t-4});というモデルにお...
***Goldfeld-Quandt test [written by [[Akihiko Noda:http:/...
RIGHT:参考：浅野・中村(2000)[pp.134]
コマンドは、
gqtest()
(パッケージ：lmtest)
です。
基本的な考え方としては、存在するデータを2つに分割して、そ...
具体的に、&mimetex(I_t=a+bY_t+Cr_{t-4});というモデルにお...
> gqtest(I ~ Y + r1)
あるいは、
> gqtest(lm(I ~ Y + r1))
> gqtest(summary((lm(I ~ Y + r1))))
でも同じ結果が返されるようです。
COLOR(RED){（注意）：パッケージの読み込みを忘れずに。}
これを実行すると、
Goldfeld-Quandt test
data: I ~ Y + r1
GQ = 0.6811, df1 = 50, df2 = 49, p-value = 0.9102
という結果が返ってきます。よって、"H0:2つの分散は等しい"...
**系列相関のチェック [#q9258075]
***図を描いて目視 [#i2020d58]
RIGHT:参考：浅野・中村(2000)[pp.146-147]
図の描き方として、二つあります。一つ目としては、横軸に残...
このプロットの仕方は、[[分散均一性のチェック:http://www.s...
> plot(lm(I ~ Y + r1)$res[1:(length(lm(I ~ Y + r1)$res)-...
と入力してみてください。ちょっとコマンドが複雑で見通しが...
> lm(I ~ Y + r1)$res->zansa
> plot(zansa[1:(length(zansa)-1)],zansa[2:length(zansa)])
と入力しても同じです。
プロットされた図は、円状ではありません。何か、右上がりの...
二つ目の図の描き方としては、横軸の時間(t),縦軸の残差をと...
> plot(c(1:length(zansa)),zansa)
とかやれば目的は果たせます。
> plot(zansa)
でも出来るみたいです。
出てきた図を見ると、ブリッジ上に反り返っているので、やは...
横軸に時間(t)をとる方法も説明しておきます。
> plot(ts(zansa,start=c(1971,1),frequency=4),type="p")
もしくは
> plot(ts(zansa,start=c(1971,1),frequency=4),type="o")
などと入力すればよいです。いま、もともとのデータの観測期...
そこで、これを統計的にちゃんと判断する方法を以下で説明し...
***Durbin-Watson [#e59ae9ad]
RIGHT:参考：浅野・中村(2000)[pp.148]
コマンドは、
dwtest()
(パッケージ：lmtest)
もしくは、
durbin.watson
(パッケージ：car)
です。
具体的に、&mimetex(I_t=a+bY_t+cr_{t-4});というモデルにお...
> dwtest(I ~ Y + r1)
あるいは、
> dwtest(lm(I ~ Y + r1))
> dwtest(summary((lm(I ~ Y + r1))))
でも同じ結果が返されるようです。
COLOR(RED){（注意）パッケージの読み込みを忘れずに。}
これを実行すると、
Durbin-Watson test
data: summary((lm(I ~ Y + r1)))
DW = 0.2924, p-value < 2.2e-16
alternative hypothesis: true autocorrelation is greater ...
という結果が返されます。したがって、"H0:系列相関はない"が...
しかし、DWテストでは正規分布のAR(1)しか想定していません。...
***Breusch-Godfrey test [#c12150d4]
RIGHT:参考：浅野・中村(2000)[pp.152]
コマンドは、
bgtest()
(パッケージ：lmtest)
です。
具体的に、&mimetex(I_t=a+bY_t+cr_{t-4});というモデルにお...
> bgtest(I ~ Y + r1)
あるいは、
> bgtest(lm(I ~ Y + r1))
> bgtest(summary((lm(I ~ Y + r1))))
でも同じ結果が返されるようです。
COLOR(RED){（注意）パッケージの読み込みを忘れずに。}
これを実行すると、
Breusch-Godfrey test for serial correlation of o...
data: lm(I ~ Y + r1)
LM test = 77.1659, df = 1, p-value < 2.2e-16
という結果が返されます。したがって、"H0:系列相関はない"が...
***Q Test(Box-Pierce, Ljung-Box test) [#e6d7a41c]
コマンドは
> Box.test()
です。[[管理人(SUGIYAMA Shunsuke):http://sugi-shun.com]]...
**RESETテスト [#j1869d35]
resettest {lmtest}
を使います．
COLOR(Fuchsia){（保留）：そのうち加筆予定。関数形が正しい...
**ここまでの回帰診断のまとめ [#f1a7ac05]
どうやら、分散不均一性は発生していないようです。しかし、...
しかし系列相関が発生してもOLS推定量は最小分散性を失うだけ...
COLOR(Fuchsia){（保留）：Hausman test, RESET についてその...
*どう工夫しても仮定が崩れるときの対処法 [#d23942c0]
[[回帰診断:http://www.sugi-shun.com/econwiki/index.php?R%...
仮定1,仮定2は、OLS推定量が一致性を持つために必要です。仮...
仮定3,仮定4は、OLS推定量が最小分散性（有効性＝効率性）を...
仮定5はなくても漸近的に何も問題ありません。直感的には、ど...
考慮すべきは、仮定1,仮定3,仮定4の三つということになります...
というわけで、古典的仮定が崩れたときの対処法をコンパクト...
|崩れる仮定|対処法|
|仮定1|2SLS(IV)|
|仮定3|HCSE(White), GLS(WLSなど)|
|仮定4|HCSE(Newey-West), GLS(Cochrane-Orcutt法など)|
|仮定3and仮定4|GLS|
|仮定1and仮定3and仮定4|GMM|
これらの対処法をRでどうやってやるのかを以下で説明したいと...
**2SLS(IV):2Step Least Squares(Instrument Variables)[#qa5...
RIGHT:参考：Wooldridge(2002)[Ch5]
日本語では二段階最小二乗法(操作変数法)のことです（IV≠TSLS...
例として、ここではWooldridge(2002)[Problem 5.4]をとりあげ...
#ref(card.csv)
です。このデータについての説明は
#ref(card.txt)
を見てください。
このデータは、[[Wooldridgeのウェブサイト:http://www.msu.e...
とりあえず、
> read.table("card.csv",header=T,sep=",")->c
> attach(c)
などとしてデータを読み込みましょう。
まず、考えたいモデルは以下であるとします。
lwage = A + B*educ + C*exper + D*expersq + E*black + F*so...
とりあえず、OLS推定してみましょう。
> summary(lm(lwage ~ educ + exper + expersq + black + so...
+ + smsa + reg661 + reg662+reg663+reg664+reg665+reg666+r...
その結果として、
Call:
lm(formula = lwage ~ educ + exper + expersq + black + so...
smsa + reg661 + reg662 + reg663 + reg664 + reg665 + ...
reg667 + reg668 + smsa66)
Residuals:
Min 1Q Median 3Q Max
-1.62326 -0.22141 0.02001 0.23932 1.33340
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.7393766 0.0715282 66.259 < 2e-16 ***
educ 0.0746933 0.0034983 21.351 < 2e-16 ***
exper 0.0848320 0.0066242 12.806 < 2e-16 ***
expersq -0.0022870 0.0003166 -7.223 6.41e-13 ***
black -0.1990123 0.0182483 -10.906 < 2e-16 ***
south -0.1479550 0.0259799 -5.695 1.35e-08 ***
smsa 0.1363846 0.0201005 6.785 1.39e-11 ***
reg661 -0.1185697 0.0388301 -3.054 0.002281 **
reg662 -0.0222026 0.0282575 -0.786 0.432092
reg663 0.0259703 0.0273644 0.949 0.342670
reg664 -0.0634942 0.0356803 -1.780 0.075254 .
reg665 0.0094551 0.0361174 0.262 0.793504
reg666 0.0219476 0.0400984 0.547 0.584182
reg667 -0.0005887 0.0393793 -0.015 0.988073
reg668 -0.1750058 0.0463394 -3.777 0.000162 ***
smsa66 0.0262417 0.0194477 1.349 0.177327
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1...
Residual standard error: 0.3723 on 2994 degrees of freedom
Multiple R-Squared: 0.2998, Adjusted R-squared: 0.29...
F-statistic: 85.48 on 15 and 2994 DF, p-value: < 2.2e-16
を得ます。
しかしこのような賃金関数の推定において、educ(教育年数)は...
このときOLS推定量は一致性を失いますから、上記の推定結果は...
そこで操作変数の登場です。ここでは、nearc4（4年制大学から...
+説明変数educとnearc4は相関するか？
+εとnearc4は無相関であるか？
の二つです（これは、厳密な言い方ではありません）。
***操作変数(nearc4)は，説明変数(educ)と相関するか？ [#g15...
まず一つ目のチェックをします。これは、より厳密には、
educ = a + b*nearc4 + c*exper + d*expersq + e*black + f*s...
というモデルにおいて、b≠0である、というのが正確です。
これをチェックしてみましょう。
> summary(lm(educ ~ nearc4 + exper + expersq + black + s...
+ + reg661 + reg662+reg663+reg664+reg665+reg666+reg667+r...
と入力し、
Call:
lm(formula = educ ~ nearc4 + exper + expersq + black + s...
smsa + reg661 + reg662 + reg663 + reg664 + reg665 + ...
reg667 + reg668 + smsa66)
Residuals:
Min 1Q Median 3Q Max
-7.54513 -1.36996 -0.09103 1.27836 6.23847
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 16.8485239 0.2111222 79.805 < 2e-16 ***
nearc4 0.3198989 0.0878638 3.641 0.000276 ***
exper -0.4125334 0.0336996 -12.241 < 2e-16 ***
expersq 0.0008686 0.0016504 0.526 0.598728
black -0.9355287 0.0937348 -9.981 < 2e-16 ***
south -0.0516126 0.1354284 -0.381 0.703152
smsa 0.4021825 0.1048112 3.837 0.000127 ***
reg661 -0.2102710 0.2024568 -1.039 0.299076
reg662 -0.2889073 0.1473395 -1.961 0.049992 *
reg663 -0.2382099 0.1426357 -1.670 0.095012 .
reg664 -0.0930890 0.1859827 -0.501 0.616742
reg665 -0.4828875 0.1881872 -2.566 0.010336 *
reg666 -0.5130857 0.2096352 -2.448 0.014442 *
reg667 -0.4270887 0.2056208 -2.077 0.037880 *
reg668 0.3136204 0.2416739 1.298 0.194490
smsa66 0.0254805 0.1057692 0.241 0.809644
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1...
Residual standard error: 1.941 on 2994 degrees of freedom
Multiple R-Squared: 0.4771, Adjusted R-squared: 0.47...
F-statistic: 182.1 on 15 and 2994 DF, p-value: < 2.2e-16
を得ます。nearc4のt値は3.641です。ちゃんと有意です。これ...
***操作変数(nearc4)は，撹乱項(ε)と関係しないか？ [#vd6434...
次にnearc4が操作変数としての条件2をクリアするか検討しまし...
educ(教育年数)⇔能力(観測されない)⇔撹乱項
なので、操作変数nearc4を使います。nearc4が撹乱項と直交条...
nearc4(大学から近いところに住んでいるか？)⇔能力(観測され...
において、nearc4と能力（の代理変数であるIQ)が相関しないか...
(能力とnearc4が相関していなくても、撹乱項とnearc4に相関が...
すなわち、IQがnearc4と関係しないということを確認しないと...
・・・確認しようと思ったら、iqに欠損値があるようで、これ...
(read.tableにはオプションにna.stringsがあって、読み込みの...
とりあえず，欠損値をすべて「空白」に置き換えたcsvファイル...
#ref(card2.csv)
空白にしておけば，Rに読み込むとこれをNA，すなわち欠損値と...
> read.table("card2.csv",header=T,sep=",")->c
> attach(c)
で改めてデータを読み込みましょう．
> summary(lm(nearc4~ iq))
とやります．結果は，以下です．
Call:
lm(formula = nearc4 ~ iq)
Residuals:
Min 1Q Median 3Q Max
-0.8021 -0.6690 0.2701 0.3016 0.4031
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.4773200 0.0671001 7.114 1.55e-12 ***
iq 0.0022555 0.0006477 3.483 0.000507 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1...
Residual standard error: 0.4534 on 2059 degrees of freedom
(949 observations deleted due to missingness)
Multiple R-Squared: 0.005856, Adjusted R-squared: 0.00...
F-statistic: 12.13 on 1 and 2059 DF, p-value: 0.0005071
nearc4 に対してiqで単回帰すると、iqのt値は有意になりまし...
> summary(lm(nearc4 ~ iq+smsa66+reg661+reg662+reg669))
とやれば，結果は，以下です．
Call:
lm(formula = nearc4 ~ iq + smsa66 + reg661 + reg662 + re...
Residuals:
Min 1Q Median 3Q Max
-0.9659 -0.3916 0.1208 0.2293 0.6265
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.3346827 0.0613386 5.456 5.45e-08 ***
iq 0.0006253 0.0005919 1.056 0.290884
smsa66 0.3747183 0.0199138 18.817 < 2e-16 ***
reg661 0.1433780 0.0414904 3.456 0.000560 ***
reg662 0.1745655 0.0241401 7.231 6.72e-13 ***
reg669 0.1029124 0.0308465 3.336 0.000864 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1...
Residual standard error: 0.4082 on 2055 degrees of freedom
(949 observations deleted due to missingness)
Multiple R-Squared: 0.1959, Adjusted R-squared: 0.19...
F-statistic: 100.1 on 5 and 2055 DF, p-value: < 2.2e-16
もともと，単回帰分析を考えていたわけではにので，このよう...
条件2もクリアした，ということが確認されました．
***2SLSする [#raaf0f44]
そこでいよいよ2SLSを行います。
使うコマンドは、
>tsls()
パッケージ:sem
です。尚，
EquationsModelling {fMultivar}
や
nlsystemfit {systemfit}
や
systemfit {systemfit}
などでは，サブコマンドでOLS,2SLS,WLS,SUR,3SLSなどが網羅的...
さて，
>tsls()
パッケージ:sem
の場合，基本的に、X1のための操作変数をZとした場合、
summary(tsls(Y ~ X1 + X2 + X3, ~ Z + X2 +X3 ))
とやります。操作変数が複数の場合は、例えばX1のための操作...
summary(tsls(Y ~ X1 + X2 + X3, ~ Z1 + Z2 + X2 +X3 ))
とやればよいです。
ここで．
#ref(card.csv)
のデータにまた戻ります．Rを再起動して，
> read.table("card2.csv",header=T,sep=",")->c
> attach(c)
で再びデータを読み込んでおきましょう．
以下のように入力しましょう。
> summary(tsls(lwage ~ educ + exper + expersq + black + ...
+ + smsa + reg661 + reg662+reg663+reg664+reg665+reg666+r...
+ +reg668+smsa66,
+ ~ nearc4 + exper + expersq + black + south + smsa + re...
+ + reg662+reg663+reg664+reg665+reg666+reg667+reg668+sms...
結果は以下のようになります。
2SLS Estimates
Model Formula: lwage ~ educ + exper + expersq + black + ...
reg662 + reg663 + reg664 + reg665 + reg666 + reg667 ...
smsa66
Instruments: ~nearc4 + exper + expersq + black + south +...
reg663 + reg664 + reg665 + reg666 + reg667 + reg668 ...
Residuals:
Min. 1st Qu. Median Mean 3rd Qu. M...
-1.83e+00 -2.41e-01 2.43e-02 -5.84e-12 2.52e-01 1.43e...
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.773966 0.9349469 4.0366 5.560e-05
educ 0.131504 0.0549637 2.3926 1.679e-02
exper 0.108271 0.0236586 4.5764 4.923e-06
expersq -0.002335 0.0003335 -7.0014 3.116e-12
black -0.146776 0.0538999 -2.7231 6.504e-03
south -0.144671 0.0272846 -5.3023 1.227e-07
smsa 0.111808 0.0316620 3.5313 4.197e-04
reg661 -0.107814 0.0418137 -2.5784 9.972e-03
reg662 -0.007046 0.0329073 -0.2141 8.305e-01
reg663 0.040445 0.0317806 1.2726 2.033e-01
reg664 -0.057917 0.0376059 -1.5401 1.236e-01
reg665 0.038458 0.0469387 0.8193 4.127e-01
reg666 0.055089 0.0526597 1.0461 2.956e-01
reg667 0.026758 0.0488287 0.5480 5.837e-01
reg668 -0.190891 0.0507113 -3.7643 1.702e-04
smsa66 0.018531 0.0216086 0.8576 3.912e-01
Residual standard error: 0.3883 on 2994 degrees of freedom
ここで、2SLS推定のほうが、OLS推定よりも標準誤差が大きいこ...
ここまでは操作変数がただ一つでした。このとき、IV推定量と2...
> summary(tsls(lwage ~ educ + exper + expersq + black + ...
+ + smsa + reg661 + reg662+reg663+reg664+reg665+reg666+...
+ +reg668+smsa66,
+ ~ nearc2+nearc4 + exper + expersq + black + south + s...
+ + reg662+reg663+reg664+reg665+reg666+reg667+reg668+sm...
とするだけです。簡単ですね。
2SLS Estimates
Model Formula: lwage ~ educ + exper + expersq + black + ...
reg662 + reg663 + reg664 + reg665 + reg666 + reg667 ...
smsa66
Instruments: ~nearc2 + nearc4 + exper + expersq + black ...
reg662 + reg663 + reg664 + reg665 + reg666 + reg667 ...
smsa66
Residuals:
Min. 1st Qu. Median Mean 3rd Qu. M...
-1.94e+00 -2.51e-01 1.93e-02 -3.00e-12 2.65e-01 1.47e...
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.3396875 0.8945377 3.733423 1.924e-04
educ 0.1570593 0.0525782 2.987155 2.839e-03
exper 0.1188149 0.0228061 5.209792 2.019e-07
expersq -0.0023565 0.0003475 -6.780907 1.433e-11
black -0.1232778 0.0521500 -2.363907 1.815e-02
south -0.1431945 0.0284448 -5.034120 5.084e-07
smsa 0.1007530 0.0315193 3.196547 1.405e-03
reg661 -0.1029760 0.0434224 -2.371497 1.778e-02
reg662 -0.0002287 0.0337943 -0.006766 9.946e-01
reg663 0.0469556 0.0326490 1.438194 1.505e-01
reg664 -0.0554084 0.0391828 -1.414098 1.574e-01
reg665 0.0515041 0.0475678 1.082752 2.790e-01
reg666 0.0699968 0.0533049 1.313139 1.892e-01
reg667 0.0390596 0.0497499 0.785119 4.324e-01
reg668 -0.1980371 0.0525350 -3.769623 1.667e-04
smsa66 0.0150626 0.0223360 0.674364 5.001e-01
Residual standard error: 0.4053 on 2994 degrees of freedom
以上で2SLSによる推定がRで出来るようになりました。
尚、上記の推定結果は、Woolridge(2002)[Problem 5.4]の出題...
***データに関する注意 [#iac47a73]
COLOR(RED){（注意）：データに関する注意．}
expersqとlwageについて、Wooldridgeが用意したデータを使う...
> exper2<-exper^2
> logwage<-log(wage)
とかやってRで加工した系列を使うと、微妙に異なる推定結果を...
具体的には，
#ref(card_3.csv)
のように，expersq, lwageが入っていない状態のデータのcsvを...
read.table("card_3.csv",header=T,sep=",")->c
attach(c)
library(sem)
exper2<-exper^2
logwage<-log(wage)
summary(tsls(logwage ~ educ + exper + exper2 + black + s...
+ smsa + reg661 + reg662+reg663+reg664+reg665+reg666+reg...
+reg668+smsa66,
~ nearc4 + exper + exper2 + black + south + smsa + reg661
+ reg662+reg663+reg664+reg665+reg666+reg667+reg668+smsa6...
とかやれば，
2SLS Estimates
Model Formula: logwage ~ educ + exper + exper2 + black +...
reg662 + reg663 + reg664 + reg665 + reg666 + reg667 ...
smsa66
Instruments: ~nearc4 + exper + exper2 + black + south + ...
reg663 + reg664 + reg665 + reg666 + reg667 + reg668 ...
Residuals:
Min. 1st Qu. Median Mean 3rd Qu. M...
-1.84e+00 -2.45e-01 2.58e-02 2.17e-11 2.54e-01 1.43e...
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.742871 0.9520939 3.9312 8.644e-05
educ 0.134309 0.0563767 2.3824 1.726e-02
exper 0.106333 0.0230529 4.6126 4.144e-06
exper2 -0.002209 0.0003321 -6.6495 3.483e-11
black -0.138156 0.0572038 -2.4152 1.579e-02
south -0.143702 0.0275365 -5.2186 1.926e-07
smsa 0.108207 0.0329625 3.2827 1.040e-03
reg661 -0.106913 0.0422370 -2.5313 1.142e-02
reg662 -0.006620 0.0332152 -0.1993 8.420e-01
reg663 0.042308 0.0322907 1.3102 1.902e-01
reg664 -0.058365 0.0378677 -1.5413 1.234e-01
reg665 0.040148 0.0476794 0.8420 3.998e-01
reg666 0.052881 0.0526310 1.0048 3.151e-01
reg667 0.026308 0.0491505 0.5352 5.925e-01
reg668 -0.191616 0.0511800 -3.7440 1.845e-04
smsa66 0.018182 0.0217943 0.8343 4.042e-01
Residual standard error: 0.3913 on 2994 degrees of freedom
となり，微妙に異なる値が得られている．
WooldridgeのウェブサイトでDLできるexpersqとlwageのデータ...
・・・このように，計量分析では，論文の追試をやると，微妙...
***補足 [#nf96aed0]
systemfit
パッケージ：systemfit
でも2SLS推定が出来るようです。2SLSだけでなく、3SLSとかも...
?systemfit
で見てください。
***過剰識別制約テスト [#u19330c8]
COLOR(Fuchsia){（保留）：操作変数の過剰識別制約テストの使...
**HCSE(Heteroskedasticity Consistene Standard Error) [#x0...
RIGHT:参考："[[Econometrics in R:http://cran.r-project.or...
***White[#i69bbbc1]
分散不均一が発生してもOLS推定量は不偏性（一致性）は満たし...
コマンドは、
vcovHC
パッケージ:sandwich
もしくは
hccm
パッケージ：car
です。どちらを使っても同じ結果が得られることをちゃんと確...
使うデータは、[[最小二乗法による推定:http://www.sugi-shun...
[[lm()が自動出力するt検定量とF検定量を自分で計算してみよ...
COLOR(RED){（注意）パッケージ：car, sandwichを読み込みま...
> unrestricted<-lm(I ~ Y + r1)
> rhs2<-c(0)
> hm2 <- rbind(c(0,1,0))
> lht(unrestricted,hm2,rhs2,test=c("Chisq"))
と入力すれば、
Linear hypothesis test
Hypothesis:
Y = 0
Model 1: I ~ Y + r1
Model 2: restricted model
Res.Df RSS Df Sum of Sq Chisq Pr(>Chisq)
1 102 3254.5
2 103 5589.0 -1 -2334.5 73.166 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1...
と出力されます。欲しい情報はカイ二乗統計量なので、
> lht(unrestricted,hm2,rhs2,test=c("Chisq"))$Chi[2]
です。t値はこれの平方根なのでsqrtをとります。
> sqrt(lht(unrestricted,hm2,rhs2,test=c("Chisq"))$Chi[2])
[1] 8.553713
そしてこれはYのt値の他ならないことをすでに見ました。
このt値は分散均一の仮定のもとで計算されたものです。この仮...
まず、修正した共分散行列を使って仮説検定する方法は
> lht(unrestricted,hm2,rhs2,vcov=vcovHC(unrestricted) )
あるいは
> lht(unrestricted,hm2,rhs2,vcov=hccm(unrestricted) )
です。欲しいのはカイ二乗統計量の平方根なので、
> lht(unrestricted,hm2,rhs2,vcov=hccm(unrestricted),test...
> sqrt(whitese)
もしくは
> lht(unrestricted,hm2,rhs2,vcov=vcovHC(unrestricted),te...
> sqrt(whitese2)
と入力します。hccm()でもvcovHC()でも
[1] 6.390953
が得られました。これがWhiteの一致標準誤差を使った場合のt...
尚、White修正にはいろいろなバージョンがあります。
詳しくは、
> ?hccm
> ?vcovHC
で見てください。ここの"type"というサブコマンドをどのよう...
hccm()を使った場合は以下です。
> sqrt(lht(unrestricted,hm2,rhs2,vcov=hccm(unrestricted,...
[1] 6.741862
> sqrt(lht(unrestricted,hm2,rhs2,vcov=hccm(unrestricted,...
[1] 6.644852
> sqrt(lht(unrestricted,hm2,rhs2,vcov=hccm(unrestricted,...
[1] 6.564739
> sqrt(lht(unrestricted,hm2,rhs2,vcov=hccm(unrestricted,...
[1] 6.390953
> sqrt(lht(unrestricted,hm2,rhs2,vcov=hccm(unrestricted,...
[1] 6.374456
vcovHC()を使った場合は以下です。
> sqrt(lht(unrestricted,hm2,rhs2,vcov=vcovHC(unrestricte...
[1] 6.741862
> sqrt(lht(unrestricted,hm2,rhs2,vcov=vcovHC(unrestricte...
[1] 6.741862
> sqrt(lht(unrestricted,hm2,rhs2,vcov=vcovHC(unrestricte...
[1] 6.644852
> sqrt(lht(unrestricted,hm2,rhs2,vcov=vcovHC(unrestricte...
[1] 6.564739
> sqrt(lht(unrestricted,hm2,rhs2,vcov=vcovHC(unrestricte...
[1] 6.390953
> sqrt(lht(unrestricted,hm2,rhs2,vcov=vcovHC(unrestricte...
[1] 6.374456
ここで、"hc0"(hccm()を使った場合),"HC"もしくは"HC0"（vcov...
尚、vcovを使った場合、typeとして"const"を指定できます。こ...
> sqrt(lht(unrestricted,hm2,rhs2,vcov=vcovHC(unrestricte...
[1] 8.553713
ここの項目は、"[[Econometrics in R:http://cran.r-project....
を参考にしました。
***Newey-West [#s56bb41a]
Newey-Westの方法は、White修正の「系列相関があっても大丈夫...
コマンドは、
NeweyWest
パッケージ：sandwich
です。
**GLS(Generalized Least Squares) [#sc894fcf]
RIGHT:参考：浅野・中村(2000)[pp.130]
日本語では一般化最小二乗法です。[[管理人(SUGIYAMA Shunsuk...
そのうち、自分でプログラムを書くなり、組み込み関数を探す...
以下は覚書です。関連しそうなのは、
lm.gls
パッケージ:MASS
gls
パッケージ:nlme
glsD
パッケージ:Design
あたりか？
***WLS:Weighted Least Squares [#m1633509]
RIGHT:参考：浅野・中村(2000)[pp.127-130]
日本語では、加重最小二乗法です。
そのうち、自分でプログラムを書くなり、組み込み関数を探す...
***Cochrane Orcutt法 [#k6400586]
RIGHT:参考：浅野・中村(2000)[pp.153]
そのうち、自分でプログラムを書くなり、組み込み関数を探す...
**GMM [#u08f350f]
[[GMM:http://www.sugi-shun.com/econwiki/index.php?R%A4%C7...
*時系列分析 [#mf5cad5a]
[[Rで時系列分析]]を見てください。
*計量ミクロ [#ced0ba2a]
[[Rで計量ミクロ]]を見てください。
*OLS以外の推定方法 [#x0baa89f]
[[RでOLS以外の推定方法]]を見てください。
*Reference [#h7afec48]
**書籍 [#h8a36313]
-[[浅野・中村(2000):http://www.amazon.co.jp/gp/product/46...
-[[Greene(2003):http://www.amazon.co.jp/gp/product/013110...
-[[Hamilton(1994):http://www.amazon.co.jp/gp/product/0691...
-[[Wooldridge(2002):http://www.amazon.co.jp/Econometric-A...
**リンク [#f4df0207]
-[[Econometrics in R:http://cran.r-project.org/doc/contri...
-[[R for economists:http://www.mayin.org/ajayshah/KB/R/R_...
-[[A Brief Guide to R for Beginners in Econometrics:http:...
-[[CRAN Task View: Computational Econometrics:http://cran...
-[[パッケージ一覧:http://cran.md.tsukuba.ac.jp/src/contri...
ページ名:
![[PukiWiki]](image/econwiki.png "[PukiWiki]")
