rddsim.R

true = mean(pmax(0.5 * rnorm(1e5), 0))

r = replicate(100, {
  n <- 2000
  p <- 5
  X <- matrix(rnorm(n * p), n, p)
  Z <- runif(n, -4, 4)
  cutoff <- 0
  W <- as.numeric(Z >= cutoff)
  tau <- pmax(0.5 * X[, 1], 0)
  Y <- tau * W  + 1 / (1 + exp(2 * Z)) + 0.2 * rnorm(n)

  # Compute the Imbens-Kalyanaraman MSE-optimal bandwidth for a local linear regression.
  if (require("rdd", quietly = TRUE)) {
    bandwidth <- IKbandwidth(Z, Y, cutoff)
  } else {
    bandwidth <- 1.1
  }
  # Compute kernel weights corresponding to a triangular kernel.
  dist <- abs((Z - cutoff) / bandwidth)
  sample.weights <- (1 - dist) * (dist <= 1) / bandwidth

  lmf <- lm_forest(X, Y, cbind(W, Z), sample.weights = sample.weights, gradient.weights = c(1, 0))
  tau.hat <- predict(lmf)$predictions[, 1, ]

  mean(tau.hat)
})

hist(r, breaks = 10)
abline(v = true, col = 'red')
summary(r)
# Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
# 0.06025 0.14157 0.17896 0.17591 0.20802 0.31772
true
# [1] 0.2005563