(IV)操作変数法を実際にRで走らせる方法をご紹介します。

まず操作変数法のアイデアを軽く復習すると、
内生性のある説明変数($x$)を、外生変数(IV)($z$)で説明できる部分とできない部分に切り分け、説明できる部分だけOLSに用いることによって先の内生変数を外生変数として扱い、OLS推定量($\hat\beta$)の一致性を確保しようということでした。

具体的な手法としては1st stageとして内生変数をIVにOLS回帰、2nd stage として 1st stage のfitted valueを元々の目的の回帰関数の説明変数として加えてOLS回帰という、一般的には Two Stage Least Squared (二段階最小二乗法、TSLS)と呼ばれる手法を用います。

つまり適切なIVが見つけられさえすれば、あとはOLSを二回実行するだけです。

＊ここでいう「内生性(endogeneity)」、「外生性(exogentity)」とは、一般的な経済学で言うものとは違い、計量経済学では説明変数がerror term($u$)と相関を持つ場合内生、持たない場合外生と言うことに注意してください。状況によっては両方の意味合いを含む場合もあります。(需要曲線、供給曲線を推定する同次方程式モデルの場合など)

###TSLSを走らせるための一般的状況と仮定(条件)を確認
簡単な単回帰の場合だと
$$y_i=\beta_0+\beta_1x_i+u_i$$

このようなPopulation Regression Function において

$$Cov(x,u)\neq0$$

つまり$x$は内生変数であるという状況で、IV($z$)の条件は、

1. Instrument Relevance

$$Cov(z,x)\neq0$$
$z$は$x$に対して説明力を持つこと。

2. Instrument Exogeneity

$$Cov(z,u)=0$$
$z$とerrorの間には相関なし。
相関を持つなら、$z$では内生変数$x$の外生的な部分だけ取り出すのは無理。

でした。数学的にどこに聞いてくるのかはのちに説明します。

3. Over Identification

$$IVの数\geq内生変数の数$$

これはIVの条件ではないがOLS推定量($\hat\beta$)を計算するための一階条件。これが成り立っていないと、条件式の数がパラメータ($\beta$)の数より少なくなってしまい、$\hat\beta$はデータから一意に決まらなくなってしまいます。

これらの条件は今回はすべて満たされるものとします。

###今回の目的
$$\log (wage_i)=\beta_0+\beta_1educ_i+u_i$$

$$wage_i:n人中i番目の個人の実質賃金,educ:n人中i番目の個人の教育年数 $$
という定式化が正しいという前提のもとで、教育年数が一年上がると実質賃金が何%(logを取っているので)上がるかを考えたいとします。実際は教育年数はどの一年を増やすかによって賃金の上昇率は上がりそうですが、簡単化のために今回はどこの一年を増やしても同じという定式化で行きます。

ここで$educ$は内生変数、つまりerror($u$)と相関がありそうです。
$u$の中には住んでいる地域、出会った学校の先生の質など、データがなかったとして、$wage$に説明力をもち、$educ$とも相関をもちそうなものが含まれていると考えられるからです。この場合この回帰式の$\beta$をOLSで推定しても一致性は満たされません。

さらにIVを父親の教育年数($fatheduc$)として、データは利用可能だったとします。
父親の教育年数は子供の教育年数と相関がありそうなのでIVの１つ目の条件($Cov(z,x)\neq0$)は満たしていて、実際微妙ですが子供の賃金に影響を与えていて説明変数に加えられていないものと相関はなかったとし、２つ目の条件($Cov(z,u)=0$)も満たされているとします。

###IV estimatorのRでの計算方法は4種類
①両辺IVとのCovarianceとる
②コマンドのivreg
③手動で1st stage とsecond stage を実行
④コマンドのtsls

Wooldridge の Introductory Econometrics の公開されているデータセットの中のmrozというデータを用います。

準備を先にしておきます。

library(foreign);library(AER);library(stargazer)
mroz<-read.dta("http://fmwww.bc.edu/ec-p/data/wooldridge/mroz.dta")
sampleset <- subset(mroz, !is.na(wage))

ではまず①から見ていきましょう。

####①両辺IVとのCovarianceとる(単回帰の時のみ)

with(sampleset, cov(log(wage),fatheduc) / cov(educ,fatheduc) )

[1] 0.05917348

$$Cov(fatheduc_i,\log (wage_i))=\beta_0Cov(fatheduc_i,1)+\beta_1Cov(fatheduc_i,educ_i)+Cov(fatheduc_i,u_i)$$

$$\beta_1=\frac{Cov(fatheduc_i,\log (wage_i))}{Cov(fatheduc_i,educ_i)}+\frac{Cov(fatheduc_i,u_i)}{Cov(fatheduc_i,educ_i)}$$

ここで条件の1が満たされていないと、第一項も第二項も分母を0で割ってしまうことになりますし、条件の２つ目が満たされていないと、そもそも第二項の分子は観測できない$u_i$を使っていることから、推定量は計算できないことになります。
よって条件が満たされていると次のようになります。
$$\beta_1=\frac{Cov(fatheduc_i,\log (wage_i))}{Cov(fatheduc_i,educ_i)}$$

0.05917348は母集団におけるこの式をそのままサンプルに当てはめた結果です。(Covarianceの推定量は$\frac{1}{n}\sum_{i=1}^{n} (fatheduc_i-\overline{fatheduc})(educ_i-\overline{educ})$などを使っているので不偏性を持ちませんが、そもそもIV推定量は構造上不偏性を持ち得ません。なので最低限の一致性を満たすような推定方法としてこの①の方法も悪くないです。)
$$\hat\beta_1=0.05917348$$

####②コマンドのivreg

$$1st~stage :educ_i=\pi_0+\pi_1fatheduc_i+v_i$$

$$2nd~stage :\log(wage_i)=\beta_0+\beta_1\hat{educ}_i+u_i$$

のIV回帰を自動でやってくれるコマンドです。一致性を持たないOLSとどれぐらい違うのか比べてみます。

reg.ols <-   lm(log(wage) ~ educ, data=sampleset)
reg.iv <- ivreg(log(wage) ~ educ | fatheduc, data=sampleset) 
stargazer(reg.ols,reg.iv, type="text")

===================================================================
                                       Dependent variable:         
                               ------------------------------------
                                            log(wage)              
                                         OLS           instrumental
                                                         variable  
                                         (1)               (2)     
-------------------------------------------------------------------
educ                                  0.109***            0.059*   
                                       (0.014)           (0.035)   
                                                                   
Constant                               -0.185             0.441    
                                       (0.185)           (0.446)   
                                                                   
-------------------------------------------------------------------
Observations                             428               428     
R2                                      0.118             0.093    
Adjusted R2                             0.116             0.091    
Residual Std. Error (df = 426)          0.680             0.689    
F Statistic                    56.929*** (df = 1; 426)             
===================================================================
Note:                                   *p<0.1; **p<0.05; ***p<0.01

$$\hat\beta_1=0.05917348$$
で①の時と同じで、0.109と推定してるOLSとはだいぶ違いますね。OLSでは、教育年数の効果を過大評価していたことになります。

####③手動で1st stage と 2nd stage を実行

③、④では、より一般的な形として、外生変数$exper$と$exper^2$を説明変数に加え、IVも$fatheduc$に加え$motheduc$(母親の教育年数)も加えます。ここでIVは$fatheduc$と$motheduc$だけでもいいのですが、composit IVとして$exper$と$exper^2$を加えたものを１つのIVとする方がefficiency(有効性)が高いのでcomposit IVを使うことにします。

### 1st stage: reduced form(誘導形)

stage1 <- lm(educ~exper+I(exper^2)+motheduc+fatheduc, data=sampleset)

### 2nd stage
(man.2SLS<-lm(log(wage)~fitted(stage1)+exper+I(exper^2), data=sampleset))

Call:
lm(formula = educ ~ exper + I(exper^2) + motheduc + fatheduc, 
    data = sampleset)

Coefficients:
(Intercept)        exper   I(exper^2)     motheduc     fatheduc  
   9.102640     0.045225    -0.001009     0.157597     0.189548  


Call:
lm(formula = log(wage) ~ fitted(stage1) + exper + I(exper^2), 
    data = sampleset)

Coefficients:
   (Intercept)  fitted(stage1)           exper      I(exper^2)  
      0.048100        0.061397        0.044170       -0.000899

####④コマンドのtsls

library(sem)
(aut.2SLS<-tsls(log(wage)~educ+exper+I(exper^2),instruments=~motheduc+fatheduc+exper+I(exper^2),data=sampleset))

Model Formula: log(wage) ~ educ + exper + I(exper^2)

Instruments: ~motheduc + fatheduc + exper + I(exper^2)

Coefficients:
  (Intercept)          educ         exper    I(exper^2) 
 0.0481002982  0.0613966289  0.0441703937 -0.0008989696

$$\hat\beta_1=0.0613966289$$

で③と同じですね。

＊notationは計量経済学の世界標準的なものを利用していますが、正確な定義をここでは記載していないのでもし気になるところがあればご質問ください。

##参考
Introductory Econometrics Chapter15