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("jtrain")
jtrain_1988 <- subset(jtrain, year == 1988)
# (i)
# Model: log(scrap) = beta_0 + beta_1 * grant + u
model1 <- lm(log(scrap) ~ grant, data = jtrain_1988)
cat("### (i) Simple Regression Model Summary:\n")
## ### (i) Simple Regression Model Summary:
summary(model1)
##
## Call:
## lm(formula = log(scrap) ~ grant, data = jtrain_1988)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.4043 -0.9536 -0.0465 0.9636 2.8103
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.4085 0.2406 1.698 0.0954 .
## grant 0.0566 0.4056 0.140 0.8895
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.423 on 52 degrees of freedom
## (103 observations deleted due to missingness)
## Multiple R-squared: 0.0003744, Adjusted R-squared: -0.01885
## F-statistic: 0.01948 on 1 and 52 DF, p-value: 0.8895
# (ii)
# Model: log(scrap) = beta_0 + beta_1 * grant + beta_2 * employ + u
model2 <- lm(log(scrap) ~ grant + employ, data = jtrain_1988)
cat("\n### (ii) Regression Model with employ:\n")
##
## ### (ii) Regression Model with employ:
summary(model2)
##
## Call:
## lm(formula = log(scrap) ~ grant + employ, data = jtrain_1988)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.4142 -0.9370 -0.0316 0.9216 2.8395
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.31395 0.30354 1.034 0.306
## grant 0.11357 0.42663 0.266 0.791
## employ 0.00123 0.00278 0.443 0.660
##
## Residual standard error: 1.451 on 49 degrees of freedom
## (105 observations deleted due to missingness)
## Multiple R-squared: 0.00615, Adjusted R-squared: -0.03442
## F-statistic: 0.1516 on 2 and 49 DF, p-value: 0.8597
grant_pval <- summary(model2)$coefficients["grant", "Pr(>|t|)"]
cat("\nP-value for grant coefficient (one-sided test):", grant_pval, "\n")
##
## P-value for grant coefficient (one-sided test): 0.7911944
# (iii) Test the null hypothesis H0: beta[employ] = 1
employ_coefficient <- coef(model2)["employ"]
employ_se <- coef(summary(model2))["employ", "Std. Error"]
t_stat <- (employ_coefficient - 1) / employ_se
p_value <- 2 * (1 - pt(abs(t_stat), df = df.residual(model2)))
cat("\n### (iii) T-statistic for H0: beta[employ] = 1:", t_stat, "\n")
##
## ### (iii) T-statistic for H0: beta[employ] = 1: -359.3336
cat("P-value for H0: beta[employ] = 1:", p_value, "\n")
## P-value for H0: beta[employ] = 1: 0
# (iv) Repeat using heteroskedasticity-robust standard errors
robust_se <- vcovHC(model2, type = "HC1")
robust_summary <- coeftest(model2, vcov = robust_se)
cat("\n### (iv) Regression Model with Robust Standard Errors:\n")
##
## ### (iv) Regression Model with Robust Standard Errors:
print(robust_summary)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3139496 0.3070107 1.0226 0.3115
## grant 0.1135737 0.3885232 0.2923 0.7713
## employ 0.0012300 0.0020959 0.5868 0.5600
employ_robust_se <- sqrt(robust_se["employ", "employ"])
t_stat_robust <- (employ_coefficient - 1) / employ_robust_se
p_value_robust <- 2 * (1 - pt(abs(t_stat_robust), df = df.residual(model2)))
cat("\nT-statistic (Robust SE) for H0: beta[employ] = 1:", t_stat_robust, "\n")
##
## T-statistic (Robust SE) for H0: beta[employ] = 1: -476.5308
cat("P-value (Robust SE) for H0: beta[employ] = 1:", p_value_robust, "\n")
## P-value (Robust SE) for H0: beta[employ] = 1: 0