library(wooldridge)
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(sandwich)
data("meap00_01")
# (i)
ols_model <- lm(math4 ~ lunch + log(enroll) + log(exppp), data = meap00_01)
summary(ols_model)
##
## Call:
## lm(formula = math4 ~ lunch + log(enroll) + log(exppp), data = meap00_01)
##
## Residuals:
## Min 1Q Median 3Q Max
## -56.326 -8.758 0.914 9.164 51.416
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 91.93240 19.96170 4.605 4.42e-06 ***
## lunch -0.44874 0.01464 -30.648 < 2e-16 ***
## log(enroll) -5.39915 0.94041 -5.741 1.11e-08 ***
## log(exppp) 3.52475 2.09785 1.680 0.0931 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.3 on 1688 degrees of freedom
## Multiple R-squared: 0.3729, Adjusted R-squared: 0.3718
## F-statistic: 334.6 on 3 and 1688 DF, p-value: < 2.2e-16
robust_se <- coeftest(ols_model, vcov = vcovHC(ols_model, type = "HC1"))
print(robust_se)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 91.932404 23.087103 3.9820 7.125e-05 ***
## lunch -0.448743 0.016584 -27.0583 < 2.2e-16 ***
## log(enroll) -5.399152 1.131175 -4.7730 1.971e-06 ***
## log(exppp) 3.524751 2.353737 1.4975 0.1344
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# (ii)
white_test <- bptest(ols_model, ~ fitted(ols_model) + I(fitted(ols_model)^2), data = meap00_01)
print(white_test)
##
## studentized Breusch-Pagan test
##
## data: ols_model
## BP = 229.78, df = 2, p-value < 2.2e-16
# (iii)
residuals_squared <- residuals(ols_model)^2
log_residuals_squared <- log(residuals_squared)
aux_model <- lm(log_residuals_squared ~ fitted(ols_model) + I(fitted(ols_model)^2), data = meap00_01)
g_hat <- fitted(aux_model)
h_hat <- exp(g_hat)
wls_model <- lm(math4 ~ lunch + log(enroll) + log(exppp), data = meap00_01, weights = 1 / h_hat)
summary(wls_model)
##
## Call:
## lm(formula = math4 ~ lunch + log(enroll) + log(exppp), data = meap00_01,
## weights = 1/h_hat)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -10.7122 -1.1914 0.1403 1.3741 3.9151
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 50.47815 16.51032 3.057 0.002268 **
## lunch -0.44859 0.01461 -30.697 < 2e-16 ***
## log(enroll) -2.64728 0.83594 -3.167 0.001569 **
## log(exppp) 6.47420 1.68588 3.840 0.000127 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.9 on 1688 degrees of freedom
## Multiple R-squared: 0.36, Adjusted R-squared: 0.3588
## F-statistic: 316.5 on 3 and 1688 DF, p-value: < 2.2e-16
# (iv)
wls_robust_se <- coeftest(wls_model, vcov = vcovHC(wls_model, type = "HC1"))
print(wls_robust_se)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 50.478152 18.925179 2.6672 0.0077206 **
## lunch -0.448592 0.014254 -31.4711 < 2.2e-16 ***
## log(enroll) -2.647282 1.054608 -2.5102 0.0121590 *
## log(exppp) 6.474200 1.812868 3.5712 0.0003652 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# (v) Focus on the coefficient of log(exppp)