R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

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 JTRAIN Dataset
data("jtrain", package = "wooldridge")

# Subset Data for the Year 1988
jtrain_1988 <- subset(jtrain, year == 1988)

# Check for missing values in scrap and lscrap_1
jtrain_1988 <- na.omit(jtrain_1988[, c("scrap", "grant", "lscrap_1")])

# STEP 1: Simple Regression Model
# Model: log(scrap) = β0 + β1 * grant + u
model1 <- lm(log(scrap) ~ grant, data = jtrain_1988)
cat("\nSimple Regression Model Summary (log(scrap) ~ grant):\n")
## 
## Simple Regression Model Summary (log(scrap) ~ grant):
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
## Multiple R-squared:  0.0003744,  Adjusted R-squared:  -0.01885 
## F-statistic: 0.01948 on 1 and 52 DF,  p-value: 0.8895
# STEP 2: Add lscrap_1 as an Explanatory Variable
# Model: log(scrap) = β0 + β1 * grant + β2 * lscrap_1 + u
model2 <- lm(log(scrap) ~ grant + lscrap_1, data = jtrain_1988)
cat("\nRegression with lscrap_1 (log(scrap87)):\n")
## 
## Regression with lscrap_1 (log(scrap87)):
summary(model2)
## 
## Call:
## lm(formula = log(scrap) ~ grant + lscrap_1, data = jtrain_1988)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9146 -0.1763  0.0057  0.2308  1.5991 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.02124    0.08910   0.238   0.8126    
## grant       -0.25397    0.14703  -1.727   0.0902 .  
## lscrap_1     0.83116    0.04444  18.701   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5127 on 51 degrees of freedom
## Multiple R-squared:  0.8728, Adjusted R-squared:  0.8678 
## F-statistic: 174.9 on 2 and 51 DF,  p-value: < 2.2e-16
# STEP 3: Hypothesis Testing
# Test H0: β_grant = 0 vs H1: β_grant < 0 (one-sided test)
coef_grant <- coef(summary(model2))["grant", "Estimate"]
se_grant <- coef(summary(model2))["grant", "Std. Error"]
t_stat <- coef_grant / se_grant
p_value <- pt(t_stat, df = df.residual(model2))  # One-sided p-value
cat("\nOne-Sided Hypothesis Test for β_grant:\n")
## 
## One-Sided Hypothesis Test for β_grant:
cat("t-statistic:", t_stat, "\n")
## t-statistic: -1.727319
cat("p-value:", p_value, "\n")
## p-value: 0.04508135
# STEP 4: Robust Standard Errors
# Heteroskedasticity-robust standard errors
robust_se <- sqrt(diag(vcovHC(model2, type = "HC1")))

# Robust p-value for grant
t_stat_robust <- coef_grant / robust_se["grant"]
p_value_robust <- pt(t_stat_robust, df = df.residual(model2))
cat("\nHypothesis Test with Robust Standard Errors:\n")
## 
## Hypothesis Test with Robust Standard Errors:
cat("t-statistic (robust):", t_stat_robust, "\n")
## t-statistic (robust): -1.735089
cat("p-value (robust):", p_value_robust, "\n")
## p-value (robust): 0.04438248
# Compare Regular vs Robust Results
cat("\nComparison of Results:\n")
## 
## Comparison of Results:
cat("Standard Errors - Regular:", se_grant, "| Robust:", robust_se["grant"], "\n")
## Standard Errors - Regular: 0.1470311 | Robust: 0.1463727
cat("p-value - Regular:", p_value, "| Robust:", p_value_robust, "\n")
## p-value - Regular: 0.04508135 | Robust: 0.04438248
''
## [1] ""
## Including Plots

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.