R語言代寫工具變量與兩階段最小二乘法

原文鏈接：http://tecdat.cn/?p=5374

 

我們要估計的模型是

y=a+bx+cd+ey=a+bx+cd+e，

其中是解釋變量，，和是我們想要估計的系數。是控制變量，是治療變量。我們特別關注我們的治療效果對。 

生成數據

首先，讓我們生成數據。

假設 的工具變量和之間的相關矩陣如下： 

                 <span style="color:#009999">0.001</span>,<span style="color:#009999">1</span>,<span style="color:#009999">0.7</span>,<span style="color:#009999">0.3</span>,
 rownames(R)<-colnames(R)<-c(<span style="color:#dd1144">"x"</span>,<span style="color:#dd1144">"d"</span>,<span style="color:#dd1144">"z"</span>,<span style="color:#dd1144">"e"</span>)
R</code></span></span>

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

具體而言，相關性表明

1. cor（d，e）= 0.3，這意味着是內生的;dd
2. cor（d，z）= 0.7，這意味着是的強大工具變量;zzdd
3. cor（z，e）= 0.001，這意味着工具變量滿足排除限制，因為它只影響到。zzyydd

現在，讓我們使用指定的相關性為，，和生成數據。xxddzzee

 nvars = dim(U)[<span style="color:#009999">1</span>]
numobs = <span style="color:#009999">1000</span>
 random.normal = matrix(rnorm(nvars*numobs,<span style="color:#009999">0</span>,<span style="color:#009999">1</span>), nrow=nvars, ncol=numobs);
X = U %*% random.normal
newX = t(X)
data = as.data.frame(newX)
<span style="color:#990000"><strong>attach</strong></span>(data)</code></span></span>

數據看起來像這樣：

<span style="color:#333333"><span style="color:#333333"><code>head(data)</code></span></span>

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

以及數據之間的相關性

<span style="color:#333333"><span style="color:#333333"><code>cor(data)</code></span></span>

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

<span style="color:#333333"><span style="color:#333333"><code>y<-<span style="color:#009999">10</span>+<span style="color:#009999">1</span>*x+<span style="color:#009999">1</span>*d+e</code></span></span>

如果我們假裝我們不知道真正的關系並使用和來解釋，我們對和正確系數應該接近到。  

OLS

如果我們只使用OLS來估計系數：

<span style="color:#333333"><span style="color:#333333"><code>ols<-lm(formula = y~x+d)
summary(ols)</code></span></span>

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

<span style="color:#333333"><span style="color:#333333"><code>tsls1<-lm(d~x+z)
summary(tsls1)</code></span></span>

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

<span style="color:#333333"><span style="color:#333333"><code>d.hat<-fitted.values(tsls1)</code></span></span>

第2階段：在和上回歸yyxxd.hatd.hat

<span style="color:#333333"><span style="color:#333333"><code>tsls2<-lm(y~x+d.hat)
summary(tsls2)</code></span></span>

## 
## 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。

 

 

如果您有任何疑問，請在下面發表評論。