# Load necessary library
library(dplyr)

# Set seed for reproducibility
set.seed(1500)

# Generate data
n <- 4000
age <- rnorm(n, mean = 25, sd = 2)
gpa <- runif(n, min = 1, max = 4)
gpa <- gpa - mean(gpa) # Center GPA

y0 <- 15000 + 0.5 * age + 1.5 * gpa + rnorm(n, mean = 1000, sd = 250)
y1 <- y0 + 1500 # ATE = 1500

# Create treatment indicator (d)
d <- rep(0, n)
d[1:1500] <- 1 # Assign treatment to the first 1500 units

# Generate earnings
earnings <- d * y1 + (1 - d) * y0

# Create data frame
data <- data.frame(age = age, gpa = gpa, y0 = y0, y1 = y1, delta = y1 - y0, d = d, earnings = earnings)

# Step 1: Auxiliary regression
aux_reg <- lm(d ~ age + gpa, data = data)
data$dhat <- predict(aux_reg)

# Step 2: Residualize dhat
data$dtilde <- data$d - data$dhat

# Step 3: Regression with residualized variable
fwl_reg <- lm(earnings ~ dtilde, data = data)
summary(fwl_reg)

# Full multivariate regression for comparison
full_reg <- lm(earnings ~ d + age + gpa, data = data)
summary(full_reg)