# install.packages("leaps") 
# install.packages("ivreg")
library("leaps")
library("ivreg")

data = read.csv("http://hedibert.org/wp-content/uploads/2024/11/card.csv",header=TRUE) 

attach(data) 

n = nrow(data) 

1 OLS regression with several exogenous variables

outcome   = lwage
treatment = educ
W         = cbind(exper,expersq,black,smsa,south,smsa66,
                  reg662,reg663,reg664,reg665,reg666,reg667,reg668,reg669)
X         = cbind(treatment,W)
reg.ols   = lm(outcome~X) 
summary(reg.ols)
## 
## Call:
## lm(formula = outcome ~ X)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.62326 -0.22141  0.02001  0.23932  1.33340 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.6208069  0.0742327  62.248  < 2e-16 ***
## Xtreatment   0.0746933  0.0034983  21.351  < 2e-16 ***
## Xexper       0.0848320  0.0066242  12.806  < 2e-16 ***
## Xexpersq    -0.0022870  0.0003166  -7.223 6.41e-13 ***
## Xblack      -0.1990123  0.0182483 -10.906  < 2e-16 ***
## Xsmsa        0.1363846  0.0201005   6.785 1.39e-11 ***
## Xsouth      -0.1479550  0.0259799  -5.695 1.35e-08 ***
## Xsmsa66      0.0262417  0.0194477   1.349  0.17733    
## Xreg662      0.0963671  0.0358979   2.684  0.00730 ** 
## Xreg663      0.1445400  0.0351244   4.115 3.97e-05 ***
## Xreg664      0.0550756  0.0416573   1.322  0.18623    
## Xreg665      0.1280248  0.0418395   3.060  0.00223 ** 
## Xreg666      0.1405174  0.0452469   3.106  0.00192 ** 
## Xreg667      0.1179810  0.0448025   2.633  0.00850 ** 
## Xreg668     -0.0564361  0.0512579  -1.101  0.27098    
## Xreg669      0.1185697  0.0388301   3.054  0.00228 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3723 on 2994 degrees of freedom
## Multiple R-squared:  0.2998, Adjusted R-squared:  0.2963 
## F-statistic: 85.48 on 15 and 2994 DF,  p-value: < 2.2e-16

2 IV regression with several exogenous variables

iv     = nearc4
Z      = cbind(iv,W)
reg.iv = ivreg(outcome ~ X | Z)
summary(reg.iv)
## 
## Call:
## ivreg(formula = outcome ~ X | Z)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.83164 -0.24075  0.02428  0.25208  1.42760 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.6661519  0.9248295   3.964 7.54e-05 ***
## Xtreatment   0.1315038  0.0549637   2.393  0.01679 *  
## Xexper       0.1082711  0.0236586   4.576 4.92e-06 ***
## Xexpersq    -0.0023349  0.0003335  -7.001 3.12e-12 ***
## Xblack      -0.1467758  0.0538999  -2.723  0.00650 ** 
## Xsmsa        0.1118084  0.0316620   3.531  0.00042 ***
## Xsouth      -0.1446715  0.0272846  -5.302 1.23e-07 ***
## Xsmsa66      0.0185311  0.0216086   0.858  0.39119    
## Xreg662      0.1007678  0.0376857   2.674  0.00754 ** 
## Xreg663      0.1482588  0.0368141   4.027 5.78e-05 ***
## Xreg664      0.0498971  0.0437398   1.141  0.25406    
## Xreg665      0.1462719  0.0470639   3.108  0.00190 ** 
## Xreg666      0.1629029  0.0519096   3.138  0.00172 ** 
## Xreg667      0.1345722  0.0494023   2.724  0.00649 ** 
## Xreg668     -0.0830770  0.0593313  -1.400  0.16155    
## Xreg669      0.1078142  0.0418137   2.578  0.00997 ** 
## 
## Diagnostic tests:
##                   df1  df2 statistic  p-value    
## Weak instruments    1 2994    13.256 0.000276 ***
## Wu-Hausman          1 2993     1.168 0.279973    
## Sargan              0   NA        NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3883 on 2994 degrees of freedom
## Multiple R-Squared: 0.2382,  Adjusted R-squared: 0.2343 
## Wald test: 51.01 on 15 and 2994 DF,  p-value: < 2.2e-16

3 OLS regression with selection of exogenous variables

X1   = X
data = data.frame(y=outcome,X=X1)
p    = ncol(X1)
regs = regsubsets(y~X,data=data,method="exhaustive",nvmax=p)
bic = summary(regs)$bic
ii  = 1:p
coef(regs,ii[bic==min(bic)])
##  (Intercept)   Xtreatment       Xexper     Xexpersq       Xblack        Xsmsa 
##  4.713907129  0.074619220  0.084431565 -0.002270828 -0.190925636  0.157572489 
##       Xsouth      Xreg663      Xreg668 
## -0.116091147  0.054617876 -0.143638697
ols = lm(outcome ~ treatment+exper+expersq+black+smsa+south+reg663+reg668)
summary(ols)
## 
## Call:
## lm(formula = outcome ~ treatment + exper + expersq + black + 
##     smsa + south + reg663 + reg668)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.58335 -0.22424  0.01746  0.24290  1.34410 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.7139071  0.0676384  69.693  < 2e-16 ***
## treatment    0.0746192  0.0034959  21.344  < 2e-16 ***
## exper        0.0844316  0.0066260  12.742  < 2e-16 ***
## expersq     -0.0022708  0.0003168  -7.169 9.46e-13 ***
## black       -0.1909256  0.0175808 -10.860  < 2e-16 ***
## smsa         0.1575725  0.0155371  10.142  < 2e-16 ***
## south       -0.1160911  0.0159501  -7.278 4.29e-13 ***
## reg663       0.0546179  0.0183605   2.975  0.00296 ** 
## reg668      -0.1436387  0.0418657  -3.431  0.00061 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3728 on 3001 degrees of freedom
## Multiple R-squared:  0.2961, Adjusted R-squared:  0.2942 
## F-statistic: 157.8 on 8 and 3001 DF,  p-value: < 2.2e-16

4 IV regression with selection of exogenous variables

data = data.frame(y=outcome,Z=Z)
p    = ncol(data)
regs = regsubsets(y~Z,data=data,method="exhaustive",nvmax=p,force.in=Z[,1])
bic = summary(regs)$bic
ii  = 1:p
coef(regs,ii[bic==min(bic)])
## (Intercept)         Ziv      Zexper    Zexpersq      Zblack       Zsmsa 
##  5.93718986  0.04594223  0.05358652 -0.00220199 -0.26286597  0.18248462 
##      Zsouth     Zreg663 
## -0.12741608  0.06092994
reg.iv = ivreg(outcome ~ treatment+exper+expersq+black+smsa+south+reg663 | 
                         iv+exper+expersq+black+smsa+south+reg663)
summary(reg.iv)
## 
## Call:
## ivreg(formula = outcome ~ treatment + exper + expersq + black + 
##     smsa + south + reg663 | iv + exper + expersq + black + smsa + 
##     south + reg663)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.83851 -0.23743  0.02038  0.25642  1.46049 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.6610748  0.8361146   4.379 1.23e-05 ***
## treatment    0.1365294  0.0495914   2.753 0.005939 ** 
## exper        0.1095925  0.0214619   5.106 3.49e-07 ***
## expersq     -0.0023032  0.0003357  -6.860 8.33e-12 ***
## black       -0.1254105  0.0532781  -2.354 0.018643 *  
## smsa         0.1271510  0.0303902   4.184 2.95e-05 ***
## south       -0.0862605  0.0241822  -3.567 0.000367 ***
## reg663       0.0661039  0.0193090   3.423 0.000627 ***
## 
## Diagnostic tests:
##                   df1  df2 statistic  p-value    
## Weak instruments    1 3002    16.621 4.68e-05 ***
## Wu-Hausman          1 3001     1.762    0.185    
## Sargan              0   NA        NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3928 on 3002 degrees of freedom
## Multiple R-Squared: 0.2184,  Adjusted R-squared: 0.2166 
## Wald test: 104.1 on 7 and 3002 DF,  p-value: < 2.2e-16
LS0tCnRpdGxlOiAiR2F1c3NpYW4gbGluZWFyIHJlZ3Jlc3Npb24gd2l0aCBJViIKc3VidGl0bGU6ICJvdXRjb21lPWxvZ3dhZ2UgLSB0cmVhdG1lbnQ9ZWR1YyAtIGl2PW5lYXJjNCIKYXV0aG9yOiAiSGVkaWJlcnQgRnJlaXRhcyBMb3BlcyIKZGF0ZTogIjE1LzExLzIwMjUiCm91dHB1dDoKICBodG1sX2RvY3VtZW50OgogICAgdG9jOiB0cnVlCiAgICB0b2NfZGVwdGg6IDMKICAgIHRvY19jb2xsYXBzZWQ6IHRydWUKICAgIGNvZGVfZG93bmxvYWQ6IHRydWUKICAgIG51bWJlcl9zZWN0aW9uczogdHJ1ZQogIHBkZl9kb2N1bWVudDoKICAgIHRvYzogdHJ1ZQogICAgdG9jX2RlcHRoOiAnMycKLS0tCgpgYGB7ciBzZXR1cCwgaW5jbHVkZT1GQUxTRX0Ka25pdHI6Om9wdHNfY2h1bmskc2V0KGVjaG8gPSBUUlVFKQpgYGAKCgpgYGB7ciBmaWcuYWxpZ249J2NlbnRlcicsIGZpZy53aWR0aD04LCBmaWcuaGVpZ2h0PTZ9CiMgaW5zdGFsbC5wYWNrYWdlcygibGVhcHMiKSAKIyBpbnN0YWxsLnBhY2thZ2VzKCJpdnJlZyIpCmxpYnJhcnkoImxlYXBzIikKbGlicmFyeSgiaXZyZWciKQoKZGF0YSA9IHJlYWQuY3N2KCJodHRwOi8vaGVkaWJlcnQub3JnL3dwLWNvbnRlbnQvdXBsb2Fkcy8yMDI0LzExL2NhcmQuY3N2IixoZWFkZXI9VFJVRSkgCgphdHRhY2goZGF0YSkgCgpuID0gbnJvdyhkYXRhKSAKYGBgCgoKIyBPTFMgcmVncmVzc2lvbiB3aXRoIHNldmVyYWwgZXhvZ2Vub3VzIHZhcmlhYmxlcwoKYGBge3J9Cm91dGNvbWUgICA9IGx3YWdlCnRyZWF0bWVudCA9IGVkdWMKVyAgICAgICAgID0gY2JpbmQoZXhwZXIsZXhwZXJzcSxibGFjayxzbXNhLHNvdXRoLHNtc2E2NiwKICAgICAgICAgICAgICAgICAgcmVnNjYyLHJlZzY2MyxyZWc2NjQscmVnNjY1LHJlZzY2NixyZWc2NjcscmVnNjY4LHJlZzY2OSkKWCAgICAgICAgID0gY2JpbmQodHJlYXRtZW50LFcpCnJlZy5vbHMgICA9IGxtKG91dGNvbWV+WCkgCnN1bW1hcnkocmVnLm9scykKYGBgCgojIElWIHJlZ3Jlc3Npb24gd2l0aCBzZXZlcmFsIGV4b2dlbm91cyB2YXJpYWJsZXMKCmBgYHtyfQppdiAgICAgPSBuZWFyYzQKWiAgICAgID0gY2JpbmQoaXYsVykKcmVnLml2ID0gaXZyZWcob3V0Y29tZSB+IFggfCBaKQpzdW1tYXJ5KHJlZy5pdikKYGBgCgoKIyBPTFMgcmVncmVzc2lvbiB3aXRoIHNlbGVjdGlvbiBvZiBleG9nZW5vdXMgdmFyaWFibGVzCgpgYGB7cn0KWDEgICA9IFgKZGF0YSA9IGRhdGEuZnJhbWUoeT1vdXRjb21lLFg9WDEpCnAgICAgPSBuY29sKFgxKQpyZWdzID0gcmVnc3Vic2V0cyh5flgsZGF0YT1kYXRhLG1ldGhvZD0iZXhoYXVzdGl2ZSIsbnZtYXg9cCkKYmljID0gc3VtbWFyeShyZWdzKSRiaWMKaWkgID0gMTpwCmNvZWYocmVncyxpaVtiaWM9PW1pbihiaWMpXSkKb2xzID0gbG0ob3V0Y29tZSB+IHRyZWF0bWVudCtleHBlcitleHBlcnNxK2JsYWNrK3Ntc2Erc291dGgrcmVnNjYzK3JlZzY2OCkKc3VtbWFyeShvbHMpCmBgYAoKIyBJViByZWdyZXNzaW9uIHdpdGggc2VsZWN0aW9uIG9mIGV4b2dlbm91cyB2YXJpYWJsZXMKCmBgYHtyfQpkYXRhID0gZGF0YS5mcmFtZSh5PW91dGNvbWUsWj1aKQpwICAgID0gbmNvbChkYXRhKQpyZWdzID0gcmVnc3Vic2V0cyh5flosZGF0YT1kYXRhLG1ldGhvZD0iZXhoYXVzdGl2ZSIsbnZtYXg9cCxmb3JjZS5pbj1aWywxXSkKYmljID0gc3VtbWFyeShyZWdzKSRiaWMKaWkgID0gMTpwCmNvZWYocmVncyxpaVtiaWM9PW1pbihiaWMpXSkKCnJlZy5pdiA9IGl2cmVnKG91dGNvbWUgfiB0cmVhdG1lbnQrZXhwZXIrZXhwZXJzcStibGFjaytzbXNhK3NvdXRoK3JlZzY2MyB8IAogICAgICAgICAgICAgICAgICAgICAgICAgaXYrZXhwZXIrZXhwZXJzcStibGFjaytzbXNhK3NvdXRoK3JlZzY2MykKc3VtbWFyeShyZWcuaXYpCmBgYAo=