ch9

# Load the JTRAIN dataset
library(wooldridge)
data("jtrain")

# Filter for the year 1988
jtrain_1988 <- jtrain[jtrain$year == 1988, ]

# (ii) Estimate the simple regression model for 1988
model1 <- lm(log(scrap) ~ grant, data = jtrain_1988)
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

summary_model1 <- summary(model1)

# Calculate the effect of receiving a grant
grant_coef <- summary_model1$coefficients["grant", "Estimate"]
grant_effect <- (exp(grant_coef) - 1) * 100
grant_pvalue <- summary_model1$coefficients["grant", "Pr(>|t|)"]

cat("Simple Regression Results:\n")

## Simple Regression Results:

cat("Coefficient on grant:", round(grant_coef, 4), "\n")

## Coefficient on grant: 0.0566

cat("Percentage effect of receiving a grant:", round(grant_effect, 2), "%\n")

## Percentage effect of receiving a grant: 5.82 %

cat("p-value:", round(grant_pvalue, 4), "\n")

## p-value: 0.8895

# (iii) Add scrap as an explanatory variable
model2 <- lm(log(scrap) ~ grant + scrap, data = jtrain_1988)
summary_model2 <- summary(model2)
summary(model2)

## 
## Call:
## lm(formula = log(scrap) ~ grant + scrap, data = jtrain_1988)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.6417 -0.4981  0.1911  0.7428  1.0521 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.36322    0.17572  -2.067   0.0438 *  
## grant        0.26797    0.25828   1.037   0.3044    
## scrap        0.18411    0.02081   8.849 7.02e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9025 on 51 degrees of freedom
##   (103 observations deleted due to missingness)
## Multiple R-squared:  0.6058, Adjusted R-squared:  0.5903 
## F-statistic: 39.18 on 2 and 51 DF,  p-value: 4.922e-11

# Calculate the effect of receiving a grant in the expanded model
grant_coef_expanded <- summary_model2$coefficients["grant", "Estimate"]
grant_effect_expanded <- (exp(grant_coef_expanded) - 1) * 100
grant_pvalue_expanded <- summary_model2$coefficients["grant", "Pr(>|t|)"]

cat("\nExpanded Regression Results:\n")

## 
## Expanded Regression Results:

cat("Coefficient on grant:", round(grant_coef_expanded, 4), "\n")

## Coefficient on grant: 0.268

cat("Percentage effect of receiving a grant:", round(grant_effect_expanded, 2), "%\n")

## Percentage effect of receiving a grant: 30.73 %

cat("p-value:", round(grant_pvalue_expanded, 4), "\n")

## p-value: 0.3044

# (iv) Test if the parameter on scrap is 1
scrap_coef <- summary_model2$coefficients["scrap", "Estimate"]
scrap_se <- summary_model2$coefficients["scrap", "Std. Error"]
scrap_tstat <- scrap_coef / scrap_se
scrap_pvalue <- 2 * (1 - pt(abs(scrap_tstat), df = summary_model2$df[2]))

cat("\nTest of scrap coefficient = 1:\n")

## 
## Test of scrap coefficient = 1:

cat("Coefficient:", round(scrap_coef, 4), "\n")

## Coefficient: 0.1841

cat("t-statistic:", round(scrap_tstat, 4), "\n")

## t-statistic: 8.8494

cat("p-value:", round(scrap_pvalue, 4), "\n")

## p-value: 0

# (v) Repeat (iii) and (iv) using heteroskedasticity-robust standard errors
library(sandwich)
library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

model2_robust <- coeftest(model2, vcov = vcovHC)

cat("\nExpanded Regression Results with Robust SEs:\n")

## 
## Expanded Regression Results with Robust SEs:

print(model2_robust)

## 
## t test of coefficients:
## 
##              Estimate Std. Error t value  Pr(>|t|)    
## (Intercept) -0.363218   0.203040 -1.7889   0.07957 .  
## grant        0.267966   0.237909  1.1263   0.26529    
## scrap        0.184112   0.026119  7.0490 4.554e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1