#0.数据读取与环境设置
options(stringsAsFactors = FALSE)
set.seed(2026)
core_packages <- c("fixest", "ggplot2", "MatchIt")
for (p in core_packages) {
  if (!requireNamespace(p, quietly = TRUE)) install.packages(p, repos = "https://cloud.r-project.org")
  if (!requireNamespace(p, quietly = TRUE)) stop(paste("Package unavailable:", p))
}
invisible(lapply(core_packages, library, character.only = TRUE))
optional_package <- function(p) {
  if (!requireNamespace(p, quietly = TRUE)) {
    try(suppressWarnings(install.packages(p, repos = "https://cloud.r-project.org")), silent = TRUE)
  }
  requireNamespace(p, quietly = TRUE)
}
data_file <- "data all.csv"
if (!file.exists(data_file)) data_file <- file.choose()
data_file <- normalizePath(data_file, mustWork = TRUE)
output_dir <- file.path(dirname(data_file), "EQT_Results")
dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
dat <- read.csv(data_file, check.names = FALSE, fileEncoding = "UTF-8-BOM")
required_vars <- c("id", "year", "DI", "Energy_Trading", "lnPGDP", "Gov", "Opening", "lnEnergy_Str", "Urban", "lnDensity", "LAC", "CER_Policy", "Emissions_Trading", "NECP", "Innovation", "Structure_Advanced", "HumanCapital", "RoadInfra", "Informatization")
if (length(setdiff(required_vars, names(dat))) > 0) stop(paste("Missing variables:", paste(setdiff(required_vars, names(dat)), collapse = ", ")))
if (anyDuplicated(dat[c("id", "year")]) > 0) stop("Duplicate province-year observations")
if (anyNA(dat[required_vars])) stop("Missing values in required variables")
dat <- dat[order(dat$id, dat$year), , drop = FALSE]
row.names(dat) <- NULL
control_vars <- c("lnPGDP", "Gov", "Opening", "lnEnergy_Str", "Urban", "lnDensity")
control_rhs <- paste(control_vars, collapse = " + ")
policy_year <- 2016L
treated_ids <- unique(dat$id[dat$Energy_Trading == 1])
dat$treat_group <- as.integer(dat$id %in% treated_ids)
if (!all(dat$Energy_Trading == as.integer(dat$treat_group == 1 & dat$year >= policy_year))) stop("Treatment timing differs from 2016 DID design")
winsorize_1pct <- function(x) {
  q <- as.numeric(quantile(x, c(0.01, 0.99), na.rm = TRUE))
  pmin(pmax(x, q[1]), q[2])
}
export_models <- function(models, filename) {
  pieces <- lapply(names(models), function(nm) {
    m <- models[[nm]]
    tab <- as.data.frame(fixest::coeftable(m))
    names(tab)[seq_len(4)] <- c("estimate", "std_error", "t_value", "p_value")
    p <- tab$p_value
    tab$stars <- ifelse(p < 0.01, "***", ifelse(p < 0.05, "**", ifelse(p < 0.10, "*", "")))
    data.frame(model = nm, term = rownames(tab), tab, n = nobs(m), r2 = as.numeric(fixest::fitstat(m, "r2")[[1]]), row.names = NULL, check.names = FALSE)
  })
  result <- do.call(rbind, pieces)
  write.csv(result, file.path(output_dir, paste0(filename, ".csv")), row.names = FALSE, fileEncoding = "UTF-8")
  writeLines(capture.output(fixest::etable(models, digits = 4, signif.code = c("***" = 0.01, "**" = 0.05, "*" = 0.10))), file.path(output_dir, paste0(filename, ".txt")), useBytes = TRUE)
  print(fixest::etable(models, digits = 4, signif.code = c("***" = 0.01, "**" = 0.05, "*" = 0.10)))
  invisible(result)
}
write.csv(dat, file.path(output_dir, "Analysis_Data.csv"), row.names = FALSE, fileEncoding = "UTF-8")
writeLines(capture.output(sessionInfo()), file.path(output_dir, "R_Session_Info.txt"))

#1.描述性统计
stat_vars <- c("lnY", "lnK", "lnC", "lnL", "Price_Carbon", "DI", "lnPGDP", "Gov", "Opening", "lnEnergy_Str", "Urban", "lnDensity", "LAC", "Innovation", "Structure_Advanced", "HumanCapital", "RoadInfra", "Informatization")
stat_vars <- intersect(stat_vars, names(dat))
descriptive <- do.call(rbind, lapply(stat_vars, function(v) {
  x <- dat[[v]]
  data.frame(variable = v, N = sum(!is.na(x)), Mean = mean(x, na.rm = TRUE), SD = sd(x, na.rm = TRUE), Min = min(x, na.rm = TRUE), Median = median(x, na.rm = TRUE), Max = max(x, na.rm = TRUE))
}))
for (v in c("Mean", "SD", "Min", "Median", "Max")) descriptive[[v]] <- sprintf("%.4f", descriptive[[v]])
write.csv(descriptive, file.path(output_dir, "01_Descriptive_Statistics.csv"), row.names = FALSE, fileEncoding = "UTF-8")
print(descriptive, row.names = FALSE)

#2.基准回归
baseline_1 <- feols(DI ~ Energy_Trading, data = dat, cluster = ~id)
baseline_2 <- feols(DI ~ Energy_Trading | id + year, data = dat, cluster = ~id)
baseline_3 <- feols(as.formula(paste("DI ~ Energy_Trading +", control_rhs)), data = dat, cluster = ~id)
baseline_4 <- feols(as.formula(paste("DI ~ Energy_Trading +", control_rhs, "| id + year")), data = dat, cluster = ~id)
baseline_models <- list("(1) No FE, No Controls" = baseline_1, "(2) TWFE" = baseline_2, "(3) Controls" = baseline_3, "(4) TWFE + Controls" = baseline_4)
export_models(baseline_models, "02_Baseline_Regression")

#3.平行趋势与动态效应
dat$event_time <- dat$year - policy_year
dat$event_time_binned <- pmax(-4, pmin(5, dat$event_time))
event_model <- feols(as.formula(paste("DI ~ i(event_time_binned, treat_group, ref = -1) +", control_rhs, "| id + year")), data = dat, cluster = ~id)
export_models(list("Event Study" = event_model), "03_Event_Study_Regression")
es_terms <- grep("^event_time_binned::", names(coef(event_model)), value = TRUE)
es_ci <- confint(event_model, level = 0.95)
es_tab <- as.data.frame(coeftable(event_model))
names(es_tab)[1:4] <- c("estimate", "std_error", "t_value", "p_value")
es_data <- data.frame(term = es_terms, event_time = as.integer(sub("^event_time_binned::(-?[0-9]+):.*$", "\\1", es_terms)), estimate = es_tab[es_terms, "estimate"], std_error = es_tab[es_terms, "std_error"], t_value = es_tab[es_terms, "t_value"], p_value = es_tab[es_terms, "p_value"], ci_low = es_ci[es_terms, 1], ci_high = es_ci[es_terms, 2], row.names = NULL)
es_data <- es_data[order(es_data$event_time), , drop = FALSE]
es_data$event_factor <- factor(es_data$event_time, levels = sort(unique(es_data$event_time)))
write.csv(es_data[, setdiff(names(es_data), "event_factor")], file.path(output_dir, "03_Event_Study_Effects.csv"), row.names = FALSE)
event_plot <- ggplot(es_data, aes(x = event_factor, y = estimate, group = 1)) + geom_hline(yintercept = 0, linetype = "dashed") + geom_errorbar(aes(ymin = ci_low, ymax = ci_high), width = 0.15, color = "#2c3e50", linewidth = 0.8) + geom_point(size = 2.6, color = "#c0392b") + geom_line(color = "#c0392b", linewidth = 0.9) + labs(x = "Years relative to EQT implementation", y = "Estimated effect on Carbon Distortion", title = "Parallel Trend Test and Dynamic Effects") + theme_classic(base_size = 12)
ggsave(file.path(output_dir, "03_Event_Study.png"), event_plot, width = 10, height = 6.5, dpi = 400)
print(event_plot)
pretrend <- wald(event_model, keep = "event_time_binned::-(4|3|2):treat_group")
writeLines(capture.output(pretrend), file.path(output_dir, "03_Pretrend_Joint_Test.txt"))
print(pretrend)

#4.合成双重差分SDID
if (!optional_package("synthdid") && optional_package("remotes")) {
  try(remotes::install_github("synth-inference/synthdid", upgrade = "never", quiet = TRUE), silent = TRUE)
}
if (requireNamespace("synthdid", quietly = TRUE)) {
  sdid_setup <- synthdid::panel.matrices(as.data.frame(dat), unit = "id", time = "year", outcome = "DI", treatment = "Energy_Trading")
  sdid_model <- synthdid::synthdid_estimate(sdid_setup$Y, sdid_setup$N0, sdid_setup$T0)
  set.seed(2026)
  sdid_se <- tryCatch(as.numeric(sqrt(vcov(sdid_model, method = "placebo", replications = 200))), error = function(e) NA_real_)
  sdid_att <- as.numeric(sdid_model)
  sdid_z <- sdid_att / sdid_se
  sdid_p <- 2 * pnorm(-abs(sdid_z))
  sdid_result <- data.frame(method = "SDID", ATT = sdid_att, SE = sdid_se, z_value = sdid_z, p_value = sdid_p, N = nrow(dat), treated_provinces = length(treated_ids))
  write.csv(sdid_result, file.path(output_dir, "04_SDID_Results.csv"), row.names = FALSE)
  pdf(file.path(output_dir, "04_SDID_Plot.pdf"), width = 10, height = 6)
  plot(sdid_model)
  dev.off()
  print(sdid_result)
} else {
  cat("SDID package unavailable; section 4 not estimated.\n")
}

#5.工具变量法
if (any(tapply(dat$LAC, dat$id, function(x) length(unique(x))) != 1)) stop("LAC must be province-invariant for the IV construction")
iv_dat <- subset(dat, !id %in% c(1, 2, 9, 21, 22))
iv_dat$IV_Trend <- iv_dat$LAC * (iv_dat$year - 2009)
win_vars_iv <- c("DI", "IV_Trend", control_vars)
for (v in win_vars_iv) iv_dat[[paste0(v, "_w")]] <- winsorize_1pct(iv_dat[[v]])
iv_ctrl <- paste(paste0(control_vars, "_w"), collapse = " + ")
iv_model <- feols(as.formula(paste("DI_w ~", iv_ctrl, "| id + year | Energy_Trading ~ IV_Trend_w")), data = iv_dat, cluster = ~id)
writeLines(capture.output(etable(iv_model, stage = 1:2, digits = 4, signif.code = c("***" = 0.01, "**" = 0.05, "*" = 0.10), fitstat = c("n", "ivf"))), file.path(output_dir, "05_IV_First_Second_Stage.txt"))
print(etable(iv_model, stage = 1:2, digits = 4, signif.code = c("***" = 0.01, "**" = 0.05, "*" = 0.10), fitstat = c("n", "ivf")))
export_models(list("IV Second Stage" = iv_model), "05_IV_Second_Stage")
write.csv(iv_dat[, c("id", "year", "LAC", "IV_Trend", "IV_Trend_w")], file.path(output_dir, "05_IV_Construction.csv"), row.names = FALSE)

#6.全组合穷举安慰剂检验
all_ids <- sort(unique(dat$id))
combinations <- combn(seq_along(all_ids), length(treated_ids))
M <- ncol(combinations)
X <- model.matrix(as.formula(paste("~ factor(id) + factor(year) +", control_rhs)), data = dat)
qr_X <- qr(X)
y_res <- as.numeric(qr.resid(qr_X, as.matrix(dat$DI)))
Fmat <- sapply(all_ids, function(k) as.numeric(dat$id == k & dat$year >= policy_year))
F_res <- qr.resid(qr_X, Fmat)
placebo_coef <- numeric(M)
placebo_p <- numeric(M)
groups <- split(seq_len(nrow(dat)), dat$id)
G <- length(groups)
correction <- (G / (G - 1)) * ((nrow(dat) - 1) / (nrow(dat) - qr_X$rank - 1))
for (i in seq_len(M)) {
  z <- rowSums(F_res[, combinations[, i], drop = FALSE])
  denom <- sum(z^2)
  b <- sum(z * y_res) / denom
  e <- y_res - b * z
  scores <- vapply(groups, function(idx) sum(z[idx] * e[idx]), numeric(1))
  se <- sqrt(correction * sum(scores^2) / denom^2)
  placebo_coef[i] <- b
  placebo_p[i] <- 2 * pt(-abs(b / se), df = G - 1)
  if (i %% 1000 == 0 || i == M) cat(sprintf("Placebo: %d / %d\n", i, M))
}
true_coef <- unname(coef(baseline_4)["Energy_Trading"])
if (abs(placebo_coef[which(vapply(seq_len(M), function(i) identical(all_ids[combinations[, i]], sort(treated_ids)), logical(1)))][1] - true_coef) > 1e-6) stop("Placebo residualization check failed")
placebo_results <- data.frame(combination = seq_len(M), coef = placebo_coef, pval = placebo_p)
p_left <- mean(placebo_coef <= true_coef)
p_two <- mean(abs(placebo_coef) >= abs(true_coef))
write.csv(placebo_results, file.path(output_dir, "06_Placebo_All_Combinations.csv"), row.names = FALSE)
placebo_summary <- data.frame(n_combinations = M, true_coefficient = true_coef, one_sided_p = p_left, two_sided_p = p_two)
write.csv(placebo_summary, file.path(output_dir, "06_Placebo_Exact_P.csv"), row.names = FALSE)
print(placebo_summary)
dens <- density(placebo_coef, adjust = 2.5, n = 4096)
dens_df <- data.frame(x = dens$x, y = dens$y)
scale_factor <- max(dens_df$y)
placebo_plot <- ggplot(placebo_results, aes(x = coef)) + geom_histogram(aes(y = after_stat(density)), bins = 45, fill = "grey80", color = "white", alpha = 0.55) + geom_point(aes(y = pval * scale_factor), shape = 1, color = "navy", size = 0.45, alpha = 0.08) + geom_line(data = dens_df, aes(x = x, y = y), inherit.aes = FALSE, color = "navy", linewidth = 1.1) + geom_vline(xintercept = 0, linetype = "dotted") + geom_vline(xintercept = true_coef, linetype = "dashed", color = "red", linewidth = 1.1) + geom_hline(yintercept = 0.05 * scale_factor, color = "red", linetype = "dotted") + scale_y_continuous(name = "Density", sec.axis = sec_axis(~ . / scale_factor, name = "P-value")) + labs(x = "Estimated Coefficients of Fake EQT Policy", subtitle = sprintf("Exhaustive permutation test (N = %s); one-sided p = %.3f; two-sided p = %.3f", format(M, big.mark = ","), p_left, p_two)) + theme_classic(base_size = 12)
ggsave(file.path(output_dir, "06_Exhaustive_Placebo.png"), placebo_plot, width = 10, height = 6, dpi = 400)
print(placebo_plot)

#7.PSM-DID与2013年虚拟政策
pre_means <- aggregate(dat[dat$year < policy_year, control_vars], by = list(id = dat$id[dat$year < policy_year]), FUN = mean)
pre_means$Treat_Prov <- as.integer(pre_means$id %in% treated_ids)
psm_formula <- as.formula(paste("Treat_Prov ~", control_rhs))
set.seed(123)
psm_match <- MatchIt::matchit(psm_formula, data = pre_means, method = "nearest", ratio = 3)
matched_prov <- MatchIt::match.data(psm_match)
matched_ids <- matched_prov$id
psm_dat <- subset(dat, id %in% matched_ids)
psm_model <- feols(as.formula(paste("DI ~ Energy_Trading +", control_rhs, "| id + year")), data = psm_dat, cluster = ~id)
pre_dat <- subset(dat, year < policy_year)
pre_dat$Fake_EQT <- as.integer(pre_dat$id %in% treated_ids & pre_dat$year >= 2013)
time_placebo_model <- feols(as.formula(paste("DI ~ Fake_EQT +", control_rhs, "| id + year")), data = pre_dat, cluster = ~id)
export_models(list("Baseline DID" = baseline_4, "PSM-DID" = psm_model, "Time Placebo (2013)" = time_placebo_model), "07_PSM_and_Time_Placebo")
write.csv(data.frame(id = matched_prov$id, match_weight = matched_prov$weights, treated = matched_prov$Treat_Prov), file.path(output_dir, "07_PSM_Matched_Provinces.csv"), row.names = FALSE)
writeLines(capture.output(summary(psm_match)), file.path(output_dir, "07_PSM_Balance.txt"))

#8.排除其他政策干扰
policy_models <- list(
  "CER Policy" = feols(as.formula(paste("DI ~ Energy_Trading + CER_Policy +", control_rhs, "| id + year")), data = dat, cluster = ~id),
  "Emissions Trading" = feols(as.formula(paste("DI ~ Energy_Trading + Emissions_Trading +", control_rhs, "| id + year")), data = dat, cluster = ~id),
  "NECP" = feols(as.formula(paste("DI ~ Energy_Trading + NECP +", control_rhs, "| id + year")), data = dat, cluster = ~id),
  "All Three Policies" = feols(as.formula(paste("DI ~ Energy_Trading + CER_Policy + Emissions_Trading + NECP +", control_rhs, "| id + year")), data = dat, cluster = ~id)
)
export_models(policy_models, "08_Concurrent_Policy_Check")

#9.剔除直辖市、极端值与疫情样本
dat_no_muni <- subset(dat, !id %in% c(1, 2, 9, 22))
dat_win <- dat
for (v in c("DI", control_vars)) dat_win[[v]] <- winsorize_1pct(dat_win[[v]])
dat_pre_covid <- subset(dat, year <= 2019)
robust_models <- list(
  "Excl. Municipalities" = feols(as.formula(paste("DI ~ Energy_Trading +", control_rhs, "| id + year")), data = dat_no_muni, cluster = ~id),
  "Winsorized 1%" = feols(as.formula(paste("DI ~ Energy_Trading +", control_rhs, "| id + year")), data = dat_win, cluster = ~id),
  "Pre-COVID Sample" = feols(as.formula(paste("DI ~ Energy_Trading +", control_rhs, "| id + year")), data = dat_pre_covid, cluster = ~id)
)
export_models(robust_models, "09_Alternative_Samples")

#10.双重机器学习
ml_packages <- c("DoubleML", "mlr3", "mlr3learners", "data.table", "glmnet", "ranger")
if (all(vapply(ml_packages, optional_package, logical(1)))) {
  dm <- fixest::demean(as.formula(paste("DI + Energy_Trading +", control_rhs, "~ id + year")), data = dat)
  dat_dm <- data.table::as.data.table(dm)
  dat_dm[, id := dat$id]
  dml_data <- DoubleML::DoubleMLClusterData$new(dat_dm, y_col = "DI", d_cols = "Energy_Trading", x_cols = control_vars, cluster_cols = "id")
  fit_dml <- function(learner, n_folds) {
    if (learner == "Random Forest") {
      ml_l <- mlr3::lrn("regr.ranger", num.trees = 500)
    } else {
      ml_l <- mlr3::lrn("regr.cv_glmnet", alpha = 1, s = "lambda.min")
    }
    ml_m <- ml_l$clone()
    fit <- DoubleML::DoubleMLPLR$new(dml_data, ml_l = ml_l, ml_m = ml_m, n_folds = n_folds, score = "partialling out")
    fit$fit()
    data.frame(Learner = learner, Folds = n_folds, Estimate = as.numeric(fit$coef[1]), SE = as.numeric(fit$se[1]), t_value = as.numeric(fit$t_stat[1]), p_value = as.numeric(fit$pval[1]), N = nrow(dat_dm), Province_Clustering = "Yes", Controls = "Yes", Province_FE = "Yes", Year_FE = "Yes")
  }
  set.seed(1234)
  dml_results <- rbind(fit_dml("Random Forest", 5), fit_dml("Lasso", 5), fit_dml("Random Forest", 3), fit_dml("Random Forest", 10))
  dml_results$stars <- ifelse(dml_results$p_value < 0.01, "***", ifelse(dml_results$p_value < 0.05, "**", ifelse(dml_results$p_value < 0.1, "*", "")))
  write.csv(dml_results, file.path(output_dir, "10_DoubleML_Clustered.csv"), row.names = FALSE)
  print(dml_results)
} else {
  cat("DoubleML packages unavailable; section 10 not estimated.\n")
}

#11.机制检验两步法
dat$ln_Innovation <- log1p(dat$Innovation)
dat$ln_Structure <- log1p(dat$Structure_Advanced)
mechanism_first <- list(
  "Innovation" = feols(as.formula(paste("ln_Innovation ~ Energy_Trading +", control_rhs, "| id + year")), data = dat, cluster = ~id),
  "Structure Advanced" = feols(as.formula(paste("ln_Structure ~ Energy_Trading +", control_rhs, "| id + year")), data = dat, cluster = ~id)
)
mechanism_second <- list(
  "Innovation and DI" = feols(as.formula(paste("DI ~ Energy_Trading + ln_Innovation +", control_rhs, "| id + year")), data = dat, cluster = ~id),
  "Structure and DI" = feols(as.formula(paste("DI ~ Energy_Trading + ln_Structure +", control_rhs, "| id + year")), data = dat, cluster = ~id),
  "Both Mediators and DI" = feols(as.formula(paste("DI ~ Energy_Trading + ln_Innovation + ln_Structure +", control_rhs, "| id + year")), data = dat, cluster = ~id)
)
export_models(mechanism_first, "11_Mechanism_First_Step")
export_models(mechanism_second, "11_Mechanism_Second_Step")

#12.调节效应检验
moderators <- c("HumanCapital", "RoadInfra", "Informatization")
moderation_models <- list()
marginal_tables <- list()
for (v in moderators) {
  centered <- paste0(v, "_c")
  dat[[centered]] <- dat[[v]] - mean(dat[[v]], na.rm = TRUE)
  f <- as.formula(paste("DI ~ Energy_Trading *", centered, "+", control_rhs, "| id + year"))
  m <- feols(f, data = dat, cluster = ~id)
  moderation_models[[v]] <- m
  interaction <- paste0("Energy_Trading:", centered)
  if (!interaction %in% names(coef(m))) stop(paste("Missing interaction coefficient:", interaction))
  vals <- seq(min(dat[[v]]), max(dat[[v]]), length.out = 100)
  cent <- vals - mean(dat[[v]])
  bb <- coef(m)
  vv <- vcov(m)
  effects <- bb["Energy_Trading"] + bb[interaction] * cent
  ses <- sqrt(vv["Energy_Trading", "Energy_Trading"] + cent^2 * vv[interaction, interaction] + 2 * cent * vv["Energy_Trading", interaction])
  crit <- qt(0.975, df = length(unique(dat$id)) - 1)
  eff_tab <- data.frame(moderator = v, moderator_value = vals, effect = as.numeric(effects), ci_low = as.numeric(effects - crit * ses), ci_high = as.numeric(effects + crit * ses))
  marginal_tables[[v]] <- eff_tab
  p <- ggplot(eff_tab, aes(x = moderator_value, y = effect)) + geom_ribbon(aes(ymin = ci_low, ymax = ci_high), fill = "#c5d8e6", alpha = 0.55) + geom_line(color = "#22577a", linewidth = 1.0) + geom_hline(yintercept = 0, linetype = "dashed") + labs(x = v, y = "Marginal effect of Energy_Trading", title = paste("Moderating Effect:", v)) + theme_classic(base_size = 12)
  ggsave(file.path(output_dir, paste0("12_Moderation_", v, ".png")), p, width = 8, height = 5.5, dpi = 400)
  print(p)
}
export_models(moderation_models, "12_Moderation_Regression")
write.csv(do.call(rbind, marginal_tables), file.path(output_dir, "12_Moderation_Marginal_Effects.csv"), row.names = FALSE)
write.csv(dat, file.path(output_dir, "Analysis_Data_With_Constructed_Variables.csv"), row.names = FALSE, fileEncoding = "UTF-8")