两阶段最小二乘法与R

文章目录

两阶段最小二乘法与R
- @[toc]
- 1、ivreg包介绍
- 2 、操作流程

1、ivreg包介绍

R语言计量包ivreg用以解决线性回归模型的内生性问题。

描述：工具变量估计的线性模型通过两阶段最小二乘(2SLS) 回归或通过稳健回归M估计(2SM)或MM估计(2SMM)。主要的ivreg()模型拟合函数旨在提供一个工作流程，尽可能类似于标准的lm()回归。大量的方法是被用来拟合ivreg模型对象，除了其他标准模型工具，还包括广泛的功能，计算和图形回归诊断。

作者：Author John Fox, Christian Kleiber, Achim Zeileis

2 、操作流程

在利用工具变量法估计线性回归模型时，往往选择Stata操作，现在介绍R的操作方法：首先，我们先安装工具变量回归安装包ivreg，并加载相关其他计量包；

setwd("D:/Allcode/Rstudy/model/IV_estimate")  # 先设置路径
install.packages("ivreg")    # 安装ivreg
install.packages("haven")    # 用于stata数据导入，默认存在，可以不安装
install.packages("lmtest")   # 用于线性回归检验
install.packages("sandwich") # 提供相关异方差稳健标准误
#加载以上所有包
library("haven")
library("ivreg")
library("lmtest")
library("sandwich")

接下来准备数据集，我选用的是陈强老师主页（陈强教授的计量经济学及Stata主页 (econometrics-stata.com)）上的数据集grilic，它是stata的dta格式，因此需要转换导入

grilic <- read_dta("grilic.dta")
names(grilic)  # 查看数据框的变量名# [1] "rns"      "rns80"    "mrt"      "mrt80"    "smsa"     "smsa80"
# [7] "med"      "iq"       "kww"      "year"     "age"      "age80"
# [13] "s"        "s80"      "expr"     "expr80"   "tenure"   "tenure80"
# [19] "lw"       "lw80"

该数据集中包括以下变量:lw(工资对数),s (受教育年限) , age(年龄) , expr(工龄) , tenure(在现单位的工作年数),q(智商), med(母亲的受教育年限),kww(在"knowledge of the World ofWork"测试中的成绩),mt(婚姻虚拟变量,已婚=1),rns (美国南方虚拟变量,住在南方=1),smsa(大城市虚拟变量,住在大城市=1),year (有数据的最早年份, 1966-1973年中的某一年)。我们选择lw为被解释变量，其余变量为解释变量。先利用OLS回归作为基准模型

fit_ols1 <- lm(lw ~ s+expr+tenure+rns+smsa,data = grilic)  # 没有加入智商iq变量# Call:
#   lm(formula = lw ~ s + expr + tenure + rns + smsa, data = grilic)
#
# Residuals:
#   Min       1Q   Median       3Q      Max
# -1.11684 -0.22626 -0.01511  0.23103  1.23738
#
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)
# (Intercept)  4.103675   0.085097  48.223  < 2e-16 ***
#   s            0.102643   0.005849  17.549  < 2e-16 ***
#   expr         0.038119   0.006327   6.025 2.65e-09 ***
#   tenure       0.035615   0.007742   4.600 4.96e-06 ***
#   rns         -0.084080   0.028797  -2.920  0.00361 **
#   smsa         0.139667   0.028082   4.974 8.15e-07 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 0.3464 on 752 degrees of freedom
# Multiple R-squared:  0.3521,  Adjusted R-squared:  0.3478
# F-statistic: 81.75 on 5 and 752 DF,  p-value: < 2.2e-16
#由于没有加入iq，存在遗漏变量问题，因此加入iqfit_ols2 <- lm(lw ~ iq+s+expr+tenure+rns+smsa,data = grilic)# Call:
#   lm(formula = lw ~ iq + s + expr + tenure + rns + smsa, data = grilic)
#
# Residuals:
#   Min       1Q   Median       3Q      Max
# -1.16056 -0.21786 -0.00622  0.22771  1.20580
#
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)
# (Intercept)  3.895172   0.109110  35.699  < 2e-16 ***
#   iq           0.003279   0.001083   3.028  0.00255 **
#   s            0.092787   0.006666  13.920  < 2e-16 ***
#   expr         0.039344   0.006306   6.239 7.33e-10 ***
#   tenure       0.034209   0.007715   4.434 1.06e-05 ***
#   rns         -0.074532   0.028815  -2.587  0.00988 **
#   smsa         0.136737   0.027948   4.893 1.22e-06 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 0.3445 on 751 degrees of freedom
# Multiple R-squared:   0.36,   Adjusted R-squared:  0.3548
# F-statistic: 70.39 on 6 and 751 DF,  p-value: < 2.2e-16
# 以上回归都是基于同方差假设条件下的结果，我们将系数转换为异方差稳健标准误；这里以fit_ols2为例coeftest(fit_ols2, vcov = vcovHC, type = "HC1") # 异方差稳健标准误# t test of coefficients:
#
#   Estimate Std. Error t value  Pr(>|t|)
# (Intercept)  3.8951718  0.1159286 33.5997 < 2.2e-16 ***
#   iq           0.0032792  0.0011321  2.8965  0.003883 **
#   s            0.0927874  0.0069763 13.3004 < 2.2e-16 ***
#   expr         0.0393443  0.0066603  5.9072 5.272e-09 ***
#   tenure       0.0342090  0.0078957  4.3326 1.674e-05 ***
#   rns         -0.0745325  0.0299772 -2.4863  0.013124 *
#   smsa         0.1367369  0.0277712  4.9237 1.045e-06 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

使用工具变量法回归，内生解释变量为iq，工具变量选择med、kww、mrt、age;其余控制变量自身视为自身的工具变量；代码如下：

fit_iv <- ivreg(lw ~ iq+s+expr+tenure+rns+smsa | s+ med+ kww+mrt+age+expr+tenure+rns+smsa ,data = grilic)# 这里用"|"分隔内生解释变量与工具变量
# 提取稳健标准误
coeftest(fit_iv, vcov = vcovHC, type = "HC0")  # 异方差稳健标准误# t test of coefficients:
#
#   Estimate Std. Error t value  Pr(>|t|)
# (Intercept)  4.8378747  0.3799432 12.7332 < 2.2e-16 ***
#   iq          -0.0115468  0.0056376 -2.0482  0.040887 *
#   s            0.1373477  0.0174989  7.8489 1.446e-14 ***
#   expr         0.0338041  0.0074844  4.5166 7.295e-06 ***
#   tenure       0.0405640  0.0095848  4.2321 2.602e-05 ***
#   rns         -0.1176984  0.0359582 -3.2732  0.001112 **
#   smsa         0.1499830  0.0322276  4.6539 3.850e-06 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1# 工具变量法回归还要进行诊断
summary(fit_iv,test = TRUE) # 诊断# Call:
#   ivreg(formula = lw ~ iq + s + expr + tenure + rns + smsa | s +
#           med + kww + mrt + age + expr + tenure + rns + smsa, data = grilic)
#
# Residuals:
#   Min         1Q     Median         3Q        Max
# -1.3825405 -0.2437078  0.0009735  0.2514625  1.4609417
#
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)
# (Intercept)  4.837875   0.346424  13.965  < 2e-16 ***
#   iq          -0.011547   0.005241  -2.203 0.027889 *
#   s            0.137348   0.017042   8.059 3.02e-15 ***
#   expr         0.033804   0.007302   4.630 4.32e-06 ***
#   tenure       0.040564   0.008896   4.560 5.98e-06 ***
#   rns         -0.117698   0.035468  -3.318 0.000949 ***
#   smsa         0.149983   0.031572   4.751 2.43e-06 ***
#
#   Diagnostic tests:
#   df1 df2 statistic  p-value
# Weak instruments   4 748     10.54 2.61e-08 ***    # 弱工具变量检验（通过）
#   Wu-Hausman         1 750     10.70  0.00112 **   # 内生性检验（通过）
#   Sargan             3  NA     61.14 3.36e-13 ***  # 过度识别检验（未通过）
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 0.3851 on 751 degrees of freedom
# Multiple R-Squared: 0.2002,   Adjusted R-squared: 0.1938
# Wald test: 55.92 on 6 and 751 DF,  p-value: < 2.2e-16

由于工具变量个数大于内生解释变量个数，且工具变量过度识别检验未通过，因此需要调整工具变量；这里怀疑age与tenure可能存在过度识别，剔除后进行ivreg回归

fit_iv2 <- ivreg(lw ~ iq+s+expr+tenure+rns+smsa | med+ kww + s + expr + tenure + rns + smsa ,data = grilic)
summary(fit_iv2,test = TRUE)
# Call:
#   ivreg(formula = lw ~ iq + s + expr + tenure + rns + smsa | med +
#           kww + s + expr + tenure + rns + smsa, data = grilic)
#
# Residuals:
#   Min         1Q     Median         3Q        Max
# -1.3025533 -0.2405658  0.0005969  0.2349962  1.2621665
#
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)
# (Intercept)  3.218043   0.384814   8.363 2.97e-16 ***
#   iq           0.013928   0.005884   2.367 0.018186 *
#   s            0.060780   0.018735   3.244 0.001230 **
#   expr         0.043324   0.007038   6.156 1.22e-09 ***
#   tenure       0.029644   0.008561   3.463 0.000565 ***
#   rns         -0.043527   0.034922  -1.246 0.213000
# smsa         0.127222   0.030137   4.221 2.72e-05 ***
#
#   Diagnostic tests:
#   df1 df2 statistic  p-value
# Weak instruments   2 750    14.906 4.49e-07 *** （通过）
#   Wu-Hausman         1 750     3.858   0.0499 *  （通过）
#   Sargan             1  NA     0.130   0.7185    （通过）
# ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 0.3661 on 751 degrees of freedom
# Multiple R-Squared: 0.2775,   Adjusted R-squared: 0.2718
# Wald test: 61.94 on 6 and 751 DF,  p-value: < 2.2e-16

-END-

参考文献

陈强.高级计量经济学[M].高等教育出版社
https://cran.r-project.org/web/packages/gmm/gmm.pdf
http://www.econometrics-stata.com/col.jsp?id=101