# Week 07: Regression discontinuity
# Synthetic scholarship cutoff data. No external data files required.
# ggplot2 is optional and used only for visualisation.

set.seed(20260907)

n <- 2500
c <- 0
X <- runif(n, min = -50, max = 50)          # running variable centred at cutoff
X_c <- X - c
Z <- as.integer(X_c >= 0)                  # eligibility at the cutoff
D <- Z                                     # sharp RDD: treatment received equals eligibility
W <- 20 + 0.08 * X_c + rnorm(n, sd = 2)    # auxiliary covariate with no true jump

Y0 <- 60 + 0.35 * X_c - 0.004 * X_c^2 + rnorm(n, sd = 4)
tau_cutoff <- 7
Y <- Y0 + tau_cutoff * D

dat <- data.frame(X, X_c, Z, D, W, Y, Y0)

estimate_jump <- function(outcome, bandwidth, data = dat) {
  sub <- subset(data, abs(X_c) <= bandwidth)
  fit <- lm(as.formula(paste(outcome, "~ Z * X_c")), data = sub)
  coef(fit)["Z"]
}

bandwidths <- c(10, 20, 35)
sharp_estimates <- data.frame(
  bandwidth = bandwidths,
  rdd_estimate = sapply(bandwidths, function(h) estimate_jump("Y", h))
)

covariate_jump_W <- estimate_jump("W", 20)

# Fuzzy RDD extension: eligibility changes treatment probability but does not fully determine D.
prob_D_fuzzy <- ifelse(X_c < 0, 0.12 + 0.001 * (X_c + 50), 0.78 + 0.001 * X_c)
prob_D_fuzzy <- pmin(pmax(prob_D_fuzzy, 0.02), 0.98)
D_fuzzy <- rbinom(n, size = 1, prob = prob_D_fuzzy)
Y_fuzzy <- Y0 + tau_cutoff * D_fuzzy
dat_fuzzy <- data.frame(X, X_c, Z, D = D_fuzzy, W, Y = Y_fuzzy)

delta_y <- estimate_jump("Y", 20, data = dat_fuzzy)
delta_d <- estimate_jump("D", 20, data = dat_fuzzy)
fuzzy_wald <- delta_y / delta_d

cat("\nWeek 07: Regression discontinuity\n")
cat("---------------------------------\n")
cat("Sharp RDD estimates by bandwidth:\n")
print(round(sharp_estimates, 2), row.names = FALSE)
cat("Covariate W jump at bandwidth 20: ",
    round(covariate_jump_W, 2), "\n", sep = "")

cat("\nFuzzy RDD at bandwidth 20:\n")
cat("Delta_Y(c): ", round(delta_y, 2), "\n", sep = "")
cat("Delta_D(c): ", round(delta_d, 3), "\n", sep = "")
cat("Wald ratio: ", round(fuzzy_wald, 2), "\n", sep = "")

cat("\nDataviz process:\n")
cat("Plot 1 shows the outcome discontinuity at the cutoff, using only observations close to the cutoff.\n")
cat("Plot 2 shows the fuzzy first stage: eligibility changes treatment probability, but not everyone complies.\n")

# ggplot2 is used only for visualisation. The RDD estimates above use base R.
if (requireNamespace("ggplot2", quietly = TRUE)) {
  library(ggplot2)

  # Focus the picture on observations near the cutoff; far-away observations are less informative for local identification.
  dat_near <- subset(dat, abs(X_c) <= 30)
  dat_near$side <- ifelse(dat_near$Z == 1, "Eligible side", "Ineligible side")

  # Sharp RDD plot: fit separate local linear trends on the two sides of the cutoff.
  plot_sharp_rdd <- ggplot(dat_near, aes(x = X_c, y = Y, colour = side)) +
    geom_point(alpha = 0.25, size = 1) +
    geom_smooth(data = subset(dat_near, X_c < 0), method = "lm",
                formula = y ~ x, se = FALSE, colour = "grey25") +
    geom_smooth(data = subset(dat_near, X_c >= 0), method = "lm",
                formula = y ~ x, se = FALSE, colour = "#a84d77") +
    geom_vline(xintercept = 0, linetype = "dashed", colour = "grey35") +
    labs(
      title = "Sharp RDD near the cutoff",
      x = "Running variable X - c",
      y = "Outcome Y",
      colour = "Cutoff side"
    ) +
    theme_minimal()

  dat_fuzzy_near <- subset(dat_fuzzy, abs(X_c) <= 30)
  dat_fuzzy_near$side <- ifelse(dat_fuzzy_near$Z == 1, "Eligible side", "Ineligible side")

  # Fuzzy first-stage plot: the jump is in treatment probability, not necessarily from 0 to 1.
  plot_fuzzy_first_stage <- ggplot(dat_fuzzy_near, aes(x = X_c, y = D, colour = side)) +
    geom_jitter(height = 0.04, width = 0, alpha = 0.18, size = 1) +
    geom_smooth(data = subset(dat_fuzzy_near, X_c < 0), method = "lm",
                formula = y ~ x, se = FALSE, colour = "grey25") +
    geom_smooth(data = subset(dat_fuzzy_near, X_c >= 0), method = "lm",
                formula = y ~ x, se = FALSE, colour = "#a84d77") +
    geom_vline(xintercept = 0, linetype = "dashed", colour = "grey35") +
    labs(
      title = "Fuzzy RDD first stage near the cutoff",
      x = "Running variable X - c",
      y = "Treatment receipt D",
      colour = "Cutoff side"
    ) +
    theme_minimal()

  if (interactive()) {
    print(plot_sharp_rdd)
    print(plot_fuzzy_first_stage)
  } else {
    cat("Plots created as plot_sharp_rdd and plot_fuzzy_first_stage. Run the script in RStudio to display them.\n")
  }
} else {
  cat("ggplot2 is not installed. To draw the figures, run install.packages(\"ggplot2\") once, then rerun this script.\n")
}

cat("\nInterpretation prompt:\n")
cat("Bandwidth stability is useful evidence, not proof of continuity or absence of manipulation.\n")