# Load necessary libraries
library(sandwich) # For heteroskedasticity-robust standard errors
## Warning: package 'sandwich' was built under R version 4.3.3
library(lmtest) # For hypothesis tests
## Warning: package 'lmtest' was built under R version 4.3.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.3.2
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
# Step 1: Simulate the dataset
set.seed(123) # For reproducibility
n <- 54 # Number of observations
grant <- rbinom(n, 1, 0.5) # Dummy variable for training grant
scrap_87 <- runif(n, 5, 20) # Simulated scrap rates for 1987
u <- rnorm(n, mean = 0, sd = 0.5) # Error term
beta_0 <- 3 # Intercept
beta_1 <- -0.5 # Effect of grant
scrap <- exp(beta_0 + beta_1 * grant + u) # Simulated scrap rates for 1988
data <- data.frame(grant, scrap_87, scrap)
# Step 2: Simple regression
model_1 <- lm(log(scrap) ~ grant, data = data)
summary(model_1)
##
## Call:
## lm(formula = log(scrap) ~ grant, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.11822 -0.27413 -0.04895 0.21853 1.13003
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.96363 0.08671 34.177 < 2e-16 ***
## grant -0.39999 0.12497 -3.201 0.00234 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4588 on 52 degrees of freedom
## Multiple R-squared: 0.1646, Adjusted R-squared: 0.1485
## F-statistic: 10.24 on 1 and 52 DF, p-value: 0.002338
# Step 3: Add log(scrap_87) as a control
model_2 <- lm(log(scrap) ~ grant + log(scrap_87), data = data)
summary(model_2)
##
## Call:
## lm(formula = log(scrap) ~ grant + log(scrap_87), data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.11454 -0.27346 -0.04818 0.22366 1.13694
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.01840 0.46002 6.561 2.68e-08 ***
## grant -0.40078 0.12634 -3.172 0.00256 **
## log(scrap_87) -0.02189 0.18051 -0.121 0.90395
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4633 on 51 degrees of freedom
## Multiple R-squared: 0.1648, Adjusted R-squared: 0.1321
## F-statistic: 5.033 on 2 and 51 DF, p-value: 0.01012
# Step 4: Hypothesis test for log(scrap_87)
test_result <- coeftest(model_2, vcov = vcovHC(model_2, type = "HC1")) # Robust SE
print(test_result)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.018399 0.472969 6.3818 5.151e-08 ***
## grant -0.400781 0.125148 -3.2025 0.002348 **
## log(scrap_87) -0.021891 0.183282 -0.1194 0.905399
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Step 5: One-sided test for grant coefficient
p_value_grant <- pt(test_result["grant", "t value"], df = model_2$df.residual, lower.tail = TRUE)
cat("P-value for one-sided test of grant (H1: β_grant < 0):", p_value_grant, "\n")
## P-value for one-sided test of grant (H1: β_grant < 0): 0.001174196
# Step 6: Discuss differences with heteroskedasticity-robust SE
# Compare regular vs robust SE
cat("Comparison of standard errors (regular vs robust):\n")
## Comparison of standard errors (regular vs robust):
print(summary(model_2)$coefficients)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.01839921 0.4600177 6.5614849 2.682274e-08
## grant -0.40078089 0.1263364 -3.1723300 2.561101e-03
## log(scrap_87) -0.02189052 0.1805099 -0.1212705 9.039533e-01
print(test_result)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.018399 0.472969 6.3818 5.151e-08 ***
## grant -0.400781 0.125148 -3.2025 0.002348 **
## log(scrap_87) -0.021891 0.183282 -0.1194 0.905399
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
