我们要估计的模型是
y=a+bx+cd+ey=a+bx+cd+e，
其中是解释变量，，和是我们想要估计的系数。是控制变量，是治疗变量。我们特别关注我们的治疗效果对。
生成数据
首先，让我们生成数据。
假设 的工具变量和之间的相关矩阵如下：
0.001,1,0.7,0.3,\n rownames(R)<-colnames(R)<-c(\"x\",\"d\",\"z\",\"e\")\nR","classes":{"has":1}}" data-cke-widget-upcasted="1" data-cke-widget-keep-attr="0" data-widget="codeSnippet" style="font-family:" courier="" font-size:="" margin:="" padding:=""> 0.001,1,0.7,0.3, rownames(R)<-colnames(R)<-c("x","d","z","e") R
## x d z e ## x 1.000 0.001 0.002 0.001 ## d 0.001 1.000 0.700 0.300 ## z 0.002 0.700 1.000 0.001 ## e 0.001 0.300 0.001 1.000
具体而言，相关性表明
cor（d，e）= 0.3，这意味着是内生的;dd
cor（d，z）= 0.7，这意味着是的强大工具变量;zzdd
cor（z，e）= 0.001，这意味着工具变量满足排除限制，因为它只影响到。zzyydd
现在，让我们使用指定的相关性为，，和生成数据。xxddzzee
1]\nnumobs = 1000\n random.normal = matrix(rnorm(nvars*numobs,0,1), nrow=nvars, ncol=numobs);\nX = U %*% random.normal\nnewX = t(X)\ndata = as.data.frame(newX)\nattach(data)","classes":{"has":1}}" data-cke-widget-upcasted="1" data-cke-widget-keep-attr="0" data-widget="codeSnippet" style="font-family:" courier="" font-size:="" margin:="" padding:=""> nvars = dim(U)[1] numobs = 1000 random.normal = matrix(rnorm(nvars*numobs,0,1), nrow=nvars, ncol=numobs); X = U %*% random.normal newX = t(X) data = as.data.frame(newX) attach(data)
数据看起来像这样：
head(data)","classes":{"has":1}}" data-cke-widget-upcasted="1" data-cke-widget-keep-attr="0" data-widget="codeSnippet" style="font-family:" courier="" font-size:="" margin:="" padding:="">head(data)
## x d z e ## 1 -0.62645381 0.1830168 -0.4694601 1.7474361 ## 2 0.32950777 -0.8201385 -0.2255741 0.2818908 ## 3 0.57578135 -0.3048125 0.8670061 -0.1795257 ## 4 -0.62124058 -2.2153200 -0.7481687 -1.0350488 ## 5 -0.01619026 0.9438195 1.2471197 0.5820200 ## 6 0.91897737 0.7830549 0.6025820 -1.5924689
以及数据之间的相关性
cor(data)","classes":{"has":1}}" data-cke-widget-upcasted="1" data-cke-widget-keep-attr="0" data-widget="codeSnippet" style="font-family:" courier="" font-size:="" margin:="" padding:="">cor(data)
## x d z e ## x 1.00000000 0.00668391 -0.012319595 0.016239235 ## d 0.00668391 1.00000000 0.680741763 0.312192680 ## z -0.01231960 0.68074176 1.000000000 0.006322354 ## e 0.01623923 0.31219268 0.006322354 1.000000000
正如我们之前指定的那样。
现在让我们指定真正的数据生成过程并生成解释变量yy
y<-10+1*x+1*d+e","classes":{"has":1}}" data-cke-widget-upcasted="1" data-cke-widget-keep-attr="0" data-widget="codeSnippet" style="font-family:" courier="" font-size:="" margin:="" padding:="">y<-10+1*x+1*d+e
如果我们假装我们不知道真正的关系并使用和来解释，我们对和正确系数应该接近到。
OLS
如果我们只使用OLS来估计系数：
ols<-lm(formula = y~x+d)\nsummary(ols)","classes":{"has":1}}" data-cke-widget-upcasted="1" data-cke-widget-keep-attr="0" data-widget="codeSnippet" style="font-family:" courier="" font-size:="" margin:="" padding:="">ols<-lm(formula = y~x+d) summary(ols)
## ## Call: ## lm(formula = y ~ x + d) ## ## Residuals: ## Min 1Q Median 3Q Max ## -3.2395 -0.5952 -0.0308 0.6617 2.7592 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 9.99495 0.03105 321.89 <2e-16 *** ## x 1.01408 0.02992 33.89 <2e-16 *** ## d 1.31356 0.03023 43.46 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.9817 on 997 degrees of freedom ## Multiple R-squared: 0.7541, Adjusted R-squared: 0.7536 ## F-statistic: 1528 on 2 and 997 DF, p-value: < 2.2e-16
b的估计系数是1.31 instread of 1. ## 2SLS ##现在我们使用2SLS来估计这种关系。我们使用z作为d的工具变量
第1阶段：在和上回归，并将d的拟合值保存为d。ddxxzz
tsls1<-lm(d~x+z)\nsummary(tsls1)","classes":{"has":1}}" data-cke-widget-upcasted="1" data-cke-widget-keep-attr="0" data-widget="codeSnippet" style="font-family:" courier="" font-size:="" margin:="" padding:="">tsls1<-lm(d~x+z) summary(tsls1)
## ## Call: ## lm(formula = d ~ x + z) ## ## Residuals: ## Min 1Q Median 3Q Max ## -2.59344 -0.52572 0.04978 0.53115 2.01555 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -0.01048 0.02383 -0.44 0.660 ## x 0.01492 0.02296 0.65 0.516 ## z 0.68594 0.02337 29.36 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.7534 on 997 degrees of freedom ## Multiple R-squared: 0.4636, Adjusted R-squared: 0.4626 ## F-statistic: 430.9 on 2 and 997 DF, p-value: < 2.2e-16
d.hat<-fitted.values(tsls1)","classes":{"has":1}}" data-cke-widget-upcasted="1" data-cke-widget-keep-attr="0" data-widget="codeSnippet" style="font-family:" courier="" font-size:="" margin:="" padding:="">d.hat<-fitted.values(tsls1)
第2阶段：在和上回归yyxxd.hatd.hat
tsls2<-lm(y~x+d.hat)\nsummary(tsls2)","classes":{"has":1}}" data-cke-widget-upcasted="1" data-cke-widget-keep-attr="0" data-widget="codeSnippet" style="font-family:" courier="" font-size:="" margin:="" padding:="">tsls2<-lm(y~x+d.hat) summary(tsls2)
## ## Call: ## lm(formula = y ~ x + d.hat) ## ## Residuals: ## Min 1Q Median 3Q Max ## -4.4531 -1.0333 0.0228 1.0657 4.0104 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 9.99507 0.04786 208.85 <2e-16 *** ## x 1.01609 0.04612 22.03 <2e-16 *** ## d.hat 1.00963 0.06842 14.76 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 1.513 on 997 degrees of freedom ## Multiple R-squared: 0.4158, Adjusted R-squared: 0.4146 ## F-statistic: 354.8 on 2 and 997 DF, p-value: < 2.2e-16
结果
b的真值：1 OLS estiamte of b：.00963 2SLS estiamte of b：1.31356
如果治疗变量是内生的，我们使用2SLS。
网友评论