# Load required libraries
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)

# Load the meap00 dataset
data(meap00_01)

# (i) Estimate the model
model <- lm(math4 ~ lunch + log(enroll) + log(exppp), data = meap00_01)
summary(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
# (ii) Apply the special case of the White test for heteroskedasticity
bgtest(model)
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  model
## LM test = 37.934, df = 1, p-value = 7.32e-10
# (iii) Obtain the fitted values and the OLS residuals
fitted_values <- fitted(model)
residuals <- resid(model)

# Use the exp() function to obtain the WLS estimates
wls_coefs <- exp(coef(model))
cat("OLS Fitted Values:\n")
## OLS Fitted Values:
print(head(fitted_values))
##        1        2        3        4        5        6 
## 71.10372 86.27345 80.95395 77.16383 77.06931 77.25098
cat("\nOLS Residuals:\n")
## 
## OLS Residuals:
print(head(residuals))
##         1         2         3         4         5         6 
## 12.196279 -0.573456 -3.653947  8.036168 14.730695 10.449013
cat("\nWLS Coefficients:\n")
## 
## WLS Coefficients:
print(wls_coefs)
##  (Intercept)        lunch  log(enroll)   log(exppp) 
## 8.428220e+39 6.384299e-01 4.520414e-03 3.394531e+01
# (iv) Obtain the standard errors for WLS that allow misspecification of the variance function
coeftest(model, vcov = vcovHC)
## 
## t test of coefficients:
## 
##              Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept) 91.932404  23.203501   3.9620 7.742e-05 ***
## lunch       -0.448743   0.016637 -26.9732 < 2.2e-16 ***
## log(enroll) -5.399152   1.134779  -4.7579 2.123e-06 ***
## log(exppp)   3.524751   2.367336   1.4889    0.1367    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# (v) Compare the precision of OLS and WLS for estimating the effect of spending on math4
summary(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