第7章 操作変数法
mydata <- read.csv("ipehd_qje2009_master.csv")実証例 7.1 （単回帰）
(p.275)
操作変数法による推定を行うためにはestimatrパッケージのiv_robust()関数が便利です．iv_robust(y1 ~ y2 + x | x + z, data)のようにして使います．ただし，y1，y2は内生変数で，xは外生変数，zは操作変数を表しています．|の右側に，操作変数だけでなく外生変数も含めることに注意してください．
library(estimatr)
ols1 <- lm_robust(f_rw ~ f_prot, data = mydata, se_type = "stata") # OLS
ols1
## Estimate Std. Error t value Pr(>|t|) CI Lower
## (Intercept) 82.37419182 1.41304363 58.295575 7.716046e-212 79.59720830
## f_prot 0.07998031 0.01611302 4.963708 9.830547e-07 0.04831421
## CI Upper DF
## (Intercept) 85.1511754 450
## f_prot 0.1116464 450
iv1 <- iv_robust(f_rw ~ f_prot | kmwittenberg, data = mydata, se_type = "stata") # IV
iv1
## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper
## (Intercept) 60.4509045 5.1130987 11.822753 2.835629e-28 50.4023889 70.4994201
## f_prot 0.4215664 0.0710932 5.929771 6.049522e-09 0.2818505 0.5612823
## DF
## (Intercept) 450
## f_prot 450実証例 7.2 （標準誤差）
(p.280)
library(estimatr)
iv1 <- iv_robust(f_rw ~ f_prot | kmwittenberg, data = mydata, se_type = "stata") # IV
summary(iv1)
##
## Call:
## iv_robust(formula = f_rw ~ f_prot | kmwittenberg, data = mydata,
## se_type = "stata")
##
## Standard error type: HC1
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
## (Intercept) 60.4509 5.11310 11.82 2.836e-28 50.4024 70.4994 450
## f_prot 0.4216 0.07109 5.93 6.050e-09 0.2819 0.5613 450
##
## Multiple R-squared: -0.9827 , Adjusted R-squared: -0.9871
## F-statistic: 35.16 on 1 and 450 DF, p-value: 6.05e-09実証例 7.3 （2SLS推定）
(pp.295-296)
library(estimatr)
iv2 <- iv_robust(f_rw ~ f_prot + f_young + f_jew + f_fem +
f_ortsgeb + f_pruss + hhsize + lnpop +
gpop + f_blind + f_deaf + f_dumb + f_miss |
kmwittenberg +
f_young + f_jew + f_fem + f_ortsgeb +
f_pruss + hhsize + lnpop + gpop +
f_blind + f_deaf + f_dumb + f_miss,
data = mydata, se_type = "stata")
summary(iv2)
##
## Call:
## iv_robust(formula = f_rw ~ f_prot + f_young + f_jew + f_fem +
## f_ortsgeb + f_pruss + hhsize + lnpop + gpop + f_blind + f_deaf +
## f_dumb + f_miss | kmwittenberg + f_young + f_jew + f_fem +
## f_ortsgeb + f_pruss + hhsize + lnpop + gpop + f_blind + f_deaf +
## f_dumb + f_miss, data = mydata, se_type = "stata")
##
## Standard error type: HC1
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
## (Intercept) 177.3665 29.49620 6.0132 3.838e-09 119.3948 235.33821 438
## f_prot 0.1885 0.02707 6.9640 1.217e-11 0.1353 0.24170 438
## f_young -1.9524 0.16773 -11.6400 1.806e-27 -2.2820 -1.62271 438
## f_jew -0.4367 0.34981 -1.2483 2.126e-01 -1.1242 0.25084 438
## f_fem -1.0726 0.33087 -3.2419 1.278e-03 -1.7229 -0.42234 438
## f_ortsgeb 0.6066 0.05177 11.7170 9.056e-28 0.5048 0.70835 438
## f_pruss -0.1807 0.14177 -1.2746 2.031e-01 -0.4593 0.09794 438
## hhsize 0.8849 1.64334 0.5385 5.905e-01 -2.3449 4.11475 438
## lnpop -1.3181 0.95739 -1.3768 1.693e-01 -3.1998 0.56351 438
## gpop 0.4099 0.10914 3.7558 1.961e-04 0.1954 0.62441 438
## f_blind -27.8647 14.98451 -1.8596 6.362e-02 -57.3152 1.58580 438
## f_deaf -52.3831 10.41183 -5.0311 7.130e-07 -72.8464 -31.91972 438
## f_dumb 7.5290 1.70627 4.4125 1.288e-05 4.1755 10.88248 438
## f_miss -0.5049 0.37778 -1.3366 1.821e-01 -1.2474 0.23755 438
##
## Multiple R-squared: 0.6886 , Adjusted R-squared: 0.6794
## F-statistic: 66.57 on 13 and 438 DF, p-value: < 2.2e-16iv3 <- iv_robust(f_rw ~ f_prot + f_young + f_jew + f_fem +
f_ortsgeb + f_pruss + hhsize + lnpop +
gpop + f_blind + f_deaf + f_dumb + f_miss |
kmwittenberg + kmwittenberg:lnpop + kmwittenberg:gpop +
f_young + f_jew + f_fem + f_ortsgeb +
f_pruss + hhsize + lnpop + gpop +
f_blind + f_deaf + f_dumb + f_miss,
data = mydata, se_type = "stata")
summary(iv3)
##
## Call:
## iv_robust(formula = f_rw ~ f_prot + f_young + f_jew + f_fem +
## f_ortsgeb + f_pruss + hhsize + lnpop + gpop + f_blind + f_deaf +
## f_dumb + f_miss | kmwittenberg + kmwittenberg:lnpop + kmwittenberg:gpop +
## f_young + f_jew + f_fem + f_ortsgeb + f_pruss + hhsize +
## lnpop + gpop + f_blind + f_deaf + f_dumb + f_miss, data = mydata,
##
## Standard error type: HC1
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
## (Intercept) 179.1507 29.13265 6.149 1.753e-09 121.8935 236.40788 438
## f_prot 0.1853 0.02572 7.205 2.553e-12 0.1348 0.23591 438
## f_young -1.9518 0.16695 -11.691 1.145e-27 -2.2799 -1.62366 438
## f_jew -0.4553 0.34569 -1.317 1.885e-01 -1.1348 0.22406 438
## f_fem -1.0800 0.33013 -3.271 1.155e-03 -1.7288 -0.43113 438
## f_ortsgeb 0.6023 0.04980 12.093 2.996e-29 0.5044 0.70016 438
## f_pruss -0.1858 0.14074 -1.320 1.876e-01 -0.4624 0.09085 438
## hhsize 0.7897 1.60849 0.491 6.237e-01 -2.3716 3.95102 438
## lnpop -1.3133 0.95155 -1.380 1.682e-01 -3.1835 0.55683 438
## gpop 0.4020 0.10650 3.775 1.823e-04 0.1927 0.61133 438
## f_blind -28.4769 14.69031 -1.938 5.321e-02 -57.3492 0.39532 438
## f_deaf -52.1996 10.37314 -5.032 7.092e-07 -72.5869 -31.81228 438
## f_dumb 7.5271 1.68652 4.463 1.029e-05 4.2124 10.84175 438
## f_miss -0.4979 0.37693 -1.321 1.872e-01 -1.2388 0.24288 438
##
## Multiple R-squared: 0.692 , Adjusted R-squared: 0.6828
## F-statistic: 66.54 on 13 and 438 DF, p-value: < 2.2e-16実証例 7.4 （操作変数の強さ）
(pp.298-299)
library(estimatr)
fs <- lm_robust(f_prot ~ kmwittenberg + f_young + f_jew + f_fem +
f_ortsgeb + f_pruss + hhsize + lnpop + gpop +
f_blind + f_deaf + f_dumb + f_miss,
data = mydata, se_type = "stata") # first-stage
beta1_hat <- round(fs$coefficients[2], 3)
se <- round(fs$std.error[2], 3)
t_value <- beta1_hat / se
f_value <- t_value ^ 2
f_value
## kmwittenberg
## 74.58678fs2 <- lm(f_prot ~ kmwittenberg + (kmwittenberg : lnpop) + (kmwittenberg : gpop) + f_young + f_jew + f_fem + f_ortsgeb + f_pruss + hhsize + lnpop + gpop + f_blind + f_deaf + f_dumb + f_miss, data = mydata) # first-stage with the interactionlibrary(car)
## Loading required package: carData
linearHypothesis(fs2, c("kmwittenberg = 0", "kmwittenberg:lnpop = 0", "kmwittenberg:gpop = 0"), white.adjust = "hc1")
## Linear hypothesis test
##
## Hypothesis:
## kmwittenberg = 0
## kmwittenberg:lnpop = 0
## kmwittenberg:gpop = 0
##
## Model 1: restricted model
## Model 2: f_prot ~ kmwittenberg + (kmwittenberg:lnpop) + (kmwittenberg:gpop) +
## f_young + f_jew + f_fem + f_ortsgeb + f_pruss + hhsize +
## lnpop + gpop + f_blind + f_deaf + f_dumb + f_miss
##
## Note: Coefficient covariance matrix supplied.
##
## Res.Df Df F Pr(>F)
## 1 439
## 2 436 3 29.794 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1実証例 7.5 （外生性の検定）
(pp.302-303)
library(estimatr)
library(car)
iv4 <- iv_robust(f_rw ~ f_prot + f_young + f_jew + f_fem + f_ortsgeb + f_pruss + hhsize + lnpop + gpop + f_blind + f_deaf + f_dumb + f_miss | kmwittenberg + f_young + f_jew + f_fem + f_ortsgeb + f_pruss + hhsize + lnpop + gpop + f_blind + f_deaf + f_dumb + f_miss, data = mydata)
mydata$uhat1 <- mydata$f_rw - iv4$fitted.values
res_reg1 <- lm_robust(uhat1 ~ kmwittenberg + f_young + f_jew + f_fem + f_ortsgeb + f_pruss + hhsize + lnpop + gpop + f_blind + f_deaf + f_dumb + f_miss, data = mydata, se_type = "classical")
summary(res_reg1)
##
## Call:
## lm_robust(formula = uhat1 ~ kmwittenberg + f_young + f_jew +
## f_fem + f_ortsgeb + f_pruss + hhsize + lnpop + gpop + f_blind +
## f_deaf + f_dumb + f_miss, data = mydata, se_type = "classical")
##
## Standard error type: classical
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
## (Intercept) 3.287e-13 29.346274 1.120e-14 1 -57.677016 57.677016 438
## kmwittenberg 2.011e-16 0.002697 7.455e-14 1 -0.005301 0.005301 438
## f_young 1.559e-15 0.172282 9.051e-15 1 -0.338603 0.338603 438
## f_jew -6.197e-16 0.304663 -2.034e-15 1 -0.598783 0.598783 438
## f_fem 1.571e-15 0.329405 4.768e-15 1 -0.647412 0.647412 438
## f_ortsgeb -2.437e-16 0.032940 -7.400e-15 1 -0.064739 0.064739 438
## f_pruss -4.301e-15 0.196687 -2.187e-14 1 -0.386567 0.386567 438
## hhsize -1.175e-14 1.417367 -8.290e-15 1 -2.785685 2.785685 438
## lnpop 6.110e-15 0.952518 6.414e-15 1 -1.872075 1.872075 438
## gpop -1.145e-15 0.099179 -1.155e-14 1 -0.194925 0.194925 438
## f_blind -1.159e-15 12.516466 -9.260e-17 1 -24.599798 24.599798 438
## f_deaf 5.180e-14 8.271965 6.262e-15 1 -16.257678 16.257678 438
## f_dumb -2.739e-14 2.068650 -1.324e-14 1 -4.065715 4.065715 438
## f_miss 1.193e-15 0.347936 3.429e-15 1 -0.683831 0.683831 438
##
## Multiple R-squared: 9.992e-16 , Adjusted R-squared: -0.02968
## F-statistic: 3.367e-14 on 13 and 438 DF, p-value: 1iv5 <- iv_robust(f_rw ~ f_prot + f_young + f_jew + f_fem +
f_ortsgeb + f_pruss + hhsize + lnpop +
gpop + f_blind + f_deaf + f_dumb + f_miss |
kmwittenberg + (kmwittenberg : lnpop) + (kmwittenberg : gpop) +
f_young + f_jew + f_fem + f_ortsgeb +
f_pruss + hhsize + lnpop + gpop +
f_blind + f_deaf + f_dumb + f_miss,
data = mydata)
mydata$uhat2 <- mydata$f_rw - iv5$fitted.values
res_reg2 <- lm_robust(uhat2 ~ kmwittenberg + (kmwittenberg : lnpop) + (kmwittenberg : gpop) + f_young + f_jew + f_fem + f_ortsgeb + f_pruss + hhsize + lnpop + gpop + f_blind + f_deaf + f_dumb + f_miss, data = mydata, se_type = "classical")
summary(res_reg2)##
## Call:
## lm_robust(formula = uhat2 ~ kmwittenberg + (kmwittenberg:lnpop) +
## (kmwittenberg:gpop) + f_young + f_jew + f_fem + f_ortsgeb +
## f_pruss + hhsize + lnpop + gpop + f_blind + f_deaf + f_dumb +
## f_miss, data = mydata, se_type = "classical")
##
## Standard error type: classical
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) CI Lower
## (Intercept) -1.193e+01 3.738e+01 -0.319277 0.7497 -85.404372
## kmwittenberg 3.506e-02 7.091e-02 0.494444 0.6212 -0.104309
## f_young 6.741e-04 1.778e-01 0.003792 0.9970 -0.348694
## f_jew -5.388e-03 3.044e-01 -0.017703 0.9859 -0.603592
## f_fem 1.195e-02 3.287e-01 0.036363 0.9710 -0.634000
## f_ortsgeb -5.907e-04 3.285e-02 -0.017983 0.9857 -0.065151
## f_pruss -1.984e-03 1.960e-01 -0.010122 0.9919 -0.387260
## hhsize 8.812e-02 1.419e+00 0.062080 0.9505 -2.701716
## lnpop 1.034e+00 2.285e+00 0.452376 0.6512 -3.457959
## gpop -5.554e-02 2.489e-01 -0.223140 0.8235 -0.544739
## f_blind 8.972e-01 1.253e+01 0.071632 0.9429 -23.720658
## f_deaf 4.100e-01 8.318e+00 0.049292 0.9607 -15.938414
## f_dumb -2.080e-02 2.062e+00 -0.010087 0.9920 -4.073170
## f_miss -1.432e-02 3.480e-01 -0.041155 0.9672 -0.698345
## kmwittenberg:lnpop -3.290e-03 6.587e-03 -0.499485 0.6177 -0.016236
## kmwittenberg:gpop 1.663e-04 6.549e-04 0.253890 0.7997 -0.001121
## CI Upper DF
## (Intercept) 61.534577 436
## kmwittenberg 0.174432 436
## f_young 0.350043 436
## f_jew 0.592816 436
## f_fem 0.657902 436
## f_ortsgeb 0.063970 436
## f_pruss 0.383291 436
## hhsize 2.877955 436
## lnpop 5.525712 436
## gpop 0.433659 436
## f_blind 25.515112 436
## f_deaf 16.758436 436
## f_dumb 4.031575 436
## f_miss 0.669699 436
## kmwittenberg:lnpop 0.009656 436
## kmwittenberg:gpop 0.001453 436
##
## Multiple R-squared: 0.0006048 , Adjusted R-squared: -0.03378
## F-statistic: 0.01759 on 15 and 436 DF, p-value: 1linearHypothesis(res_reg2, c("kmwittenberg", "kmwittenberg:lnpop = 0", "kmwittenberg:gpop = 0"))## Linear hypothesis test
##
## Hypothesis:
## kmwittenberg = 0
## kmwittenberg:lnpop = 0
## kmwittenberg:gpop = 0
##
## Model 1: restricted model
## Model 2: uhat2 ~ kmwittenberg + (kmwittenberg:lnpop) + (kmwittenberg:gpop) +
## f_young + f_jew + f_fem + f_ortsgeb + f_pruss + hhsize +
## lnpop + gpop + f_blind + f_deaf + f_dumb + f_miss
##
## Res.Df Df Chisq Pr(>Chisq)
## 1 439
## 2 436 3 0.2638 0.9667