Rを使って計量経済分析

概要:このページは、統計言語Rで計量経済分析する方法についての解説ページです。暇をみつけて更新していきます。

関連するサイトは、リンクを見てください。

親ページ:このページの親ページはRです。

---

古典的線形回帰モデル †

最小二乗法による推定 †

とりあえず、RでOLSする方法を学びます。

例として、こんな問題を考えて見ましょう。

いや、あくまで例ですよ。

以下のデータを使います。

このデータは、日本の企業設備投資と国内総支出と長期貸し出しレートのデータです。頻度は四半期データで、投資と支出は実質表示、観測期間は1970年第1四半期-1997:第1四半期です。

見出しの意味は、

です。

とりあえず、

> 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にこのデータを読み込ませましょう。（注意）：ディレクトリの変更を忘れずに。

次に、

> attach(i)

と入力しましょう。こうすると、見出しの変数名を読み込んでくれます。たとえば、

> INV

などと打つと、企業設備投資のデータが出力されることを確認しましょう。

単回帰 †

とりあえず、単回帰してみましょう。

投資(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

と入力しておきましょう。これで長さがそろました。

（注意）：Rでは、大文字と小文字を区別します。だから、iとIは別物とRは認識してくれます。

次に、単回帰式, を推定してみましょう。

まず、を推定してみましょう。

> lm(I ~ Y)

と入力してみましょう。すると、

Call:
lm(formula = I ~ Y)

Coefficients:
(Intercept)            Y  
     -2.894        1.947 

と出力されます。

このコマンドの意味は、

> lm(被説明変数　~　説明変数)

です。結果を見ればわかるとおり、定数項ありモデルを推定しています。定数項なしモデルを推定したい場合は、

> lm(I ~ 0 + Y)

と入力します。

この結果は、と推定されたことを意味します。計量経済分析では、係数の推定値だけではなく、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 

という結果が返ってきます。

では次に、を推定してみましょう。

> 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 

と出力されます。

尚，どうして伸び率データを使うかというと，いま使っている原データは単位根を持つことが広く知られているからです．

重回帰分析 †

単回帰のやり方を学んだので、次に重回帰のやり方を覚えましょう。といっても、単回帰とほとんど一緒です。

たとえば、を推定することを考えましょう。

以下のように入力してやりましょう。

> 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パッケージにありますから、いちいちパッケージの読み込みなどしなくてもいいです。

予測 †

演習として、予測もやってみましょう。

> 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推定された直線を表してくれることがわかります。

ここまでの推定結果の解釈 †

1. まず、投資関数とかいう名前で習うようなというモデルにおいて、のt値は有意ではありません。われわれは普通、「投資は利子率の減少関数」と習いますが、実際のデータではこの関係は認めらないということになります。
2. もう一つの単回帰式については、のt値は有意です。符号は正です．GDPが上昇するとき、投資も上昇するよ、ということが観察されました。
3. 重回帰式の推定結果を見ると、,の両方についてt値が有意になりました。これは、「利子率は、ほかの事情一定としないで見ると、投資に対して有意な負の影響を持たない。しかし、GDPという要因を一定にしたうえで見ると、投資に対して有意な負の影響を持つ」ということです。
4. さらに、の係数の大きさを見てみると、重回帰のときは-0.7051で、単回帰のときは-0.3836であり、単回帰は重回帰よりも大きくなっている。ここで一般に、利子率が高いときは、景気がいいときであることを思い出しましょう。すると、「r↑のときはY↑という経路からくる影響（Y↑のときはI↑となる影響）のため，重回帰のrの係数より単回帰のrの係数が大きくなっている」と言えそうです。

より良い計量モデルを考えてみる †

ここまでの推定結果の解釈の1.として述べた点について、次のような反論が考えられます。

「企業が今年実行する設備投資計画は、去年、計画したものだ。だから、t期の投資は、t期より以前の利子率に影響されるだろう」

というわけで、ではなく、のようなモデリングをするべきだ、という考え方があります。実際にこれを推定してみましょう。

まず、""に対応する変数を定義しないといけませんね。

> 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も回帰式に入れてみたら、「よりよさそう」な推定結果が得られていました。だから、というモデルを考えて見ましょう。

> 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値はさらに有意になりました。符合条件も、理論通りです。決定係数も、ここまでの試したモデルの中で、一番高いです。*1どうやらこれが現実をよく捉えた式なのではないか、という風に考えてよさそうです。

仮説検定 †

ここでは、LSE(Linear Squares Estimation)における仮説検定について記述します。Wald test, LM test, 制約付残差に基づくF検定が有名です。また、これら三つは、結果的に同じ検定量となることも有名です。

lm()が自動出力するt検定量とF検定量を自分で計算してみよう †

lm()が自動的にやってくれる仮説検定は、係数が有意に0になるかどうかのt検定とF検定です。これを例えばRの自動出力にまかせずに自分で計算する方法を説明します。

waldtest 
パッケージ：lmtest

というコマンドがあります。Wald検定はこれを使うはずなんですが、管理人(SUGIYAMA Shunsuke)はこのコマンドを実践で使ったことがありません。使い方がわかる人は教えてください。

（保留）：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の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 

となります。ここでわれわれが興味があるのは、分布に従う統計量の73.166だけです。これだけを取り出すのには、

> lht(unrestricted,hm2,rhs2,test=c("Chisq"))$Chi[2]

とやればよいです。これは自由度が1の分布する検定量なわけですから、これの平方根をとる統計量は標準正規分布します。

平方根を計算してやると、

> sqrt(lht(unrestricted,hm2,rhs2,test=c("Chisq"))$Chi[2])
[1] 8.553713

を得ますが、これはYの係数のt値にほかなりません。

ここまでやってきたことで、実は次のようなことがわかりました。すなわち、いつも見ているt検定量やF検定量も、実はWald検定の特殊ケースでしかないということです。次にWald検定について説明します。

WALD検定 †

さて、すでに述べたとおり、Wald test, LM test, 制約付残差に基づくF検定の三つの仮説検定量は、それぞれ違う考え方に基づいて考案された仮説検定量ですが、結果的に同じ計算式になることがわかっています。そして、統計学の知識を使えば、F分布に従う統計量と、分布に従う統計量の二つのタイプの仮説検定量を計算することが出来ます。

結果のみ示すと、

分布に従う統計量 = 制約の数 × 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ですから、確かに分布に従う統計量(86.357)は、制約の数(2) × F分布に従う統計量(43.178)が成立しています。

回帰診断 †

OLS推定したら、「この回帰分析は成功なのか、失敗なのか」ということを考えなければなりません。それは、「古典的仮定は満たされているの？満たされていない場合、どうしたらいいの？」ということを考えることに等しいです。これが回帰診断です。

以下の表の「判断方法」に記された検定などの結果、全ての仮定が満たされていると判断された場合、その回帰分析は成功であると考えられます。満たされないものがあった場合、満たされるように、工夫しましょう。

ここで工夫とは、「古典的仮定が満たされていなかった場合、変数を追加したり減らしたり、あるいはデータ加工をするなどして、古典的仮定がすべで満たされているような推定式を探す」ということです。*2

その上で、どう工夫しても改善できなかった場合、表の一番右に記した「対処法」をとりましょう。

表では、古典的仮定については、浅野・中村(2000)[pp.23]に基づきました。

なお、仮定１～仮定４のすべてが崩れた場合は、GMMという最終手段が待っています。

それでは、実際にこれらのチェックをしてみましょう。 とりあえず、というモデリングがいまのところの暫定チャンピオンです。以下では、このモデルをベースに議論します。

説明変数の外生性のチェック †

hausman.systemfit {systemfit}

を使います． （保留）：（説明変数と撹乱項の直交条件）。そのうち加筆予定。Hausman testについて説明します。

分散均一性のチェック †

図を描いて目視 †

横軸に説明変数、縦軸に残差をプロットしてみて、目でまずは確認してみましょう。

で、残差をプロットする方法を説明します。その前に、残差系列を取り出す方法を知らないといけません。

> ?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で見てみたら、このページの"[質問] lm() 関数の返り値のF統計量のp-value <= なんでも掲示板から移動"に解決法が書かれていました。

さて、話をもとに戻します。横軸に説明変数、縦軸に残差をプロットするには、

> plot(r1,lm(I ~ Y + r1)$res)

あるいは、

> plot(Y,lm(I ~ Y + r1)$res)

とすればよいです。この図の結果を見てみると、まぁ微妙。分散は均一のような、不均一のような、目視では微妙です。管理人には、あるように見えます。やっぱりないかな？わかりません。そこで、ちゃんと統計的にこの判断をする方法を以下に記します。

Breush Pagan test †

コマンドは、

bptest()
(パッケージ：lmtest)

です。

具体的に、というモデルにおけるBP検定を行う方法は、

> bptest(I ~ Y + r1)

あるいは、

> bptest(lm(I ~ Y + r1))
> bptest(summary((lm(I ~ Y + r1))))

でも同じ結果が返されるようです。

（注意）：パッケージの読み込みを忘れずに。

これを実行すると、

        studentized Breusch-Pagan test

data:  summary((lm(I ~ Y + r1))) 
BP = 1.8808, df = 2, p-value = 0.3905

という結果が返ってきます。よって、"H0:分散は均一である"は棄却されず、これを受容し、分散は均一であると判断されました。

（保留）：studentizedってなんですか？疑問です．

White Test †

コマンドは、

white.test()
(パッケージ：tseries)

です。

具体的に、というモデルにおけるwhite.test検定を行う方法は、わかりません。管理人(SUGIYAMA Shunsuke)は、Rでwhite.testの実践経験がないので、説明できません。

Goldfeld-Quandt test [written by Akihiko Noda] †

コマンドは、

gqtest()
(パッケージ：lmtest)

です。

基本的な考え方としては、存在するデータを2つに分割して、それらを用いて同じモデルを回帰した場合に誤差項の分散が等しくなるかどうかを検定しています。

具体的に、というモデルにおけるGQ検定を行う方法は、

> gqtest(I ~ Y + r1)

あるいは、

> gqtest(lm(I ~ Y + r1))
> gqtest(summary((lm(I ~ Y + r1))))

でも同じ結果が返されるようです。

（注意）：パッケージの読み込みを忘れずに。

これを実行すると、

       Goldfeld-Quandt test

data:  I ~ Y + r1 
GQ = 0.6811, df1 = 50, df2 = 49, p-value = 0.9102

という結果が返ってきます。よって、"H0:2つの分散は等しい"は棄却されず、これを受容し、分散は均一であると判断されました。

系列相関のチェック †

図を描いて目視 †

図の描き方として、二つあります。一つ目としては、横軸に残差、縦軸に1期ずれの残差をプロットしてみましょう。系列相関がなければ、このプロット図は、円状になるはずです。3次元プロットでは、球状になります。

このプロットの仕方は、分散均一性のチェックのところで説明した知識を使えば、わかるはずです。

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

コマンドは、

dwtest()
(パッケージ：lmtest)

もしくは、

durbin.watson
(パッケージ：car)

です。

具体的に、というモデルにおけるDW検定を行う方法は、

> dwtest(I ~ Y + r1)

あるいは、

> dwtest(lm(I ~ Y + r1))
> dwtest(summary((lm(I ~ Y + r1))))

でも同じ結果が返されるようです。

（注意）パッケージの読み込みを忘れずに。

これを実行すると、

        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 †

コマンドは、

bgtest()
(パッケージ：lmtest)

です。

具体的に、というモデルにおけるDW検定を行う方法は、

> bgtest(I ~ Y + r1)

あるいは、

> bgtest(lm(I ~ Y + r1))
> bgtest(summary((lm(I ~ Y + r1))))

でも同じ結果が返されるようです。

（注意）パッケージの読み込みを忘れずに。

これを実行すると、

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

コマンドは

> Box.test()

です。管理人(SUGIYAMA Shunsuke)は、RでBox.testを使った分散均一性のチェックの実践経験がないので、説明できません。

RESETテスト †

resettest {lmtest} 

を使います．

（保留）：そのうち加筆予定。関数形が正しいとdefendするために行う検定。

ここまでの回帰診断のまとめ †

どうやら、分散不均一性は発生していないようです。しかし、系列相関が発生してしまっているようです。したがって、モデリングを工夫して、「系列相関を消す」ことを考えるのが一つの戦略です。いくら工夫しても系列相関が消せない場合は、最終の対処法を取ります。この場合は、GLS(Cochrane-Orcutt法など)という推定方法を使うことになります。GLSについては、GLSを見てください。

しかし系列相関が発生してもOLS推定量は最小分散性を失うだけです。一致性は失われないので問題なし、とここでは結論付けることにします。

（保留）：Hausman test, RESET についてそのうち言及予定。

どう工夫しても仮定が崩れるときの対処法 †

回帰診断で記した表をおさらいしましょう。ここで、5個の仮定は重要な順に並んでいると考えてよいです。

仮定1,仮定2は、OLS推定量が一致性を持つために必要です。仮定2はいつも満たされると考えてしまうので、実際には仮定1の妥当性を考えればよいです。しかしこれには計量経済学者はとても敏感です。

仮定3,仮定4は、OLS推定量が最小分散性（有効性＝効率性）を持つために必要です。これらの仮定はしばしば崩れますが、一致性を確保できるならば問題ない、という立場もありますから、それほど敏感にならなくてもよいかもしれません。

仮定5はなくても漸近的に何も問題ありません。直感的には、どうせCLT(中心極限定理)によってデータ数が増えれば正規分布に分布収束するからです。

考慮すべきは、仮定1,仮定3,仮定4の三つということになります（もう少しすっきり表現すると、仮定3と仮定4は合体させて、撹乱項の共分散行列が(ここでは単位行列)である、と記述できます）。

というわけで、古典的仮定が崩れたときの対処法をコンパクトにまとめると、こうなります。

これらの対処法をRでどうやってやるのかを以下で説明したいと思います。

2SLS(IV):2Step Least Squares(Instrument Variables) †

日本語では二段階最小二乗法(操作変数法)のことです（IV≠TSLS＝2SLS。結果が同じになるだけ(･∀･)ﾆﾔﾆﾔ。参考：(；ﾟдﾟ) R言語で同時方程式を推計せよ。）。

例として、ここではWooldridge(2002)[Problem 5.4]をとりあげます。使うデータは

です。このデータについての説明は

を見てください。

このデータは、Wooldridgeのウェブサイトからとりました。

とりあえず、

> 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の外生性については、説明変数の外生性のチェックで実行したHausman検定の結果を思い出しましょう）。

このときOLS推定量は一致性を失いますから、上記の推定結果はぜんぜん信用できない、ということになります。

そこで操作変数の登場です。ここでは、nearc4（4年制大学から近いかどうか）という変数を、操作変数として考えます。Wooldridgeがnearc4を操作変数として考えなさい、といってくれています。まぁ、演習問題ですから。nearc4が操作変数としてふさわしい条件を満たすかどうか、チェックしないといけません。このチェックすべき条件とは、

1. 説明変数educとnearc4は相関するか？
2. εとnearc4は無相関であるか？

の二つです（これは、厳密な言い方ではありません）。

操作変数(nearc4)は，説明変数(educ)と相関するか？ †

まず一つ目のチェックをします。これは、より厳密には、

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)は，撹乱項(ε)と関係しないか？ †

次にnearc4が操作変数としての条件2をクリアするか検討しましょう。nearc4がIQと関係してしまうと、IQは「個人の能力」の代理変数とみなせますから、この経路を通して撹乱項εと相関してしまいます。まとめると、

educ(教育年数)⇔能力(観測されない)⇔撹乱項

なので、操作変数nearc4を使います。nearc4が撹乱項と直交条件を満たすために、

nearc4(大学から近いところに住んでいるか？)⇔能力(観測されない)⇔撹乱項

において、nearc4と能力（の代理変数であるIQ)が相関しないか、チェックしないといけません。 (能力とnearc4が相関していなくても、撹乱項とnearc4に相関がないということにはならないよ)

すなわち、IQがnearc4と関係しないということを確認しないといけません。

・・・確認しようと思ったら、iqに欠損値があるようで、これをとりのぞかねばなりません。欠損値は、Woolridgeのデータで"."となっているようなので、これをすっ飛ばしてデータフレームを再定義するような関数をつくらないといけません。 (read.tableにはオプションにna.stringsがあって、読み込みの際にna.strings="."をすれば欠損値として読み込んでくれる)

とりあえず，欠損値をすべて「空白」に置き換えた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する †

そこでいよいよ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 ))

とやればよいです。

ここで．

のデータにまた戻ります．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)は確認済みです。

データに関する注意 †

（注意）：データに関する注意．

expersqとlwageについて、Wooldridgeが用意したデータを使うのではなく、Rに読み込んでから自分で

> exper2<-exper^2
> logwage<-log(wage)

とかやってRで加工した系列を使うと、微妙に異なる推定結果を得る。Limdepで同様の作業をしても同じ結果を得る。

具体的には，

のように，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ヶ月前にやったときは，こういう結果が得られていたんだ，と主張しても，誰も耳を貸しません．なるべく詳細に，自分が行った計量分析の手順を記録すべきです．

補足 †

systemfit 
パッケージ：systemfit

でも2SLS推定が出来るようです。2SLSだけでなく、3SLSとかも出来るようです。詳しくは、

?systemfit 

で見てください。

過剰識別制約テスト †

（保留）：操作変数の過剰識別制約テストの使い方はそのうち更新します。Hausman原理でも出来るし，単にt値を見る方法もある．

HCSE(Heteroskedasticity Consistene Standard Error) †

White †

分散不均一が発生してもOLS推定量は不偏性（一致性）は満たします。失うのは最小分散性だけです。その理由は、OLS推定量の共分散行列の推定が間違うからです。間違うならば、正せばよい、というのがWhite修正された共分散行列で、これは一致性をもつことが証明されています。いわゆるWhiteの一致標準誤差とか呼ばれるものです。

コマンドは、

vcovHC 
パッケージ:sandwich

もしくは

hccm 
パッケージ：car

です。どちらを使っても同じ結果が得られることをちゃんと確認してみましょう。

使うデータは、最小二乗法による推定で使った投資関数のデータです。考えるモデルは、としておきましょう。本当は、分散不均一はこのモデルでは発生していなかったので、white修正する必要はないわけですが、ここでは演習と割り切りましょう。

lm()が自動出力するt検定量とF検定量を自分で計算してみよう では、Yのt値を自分で計算してみました。これをもう一度やりましょう。

（注意）パッケージ：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の3.3 Heteroskedasticity and Friendsと3.4 Linear Hypothesis Testing (Wald and F)" を参考にしました。

Newey-West †

Newey-Westの方法は、White修正の「系列相関があっても大丈夫だよ」バージョンです。

コマンドは、

NeweyWest
パッケージ：sandwich

です。

GLS(Generalized Least Squares) †

日本語では一般化最小二乗法です。管理人(SUGIYAMA Shunsuke)はRでGLSの実戦経験がないので、説明できません。

そのうち、自分でプログラムを書くなり、組み込み関数を探すなりして、加筆修正します。

以下は覚書です。関連しそうなのは、

lm.gls 
パッケージ:MASS

gls 
パッケージ:nlme

glsD 
パッケージ:Design

あたりか？

WLS:Weighted Least Squares †

日本語では、加重最小二乗法です。

そのうち、自分でプログラムを書くなり、組み込み関数を探すなりして、加筆修正します。

Cochrane Orcutt法 †

そのうち、自分でプログラムを書くなり、組み込み関数を探すなりして、加筆修正します。

GMM †

GMMを見てください。

時系列分析 †

Rで時系列分析を見てください。

計量ミクロ †

Rで計量ミクロを見てください。

OLS以外の推定方法 †

RでOLS以外の推定方法を見てください。

Reference †

書籍 †

• 浅野・中村(2000)

• Greene(2003)

• Hamilton(1994)

• Wooldridge(2002)

リンク †

• Econometrics in R

• R for economists

• A Brief Guide to R for Beginners in Econometrics

• CRAN Task View: Computational Econometrics

• パッケージ一覧