![[PukiWiki]](image/econwiki.png "[PukiWiki]")
|このページについて| ''概要'':このページは、統計言語Rで計量経済分析する方法についての、[[管理人(SUGIYAMA Shunsuke):http://sugi-shun.com]]によるオリジナルな解説ページです。このページは、2006-08-10につくりました。暇をみつけて更新していきます。いつの日か、「Rで計量経済分析をするならここを見ろ」となるくらいのレベルを目指します。ページサイズがでかくなったら、適宜、テーマごとにページを分割していきます。 ''概要'':このページは、統計言語Rで計量経済分析する方法についての解説ページです。暇をみつけて更新していきます。 関連するサイトは、[[リンク:http://www.sugi-shun.com/econwiki/index.php?R%A4%F2%BB%C8%A4%C3%A4%C6%B7%D7%CE%CC%B7%D0%BA%D1%CA%AC%C0%CF#f4df0207]]を見てください。 ''親ページ'':このページの親ページは[[R]]です。 |目次| #contents ---- *古典的線形回帰モデル [#yed9bf5c] **最小二乗法による推定 [#u30dda11] とりあえず、RでOLSする方法を学びます。 例として、こんな問題を考えて見ましょう。 |投資関数は利子率の減少関数って習ったけど、これって本当？| いや、あくまで例ですよ。 以下のデータを使います。 #ref(inv.csv) このデータは、日本の企業設備投資と国内総支出と長期貸し出しレートのデータです。頻度は四半期データで、投資と支出は実質表示、観測期間は1970年第1四半期-1997:第1四半期です。 見出しの意味は、 |企業設備投資|国内総支出|Long-term Prime Lending Rate| |INV | GDE |R| です。 とりあえず、 > read.table("http://sugi-shun.com/econwiki/index.php?plugin=attach&refer=R%A4%F2%BB%C8%A4%C3%A4%C6%B7%D7%CE%CC%B7%D0%BA%D1%CA%AC%C0%CF&openfile=inv.csv",header=T,sep=",")->i と入力して、Rにこのデータを読み込ませましょう。COLOR(RED){（注意）：ディレクトリの変更を忘れずに。} 次に、 > attach(i) と入力しましょう。こうすると、見出しの変数名を読み込んでくれます。たとえば、 > INV などと打つと、企業設備投資のデータが出力されることを確認しましょう。 ***単回帰 [#b8aa6a65] とりあえず、単回帰してみましょう。 投資(Inv)と国内支出(GDE)の対前年同期比伸び率をまず求め、それぞれI, Yと定義しましょう。 |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では、大文字と小文字を区別します。だから、iとIは別物と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);と推定されたことを意味します。計量経済分析では、係数の推定値だけではなく、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 ' ' 1 Residual standard error: 6.195 on 103 degrees of freedom Multiple R-Squared: 0.3422, Adjusted R-squared: 0.3359 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 ' ' 1 Residual standard error: 7.604 on 103 degrees of freedom Multiple R-Squared: 0.009055, Adjusted R-squared: -0.0005663 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 ' ' 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)[1,] fit lwr upr 5.646562 4.104640 7.188485 として出てきます。 単回帰の場合の予測をしてみましょう。 predict(lm(I ~ Y ))->yosoku1 plot(Y,yosoku1) とやれば、ちゃんと理論値はOLS推定された直線を表してくれることがわかります。 **ここまでの推定結果の解釈 [#f089e914] +まず、投資関数とかいう名前で習うような&mimetex(I_t=a+br_t);というモデルにおいて、&mimetex(r_t);のt値は有意ではありません。われわれは普通、「投資は利子率の減少関数」と習いますが、実際のデータではこの関係は認めらないということになります。 +もう一つの単回帰式&mimetex(I_t=a+bY_t);については、&mimetex(Y_t);のt値は有意です。符号は正です．GDPが上昇するとき、投資も上昇するよ、ということが観察されました。 +重回帰式&mimetex(I_t=a+bY_t+cr_t);の推定結果を見ると、&mimetex(Y_t);,&mimetex(r_t);の両方についてt値が有意になりました。これは、「利子率は、ほかの事情一定としないで見ると、投資に対して有意な負の影響を持たない。しかし、GDPという要因を一定にしたうえで見ると、投資に対して有意な負の影響を持つ」ということです。 +さらに、&mimetex(r_t);の係数の大きさを見てみると、重回帰のときは-0.7051で、単回帰のときは-0.3836であり、単回帰は重回帰よりも大きくなっている。ここで一般に、利子率が高いときは、景気がいいときであることを思い出しましょう。すると、「r↑のときはY↑という経路からくる影響（Y↑のときはI↑となる影響）のため，重回帰のrの係数より単回帰のrの係数が大きくなっている」と言えそうです。 **より良い計量モデルを考えてみる [#m8485961] [[ここまでの推定結果の解釈:http://www.sugi-shun.com/econwiki/index.php?R%A4%F2%BB%C8%A4%C3%A4%C6%B7%D7%CE%CC%B7%D0%BA%D1%CA%AC%C0%CF#f089e914]]の1.として述べた点について、次のような反論が考えられます。 「企業が今年実行する設備投資計画は、去年、計画したものだ。だから、t期の投資は、t期より以前の利子率に影響されるだろう」 というわけで、&mimetex(I_t=a+br_t);ではなく、&mimetex(I_t=a+br_{t-4});のようなモデリングをするべきだ、という考え方があります。実際にこれを推定してみましょう。 まず、"&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 ' ' 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も回帰式に入れてみたら、「よりよさそう」な推定結果が得られていました。だから、&mimetex(I_t=a+bY_t+cr_{t-4});というモデルを考えて見ましょう。 > 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 ' ' 1 Residual standard error: 5.649 on 102 degrees of freedom Multiple R-Squared: 0.4585, Adjusted R-squared: 0.4479 F-statistic: 43.18 on 2 and 102 DF, p-value: 2.597e-14 ほら、t値はさらに有意になりました。符合条件も、理論通りです。決定係数も、ここまでの試したモデルの中で、一番高いです。((実際に世の中にある論文を読んでみると，決定係数はモデル選択基準としてまったく使われていません．ほとんど無視されています．決定係数が低かろうが高かろうがお構いなし，というのが現実の研究者の相場観です．異なる説明変数を使った式どうしで決定係数を比較することは無意味であることが一因です．))どうやらこれが現実をよく捉えた式なのではないか、という風に考えてよさそうです。 **仮説検定 [#i997fb8b] RIGHT:参考：浅野・中村(2000)[pp.78-80]、"[[Econometrics in R:http://cran.r-project.org/doc/contrib/Farnsworth-EconometricsInR.pdf]]の3.4 Linear Hypothesis Testing (Wald and F)" ここでは、LSE(Linear Squares Estimation)における仮説検定について記述します。Wald test, LM test, 制約付残差に基づくF検定が有名です。また、これら三つは、結果的に同じ検定量となることも有名です。 ***lm()が自動出力するt検定量とF検定量を自分で計算してみよう [#ke1db65d] lm()が自動的にやってくれる仮説検定は、係数が有意に0になるかどうかのt検定とF検定です。これを例えばRの自動出力にまかせずに自分で計算する方法を説明します。 waldtest パッケージ：lmtest というコマンドがあります。Wald検定はこれを使うはずなんですが、[[管理人(SUGIYAMA Shunsuke):http://sugi-shun.com]]はこのコマンドを実践で使ったことがありません。使い方がわかる人は教えてください。 COLOR(Fuchsia){（保留）：waldtest（パッケージ：lmtest）の使い方はそのうち更新します。} 代わりのコマンドを紹介します。 linear.hypothesis() パッケージ：car もしくは、 lht() パッケージ：car を使う方法を記述します。 > summary(lm(I ~ Y + r1)) を例にとりましょう。すでに、"H0:回帰係数はすべて0"が真のもとでのF統計量は F-statistic: 43.18 on 2 and 102 DF, p-value: 2.597e-14 という結果を得ています。これと同じ結果を、自分で計算しても得られるかどうか確認してみましょう。 ここでは、"[[Econometrics in R:http://cran.r-project.org/doc/contrib/Farnsworth-EconometricsInR.pdf]]の3.4 Linear Hypothesis Testing (Wald and F)"に基づいて説明します。 次のように入力します。 > 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''は(P*K)の既知行列、''β''は(K*1)の係数ベクトル、''r''は(P*1)の制約ベクトルです。Pが制約の数、Kは定数項を含む説明変数の数です。 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 ' ' 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の係数が0"ということになります。出力結果は 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 ' ' 1 となります。ここでわれわれが興味があるのは、&mimetex(\chi^2);分布に従う統計量の73.166だけです。これだけを取り出すのには、 > 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値にほかなりません。 ここまでやってきたことで、実は次のようなことがわかりました。すなわち、いつも見ているt検定量やF検定量も、実はWald検定の特殊ケースでしかないということです。次にWald検定について説明します。 ***WALD検定 [#n7d120c5] さて、すでに述べたとおり、Wald test, LM test, 制約付残差に基づくF検定の三つの仮説検定量は、それぞれ違う考え方に基づいて考案された仮説検定量ですが、結果的に同じ計算式になることがわかっています。そして、統計学の知識を使えば、F分布に従う統計量と、&mimetex(\chi^2);分布に従う統計量の二つのタイプの仮説検定量を計算することが出来ます。 RIGHT:参考：Greene(2003)[pp.853-854] 結果のみ示すと、 &mimetex(\chi^2);分布に従う統計量 = 制約の数 × F分布に従う統計量 となります。 先ほどは、何も考えずにF検定量が出ました。これは、RコマンドでF統計量を出すのがデフォルトになっているためです。χ^2統計量を出力するためには、 > 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 ' ' 1 と出力されます。いま制約の数は2ですから、確かに&mimetex(\chi^2);分布に従う統計量(86.357)は、制約の数(2) × F分布に従う統計量(43.178)が成立しています。 *回帰診断 [#j659cc8a] OLS推定したら、「この回帰分析は成功なのか、失敗なのか」ということを考えなければなりません。それは、「古典的仮定は満たされているの？満たされていない場合、どうしたらいいの？」ということを考えることに等しいです。これが回帰診断です。 以下の表の「判断方法」に記された検定などの結果、全ての仮定が満たされていると判断された場合、その回帰分析は成功であると考えられます。満たされないものがあった場合、満たされるように、工夫しましょう。 ここで工夫とは、「古典的仮定が満たされていなかった場合、変数を追加したり減らしたり、あるいはデータ加工をするなどして、古典的仮定がすべで満たされているような推定式を探す」ということです。((これには反対意見もあります。理論経済学者からしたら、「そんなデータ加工を施したモデルで理論を考えたわけではない」と反発されるでしょう。というか、本当は管理人もそう思うのですが、それでは議論が進まないので、立ち止まるよりは前進すべき、という考え方から、このように記しました。)) その上で、どう工夫しても改善できなかった場合、表の一番右に記した「対処法」をとりましょう。 |仮定 |内容 |判断方法|対処法| |仮定1|説明変数は非確率的である(or,説明変数と撹乱項は直交する)|(1)理論から先見的に判断する。(2)Hausman検定（操作変数が天国から降ってきているという大前提のもとであることに注意）。|2SLS(IV)| |仮定2|撹乱項の期待値は0である。 |残差の期待値はいつも0であるから、この仮定はいつも満たされると考えてしまうのが普通（よく考えると、そう考えることに根拠はないです）|| |仮定3|撹乱項の分散は一定である。 |Breush Pagan test, White Test, Goldfeld-Quandt test |HCSE(White),GLS(WLSなど)| |仮定4|撹乱項は自己相関しない。 |Durbin-Watson, Breusch-Godfrey test ,Q test |HCSE(Newey-West),GLS(Cochrane Orcutt法など)| |仮定5|撹乱項は正規分布に従う。 |残差に対する各種の正規性検定|この仮定はなくてもOLSはBLUEになる。しかも漸近理論によれば、この仮定5はなくてもまったく問題ない。| 表では、古典的仮定については、[[浅野・中村(2000):http://www.amazon.co.jp/gp/product/4641160805/sr=8-1/qid=1155206048/ref=sr_1_1/250-7742157-1058652?ie=UTF8&s=gateway]][pp.23]に基づきました。 なお、仮定１～仮定４のすべてが崩れた場合は、GMMという最終手段が待っています。 それでは、実際にこれらのチェックをしてみましょう。 とりあえず、&mimetex(I_t=a+bY_t+Cr_{t-4});というモデリングがいまのところの暫定チャンピオンです。以下では、このモデルをベースに議論します。 **説明変数の外生性のチェック [#fd0a5ff1] hausman.systemfit {systemfit} を使います． COLOR(Fuchsia){（保留）：（説明変数と撹乱項の直交条件）。そのうち加筆予定。Hausman testについて説明します。} **分散均一性のチェック [#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は認識しているようです。だから、これの3行3列にある、r1のt値だけ取り出すには、 > 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で見てみたら、[[このページ:http://www.okada.jp.org/RWiki/index.php?cmd=read&page=%A3%D1%A1%F5%A3%C1%28%B5%EC1%29&word=lm%28%29%20%B7%E8%C4%EA%B7%B8%BF%F4#content_1_46]]の"[質問] lm() 関数の返り値のF統計量のp-value <= なんでも掲示板から移動"に解決法が書かれていました。 さて、話をもとに戻します。横軸に説明変数、縦軸に残差をプロットするには、 > 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});というモデルにおけるBP検定を行う方法は、 > 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});というモデルにおけるwhite.test検定を行う方法は、わかりません。[[管理人(SUGIYAMA Shunsuke):http://sugi-shun.com]]は、Rでwhite.testの実践経験がないので、説明できません。 ***Goldfeld-Quandt test [written by [[Akihiko Noda:http://at-noda.com]]] [#c7c67892] RIGHT:参考：浅野・中村(2000)[pp.134] コマンドは、 gqtest() (パッケージ：lmtest) です。 基本的な考え方としては、存在するデータを2つに分割して、それらを用いて同じモデルを回帰した場合に誤差項の分散が等しくなるかどうかを検定しています。 具体的に、&mimetex(I_t=a+bY_t+Cr_{t-4});というモデルにおけるGQ検定を行う方法は、 > 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] 図の描き方として、二つあります。一つ目としては、横軸に残差、縦軸に1期ずれの残差をプロットしてみましょう。系列相関がなければ、このプロット図は、円状になるはずです。3次元プロットでは、球状になります。 このプロットの仕方は、[[分散均一性のチェック:http://www.sugi-shun.com/econwiki/index.php?R%A4%F2%BB%C8%A4%C3%A4%C6%B7%D7%CE%CC%B7%D0%BA%D1%CA%AC%C0%CF#t8c26c13]]のところで説明した知識を使えば、わかるはずです。 > plot(lm(I ~ Y + r1)$res[1:(length(lm(I ~ Y + r1)$res)-1)],lm(I ~ Y + r1)$res[2:length(lm(I ~ Y + r1)$res)]) と入力してみてください。ちょっとコマンドが複雑で見通しが悪いですね。 > lm(I ~ Y + r1)$res->zansa > plot(zansa[1:(length(zansa)-1)],zansa[2:length(zansa)]) と入力しても同じです。 プロットされた図は、円状ではありません。何か、右上がりの関係がありそうです。すなわち、正の系列相関が強く疑われます。 二つ目の図の描き方としては、横軸の時間(t),縦軸の残差をとってプロットしてみましょう。横軸に時間(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") などと入力すればよいです。いま、もともとのデータの観測期間は1970年第1四半期からでしたが、対前年同期比伸び率で考えていますから、時系列データのスタート時点は、1971年第1四半期であることに注意です「type="p"」とかは、散布図にするためのサブコマンドです。 そこで、これを統計的にちゃんと判断する方法を以下で説明します。 ***Durbin-Watson [#e59ae9ad] RIGHT:参考：浅野・中村(2000)[pp.148] コマンドは、 dwtest() (パッケージ：lmtest) もしくは、 durbin.watson (パッケージ：car) です。 具体的に、&mimetex(I_t=a+bY_t+cr_{t-4});というモデルにおけるDW検定を行う方法は、 > 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 than 0 という結果が返されます。したがって、"H0:系列相関はない"が棄却され、系列相関はあるとわかりました。 しかし、DWテストでは正規分布のAR(1)しか想定していません。高階に渡る自己相関をチェックするために、以下のBreusch-Godfrey testを説明します。 ***Breusch-Godfrey test [#c12150d4] RIGHT:参考：浅野・中村(2000)[pp.152] コマンドは、 bgtest() (パッケージ：lmtest) です。 具体的に、&mimetex(I_t=a+bY_t+cr_{t-4});というモデルにおけるDW検定を行う方法は、 > bgtest(I ~ Y + r1) あるいは、 > bgtest(lm(I ~ Y + r1)) > bgtest(summary((lm(I ~ Y + r1)))) でも同じ結果が返されるようです。 COLOR(RED){（注意）パッケージの読み込みを忘れずに。} これを実行すると、 Breusch-Godfrey test for serial correlation of order 1 data: lm(I ~ Y + r1) LM test = 77.1659, df = 1, p-value < 2.2e-16 という結果が返されます。したがって、"H0:系列相関はない"が棄却され、系列相関はあるとわかりました。目視によって、すでにこの結果は強く強く予想されていましたが、実際p値がほとんど0に近いという結果を得ていますね。 ***Q Test(Box-Pierce, Ljung-Box test) [#e6d7a41c] コマンドは > Box.test() です。[[管理人(SUGIYAMA Shunsuke):http://sugi-shun.com]]は、RでBox.testを使った分散均一性のチェックの実践経験がないので、説明できません。 **RESETテスト [#j1869d35] resettest {lmtest} を使います． COLOR(Fuchsia){（保留）：そのうち加筆予定。関数形が正しいとdefendするために行う検定。} **ここまでの回帰診断のまとめ [#f1a7ac05] どうやら、分散不均一性は発生していないようです。しかし、系列相関が発生してしまっているようです。したがって、モデリングを工夫して、「系列相関を消す」ことを考えるのが一つの戦略です。いくら工夫しても系列相関が消せない場合は、最終の対処法を取ります。この場合は、GLS(Cochrane-Orcutt法など)という推定方法を使うことになります。GLSについては、[[GLS:http://www.sugi-shun.com/econwiki/index.php?R%A4%F2%BB%C8%A4%C3%A4%C6%B7%D7%CE%CC%B7%D0%BA%D1%CA%AC%C0%CF#sc894fcf]]を見てください。 しかし系列相関が発生してもOLS推定量は最小分散性を失うだけです。一致性は失われないので問題なし、とここでは結論付けることにします。 COLOR(Fuchsia){（保留）：Hausman test, RESET についてそのうち言及予定。} *どう工夫しても仮定が崩れるときの対処法 [#d23942c0] [[回帰診断:http://www.sugi-shun.com/econwiki/index.php?R%A4%F2%BB%C8%A4%C3%A4%C6%B7%D7%CE%CC%B7%D0%BA%D1%CA%AC%C0%CF#j659cc8a]]で記した表をおさらいしましょう。ここで、5個の仮定は重要な順に並んでいると考えてよいです。 仮定1,仮定2は、OLS推定量が一致性を持つために必要です。仮定2はいつも満たされると考えてしまうので、実際には仮定1の妥当性を考えればよいです。しかしこれには計量経済学者はとても敏感です。 仮定3,仮定4は、OLS推定量が最小分散性（有効性＝効率性）を持つために必要です。これらの仮定はしばしば崩れますが、一致性を確保できるならば問題ない、という立場もありますから、それほど敏感にならなくてもよいかもしれません。 仮定5はなくても漸近的に何も問題ありません。直感的には、どうせCLT(中心極限定理)によってデータ数が増えれば正規分布に分布収束するからです。 考慮すべきは、仮定1,仮定3,仮定4の三つということになります（もう少しすっきり表現すると、仮定3と仮定4は合体させて、撹乱項の共分散行列が&mimetex(\sigma^2 I);(ここで&mimetex(I);は単位行列)である、と記述できます）。 というわけで、古典的仮定が崩れたときの対処法をコンパクトにまとめると、こうなります。 |崩れる仮定|対処法| |仮定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)[#qa58e42e] RIGHT:参考：Wooldridge(2002)[Ch5] 日本語では二段階最小二乗法(操作変数法)のことです（IV≠TSLS＝2SLS。結果が同じになるだけ(･∀･)ﾆﾔﾆﾔ。参考：[[(；ﾟдﾟ) R言語で同時方程式を推計せよ。:http://uncorrelated.no-ip.com/cgi-bin/view.cgi/20090213/T]]）。 例として、ここではWooldridge(2002)[Problem 5.4]をとりあげます。使うデータは #ref(card.csv) です。このデータについての説明は #ref(card.txt) を見てください。 このデータは、[[Wooldridgeのウェブサイト:http://www.msu.edu/user/ec/faculty/wooldridge/book2.htm]]からとりました。 とりあえず、 > read.table("card.csv",header=T,sep=",")->c > attach(c) などとしてデータを読み込みましょう。 まず、考えたいモデルは以下であるとします。 lwage = A + B*educ + C*exper + D*expersq + E*black + F*south + G*smsa + H*reg661 + I*reg662 + J*reg663 + K*reg664 + L*reg665 + M*reg666 + N*reg667 + O*reg668 + P*smsa66 + ε とりあえず、OLS推定してみましょう。 > summary(lm(lwage ~ educ + exper + expersq + black + south + + smsa + reg661 + reg662+reg663+reg664+reg665+reg666+reg667+reg668+smsa66)) その結果として、 Call: lm(formula = lwage ~ educ + exper + expersq + black + south + smsa + reg661 + reg662 + reg663 + reg664 + reg665 + reg666 + 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 ' ' 1 Residual standard error: 0.3723 on 2994 degrees of freedom Multiple R-Squared: 0.2998, Adjusted R-squared: 0.2963 F-statistic: 85.48 on 15 and 2994 DF, p-value: < 2.2e-16 を得ます。 しかしこのような賃金関数の推定において、educ(教育年数)は撹乱項と相関すると理論的に考えられます。その理由は、「個人の能力」という変数を観測することがとても難しいからです。「個人の能力」が計測できない以上、モデルにも「個人の能力」は現れません。したがって、「個人の能力」は撹乱項としてモデルに現れることになります。ところが、「個人の能力」とeduc（教育年数）は明らかに相関すると考えられます。教育年数が高ければ、能力も高まると考えられるからです。したがって、このモデリングではeduc(教育年数)は撹乱項と相関すると考えられるわけです（educの外生性については、[[説明変数の外生性のチェック:http://www.sugi-shun.com/econwiki/index.php?R%A4%F2%BB%C8%A4%C3%A4%C6%B7%D7%CE%CC%B7%D0%BA%D1%CA%AC%C0%CF]]で実行したHausman検定の結果を思い出しましょう）。 このときOLS推定量は一致性を失いますから、上記の推定結果はぜんぜん信用できない、ということになります。 そこで操作変数の登場です。ここでは、nearc4（4年制大学から近いかどうか）という変数を、操作変数として考えます。Wooldridgeがnearc4を操作変数として考えなさい、といってくれています。まぁ、演習問題ですから。nearc4が操作変数としてふさわしい条件を満たすかどうか、チェックしないといけません。このチェックすべき条件とは、 +説明変数educとnearc4は相関するか？ +εとnearc4は無相関であるか？ の二つです（これは、厳密な言い方ではありません）。 ***操作変数(nearc4)は，説明変数(educ)と相関するか？ [#g1553652] まず一つ目のチェックをします。これは、より厳密には、 educ = a + b*nearc4 + c*exper + d*expersq + e*black + f*south + g*smsa + h*reg661 + i*reg662 + j*reg663 + k*reg664 + l*reg665 + m*reg666 + n*reg667 + o*reg668 + p*smsa66 + υ というモデルにおいて、b≠0である、というのが正確です。 これをチェックしてみましょう。 > summary(lm(educ ~ nearc4 + exper + expersq + black + south + smsa + + reg661 + reg662+reg663+reg664+reg665+reg666+reg667+reg668+smsa66)) と入力し、 Call: lm(formula = educ ~ nearc4 + exper + expersq + black + south + smsa + reg661 + reg662 + reg663 + reg664 + reg665 + reg666 + 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 ' ' 1 Residual standard error: 1.941 on 2994 degrees of freedom Multiple R-Squared: 0.4771, Adjusted R-squared: 0.4745 F-statistic: 182.1 on 15 and 2994 DF, p-value: < 2.2e-16 を得ます。nearc4のt値は3.641です。ちゃんと有意です。これで条件1がクリアです。 ***操作変数(nearc4)は，撹乱項(ε)と関係しないか？ [#vd64342c] 次にnearc4が操作変数としての条件2をクリアするか検討しましょう。nearc4がIQと関係してしまうと、IQは「個人の能力」の代理変数とみなせますから、この経路を通して撹乱項εと相関してしまいます。まとめると、 educ(教育年数)⇔能力(観測されない)⇔撹乱項 なので、操作変数nearc4を使います。nearc4が撹乱項と直交条件を満たすために、 nearc4(大学から近いところに住んでいるか？)⇔能力(観測されない)⇔撹乱項 において、nearc4と能力（の代理変数であるIQ)が相関しないか、チェックしないといけません。 (能力とnearc4が相関していなくても、撹乱項とnearc4に相関がないということにはならないよ) すなわち、IQがnearc4と関係しないということを確認しないといけません。 ・・・確認しようと思ったら、iqに欠損値があるようで、これをとりのぞかねばなりません。欠損値は、Woolridgeのデータで"."となっているようなので、これをすっ飛ばしてデータフレームを再定義するような関数をつくらないといけません。 (read.tableにはオプションにna.stringsがあって、読み込みの際に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 ' ' 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.005373 F-statistic: 12.13 on 1 and 2059 DF, p-value: 0.0005071 nearc4 に対してiqで単回帰すると、iqのt値は有意になりました。一見，条件2はクリアされていないかに見えます．ところが、nearc4に対してiq,smsa66,reg661,reg662,reg669で重回帰すれば、iqのt値は有意ではなくなります． > summary(lm(nearc4 ~ iq+smsa66+reg661+reg662+reg669)) とやれば，結果は，以下です． Call: lm(formula = nearc4 ~ iq + smsa66 + reg661 + reg662 + reg669) 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 ' ' 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.1939 F-statistic: 100.1 on 5 and 2055 DF, p-value: < 2.2e-16 もともと，単回帰分析を考えていたわけではにので，このような重回帰におけるiqのt値が有意ではないこと，というのが正しい条件です． 条件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のための操作変数をZ1,Z2とした場合、 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 + south + + smsa + reg661 + reg662+reg663+reg664+reg665+reg666+reg667 + +reg668+smsa66, + ~ nearc4 + exper + expersq + black + south + smsa + reg661 + + reg662+reg663+reg664+reg665+reg666+reg667+reg668+smsa66)) 結果は以下のようになります。 2SLS Estimates Model Formula: lwage ~ educ + exper + expersq + black + south + smsa + reg661 + reg662 + reg663 + reg664 + reg665 + reg666 + reg667 + reg668 + smsa66 Instruments: ~nearc4 + exper + expersq + black + south + smsa + reg661 + reg662 + reg663 + reg664 + reg665 + reg666 + reg667 + reg668 + smsa66 Residuals: Min. 1st Qu. Median Mean 3rd Qu. Max. -1.83e+00 -2.41e-01 2.43e-02 -5.84e-12 2.52e-01 1.43e+00 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推定量と2SLS推定量は結果的に同じになることが知られています。次に、操作変数としてnearc4とnearc2の二つを使って2SLS推定してみましょう。 > summary(tsls(lwage ~ educ + exper + expersq + black + south + + smsa + reg661 + reg662+reg663+reg664+reg665+reg666+reg667 + +reg668+smsa66, + ~ nearc2+nearc4 + exper + expersq + black + south + smsa + reg661 + + reg662+reg663+reg664+reg665+reg666+reg667+reg668+smsa66)) とするだけです。簡単ですね。 2SLS Estimates Model Formula: lwage ~ educ + exper + expersq + black + south + smsa + reg661 + reg662 + reg663 + reg664 + reg665 + reg666 + reg667 + reg668 + smsa66 Instruments: ~nearc2 + nearc4 + exper + expersq + black + south + smsa + reg661 + reg662 + reg663 + reg664 + reg665 + reg666 + reg667 + reg668 + smsa66 Residuals: Min. 1st Qu. Median Mean 3rd Qu. Max. -1.94e+00 -2.51e-01 1.93e-02 -3.00e-12 2.65e-01 1.47e+00 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]の出題元となっているCard(1995)とちゃんと同じになっていることを[[管理人(SUGIYAMA Shunsuke):http://www.sugi-shun.com/]]は確認済みです。 ***データに関する注意 [#iac47a73] COLOR(RED){（注意）：データに関する注意．} expersqとlwageについて、Wooldridgeが用意したデータを使うのではなく、Rに読み込んでから自分で > exper2<-exper^2 > logwage<-log(wage) とかやってRで加工した系列を使うと、微妙に異なる推定結果を得る。Limdepで同様の作業をしても同じ結果を得る。 具体的には， #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 + south + smsa + reg661 + reg662+reg663+reg664+reg665+reg666+reg667 +reg668+smsa66, ~ nearc4 + exper + exper2 + black + south + smsa + reg661 + reg662+reg663+reg664+reg665+reg666+reg667+reg668+smsa66)) とかやれば， 2SLS Estimates Model Formula: logwage ~ educ + exper + exper2 + black + south + smsa + reg661 + reg662 + reg663 + reg664 + reg665 + reg666 + reg667 + reg668 + smsa66 Instruments: ~nearc4 + exper + exper2 + black + south + smsa + reg661 + reg662 + reg663 + reg664 + reg665 + reg666 + reg667 + reg668 + smsa66 Residuals: Min. 1st Qu. Median Mean 3rd Qu. Max. -1.84e+00 -2.45e-01 2.58e-02 2.17e-11 2.54e-01 1.43e+00 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のデータそのものが間違っていると考えざると得ないのだが、これを使うとCardと同じ結果を得るというのは、一体どいういうことか？Cardから提供されたデータをWooldridgeがそのまま配布しているとしたら、Cardがもともと対数をとったり二乗したりするデータ加工を少し間違えたということか？例えば、CardがExcelを使ったが、そのExcelがバカなのか？ ・・・このように，計量分析では，論文の追試をやると，微妙に異なる結果が得られることが多々あります．それは，使ったデータ・使った計量ソフト・データ整備で間違えた，などなどいろいろな理由が考えられます．基本的に，分析者のITリテラシー不足or性格がいい加減（ヒューマンエラーを犯しやすい）のどちらかでしょう．そうはいってもいくら気をつけていても，ピッタリ一致することのほうが少ないかもしれません．科学は再現性が大事なのに，微妙に計量結果が異なってしまうことがあるなんて，計量経済学はなんていい加減なんだ，という批判は真摯に受け止めるべきだと思います．この批判を回避するためには，用いたデータ，データの出所，用いたプログラムなどを，出せといわれたらいつでも出せるように準備しておくことが良いと思います．とても細かいことですが，このPCで1ヶ月前にやったときは，こういう結果が得られていたんだ，と主張しても，誰も耳を貸しません．なるべく詳細に，自分が行った計量分析の手順を記録すべきです． ***補足 [#nf96aed0] systemfit パッケージ：systemfit でも2SLS推定が出来るようです。2SLSだけでなく、3SLSとかも出来るようです。詳しくは、 ?systemfit で見てください。 ***過剰識別制約テスト [#u19330c8] COLOR(Fuchsia){（保留）：操作変数の過剰識別制約テストの使い方はそのうち更新します。Hausman原理でも出来るし，単にt値を見る方法もある．} **HCSE(Heteroskedasticity Consistene Standard Error) [#x011ee85] RIGHT:参考："[[Econometrics in R:http://cran.r-project.org/doc/contrib/Farnsworth-EconometricsInR.pdf]]の3.3 Heteroskedasticity and Friendsと3.4 Linear Hypothesis Testing (Wald and F)" ***White[#i69bbbc1] 分散不均一が発生してもOLS推定量は不偏性（一致性）は満たします。失うのは最小分散性だけです。その理由は、OLS推定量の共分散行列の推定が間違うからです。間違うならば、正せばよい、というのがWhite修正された共分散行列で、これは一致性をもつことが証明されています。いわゆるWhiteの一致標準誤差とか呼ばれるものです。 コマンドは、 vcovHC パッケージ:sandwich もしくは hccm パッケージ：car です。どちらを使っても同じ結果が得られることをちゃんと確認してみましょう。 使うデータは、[[最小二乗法による推定:http://www.sugi-shun.com/econwiki/index.php?R%A4%F2%BB%C8%A4%C3%A4%C6%B7%D7%CE%CC%B7%D0%BA%D1%CA%AC%C0%CF#u30dda11]]で使った投資関数のデータです。考えるモデルは、&mimetex(I_t=a+bY_t+cr_{t-4});としておきましょう。本当は、分散不均一はこのモデルでは発生していなかったので、white修正する必要はないわけですが、ここでは演習と割り切りましょう。 [[lm()が自動出力するt検定量とF検定量を自分で計算してみよう :http://www.sugi-shun.com/econwiki/index.php?R%A4%F2%BB%C8%A4%C3%A4%C6%B7%D7%CE%CC%B7%D0%BA%D1%CA%AC%C0%CF#ke1db65d]]では、Yのt値を自分で計算してみました。これをもう一度やりましょう。 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 ' ' 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値は分散均一の仮定のもとで計算されたものです。この仮定をとっぱらったとき、すなわちWhite修正された共分散行列を使うときのt値を計算してみましょう。 まず、修正した共分散行列を使って仮説検定する方法は > lht(unrestricted,hm2,rhs2,vcov=vcovHC(unrestricted) ) あるいは > lht(unrestricted,hm2,rhs2,vcov=hccm(unrestricted) ) です。欲しいのはカイ二乗統計量の平方根なので、 > lht(unrestricted,hm2,rhs2,vcov=hccm(unrestricted),test=c("Chisq"))$Chi[2]->whitese > sqrt(whitese) もしくは > lht(unrestricted,hm2,rhs2,vcov=vcovHC(unrestricted),test=c("Chisq"))$Chi[2]->whitese2 > sqrt(whitese2) と入力します。hccm()でもvcovHC()でも [1] 6.390953 が得られました。これがWhiteの一致標準誤差を使った場合のt値です。 尚、White修正にはいろいろなバージョンがあります。 詳しくは、 > ?hccm > ?vcovHC で見てください。ここの"type"というサブコマンドをどのように指定するかで、White修正のどのバージョンを使うかを指定できます。 hccm()を使った場合は以下です。 > sqrt(lht(unrestricted,hm2,rhs2,vcov=hccm(unrestricted,type=c("hc0")),test=c("Chisq"))$Ch[2]) [1] 6.741862 > sqrt(lht(unrestricted,hm2,rhs2,vcov=hccm(unrestricted,type=c("hc1")),test=c("Chisq"))$Ch[2]) [1] 6.644852 > sqrt(lht(unrestricted,hm2,rhs2,vcov=hccm(unrestricted,type=c("hc2")),test=c("Chisq"))$Ch[2]) [1] 6.564739 > sqrt(lht(unrestricted,hm2,rhs2,vcov=hccm(unrestricted,type=c("hc3")),test=c("Chisq"))$Ch[2]) [1] 6.390953 > sqrt(lht(unrestricted,hm2,rhs2,vcov=hccm(unrestricted,type=c("hc4")),test=c("Chisq"))$Ch[2]) [1] 6.374456 vcovHC()を使った場合は以下です。 > sqrt(lht(unrestricted,hm2,rhs2,vcov=vcovHC(unrestricted,type="HC"),test=c("Chisq"))$Ch[2]) [1] 6.741862 > sqrt(lht(unrestricted,hm2,rhs2,vcov=vcovHC(unrestricted,type="HC0"),test=c("Chisq"))$Ch[2]) [1] 6.741862 > sqrt(lht(unrestricted,hm2,rhs2,vcov=vcovHC(unrestricted,type="HC1"),test=c("Chisq"))$Ch[2]) [1] 6.644852 > sqrt(lht(unrestricted,hm2,rhs2,vcov=vcovHC(unrestricted,type="HC2"),test=c("Chisq"))$Ch[2]) [1] 6.564739 > sqrt(lht(unrestricted,hm2,rhs2,vcov=vcovHC(unrestricted,type="HC3"),test=c("Chisq"))$Ch[2]) [1] 6.390953 > sqrt(lht(unrestricted,hm2,rhs2,vcov=vcovHC(unrestricted,type="HC4"),test=c("Chisq"))$Ch[2]) [1] 6.374456 ここで、"hc0"(hccm()を使った場合),"HC"もしくは"HC0"（vcovHC()を使った場合）がWhiteが原論文（White(1980) in Econometrica）で提唱した修正共分散行列です。その他のTypeは、なんかその後いろいろなWhite修正を微修正すべきと再提案されたものです。（See Greene(2003)[pp.220]） 尚、vcovを使った場合、typeとして"const"を指定できます。これはconstantの略で、分散均一の仮定を置くということになります。だから、t値はもともとの値と同じになります。 > sqrt(lht(unrestricted,hm2,rhs2,vcov=vcovHC(unrestricted,type="const"),test=c("Chisq"))$Ch[2]) [1] 8.553713 ここの項目は、"[[Econometrics in R:http://cran.r-project.org/doc/contrib/Farnsworth-EconometricsInR.pdf]]の3.3 Heteroskedasticity and Friendsと3.4 Linear Hypothesis Testing (Wald and F)" を参考にしました。 ***Newey-West [#s56bb41a] Newey-Westの方法は、White修正の「系列相関があっても大丈夫だよ」バージョンです。 コマンドは、 NeweyWest パッケージ：sandwich です。 **GLS(Generalized Least Squares) [#sc894fcf] RIGHT:参考：浅野・中村(2000)[pp.130] 日本語では一般化最小二乗法です。[[管理人(SUGIYAMA Shunsuke):http://sugi-shun.com]]はRでGLSの実戦経験がないので、説明できません。 そのうち、自分でプログラムを書くなり、組み込み関数を探すなりして、加筆修正します。 以下は覚書です。関連しそうなのは、 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%C7OLS%B0%CA%B3%B0%A4%CE%BF%E4%C4%EA%CA%FD%CB%A1#ebd958a0]]を見てください。 *時系列分析 [#mf5cad5a] [[Rで時系列分析]]を見てください。 *計量ミクロ [#ced0ba2a] [[Rで計量ミクロ]]を見てください。 *OLS以外の推定方法 [#x0baa89f] [[RでOLS以外の推定方法]]を見てください。 *Reference [#h7afec48] **書籍 [#h8a36313] -[[浅野・中村(2000):http://www.amazon.co.jp/gp/product/4641160805/sr=8-1/qid=1155206048/ref=sr_1_1/250-7742157-1058652?ie=UTF8&s=gateway]] -[[Greene(2003):http://www.amazon.co.jp/gp/product/0131108492/sr=8-1/qid=1155281794/ref=pd_bowtega_1/250-7742157-1058652?ie=UTF8&s=gateway]] -[[Hamilton(1994):http://www.amazon.co.jp/gp/product/0691042896/ref=pd_sim_gw_4/250-7742157-1058652?ie=UTF8]] -[[Wooldridge(2002):http://www.amazon.co.jp/Econometric-Analysis-Cross-Section-Panel-Data/dp/0262232197/ref=pd_bxgy_b_img_b/503-7013025-2018353?ie=UTF8]] **リンク [#f4df0207] -[[Econometrics in R:http://cran.r-project.org/doc/contrib/Farnsworth-EconometricsInR.pdf]] -[[R for economists:http://www.mayin.org/ajayshah/KB/R/R_for_economists.html]] -[[A Brief Guide to R for Beginners in Econometrics:http://people.su.se/~ma/R_intro/]] -[[CRAN Task View: Computational Econometrics:http://cran.md.tsukuba.ac.jp/src/contrib/Views/Econometrics.html]] -[[パッケージ一覧:http://cran.md.tsukuba.ac.jp/src/contrib/PACKAGES.html]]
