R Code for Kline (2011) - KOB as a reweighting estimator.

OB_ATT.R

# %% # obs lalonde data from Kline paper - init housekeeping
libreq(data.table, fixest, rio)
cps3 = import("cps3re74.dta") %>% setDT() %>% na.omit()
setnames(cps3, c("re78", "treat"), c("y", "W"));
xs = setdiff(colnames(cps3), c("y", 'W'))
W = cps3$W  %>% as.matrix(); Y = cps3$y  %>% as.matrix()
X = cbind(1, cps3[, ..xs]) %>% as.matrix()
X1 = X[W==1,]; X0 = X[W==0,]
N = length(W); N_t = sum(W)
# %% first way - KOB / kline - page 1
β_0 = solve(t(X0) %*% X0) %*% t(X0) %*% Y[W==0]
# X̄_t β_0
mean(X[W==1,] %*% β_0) # E[Y^0 | W == 1]
mean(Y[W==1]) - mean(X[W==1,] %*% β_0) # ATT - matches paper
# %% second way kline paper - N weights
S = matrix(0, nrow = N, ncol = N)
diag(S) = 1 - W
H0 = solve(t(X0) %*% X0)  %*% t(X) %*% S
H = X %*% H0
obW = (1/N_t) * t(W) %*% H
obW %*% Y %>% as.numeric # E[Y^0 | W == 1]
mean(Y[W==1]) - obW %*% Y %>% as.numeric
# %% third way - simpler (?!?!) - use selectors W and S
H00 = solve(t(X0) %*% X0) %*% t(X0)
wts = (1/N_t)* colSums(X[W==1,] %*% H00)
weighted.mean(Y[W==0], wts)
mean(Y[W==1]) - weighted.mean(Y[W==0], wts)

# %% 1st way take 2 - with SEs from Kline paper
Kline_OB = function(yn, wn, xs, df){
  # outcome model in treatment
  omod = feols(.[yn] ~ .[xs], data = df[get(wn) == 0], vcov = "hc1")
  b = omod$coefficients; V0 = vcov(omod)
  # point estimate
  df$Y0 = predict(omod, df)
  df[, dif := y - Y0]
  M1 = df[get(wn) == 1, mean(dif)]
  # standard error
  X = cbind(1, df[, ..xs])  %>% as.matrix()
  W = df[[wn]]  %>% as.matrix()
  N1 = sum(W)
  V1 = df[get(wn) == 1, .(var(y), .N)]  %>% as.numeric()
  V1 = V1[1]/V1[2]
  Vdif = V1 + (t(W)  %*% X %*% V0  %*% t(X)  %*% (W / (N1^2)) )
  SE = sqrt(Vdif)  %>% as.numeric()
  # return
  c(est = M1, se = SE)  %>% round(3)
}

Kline_OB('y', 'W', xs, cps3)
# %% 5th way - interacted regression
reg_adjust = function(yn, wn, xs, df, estimand = c("ATT", "ATE")) {
  estimand = match.arg(estimand)
  # for ATE, sweep out overall means, else sweep out treatment means
  # dd = ifelse(estimand == "ATE", df[, ..xs], df[get(wn) == 1, ..xs])
  dd = if(estimand == "ATE") df[, ..xs] else df[get(wn) == 1, ..xs]
  # take out means
  Xbar = apply(dd[, ..xs], 2, mean); Xdemeaned = sweep(df[, ..xs], 2, Xbar)
  colnames(Xdemeaned) = paste0(colnames(Xdemeaned), "_dem")
  # concat columns
  regdf = cbind(y = df[[yn]], w = df[[wn]], df[, ..xs], Xdemeaned)
  # interacted regression
  f = as.formula(paste0(
    "y ~ w +",
    paste(xs, collapse = "+"), "+", # controls effects
    "w", ":(", paste(colnames(Xdemeaned), collapse = "+"), ")" # treat effects
  ))
  m = fixest::feols(f, data = regdf, vcov = 'hc1')
  m
}

reg_adjust(yn, wn, xs, cps3, "ATT")
# %%