## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = requireNamespace("foreign", quietly = TRUE)
)

## ----load-data----------------------------------------------------------------
library(Rfuzzydid)

# Resolve package root when rendering inside source tree
find_local_pkg_dir <- function() {
  roots <- c(".", "..", "../..")

  for (root in roots) {
    pkg_desc <- file.path(root, "pkg", "DESCRIPTION")
    root_desc <- file.path(root, "DESCRIPTION")

    if (file.exists(pkg_desc)) {
      return(file.path(root, "pkg"))
    }

    if (file.exists(root_desc)) {
      desc <- tryCatch(read.dcf(root_desc), error = function(e) NULL)
      if (!is.null(desc) && "Package" %in% colnames(desc) && identical(desc[1, "Package"], "Rfuzzydid")) {
        return(root)
      }
    }
  }

  NULL
}

# First try installed package location, then source-tree fallback
data_path <- system.file("extdata", "inpresdata.dta", package = "Rfuzzydid")
if (identical(data_path, "")) {
  local_pkg_dir <- find_local_pkg_dir()
  local_path <- if (!is.null(local_pkg_dir)) {
    file.path(local_pkg_dir, "inst", "extdata", "inpresdata.dta")
  } else {
    ""
  }

  if (file.exists(local_path)) {
    data_path <- local_path
  } else {
    stop("Bundled dataset 'inpresdata.dta' was not found in installed package or source-tree extdata.")
  }
}
raw <- foreign::read.dta(data_path)

## ----sample-construction------------------------------------------------------
# Define cohort groups
cohort0 <- raw$p504thn >= 57 & raw$p504thn <= 62
cohort1 <- raw$p504thn >= 68 & raw$p504thn <= 72

# Keep observations with non-missing outcome and treatment
keep <- !is.na(raw$lhwage) & !is.na(raw$yeduc) & (cohort0 | cohort1)

dat <- raw[keep, c("birthpl", "lhwage", "yeduc", "p504thn")]

# Time indicator: t = 1 for later cohort
dat$t <- as.integer(dat$p504thn >= 68 & dat$p504thn <= 72)

## ----group-construction-------------------------------------------------------
pval_threshold <- 0.5
districts <- sort(unique(dat$birthpl))
gstar_map <- setNames(integer(length(districts)), as.character(districts))

for (g in districts) {
  i <- dat$birthpl == g
  sub <- dat[i, ]
  tab <- table(sub$yeduc, sub$t)
  pval <- suppressWarnings(chisq.test(tab)$p.value)
  evol <- mean(sub$yeduc[sub$t == 1]) - mean(sub$yeduc[sub$t == 0])
  if (!is.na(pval) && pval < pval_threshold) {
    gstar_map[[as.character(g)]] <- if (evol > 0) 1L else if (evol < 0) -1L else 0L
  } else {
    gstar_map[[as.character(g)]] <- 0L
  }
}

dat$gstar <- unname(gstar_map[as.character(dat$birthpl)])

## ----arm-inc------------------------------------------------------------------
arm_inc_data <- dat[dat$gstar >= 0, ]
arm_inc_data$G <- as.integer(arm_inc_data$gstar == 1)

fit_inc <- fuzzydid(
  data = arm_inc_data,
  formula = lhwage ~ yeduc,
  group = "G",
  time = "t",
  did = TRUE,
  tc = TRUE,
  cic = TRUE,
  newcateg = c(5, 8, 11, 14, 1000),
  nose = TRUE,
  backend = "native"
)

knitr::kable(fit_inc$late)

## ----arm-dec------------------------------------------------------------------
arm_dec_data <- dat[dat$gstar <= 0, ]
arm_dec_data$G <- as.integer(arm_dec_data$gstar == -1)

fit_dec <- fuzzydid(
  data = arm_dec_data,
  formula = lhwage ~ yeduc,
  group = "G",
  time = "t",
  did = TRUE,
  tc = TRUE,
  cic = TRUE,
  newcateg = c(5, 8, 11, 14, 1000),
  nose = TRUE,
  backend = "native"
)

knitr::kable(fit_dec$late)

## ----aggregate----------------------------------------------------------------
# Treatment intensity weights
dD_inc <- mean(arm_inc_data$yeduc[arm_inc_data$G == 1 & arm_inc_data$t == 1]) -
          mean(arm_inc_data$yeduc[arm_inc_data$G == 1 & arm_inc_data$t == 0]) -
          mean(arm_inc_data$yeduc[arm_inc_data$G == 0 & arm_inc_data$t == 1]) +
          mean(arm_inc_data$yeduc[arm_inc_data$G == 0 & arm_inc_data$t == 0])

dD_dec <- mean(arm_dec_data$yeduc[arm_dec_data$G == 1 & arm_dec_data$t == 1]) -
          mean(arm_dec_data$yeduc[arm_dec_data$G == 1 & arm_dec_data$t == 0]) -
          mean(arm_dec_data$yeduc[arm_dec_data$G == 0 & arm_dec_data$t == 1]) +
          mean(arm_dec_data$yeduc[arm_dec_data$G == 0 & arm_dec_data$t == 0])

p_inc <- mean(dat$gstar == 1)
p_dec <- mean(dat$gstar == -1)

w_inc <- p_inc * dD_inc
w_dec <- p_dec * (-dD_dec)

# Extract point estimates
did_inc <- fit_inc$late$estimate[fit_inc$late$estimator == "W_DID"]
tc_inc  <- fit_inc$late$estimate[fit_inc$late$estimator == "W_TC"]
cic_inc <- fit_inc$late$estimate[fit_inc$late$estimator == "W_CIC"]

did_dec <- fit_dec$late$estimate[fit_dec$late$estimator == "W_DID"]
tc_dec  <- fit_dec$late$estimate[fit_dec$late$estimator == "W_TC"]
cic_dec <- fit_dec$late$estimate[fit_dec$late$estimator == "W_CIC"]

# Weighted aggregation
did_agg <- (did_inc * w_inc + did_dec * w_dec) / (w_inc + w_dec)
tc_agg  <- (tc_inc  * w_inc + tc_dec  * w_dec) / (w_inc + w_dec)
cic_agg <- (cic_inc * w_inc + cic_dec * w_dec) / (w_inc + w_dec)

cat("Aggregate Estimates:\n")
cat(sprintf("DID = %.6f\n", did_agg))
cat(sprintf("TC  = %.6f\n", tc_agg))
cat(sprintf("CIC = %.6f\n", cic_agg))

## ----stata-parity-check-------------------------------------------------------
stata_vals <- c(
  did_inc = 0.123679191,
  tc_inc = 0.079868219,
  cic_inc = 0.074316049,
  did_dec = 0.117024113,
  tc_dec = 0.109567707,
  cic_dec = 0.114301855,
  did_agg = 0.122244473,
  tc_agg = 0.086270906,
  cic_agg = 0.082936285
)

r_vals <- c(
  did_inc = did_inc,
  tc_inc = tc_inc,
  cic_inc = cic_inc,
  did_dec = did_dec,
  tc_dec = tc_dec,
  cic_dec = cic_dec,
  did_agg = did_agg,
  tc_agg = tc_agg,
  cic_agg = cic_agg
)

cmp <- data.frame(
  estimate = names(stata_vals),
  R = unname(r_vals),
  Stata = unname(stata_vals),
  abs_diff = abs(unname(r_vals) - unname(stata_vals))
)

knitr::kable(cmp, digits = 6)