Learning Resources

# Lab 1: Loading Data in R
# ─────────────────────────────────────────────────

# Step 0 — Set your working directory
setwd("C:/Users/YourName/Downloads")

# Step 1 — Import the CSV
data <- read.csv("lalonde.csv")

# Step 2 — Examine the dataset structure
str(data)
head(data)

# Step 3 — Summary statistics for all variables
summary(data)

# Step 4 — Mean earnings in 1978 by treatment status
tapply(data$re78, data$treat,
       function(x) c(mean = mean(x), sd = sd(x), n = length(x)))

# Lab 1: Loading Data in Python
# ─────────────────────────────────────────────────
import pandas as pd
import os

# Step 0 — Set your working directory
os.chdir("C:/Users/YourName/Downloads")

# Step 1 — Import the CSV
data = pd.read_csv("lalonde.csv")

# Step 2 — Examine the dataset structure
print(data.dtypes)
print(data.head())

# Step 3 — Summary statistics for all variables
print(data.describe())

# Step 4 — Mean earnings in 1978 by treatment status
print(data.groupby("treat")["re78"].agg(["mean", "std", "count"]))

# Lab 2: Cleaning the Data in R
# ─────────────────────────────────────────────────

# Step 1 — Check for missing values
colSums(is.na(data))

# Step 2 — Verify binary variables
lapply(data[, c("treat", "married", "nodegree")], table)

# Step 3 — Create dummy variables from race string
data$black  <- as.integer(data$race == "black")
data$hispan <- as.integer(data$race == "hispan")
data$race   <- NULL

# Step 4 — Inspect the distribution of re78
summary(data$re78)
hist(data$re78, breaks = 30, main = "Distribution of 1978 Earnings",
     xlab = "re78")

# Lab 2: Cleaning the Data in Python
# ─────────────────────────────────────────────────
import pandas as pd
import matplotlib.pyplot as plt

data = pd.read_csv("lalonde.csv")

# Step 1 — Check for missing values
print(data.isnull().sum())

# Step 2 — Verify binary variables
for col in ["treat", "married", "nodegree"]:
    print(col, data[col].value_counts().to_dict())

# Step 3 — Create dummy variables from race string
data["black"]  = (data["race"] == "black").astype(int)
data["hispan"] = (data["race"] == "hispan").astype(int)
data = data.drop(columns=["race"])

# Step 4 — Inspect the distribution of re78
print(data["re78"].describe())
data["re78"].hist(bins=30)
plt.title("Distribution of 1978 Earnings")
plt.xlabel("re78")
plt.show()

# Lab 3: Summary Statistics in R
# ─────────────────────────────────────────────────

# Step 0 — Load and prepare data
data <- read.csv("lalonde.csv")
data$black  <- as.integer(data$race == "black")
data$hispan <- as.integer(data$race == "hispan")
data$race   <- NULL

vars <- c("age","educ","black","hispan","married","nodegree","re74","re75","re78")

# Step 1 — Full-sample summary
summary(data[, vars])

# Step 2 — Means by treatment status
do.call(rbind, lapply(vars, function(v) {
  c(variable = v,
    control  = round(mean(data[data$treat == 0, v]), 2),
    treated  = round(mean(data[data$treat == 1, v]), 2),
    full     = round(mean(data[, v]), 2))
}))

# Step 3 — t-tests for balance
lapply(vars, function(v) {
  t.test(data[data$treat == 1, v], data[data$treat == 0, v])
})

# Lab 3: Summary Statistics in Python
# ─────────────────────────────────────────────────
import pandas as pd
from scipy import stats

data = pd.read_csv("lalonde.csv")
data["black"]  = (data["race"] == "black").astype(int)
data["hispan"] = (data["race"] == "hispan").astype(int)
data = data.drop(columns=["race"])

vars = ["age","educ","black","hispan","married","nodegree","re74","re75","re78"]

# Step 1 — Full-sample summary
print(data[vars].describe().round(2))

# Step 2 — Means by treatment status
table = data.groupby("treat")[vars].mean().T.round(2)
table.columns = ["Control", "Treated"]
table["Full"] = data[vars].mean().round(2)
print(table)

# Step 3 — t-tests for balance
for v in vars:
    t, p = stats.ttest_ind(data.loc[data.treat==1, v],
                           data.loc[data.treat==0, v])
    print(f"{v}: t = {t:.3f}, p = {p:.3f}")

library(sandwich)
library(lmtest)

# Step 1. True ATT from experimental data
exp_data <- read.csv("lalonde_experimental.csv")
means <- tapply(exp_data$re78, exp_data$treat, mean)
att <- means["1"] - means["0"]
cat("True ATT (benchmark):", round(att, 2), "\n")
# Record this. It is your benchmark for Labs 5, 6, and 7.

# Steps 2-4. OLS on non-experimental data
data <- read.csv("lalonde.csv")
data$black  <- as.integer(data$race == "black")
data$hispan <- as.integer(data$race == "hispan")
data$race   <- NULL

# Step 2. OLS with treatment only
m1 <- lm(re78 ~ treat, data = data)
coeftest(m1, vcov = vcovHC(m1, type = "HC1"))
# How does this compare to your Step 1 benchmark?

# Step 3. Add demographic controls
m2 <- lm(re78 ~ treat + age + educ + black + hispan + married + nodegree,
          data = data)
coeftest(m2, vcov = vcovHC(m2, type = "HC1"))

# Step 4. Add pre-program earnings
m3 <- lm(re78 ~ treat + age + educ + black + hispan + married + nodegree +
               re74 + re75, data = data)
coeftest(m3, vcov = vcovHC(m3, type = "HC1"))

# Step 5. How far is the fully controlled estimate from your Step 1 benchmark?
cat("OLS estimate (fully controlled):", round(coef(m3)["treat"], 2), "\n")
cat("Benchmark:", round(att, 2), "\n")
cat("Gap:", round(coef(m3)["treat"] - att, 2), "\n")

import pandas as pd
import statsmodels.formula.api as smf

# Step 1. True ATT from experimental data
exp = pd.read_csv("lalonde_experimental.csv")
att = exp[exp["treat"]==1]["re78"].mean() - exp[exp["treat"]==0]["re78"].mean()
print(f"True ATT (benchmark): ${att:,.2f}")
# Record this. It is your benchmark for Labs 5, 6, and 7.

# Steps 2-4. OLS on non-experimental data
data = pd.read_csv("lalonde.csv")
data["black"]  = (data["race"] == "black").astype(int)
data["hispan"] = (data["race"] == "hispan").astype(int)
data = data.drop(columns=["race"])

# Step 2. OLS with treatment only
m1 = smf.ols("re78 ~ treat", data=data).fit(cov_type="HC1")
print(m1.summary())
# How does this compare to your Step 1 benchmark?

# Step 3. Add demographic controls
m2 = smf.ols("re78 ~ treat + age + educ + black + hispan + married + nodegree",
              data=data).fit(cov_type="HC1")
print(m2.summary())

# Step 4. Add pre-program earnings
m3 = smf.ols("re78 ~ treat + age + educ + black + hispan + married + nodegree"
              " + re74 + re75", data=data).fit(cov_type="HC1")
print(m3.summary())

# Step 5. How far is the fully controlled estimate from your Step 1 benchmark?
ols_est = m3.params["treat"]
print(f"OLS estimate (fully controlled): ${ols_est:,.2f}")
print(f"Benchmark:                       ${att:,.2f}")
print(f"Gap:                             ${ols_est - att:,.2f}")

library(sandwich)
library(lmtest)
library(margins)
data <- read.csv("lalonde.csv")
data$black  <- as.integer(data$race == "black")
data$hispan <- as.integer(data$race == "hispan")
data$race   <- NULL

# Step 1. Probit first stage
probit <- glm(treat ~ age + educ + black + hispan + married + nodegree + re74 + re75,
               data = data, family = binomial(link = "probit"))
summary(margins(probit))

# Step 2. Predicted probabilities (propensity scores)
data$phat <- predict(probit, type = "response")
tapply(data$phat, data$treat, summary)

# Step 3. Inverse Mills ratio
xb <- predict(probit, type = "link")
data$imr <- ifelse(data$treat == 1, dnorm(xb) / pnorm(xb), -dnorm(xb) / (1 - pnorm(xb)))

# Step 4. Second-stage OLS
heck <- lm(re78 ~ treat + age + educ + black + hispan + married + nodegree +
                   re74 + re75 + imr, data = data)
coeftest(heck, vcov = vcovHC(heck, type = "HC1"))

# Step 5. Compare to benchmarks
cat("Heckman:", round(coef(heck)["treat"], 2), "\n")
# Compare to your Lab 4 Step 1 experimental ATT.
# Record this and compare to Lab 6 DiD when you get there.

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from scipy.stats import norm

data = pd.read_csv("lalonde.csv")
data["black"]  = (data["race"] == "black").astype(int)
data["hispan"] = (data["race"] == "hispan").astype(int)
data = data.drop(columns=["race"])

# Step 1. Probit first stage
probit = smf.probit(
    "treat ~ age + educ + black + hispan + married + nodegree + re74 + re75",
    data=data
).fit()
print(probit.get_margeff().summary())

# Step 2. Predicted probabilities (propensity scores)
data["phat"] = probit.predict()
print(data.groupby("treat")["phat"].describe())

# Step 3. Inverse Mills ratio
xb = probit.predict(linear=True)
data["imr"] = np.where(data["treat"] == 1, norm.pdf(xb) / norm.cdf(xb), -norm.pdf(xb) / (1 - norm.cdf(xb)))

# Step 4. Second-stage OLS
heck = smf.ols(
    "re78 ~ treat + age + educ + black + hispan + married + nodegree + re74 + re75 + imr",
    data=data
).fit(cov_type="HC1")
print(heck.summary())

# Step 5. Compare to benchmarks
print(f"Heckman: ${heck.params['treat']:,.2f}")
# Compare to your Lab 4 Step 1 experimental ATT.
# Record this and compare to Lab 6 DiD when you get there.

library(dplyr)
library(tidyr)
library(sandwich)
library(lmtest)
data <- read.csv("lalonde.csv")
data$black  <- as.integer(data$race == "black")
data$hispan <- as.integer(data$race == "hispan")
data$race   <- NULL

# Step 1. Four group means and manual DiD
data %>% group_by(treat) %>%
  summarise(mean_re74 = mean(re74), mean_re75 = mean(re75), mean_re78 = mean(re78))

did_manual <- (mean(data$re78[data$treat==1]) - mean(data$re75[data$treat==1])) -
               (mean(data$re78[data$treat==0]) - mean(data$re75[data$treat==0]))
cat("Manual DiD:", did_manual, "\n")

# Step 2. Reshape and DiD regression
data_long <- data %>%
  mutate(id = row_number()) %>%
  pivot_longer(cols = c(re75, re78), names_to = "year", values_to = "earnings") %>%
  mutate(post = as.integer(year == "re78"))

did <- lm(earnings ~ post * treat, data = data_long)
coeftest(did, vcov = vcovCL(did, cluster = ~id))

# Step 3. Add controls
did_c <- lm(earnings ~ post * treat + age + educ + black + hispan + married + nodegree,
             data = data_long)
coeftest(did_c, vcov = vcovCL(did_c, cluster = ~id))

# Step 4. Pre-period trend check
data %>% group_by(treat) %>%
  summarise(change_74_75 = mean(re75) - mean(re74))
# Is the change the same for both groups?

# Step 5. Compare to benchmark
cat("DiD estimate:", round(coef(did)["post:treat"], 2), "\n")
# Compare to your Lab 4 Step 1 experimental ATT.
# How large is the gap? What does the pre-period pattern suggest about why?

import pandas as pd
import statsmodels.formula.api as smf

data = pd.read_csv("lalonde.csv")
data["black"]  = (data["race"] == "black").astype(int)
data["hispan"] = (data["race"] == "hispan").astype(int)
data = data.drop(columns=["race"])

# Step 1. Four group means and manual DiD
print(data.groupby("treat")[["re74","re75","re78"]].mean())
did_manual = (
    (data.loc[data.treat==1,"re78"].mean() - data.loc[data.treat==1,"re75"].mean()) -
    (data.loc[data.treat==0,"re78"].mean() - data.loc[data.treat==0,"re75"].mean())
)
print("Manual DiD:", did_manual)

# Step 2. Reshape and DiD regression
data_long = pd.melt(
    data.assign(id=range(len(data))),
    id_vars=["id","treat","age","educ","black","hispan","married","nodegree","re74"],
    value_vars=["re75","re78"], var_name="year", value_name="earnings"
)
data_long["post"] = (data_long["year"] == "re78").astype(int)

did = smf.ols("earnings ~ post * treat", data=data_long).fit(cov_type="cluster", cov_kwds={"groups": data_long["id"]})
print(did.summary())

# Step 3. Add controls
did_c = smf.ols(
    "earnings ~ post * treat + age + educ + black + hispan + married + nodegree",
    data=data_long
).fit(cov_type="cluster", cov_kwds={"groups": data_long["id"]})
print(did_c.summary())

# Step 4. Pre-period trend check
trends = data.groupby("treat")[["re74","re75"]].mean()
trends["change_74_75"] = trends["re75"] - trends["re74"]
print(trends)
# Is the change the same for both groups?

# Step 5. Compare to benchmark
print(f"DiD estimate: ${did.params['post:treat']:,.2f}")
# Compare to your Lab 4 Step 1 experimental ATT.
# How large is the gap? What does the pre-period pattern suggest about why?

library(MatchIt)
library(sandwich)
library(lmtest)
data <- read.csv("lalonde.csv")
data$black  <- as.integer(data$race == "black")
data$hispan <- as.integer(data$race == "hispan")
data$race   <- NULL

# Step 1. Nearest-neighbor matching on propensity score (probit)
m_out <- matchit(treat ~ age + educ + black + hispan + married + nodegree + re74 + re75,
                 data = data, method = "nearest", distance = "probit", replace = FALSE)

# Step 2. Balance check
summary(m_out, standardize = TRUE)
# Look at standardized mean differences before and after matching.

# Step 3. Treatment effect on matched sample
m_data <- match.data(m_out)
fit <- lm(re78 ~ treat + age + educ + black + hispan + married + nodegree + re74 + re75,
          data = m_data, weights = weights)
coeftest(fit, vcov = vcovCL(fit, cluster = ~subclass))
cat("Matching ATT:", round(coef(fit)["treat"], 2), "\n")
# Compare to your Lab 4 Step 1 experimental ATT.

# Step 4. Caliper matching (Getting It Right)
m_cal <- matchit(treat ~ age + educ + black + hispan + married + nodegree + re74 + re75,
                 data = data, method = "nearest", distance = "probit", caliper = 0.01, replace = FALSE)
summary(m_cal)
m_cal_data <- match.data(m_cal)
fit_cal <- lm(re78 ~ treat + age + educ + black + hispan + married + nodegree + re74 + re75,
              data = m_cal_data, weights = weights)
coeftest(fit_cal, vcov = vcovCL(fit_cal, cluster = ~subclass))
cat("Caliper ATT:", round(coef(fit_cal)["treat"], 2), "\n")
cat("Units dropped:", nrow(data) - nrow(m_cal_data), "\n")
# How does the estimate change? Does restricting to better matches help?

import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
from sklearn.neighbors import NearestNeighbors

data = pd.read_csv("lalonde.csv")
data["black"]  = (data["race"] == "black").astype(int)
data["hispan"] = (data["race"] == "hispan").astype(int)
data = data.drop(columns=["race"])

# Step 1. Propensity scores
probit = smf.probit(
    "treat ~ age + educ + black + hispan + married + nodegree + re74 + re75",
    data=data
).fit()
data["pscore"] = probit.predict()

# Step 2. Nearest-neighbor matching
treated = data[data.treat == 1].reset_index(drop=True)
control = data[data.treat == 0].reset_index(drop=True)

nn = NearestNeighbors(n_neighbors=1)
nn.fit(control[["pscore"]])
distances, indices = nn.kneighbors(treated[["pscore"]])
matched_control = control.iloc[indices.flatten()].reset_index(drop=True)

# Step 3. Balance check
covs = ["age", "educ", "black", "hispan", "married", "nodegree", "re74", "re75"]
before = data.groupby("treat")[covs].mean().T
before.columns = ["control_before", "treated_before"]
after = pd.DataFrame({
    "control_after":  matched_control[covs].mean(),
    "treated_after":  treated[covs].mean()
})
print(pd.concat([before, after], axis=1).round(2))

# Step 4. ATT estimate
att = treated["re78"].mean() - matched_control["re78"].mean()
print(f"Matching ATT: ${att:,.2f}")
# Compare to your Lab 4 Step 1 experimental ATT.

# Getting It Right: caliper matching
caliper = 0.01
mask = distances.flatten() <= caliper
att_cal = treated.loc[mask, "re78"].mean() - matched_control.loc[mask, "re78"].mean()
print("Caliper ATT (0.01):", att_cal)
print("Treated units dropped:", (~mask).sum())