# Paper-distribution entry point for the replication package.
# Run this file from the repository root with either:
#   source("replication_package/run_replication_package.R")
# or:
#   Rscript replication_package/run_replication_package.R
# The script reads replication_package/Data.xlsx and rewrites the package-ready
# CSV/PDF outputs under replication_package/tables/ and replication_package/figures/.
#
# Project-owned workflow logic is inlined here so this package is self-contained.
# External-method implementations are bundled separately under:
#   - stem_method/
#   - andrews_kasy/
# and are sourced only when needed.

# Path helpers -----------------------------------------------------------------
# These utilities make the script robust to being run via `source()` or `Rscript`.
script_path <- local({
  is_path_string <- function(value) {
    is.character(value) && length(value) == 1 && !is.na(value) && nzchar(value)
  }

  normalize_dir <- function(path) {
    normalizePath(dirname(path), winslash = "/", mustWork = TRUE)
  }

  file_arg <- grep("^--file=", commandArgs(trailingOnly = FALSE), value = TRUE)
  if (length(file_arg) > 0) {
    normalize_dir(sub("^--file=", "", file_arg[1]))
  } else {
    source_files <- rev(Filter(is_path_string, lapply(sys.frames(), function(frame) {
      if (exists("ofile", envir = frame, inherits = FALSE)) {
        get("ofile", envir = frame, inherits = FALSE)
      } else {
        NULL
      }
    })))

    if (length(source_files) > 0) {
      normalize_dir(source_files[[1]])
    } else {
      warning(
        "Could not infer the script location from R. Falling back to the current working directory: ",
        normalizePath(getwd(), winslash = "/", mustWork = TRUE),
        call. = FALSE
      )
      normalizePath(getwd(), winslash = "/", mustWork = TRUE)
    }
  }
})

package_path <- function(...) {
  file.path(script_path, ...)
}

# Output and dependency setup ---------------------------------------------------
ensure_output_dirs <- function() {
  for (dir_path in c(package_path("tables"), package_path("figures"))) {
    if (!dir.exists(dir_path)) {
      dir.create(dir_path, recursive = TRUE, showWarnings = FALSE)
    }
  }
}

require_packages <- function(packages) {
  # Make the package runnable on a fresh machine: install any missing CRAN
  # dependencies, then verify they load. If installation fails (e.g. offline),
  # stop with the exact packages to install manually.
  missing <- packages[!vapply(packages, requireNamespace, logical(1), quietly = TRUE)]
  if (length(missing) > 0) {
    message("Installing missing required packages: ", paste(missing, collapse = ", "))
    utils::install.packages(missing, repos = "https://cloud.r-project.org")
    still_missing <- missing[!vapply(missing, requireNamespace, logical(1), quietly = TRUE)]
    if (length(still_missing) > 0) {
      stop(
        "Could not install required packages: ",
        paste(still_missing, collapse = ", "),
        ". Please install them manually with install.packages().",
        call. = FALSE
      )
    }
  }
}

# Data loading and column normalization ----------------------------------------
# The distributed workbook has a few historical naming variants. These helpers
# map them onto the canonical names expected by the package code below.
canonicalize_names <- function(data) {
  names(data) <- tolower(names(data))
  names(data) <- gsub("[^a-z0-9]+", "_", names(data))
  names(data) <- gsub("_+", "_", names(data))
  names(data) <- gsub("^_|_$", "", names(data))

  alias_map <- list(
    dv_log_mean_wage = c("dv_log_mean_wage", "dv_log_mean_w", "dv_logmeanwage"),
    time_fe = c("time_fe"),
    person_fe = c("person_fe"),
    skill_fe = c("skill_fe"),
    lpub_year = c("lpub_year"),
    lcit_per_year = c("lcit_per_year"),
    topfive = c("topfive"),
    toptwenty = c("toptwenty"),
    ols_method2 = c("ols_method2")
  )

  for (target in names(alias_map)) {
    hits <- names(data) %in% alias_map[[target]]
    if (any(hits)) {
      names(data)[hits] <- target
    }
  }

  data
}

load_master_data <- function() {
  raw_path <- package_path("Data.xlsx")
  if (!file.exists(raw_path)) {
    stop("Expected local workbook not found: ", raw_path, call. = FALSE)
  }

  data <- readxl::read_excel(raw_path)
  data <- canonicalize_names(data)

  for (col in intersect(c("author", "journal", "country"), names(data))) {
    data[[col]] <- as.character(data[[col]])
  }

  data
}

# Lightweight formatting helpers used by multiple tables -----------------------
write_table_csv <- function(data, filename) {
  utils::write.csv(
    data,
    file = package_path("tables", filename),
    row.names = FALSE,
    na = ""
  )
}

# Regression result extractors --------------------------------------------------
cluster_vcov <- function(model, cluster) {
  sandwich::vcovCL(model, cluster = cluster, type = "HC1")
}

clustered_lm_row <- function(model, coef_name, cluster) {
  vcov_mat <- cluster_vcov(model, cluster = cluster)
  se <- sqrt(diag(vcov_mat))[coef_name]
  estimate <- stats::coef(model)[coef_name]
  statistic <- estimate / se
  p_value <- 2 * stats::pnorm(abs(statistic), lower.tail = FALSE)

  data.frame(
    estimate = unname(estimate),
    std_error = unname(se),
    p_value = unname(p_value)
  )
}

lm_summary_row <- function(model, coef_name) {
  coef_table <- summary(model)$coefficients
  estimate <- coef_table[coef_name, "Estimate"]
  se <- coef_table[coef_name, "Std. Error"]
  p_value <- coef_table[coef_name, "Pr(>|t|)"]

  data.frame(
    estimate = unname(estimate),
    std_error = unname(se),
    p_value = unname(p_value)
  )
}

clustered_fixest_row <- function(model, coef_name) {
  ct <- fixest::coeftable(model)
  estimate <- ct[coef_name, "Estimate"]
  se <- ct[coef_name, "Std. Error"]
  p_value <- ct[coef_name, "Pr(>|t|)"]

  data.frame(
    estimate = unname(estimate),
    std_error = unname(se),
    p_value = unname(p_value)
  )
}

xtreg_fe_intercept_row <- function(data, outcome, predictor, slope_row) {
  estimate <- mean(data[[outcome]], na.rm = TRUE) - mean(data[[predictor]], na.rm = TRUE) * slope_row$estimate
  se <- abs(mean(data[[predictor]], na.rm = TRUE)) * slope_row$std_error
  p_value <- 2 * stats::pnorm(abs(estimate / se), lower.tail = FALSE)

  data.frame(
    estimate = unname(estimate),
    std_error = unname(se),
    p_value = unname(p_value)
  )
}

# External-method wrappers ------------------------------------------------------
selection_model_env <- local({
  cached_env <- NULL

  function(reset = FALSE) {
    if (is.null(cached_env) || reset) {
      cached_env <<- new.env(parent = globalenv())
      method_dir <- package_path("andrews_kasy")
      old_wd <- getwd()
      on.exit(setwd(old_wd), add = TRUE)
      setwd(method_dir)
      sys.source("metastudiesfunctions.r", envir = cached_env)
    }

    cached_env
  }
})

run_selection_model <- function(effect, standard_error, cutoff = 1.96, symmetric = TRUE, model = "t") {
  keep <- !is.na(effect) & !is.na(standard_error)
  if (!any(keep)) {
    stop("Selection model requires non-missing effect and standard-error inputs.", call. = FALSE)
  }

  env <- selection_model_env()
  fit <- NULL
  utils::capture.output({
    fit <- env$metastudies_estimation(
      X = as.numeric(effect[keep]),
      sigma = as.numeric(standard_error[keep]),
      cutoffs = cutoff,
      symmetric = symmetric,
      model = model
    )
  })

  terms <- if (identical(model, "t")) {
    c("mu", "tau", "df", "selection_probability")
  } else {
    c("mu", "tau", "selection_probability")
  }
  estimates <- as.numeric(fit$Psihat)
  standard_errors <- as.numeric(fit$SE)
  p_values <- ifelse(
    is.na(standard_errors) | standard_errors == 0,
    NA_real_,
    2 * stats::pnorm(abs(estimates / standard_errors), lower.tail = FALSE)
  )

  data.frame(
    term = terms,
    estimate = estimates,
    std_error = standard_errors,
    p_value = p_values,
    stringsAsFactors = FALSE
  )
}

selection_model_row <- function(results, term) {
  row <- results[results$term == term, c("estimate", "std_error", "p_value"), drop = FALSE]
  if (nrow(row) != 1) {
    stop("Selection-model term not found: ", term, call. = FALSE)
  }

  data.frame(
    estimate = row$estimate[1],
    std_error = row$std_error[1],
    p_value = row$p_value[1]
  )
}

sig_stars <- function(p_value) {
  if (is.na(p_value)) {
    return("")
  }
  if (p_value < 0.01) {
    return("***")
  }
  if (p_value < 0.05) {
    return("**")
  }
  if (p_value < 0.10) {
    return("*")
  }
  ""
}

fmt_num <- function(x, digits = 3) {
  ifelse(is.na(x), "", trimws(formatC(x, digits = digits, format = "f")))
}

safe_sigma <- function(beta) {
  ifelse(is.na(beta) | beta == 0, NA_real_, 1 / abs(beta))
}

# Table-specific formatting helpers --------------------------------------------
publication_bias_table <- function(rows) {
  out <- data.frame(
    column = integer(),
    method = character(),
    panel = character(),
    publication_bias = character(),
    publication_bias_se = character(),
    effect_beyond_bias = character(),
    effect_beyond_bias_se = character(),
    implied_sigma = character(),
    notes = character(),
    stringsAsFactors = FALSE
  )

  for (row in rows) {
    out[nrow(out) + 1, ] <- list(
      row$column,
      row$method,
      row$panel,
      row$publication_bias,
      row$publication_bias_se,
      row$effect_beyond_bias,
      row$effect_beyond_bias_se,
      row$implied_sigma,
      row$notes
    )
  }

  out
}

bma_labels <- c(
  `(Intercept)` = "Constant",
  se = "SE",
  lcells = "Number of cells (log)",
  annual_freq = "Annual frequency data",
  all_workers = "All workers",
  high_exp = "High level of experience",
  low_exp = "Low level of experience",
  high_educ = "High level of education",
  low_educ = "Low level of education",
  top6_lang = "Top 6 languages",
  male = "Male",
  female = "Female",
  immi = "Immigrant population",
  dv_log_mean_wage = "DV: log mean wage",
  annual = "DV_annual",
  monthly = "DV_monthly",
  daily = "DV_daily",
  hourly = "DV_hourly",
  national = "National",
  ols_method2 = "OLS",
  time_fe = "Time FE",
  person_fe = "Personal FE",
  skill_fe = "Skill FE",
  impact_factor = "Impact factor",
  lcit_per_year = "Citations/year (log)",
  published = "Published study",
  topfive = "Top5 journal",
  lpub_year = "Publication year (log)"
)

ordered_heterogeneity_terms <- unname(bma_labels)

as_numeric_data_frame <- function(data) {
  as.data.frame(lapply(data, as.numeric), stringsAsFactors = FALSE)
}

extract_bma_df <- function(model) {
  expected_variable_count <- length(ordered_heterogeneity_terms) - 1
  coef_df <- as.data.frame(
    stats::coef(
      model,
      std.coefs = FALSE,
      order.by.pip = FALSE,
      exact = TRUE,
      include.constant = TRUE
    )
  )
  raw_terms <- rownames(coef_df)
  names(coef_df) <- c("pip", "post_mean", "post_sd", "cond_pos_sign", "idx")
  coef_df$idx <- as.integer(coef_df$idx)
  coef_df$raw_term <- raw_terms

  constant_rows <- which(coef_df$idx == 0 | raw_terms %in% c("(Intercept)", "Constant"))
  if (length(constant_rows) != 1) {
    stop("Expected exactly one BMA constant row; found ", length(constant_rows), call. = FALSE)
  }

  variable_rows <- setdiff(seq_len(nrow(coef_df)), constant_rows)
  expected_indices <- seq_len(expected_variable_count)
  actual_indices <- sort(coef_df$idx[variable_rows])
  if (!identical(actual_indices, expected_indices)) {
    missing_indices <- setdiff(expected_indices, actual_indices)
    extra_indices <- setdiff(actual_indices, expected_indices)
    stop(
      "Unexpected BMA variable indices. Missing: ",
      paste(missing_indices, collapse = ", "),
      "; extra: ",
      paste(extra_indices, collapse = ", "),
      call. = FALSE
    )
  }

  coef_df$term <- NA_character_
  coef_df$term[constant_rows] <- "Constant"
  coef_df$term[variable_rows] <- ordered_heterogeneity_terms[-1][coef_df$idx[variable_rows]]
  coef_df <- coef_df[match(ordered_heterogeneity_terms, coef_df$term), ]
  rownames(coef_df) <- NULL
  coef_df
}

build_bma_robustness_table <- function(models_by_prior) {
  term_order <- ordered_heterogeneity_terms
  out <- data.frame(term = term_order, stringsAsFactors = FALSE)

  for (prior_name in names(models_by_prior)) {
    bma_df <- extract_bma_df(models_by_prior[[prior_name]])
    match_idx <- match(term_order, bma_df$term)
    if (anyNA(match_idx)) {
      stop("Missing BMA rows for prior: ", prior_name, call. = FALSE)
    }
    out[[paste0(prior_name, "_pip")]] <- bma_df$pip[match_idx]
    out[[paste0(prior_name, "_post_mean")]] <- bma_df$post_mean[match_idx]
    out[[paste0(prior_name, "_post_sd")]] <- bma_df$post_sd[match_idx]
  }

  out
}

source_stem_method <- function() {
  stem_env <- new.env(parent = baseenv())
  sys.source(package_path("stem_method", "stem_method.R"), envir = stem_env)
  stem_env
}

study_weight <- function(idstudy) {
  1 / ave(idstudy, idstudy, FUN = length)
}

weighted_mean_safe <- function(x, w) {
  keep <- !is.na(x) & !is.na(w)
  if (!any(keep)) {
    return(NA_real_)
  }
  stats::weighted.mean(x[keep], w[keep])
}

weighted_variance_safe <- function(x, w) {
  keep <- !is.na(x) & !is.na(w)
  if (!any(keep)) {
    return(NA_real_)
  }
  x <- x[keep]
  w <- w[keep]
  m <- stats::weighted.mean(x, w)
  sum(w * (x - m)^2) / sum(w)
}

stata_aweight_mean_ci <- function(x, w, conf_level = 0.95) {
  keep <- !is.na(x) & !is.na(w)
  if (!any(keep)) {
    return(list(mean = NA_real_, se = NA_real_, ci_lower = NA_real_, ci_upper = NA_real_))
  }
  x <- x[keep]
  w <- w[keep]
  n <- length(x)
  if (n <= 1) {
    m <- weighted_mean_safe(x, w)
    return(list(mean = m, se = NA_real_, ci_lower = NA_real_, ci_upper = NA_real_))
  }

  w_sum <- sum(w)
  mean_value <- weighted_mean_safe(x, w)
  # Match Stata's `mean ... [aweight=...]` convention.
  s2 <- (n / (w_sum * (n - 1))) * sum(w * (x - mean_value)^2)
  se_mean <- sqrt(s2 / n)
  alpha <- 1 - conf_level
  crit <- stats::qt(1 - alpha / 2, df = n - 1)

  list(
    mean = mean_value,
    se = se_mean,
    ci_lower = mean_value - crit * se_mean,
    ci_upper = mean_value + crit * se_mean
  )
}

term_display_to_code <- c(
  "Elasticity" = "elasticity",
  "Standard error (SE)" = "se",
  "Number of cells" = "lcells",
  "Annual frequency data" = "annual_freq",
  "All workers" = "all_workers",
  "High level of experience" = "high_exp",
  "Low level of experience" = "low_exp",
  "High level of education" = "high_educ",
  "Low level of education" = "low_educ",
  "English" = "en",
  "Bilingual" = "biling",
  "Top 6 languages" = "top6_lang",
  "Male" = "male",
  "Female" = "female",
  "Both" = "both",
  "Farmers" = "farmers",
  "North America" = "north_america",
  "Immigrants" = "immi",
  "DV_log mean wage" = "dv_log_mean_wage",
  "Annual wage" = "annual",
  "Monthly wage" = "monthly",
  "Weekly wage" = "weekly",
  "Daily wage" = "daily",
  "Hourly wage" = "hourly",
  "National" = "national",
  "OLS" = "ols_method2",
  "Time" = "time_fe",
  "Person" = "person_fe",
  "Skill" = "skill_fe",
  "Impact factor" = "impact_factor",
  "Citations" = "lcit_per_year",
  "Published" = "published",
  "Unpublished old" = "unpublished_old",
  "Top 5 journals" = "topfive",
  "Top 20 journals" = "toptwenty",
  "Publication year" = "lpub_year"
)

variable_description_specs <- list(
  list(type = "variable", section = "", label = "Elasticity", variable = "elasticity", used = FALSE,
       description = "Estimate of the negative of the inverse elasticity of substitution between the immigrant and native labor (response variable)."),
  list(type = "variable", section = "", label = "Standard error (SE)", variable = "se", used = TRUE,
       description = "Standard error of the estimated inverse elasticity. The variable is important for gauging publication bias."),
  list(type = "section", section = "Data characteristics", label = "Data characteristics", variable = NA_character_, used = FALSE, description = ""),
  list(type = "variable", section = "Data characteristics", label = "Number of cells", variable = "lcells", used = TRUE,
       description = "The logarithm of the total number of the multidimensional partitions of the labor market usually based on education, experience, region, gender, and time."),
  list(type = "variable", section = "Data characteristics", label = "Annual frequency data", variable = "annual_freq", used = TRUE,
       description = "= 1 if annual frequency of the data is used in the estimation, reference category; = 0 if decennial or 5 year span."),
  list(type = "section", section = "Structural variation", label = "Structural variation", variable = NA_character_, used = FALSE, description = ""),
  list(type = "variable", section = "Structural variation", label = "All workers", variable = "all_workers", used = TRUE,
       description = "= 1 if all workers included (including self-employed, in-school, or part-time workers), reference category; = 0 if only full-time workers included."),
  list(type = "variable", section = "Structural variation", label = "High level of experience", variable = "high_exp", used = TRUE,
       description = "= 1 if only high experienced workers included in the estimation (years of experience equal to or higher than 20)."),
  list(type = "variable", section = "Structural variation", label = "Low level of experience", variable = "low_exp", used = TRUE,
       description = "= 1 if only low experienced workers included in the estimation (years of experience lower than 20)."),
  list(type = "variable", section = "Structural variation", label = "High level of education", variable = "high_educ", used = TRUE,
       description = "= 1 if only high educated workers included in the estimation (worker has some college degree)."),
  list(type = "variable", section = "Structural variation", label = "Low level of education", variable = "low_educ", used = TRUE,
       description = "= 1 if only low educated workers included in the estimation (worker does not have any college degree, i.e. high school graduates as well as high school dropouts)."),
  list(type = "variable", section = "Structural variation", label = "English", variable = "en", used = FALSE,
       description = "= 1 if English is official language in the country; = 0 otherwise."),
  list(type = "variable", section = "Structural variation", label = "Bilingual", variable = "biling", used = FALSE,
       description = "= 1 if the country for which the elasticity is estimated has two or more official languages (multilingualism); = 0 otherwise."),
  list(type = "variable", section = "Structural variation", label = "Top 6 languages", variable = "top6_lang", used = TRUE,
       description = "= 1 if one of the official languages of the country is in top 6 languages by total number of speakers (English, Mandarin Chinese, Hindi, Spanish, Standard Arabic, French), reference category; = 0 otherwise."),
  list(type = "variable", section = "Structural variation", label = "Male", variable = "male", used = TRUE,
       description = "= 1 if only men included in the estimation."),
  list(type = "variable", section = "Structural variation", label = "Female", variable = "female", used = TRUE,
       description = "= 1 if only women included in the estimation."),
  list(type = "variable", section = "Structural variation", label = "Both", variable = "both", used = FALSE,
       description = "= 1 if women and men (pooled) included in the estimation."),
  list(type = "variable", section = "Structural variation", label = "Farmers", variable = "farmers", used = FALSE,
       description = "= 1 if only farmers included in the estimation."),
  list(type = "variable", section = "Structural variation", label = "North America", variable = "north_america", used = FALSE,
       description = "= 1 if the country for which the elasticity is estimated is the United States or Canada; = 0 otherwise."),
  list(type = "variable", section = "Structural variation", label = "Immigrants", variable = "immi", used = TRUE,
       description = "The percentage of foreign-born population."),
  list(type = "section", section = "Estimation characteristics", label = "Estimation characteristics", variable = NA_character_, used = FALSE, description = ""),
  list(type = "variable", section = "Estimation characteristics", label = "DV_log mean wage", variable = "dv_log_mean_wage", used = TRUE,
       description = "= 1 if log mean wage used as dependent variable, reference category; = 0 if mean log wage."),
  list(type = "variable", section = "Estimation characteristics", label = "Annual wage", variable = "annual", used = TRUE,
       description = "= 1 if the study uses annual wages to examine the elasticity of substitution."),
  list(type = "variable", section = "Estimation characteristics", label = "Monthly wage", variable = "monthly", used = TRUE,
       description = "= 1 if the study uses monthly wages to examine the elasticity of substitution."),
  list(type = "variable", section = "Estimation characteristics", label = "Weekly wage", variable = "weekly", used = FALSE,
       description = "= 1 if the study uses weekly wages to examine the elasticity of substitution."),
  list(type = "variable", section = "Estimation characteristics", label = "Daily wage", variable = "daily", used = TRUE,
       description = "= 1 if the study uses daily wages to examine the elasticity of substitution."),
  list(type = "variable", section = "Estimation characteristics", label = "Hourly wage", variable = "hourly", used = TRUE,
       description = "= 1 if the study uses hourly wages to examine the elasticity of substitution."),
  list(type = "variable", section = "Estimation characteristics", label = "National", variable = "national", used = TRUE,
       description = "= 1 if national approach, reference category; = 0 if regional approach."),
  list(type = "variable", section = "Estimation characteristics", label = "OLS", variable = "ols_method2", used = TRUE,
       description = "=1 if the ordinary least squares method, its variations (LS, DOLS, WLS, GLS) or LASSO estimates are used for estimation, reference category; = 0 if instrumental variables are used, including 2SLS, 3SLS, and GMM."),
  list(type = "section", section = "Fixed effects", label = "Fixed effects", variable = NA_character_, used = FALSE, description = ""),
  list(type = "variable", section = "Fixed effects", label = "Time", variable = "time_fe", used = TRUE,
       description = "= 1 if time fixed effects are included in the model."),
  list(type = "variable", section = "Fixed effects", label = "Person", variable = "person_fe", used = TRUE,
       description = "= 1 if person fixed effects are included in the model (age, location, gender, ..)."),
  list(type = "variable", section = "Fixed effects", label = "Skill", variable = "skill_fe", used = TRUE,
       description = "= 1 if skill fixed effects are included in the model (education, experience)."),
  list(type = "section", section = "Publication characteristics", label = "Publication characteristics", variable = NA_character_, used = FALSE, description = ""),
  list(type = "variable", section = "Publication characteristics", label = "Impact factor", variable = "impact_factor", used = TRUE,
       description = "The discounted recursive RePEc impact factor of the outlet."),
  list(type = "variable", section = "Publication characteristics", label = "Citations", variable = "lcit_per_year", used = TRUE,
       description = "The logarithm of the number of per-year citations of the study in Google Scholar."),
  list(type = "variable", section = "Publication characteristics", label = "Published", variable = "published", used = TRUE,
       description = "=1 if a study is published in a peer-reviewed journal, reference category; = 0 if a study is in a form of working paper or PhD Thesis."),
  list(type = "variable", section = "Publication characteristics", label = "Unpublished old", variable = "unpublished_old", used = FALSE,
       description = "= 1 if WP was not published more than 10 years."),
  list(type = "variable", section = "Publication characteristics", label = "Top 5 journals", variable = "topfive", used = TRUE,
       description = "= 1 if a study published in the top 5 economic journals."),
  list(type = "variable", section = "Publication characteristics", label = "Top 20 journals", variable = "toptwenty", used = FALSE,
       description = "= 1 if a study published in the top 20 economic journals."),
  list(type = "variable", section = "Publication characteristics", label = "Publication year", variable = "lpub_year", used = TRUE,
       description = "The logarithm of the study's publication year.")
)

build_variable_description_table <- function(het_base) {
  weights <- study_weight(het_base$idstudy)

  rows <- lapply(variable_description_specs, function(spec) {
    if (identical(spec$type, "section")) {
      return(data.frame(
        row_type = "section",
        section = spec$section,
        variable = "",
        description = "",
        mean = NA_real_,
        sd = NA_real_,
        weighted_mean = NA_real_,
        used_in_heterogeneity = FALSE,
        stringsAsFactors = FALSE
      ))
    }

    values <- het_base[[spec$variable]]
    data.frame(
      row_type = "variable",
      section = spec$section,
      variable = spec$label,
      description = spec$description,
      mean = mean(values, na.rm = TRUE),
      sd = stats::sd(values, na.rm = TRUE),
      weighted_mean = weighted_mean_safe(values, weights),
      used_in_heterogeneity = spec$used,
      stringsAsFactors = FALSE
    )
  })

  out <- do.call(rbind, rows)
  out$observations <- ""
  out$studies <- ""
  out$observations[out$row_type == "variable"][1] <- nrow(het_base)
  out$studies[out$row_type == "variable"][1] <- dplyr::n_distinct(het_base$idstudy)
  out
}

run_bma_model <- function(het_bma, g_prior, model_prior) {
  clean <- as_numeric_data_frame(het_bma)
  colnames(clean) <- NULL

  BMS::bms(
    X.data = clean,
    burn = 1e6,
    iter = 3e6,
    g = g_prior,
    mprior = model_prior,
    nmodel = 10000,
    mcmc = "bd",
    user.int = FALSE
  )
}

plot_bma_assets <- function(model, image_file, pmp_file) {
  grDevices::pdf(package_path("figures", image_file), width = 10, height = 7)
  image(model, cex = 0.7, xlab = "", main = "")
  grDevices::dev.off()

  grDevices::pdf(package_path("figures", pmp_file), width = 10, height = 7)
  plot(model)
  grDevices::dev.off()
}

format_bma_out <- function(bma_df) {
  bma_out <- data.frame(
    label = ordered_heterogeneity_terms,
    pip = NA_real_,
    post_mean = NA_real_,
    post_sd = NA_real_,
    stringsAsFactors = FALSE
  )

  for (i in seq_len(nrow(bma_df))) {
    label <- bma_df$term[i]
    hit <- bma_out$label == label
    if (!any(hit)) {
      stop("Unexpected BMA label: ", label, call. = FALSE)
    }
    bma_out$pip[hit] <- bma_df$pip[i]
    bma_out$post_mean[hit] <- bma_df$post_mean[i]
    bma_out$post_sd[hit] <- bma_df$post_sd[i]
  }

  if (anyNA(bma_out$pip) || anyNA(bma_out$post_mean)) {
    stop("BMA output has missing posterior values after relabeling.", call. = FALSE)
  }

  bma_out
}

main <- function() {
  # 0. Environment checks and common derived variables -------------------------
  # This block ensures the distributed package has everything it needs before
  # any computationally expensive estimators are run.
  require_packages(c(
    "readxl",
    "dplyr",
    "BMS",
    "multcomp",
    "sandwich",
    "fixest",
    "quadprog"
  ))
  ensure_output_dirs()

  master <- load_master_data()
  master$prec2 <- 1 / (master$se^2)
  master$se2 <- master$se^2
  master$weight <- if ("weight" %in% names(master)) master$weight else study_weight(master$idstudy)
  master$sqrtweight <- sqrt(master$weight)

  stopifnot(nrow(master) == 1091)
  stopifnot(dplyr::n_distinct(master$idstudy) == 41)

  het_base <- master[!is.na(master$immi), ]
  stopifnot(nrow(het_base) == 1087)
  stopifnot(dplyr::n_distinct(het_base$idstudy) == 40)

  # 1. Table 1: summary statistics --------------------------------------------
  # `summary_specs` mirrors the row ordering used in the manuscript table.
  summary_specs <- data.frame(
    section = c(
      "overall",
      rep("data_characteristics", 4),
      rep("structural_variation", 19),
      rep("estimation_characteristics", 14),
      rep("publication_characteristics", 5)
    ),
    variable = c(
      "overall",
      "high_cell", "high_cell", "annual_freq", "annual_freq",
      "all_workers", "all_workers", "high_exp", "low_exp", "high_educ", "low_educ",
      "en", "en", "biling", "biling", "top5_lang", "top6_lang",
      "male", "female", "both", "farmers", "farmers", "north_america", "north_america",
      "dv_log_mean_wage", "dv_log_mean_wage", "annual", "monthly", "weekly", "daily",
      "hourly", "national", "national", "ols_method2", "ols_method2", "time_fe",
      "person_fe", "skill_fe",
      "published", "published", "unpublished_old", "topfive", "toptwenty"
    ),
    level = c(
      "all",
      "1", "0", "1", "0",
      "1", "0", "1", "1", "1", "1",
      "1", "0", "1", "0", "1", "1",
      "1", "1", "1", "1", "0", "1", "0",
      "1", "0", "1", "1", "1", "1",
      "1", "1", "0", "1", "0", "1",
      "1", "1",
      "1", "0", "1", "1", "1"
    ),
    stringsAsFactors = FALSE
  )

  compute_summary_row <- function(data, variable, level) {
    subset <- data
    label <- "all"
    if (variable != "overall") {
      keep <- data[[variable]] == as.numeric(level)
      keep[is.na(keep)] <- FALSE
      subset <- data[keep, ]
      label <- paste0(variable, "=", level)
    }
    mean_value <- mean(subset$elasticity, na.rm = TRUE)
    sd_value <- stats::sd(subset$elasticity, na.rm = TRUE)
    n_value <- sum(!is.na(subset$elasticity))
    se_mean <- sd_value / sqrt(n_value)

    weighted_stats <- stata_aweight_mean_ci(subset$elasticity, subset$weight)
    weighted_mean_value <- weighted_stats$mean

    data.frame(
      group = label,
      n = nrow(subset),
      mean = mean_value,
      sigma = safe_sigma(mean_value),
      ci_lower = mean_value - 1.96 * se_mean,
      ci_upper = mean_value + 1.96 * se_mean,
      weighted_mean = weighted_mean_value,
      weighted_sigma = safe_sigma(weighted_mean_value),
      weighted_ci_lower = weighted_stats$ci_lower,
      weighted_ci_upper = weighted_stats$ci_upper
    )
  }

  summary_rows <- Map(
    compute_summary_row,
    MoreArgs = list(data = master),
    variable = summary_specs$variable,
    level = summary_specs$level
  )
  summary_table <- cbind(summary_specs["section"], do.call(rbind, summary_rows))
  write_table_csv(summary_table, "table_1_summary_statistics.csv")

  # 2. Table A2: country-level summary statistics ------------------------------
  country_ids <- sort(unique(master$idcountry))
  build_country_row <- function(subset, country_label, studies_override = NULL) {
    mean_value <- mean(subset$elasticity, na.rm = TRUE)
    sd_value <- stats::sd(subset$elasticity, na.rm = TRUE)
    n_value <- sum(!is.na(subset$elasticity))
    se_mean <- sd_value / sqrt(n_value)
    weighted_stats <- stata_aweight_mean_ci(subset$elasticity, subset$weight)

    data.frame(
      country = country_label,
      studies = if (is.null(studies_override)) dplyr::n_distinct(subset$idstudy) else studies_override,
      n = nrow(subset),
      mean = mean_value,
      sigma = safe_sigma(mean_value),
      ci_lower = mean_value - 1.96 * se_mean,
      ci_upper = mean_value + 1.96 * se_mean,
      weighted_mean = weighted_stats$mean,
      weighted_sigma = safe_sigma(weighted_stats$mean),
      weighted_ci_lower = weighted_stats$ci_lower,
      weighted_ci_upper = weighted_stats$ci_upper,
      stringsAsFactors = FALSE
    )
  }

  country_table <- do.call(
    rbind,
    lapply(country_ids, function(country_id) {
      subset <- master[master$idcountry == country_id, ]
      build_country_row(subset, subset$country[1])
    })
  )
  us_without_farmers <- master[master$idcountry == 10 & master$farmers == 0, ]
  country_table <- rbind(
    country_table,
    build_country_row(us_without_farmers, "US without farmers")
  )
  country_table <- rbind(
    country_table,
    data.frame(
      country = "",
      studies = dplyr::n_distinct(master$idstudy),
      n = nrow(master),
      mean = NA_real_,
      sigma = NA_real_,
      ci_lower = NA_real_,
      ci_upper = NA_real_,
      weighted_mean = NA_real_,
      weighted_sigma = NA_real_,
      weighted_ci_lower = NA_real_,
      weighted_ci_upper = NA_real_,
      stringsAsFactors = FALSE
    )
  )
  write_table_csv(country_table, "table_a2_summary_by_country.csv")

  # 3. Figure 1: boxplots of estimates by study --------------------------------
  # The study order matches the paper: reverse chronological publication order.
  study_medians <- stats::aggregate(
    cbind(elasticity, midyear) ~ idstudy,
    data = master,
    FUN = stats::median
  )

  study_order <- unique(master[, c("idstudy", "author", "pub_year")])
  study_order <- study_order[order(study_order$pub_year, study_order$idstudy), ]
  ordered_ids <- rev(study_order$idstudy)
  ordered_labels <- rev(study_order$author)

  master$study_factor <- factor(master$idstudy, levels = ordered_ids, labels = ordered_labels)

  study_stats <- boxplot(split(master$elasticity, master$study_factor), plot = FALSE)

  grDevices::pdf(package_path("figures", "figure_1_estimates_by_study.pdf"), width = 7.5, height = 11)
  old_par <- par(no.readonly = TRUE)
  on.exit(par(old_par), add = TRUE)
  par(
    mar = c(4.5, 12, 1.5, 1),
    bg = "white"
  )
  plot.new()
  plot.window(xlim = c(-0.9, 0.2), ylim = c(0.5, length(study_stats$names) + 0.5))
  axis(side = 1, at = seq(-0.8, 0.2, by = 0.2))
  axis(side = 2, at = seq_along(study_stats$names), labels = study_stats$names, las = 1, cex.axis = 0.8)
  abline(h = seq_along(study_stats$names), col = "gray90", lty = "solid")
  bxp(
    study_stats,
    horizontal = TRUE,
    add = TRUE,
    axes = FALSE,
    border = "black",
    boxfill = "white",
    outpch = 1
  )
  abline(v = mean(master$elasticity), col = "red", lty = 2)
  box()
  title(xlab = expression("Negative inverse elasticity of substitution (" * -1/sigma * ")"))
  grDevices::dev.off()

  # 4. Figure 2: funnel plot ---------------------------------------------------
  # The visual follows the manuscript convention of trimming the far-left tail
  # to keep the precision/effect cloud readable.
  grDevices::pdf(package_path("figures", "figure_2_funnel_plot.pdf"), width = 7, height = 5)
  funnel_subset <- master[master$elasticity >= -0.5, ]
  uwls_mean <- stats::weighted.mean(master$elasticity, master$prec2)
  simple_mean <- mean(master$elasticity)
  old_par <- par(no.readonly = TRUE)
  on.exit(par(old_par), add = TRUE)
  par(bg = "white")
  plot(
    funnel_subset$elasticity,
    funnel_subset$invse,
    xlab = expression("Negative inverse elasticity of substitution (" * -1/sigma * ")"),
    ylab = "Precision (1/SE)",
    pch = 1,
    col = "gray45",
    xaxt = "n",
    yaxt = "n"
  )
  axis(side = 1, at = seq(-0.5, 0.1, by = 0.1))
  axis(side = 2, at = seq(0, 800, by = 200), las = 1)
  abline(h = seq(0, 800, by = 100), col = "gray90", lty = "solid")
  abline(v = seq(-0.5, 0.1, by = 0.1), col = "gray95", lty = "solid")
  points(funnel_subset$elasticity, funnel_subset$invse, pch = 1, col = "gray45")
  abline(v = uwls_mean, lty = 2, col = "black")
  abline(v = simple_mean, lty = 3, col = "gray45")
  legend(
    "topleft",
    legend = c(paste0("UWLS (", fmt_num(uwls_mean, 2), ")"), paste0("Simple Mean (", fmt_num(simple_mean, 2), ")")),
    lty = c(2, 3),
    col = c("black", "gray45"),
    bty = "n",
    inset = c(0.02, 0.02)
  )
  grDevices::dev.off()

  # 5. Supplementary figure (exploratory; not numbered in the manuscript): time trend in study medians
  grDevices::pdf(package_path("figures", "figure_supplementary_trend.pdf"), width = 7, height = 5)
  old_par <- par(no.readonly = TRUE)
  on.exit(par(old_par), add = TRUE)
  par(bg = "white")
  plot(
    study_medians$midyear,
    study_medians$elasticity,
    pch = 1,
    xlab = "Median year of data used by a study",
    ylab = expression("Median negative inverse elasticity (" * -1/sigma * ")"),
    xaxt = "n",
    yaxt = "n"
  )
  axis(side = 1)
  axis(side = 2, las = 1)
  abline(v = pretty(study_medians$midyear), col = "gray95", lty = "solid")
  abline(h = pretty(study_medians$elasticity), col = "gray90", lty = "solid")
  points(study_medians$midyear, study_medians$elasticity, pch = 1)
  abline(stats::lm(elasticity ~ midyear, data = study_medians), lty = 2)
  grDevices::dev.off()

  # 6. Figure A2: boxplots of estimates by country -----------------------------
  country_subset <- master[master$idcountry != 1, ]
  country_order <- sort(unique(country_subset$country))
  country_subset$country_factor <- factor(country_subset$country, levels = rev(country_order))
  country_stats <- boxplot(split(country_subset$elasticity, country_subset$country_factor), plot = FALSE)

  grDevices::pdf(package_path("figures", "figure_a2_estimates_by_country.pdf"), width = 8, height = 6)
  old_par <- par(no.readonly = TRUE)
  on.exit(par(old_par), add = TRUE)
  par(
    mar = c(4.5, 8, 1.5, 1),
    bg = "white"
  )
  plot.new()
  plot.window(xlim = c(-0.8, 0.2), ylim = c(0.5, length(country_stats$names) + 0.5))
  axis(side = 1, at = seq(-0.8, 0.2, by = 0.2))
  axis(side = 2, at = seq_along(country_stats$names), labels = country_stats$names, las = 1, cex.axis = 0.9)
  abline(h = seq_along(country_stats$names), col = "gray90", lty = "solid")
  bxp(
    country_stats,
    horizontal = TRUE,
    add = TRUE,
    axes = FALSE,
    border = "black",
    boxfill = "white",
    outpch = 1
  )
  abline(v = mean(master$elasticity), col = "red", lty = 2)
  box()
  title(xlab = expression("Negative inverse elasticity of substitution (" * -1/sigma * ")"))
  grDevices::dev.off()

  # 7. Table 2: publication-bias corrections ----------------------------------
  # Each estimator is fit explicitly so readers can map table columns back to
  # the underlying regression call without tracing through wrapper code.
  master$inst2 <- master$inst^2
  master$k <- ave(master$idstudy, master$idstudy, FUN = length)
  master$se_weight <- 1 / (master$se * master$k)

  pet_model <- stats::lm(tstat ~ invse, data = master)
  pet_pub <- clustered_lm_row(pet_model, "(Intercept)", master$idstudy)
  pet_eff <- clustered_lm_row(pet_model, "invse", master$idstudy)

  peese_model <- stats::lm(elasticity ~ se2, data = master, weights = prec2)
  peese_pub <- clustered_lm_row(peese_model, "se2", master$idstudy)
  peese_eff <- clustered_lm_row(peese_model, "(Intercept)", master$idstudy)

  study_weight_model <- stats::lm(elasticity ~ se, data = master, weights = se_weight)
  study_pub <- clustered_lm_row(study_weight_model, "se", master$idstudy)
  study_eff <- clustered_lm_row(study_weight_model, "(Intercept)", master$idstudy)

  fe_model <- fixest::feols(
    tstat ~ invse | idstudy,
    data = master,
    cluster = ~idstudy
  )
  fe_eff <- clustered_fixest_row(fe_model, "invse")
  fe_pub <- xtreg_fe_intercept_row(master, "tstat", "invse", fe_eff)

  waap_bound <- abs(stats::weighted.mean(master$elasticity, master$prec2)) / 2.8
  waap_subset <- master[master$se < waap_bound, ]
  waap_model <- stats::lm(tstat ~ 0 + invse, data = waap_subset)
  waap_eff <- lm_summary_row(waap_model, "invse")

  stem_env <- source_stem_method()
  stem_result <- stem_env$stem(master$elasticity, master$se, stem_env$param)
  stem_estimates <- as.data.frame(stem_result$estimates)
  stem_eff_estimate <- stem_estimates$estimate[1]
  stem_eff_se <- stem_estimates$se[1]
  stem_eff_p <- 2 * stats::pnorm(abs(stem_eff_estimate / stem_eff_se), lower.tail = FALSE)

  # Endogenous Kink (EK) implementation follows Bom & Rachinger (2020), with Stata code downloaded from
  # https://sites.google.com/site/heikorachinger/codes
  ek_data <- master[, c("elasticity", "se")]
  names(ek_data) <- c("bs", "sebs")
  M <- nrow(ek_data)
  sebs2 <- ek_data$sebs^2
  wis <- 1 / sebs2
  bs_sebs <- ek_data$bs / ek_data$sebs
  ones_sebs <- 1 / ek_data$sebs
  wis_sum <- sum(wis)

  ek_fit_data <- data.frame(
    bs_sebs = bs_sebs,
    ones_sebs = ones_sebs,
    sebs = ek_data$sebs,
    constant_term = 1
  )

  lin_model <- stats::lm(bs_sebs ~ 0 + ones_sebs + constant_term, data = ek_fit_data)
  pet_coef <- stats::coef(lin_model)[["ones_sebs"]]
  t1_linreg <- pet_coef / summary(lin_model)$coefficients["ones_sebs", "Std. Error"]
  q1_lin <- sum(stats::residuals(lin_model)^2)

  sq_model <- stats::lm(bs_sebs ~ 0 + ones_sebs + sebs, data = ek_fit_data)
  peese_coef <- stats::coef(sq_model)[["ones_sebs"]]
  q1_sq <- sum(stats::residuals(sq_model)^2)

  combreg <- if (abs(t1_linreg) > stats::qt(0.975, df = M - 2)) peese_coef else pet_coef
  q1 <- if (abs(t1_linreg) > stats::qt(0.975, df = M - 2)) q1_sq else q1_lin
  sigh2hat <- max(0, M * (((q1 / (M - 2)) - 1) / wis_sum))
  sighhat <- sqrt(sigh2hat)

  a1 <- if (combreg > 1.96 * sighhat) {
    (combreg - 1.96 * sighhat) * (combreg + 1.96 * sighhat) / (2 * 1.96 * combreg)
  } else {
    0
  }

  ek_reg_data <- data.frame(bs = bs_sebs, constant = ones_sebs, pub_bias = 1)
  if (a1 > min(ek_data$sebs) && a1 < max(ek_data$sebs)) {
    sebs_a1 <- ifelse(ek_data$sebs > a1, ek_data$sebs - a1, 0)
    ek_reg_data$pub_bias <- sebs_a1 / ek_data$sebs
    ek_model <- stats::lm(bs ~ 0 + constant + pub_bias, data = ek_reg_data)
    ek_bias_name <- "pub_bias"
  } else if (a1 < min(ek_data$sebs)) {
    ek_model <- stats::lm(bs ~ 0 + constant + pub_bias, data = ek_reg_data)
    ek_bias_name <- "pub_bias"
  } else {
    ek_model <- stats::lm(bs ~ 0 + constant, data = ek_reg_data)
    ek_bias_name <- NA_character_
  }

  ek_eff <- lm_summary_row(ek_model, "constant")
  ek_pub <- if (is.na(ek_bias_name)) {
    data.frame(estimate = NA_real_, std_error = NA_real_, p_value = NA_real_)
  } else {
    lm_summary_row(ek_model, ek_bias_name)
  }

  selection_fit <- run_selection_model(
    effect = master$elasticity,
    standard_error = master$se,
    cutoff = 1.96,
    symmetric = TRUE,
    model = "t"
  )
  selection_eff <- selection_model_row(selection_fit, "mu")
  selection_pub <- selection_model_row(selection_fit, "selection_probability")

  table_2 <- publication_bias_table(list(
    list(
      column = 1, method = "PET", panel = "A",
      publication_bias = paste0(fmt_num(pet_pub$estimate), sig_stars(pet_pub$p_value)),
      publication_bias_se = fmt_num(pet_pub$std_error),
      effect_beyond_bias = paste0(fmt_num(pet_eff$estimate), sig_stars(pet_eff$p_value)),
      effect_beyond_bias_se = fmt_num(pet_eff$std_error),
      implied_sigma = fmt_num(safe_sigma(pet_eff$estimate), 1),
      notes = ""
    ),
    list(
      column = 2, method = "PEESE", panel = "A",
      publication_bias = paste0(fmt_num(peese_pub$estimate), sig_stars(peese_pub$p_value)),
      publication_bias_se = fmt_num(peese_pub$std_error),
      effect_beyond_bias = paste0(fmt_num(peese_eff$estimate), sig_stars(peese_eff$p_value)),
      effect_beyond_bias_se = fmt_num(peese_eff$std_error),
      implied_sigma = fmt_num(safe_sigma(peese_eff$estimate), 1),
      notes = ""
    ),
    list(
      column = 3, method = "Study-level", panel = "A",
      publication_bias = paste0(fmt_num(study_pub$estimate), sig_stars(study_pub$p_value)),
      publication_bias_se = fmt_num(study_pub$std_error),
      effect_beyond_bias = paste0(fmt_num(study_eff$estimate), sig_stars(study_eff$p_value)),
      effect_beyond_bias_se = fmt_num(study_eff$std_error),
      implied_sigma = fmt_num(safe_sigma(study_eff$estimate), 1),
      notes = ""
    ),
    list(
      column = 4, method = "Fixed Effects", panel = "A",
      publication_bias = paste0(fmt_num(fe_pub$estimate), sig_stars(fe_pub$p_value)),
      publication_bias_se = fmt_num(fe_pub$std_error),
      effect_beyond_bias = paste0(fmt_num(fe_eff$estimate), sig_stars(fe_eff$p_value)),
      effect_beyond_bias_se = fmt_num(fe_eff$std_error),
      implied_sigma = fmt_num(safe_sigma(fe_eff$estimate), 1),
      notes = ""
    ),
    list(
      column = 5, method = "WAAP", panel = "B",
      publication_bias = "\u2014",
      publication_bias_se = "",
      effect_beyond_bias = paste0(fmt_num(waap_eff$estimate), sig_stars(waap_eff$p_value)),
      effect_beyond_bias_se = fmt_num(waap_eff$std_error),
      implied_sigma = fmt_num(safe_sigma(waap_eff$estimate), 1),
      notes = paste("Adequately powered subset n =", nrow(waap_subset))
    ),
    list(
      column = 6, method = "Stem Method", panel = "B",
      publication_bias = "\u2014",
      publication_bias_se = "",
      effect_beyond_bias = paste0(fmt_num(stem_eff_estimate), sig_stars(stem_eff_p)),
      effect_beyond_bias_se = fmt_num(stem_eff_se),
      implied_sigma = fmt_num(safe_sigma(stem_eff_estimate), 1),
      notes = paste("n_stem =", stem_estimates$n_stem[1])
    ),
    list(
      column = 7, method = "EK", panel = "B",
      publication_bias = if (is.na(ek_pub$estimate)) "\u2014" else paste0(fmt_num(ek_pub$estimate), sig_stars(ek_pub$p_value)),
      publication_bias_se = if (is.na(ek_pub$std_error)) "" else fmt_num(ek_pub$std_error),
      effect_beyond_bias = paste0(fmt_num(ek_eff$estimate), sig_stars(ek_eff$p_value)),
      effect_beyond_bias_se = fmt_num(ek_eff$std_error),
      implied_sigma = fmt_num(safe_sigma(ek_eff$estimate), 1),
      notes = ""
    ),
    list(
      column = 8, method = "Selection Model", panel = "B",
      publication_bias = paste0("P = ", fmt_num(selection_pub$estimate)),
      publication_bias_se = fmt_num(selection_pub$std_error),
      effect_beyond_bias = paste0(fmt_num(selection_eff$estimate), sig_stars(selection_eff$p_value)),
      effect_beyond_bias_se = fmt_num(selection_eff$std_error),
      implied_sigma = fmt_num(safe_sigma(selection_eff$estimate), 1),
      notes = ""
    )
  ))
  write_table_csv(table_2, "table_2_publication_bias.csv")

  # 8. Table B1: variable descriptions for the heterogeneity sample ------------
  write_table_csv(build_variable_description_table(het_base), "table_b1_variable_description.csv")

  # 9. Table 4 and Figure 3: heterogeneity analysis ----------------------------
  het4 <- het_base[, c(
    "elasticity", "se", "lcells", "annual_freq", "all_workers", "high_exp", "low_exp",
    "high_educ", "low_educ", "top6_lang", "male", "female", "immi", "dv_log_mean_wage",
    "annual", "monthly", "daily", "hourly", "national", "ols_method2",
    "time_fe", "person_fe", "skill_fe", "impact_factor", "lcit_per_year", "published",
    "topfive", "lpub_year"
  )]
  het4 <- stats::na.omit(het4)
  stopifnot(nrow(het4) == 1087)

  # Figure 3 uses the baseline BMA specification shown in the manuscript.
  baseline_bma <- run_bma_model(het4, g_prior = "UIP", model_prior = "uniform")
  baseline_bma_df <- extract_bma_df(baseline_bma)
  plot_bma_assets(baseline_bma, "figure_3_bma_inclusion.pdf", "pmp_uip_uniform.pdf")

  # The BMA robustness runs below are intentionally disabled in the distributed
  # replication package because they take a long time to finish. If you want to
  # reproduce the robustness figures and Tables B2 and B3, uncomment this block.
  # random_models <- list(
  #   uip = run_bma_model(het4, g_prior = "UIP", model_prior = "random"),
  #   bric = run_bma_model(het4, g_prior = "BRIC", model_prior = "random"),
  #   hyper = run_bma_model(het4, g_prior = "hyper=3", model_prior = "random"),
  #   ebl = run_bma_model(het4, g_prior = "EBL", model_prior = "random")
  # )
  # uniform_models <- list(
  #   uip = baseline_bma,
  #   bric = run_bma_model(het4, g_prior = "BRIC", model_prior = "uniform"),
  #   hyper = run_bma_model(het4, g_prior = "hyper=3", model_prior = "uniform"),
  #   ebl = run_bma_model(het4, g_prior = "EBL", model_prior = "uniform")
  # )
  #
  # plot_bma_assets(random_models$uip, "bma_uip_random.pdf", "pmp_uip_random.pdf")
  # plot_bma_assets(uniform_models$uip, "bma_uip_uniform.pdf", "pmp_uip_uniform.pdf")
  # plot_bma_assets(random_models$hyper, "bma_hyper_random.pdf", "pmp_hyper_random.pdf")
  # plot_bma_assets(uniform_models$hyper, "bma_hyper_uniform.pdf", "pmp_hyper_uniform.pdf")
  # plot_bma_assets(random_models$ebl, "bma_ebl_random.pdf", "pmp_ebl_random.pdf")
  # plot_bma_assets(uniform_models$ebl, "bma_ebl_uniform.pdf", "pmp_ebl_uniform.pdf")
  # plot_bma_assets(random_models$bric, "bma_bric_random.pdf", "pmp_bric_random.pdf")
  # plot_bma_assets(uniform_models$bric, "bma_bric_uniform.pdf", "pmp_bric_uniform.pdf")
  #
  # write_table_csv(build_bma_robustness_table(random_models), "table_b2_bma_robustness_random_priors.csv")
  # write_table_csv(build_bma_robustness_table(uniform_models), "table_b3_bma_robustness_uniform_priors.csv")

  # Frequentist model averaging (FMA) is built step by step to keep the
  # quadratic-program setup inspectable in the distributed package.
  x_data <- het4[, -1]
  x_data <- cbind(constant = 1, x_data)
  x_scaled <- sapply(seq_len(ncol(x_data)), function(i) x_data[, i] / max(x_data[, i]))
  scale_vector <- as.matrix(sapply(seq_len(ncol(x_data)), function(i) max(x_data[, i])))
  y <- as.matrix(het4[, 1])
  output_colnames <- colnames(x_data)
  full_fit <- stats::lm(y ~ x_scaled - 1)
  beta_full <- as.matrix(stats::coef(full_fit))
  M <- ncol(x_scaled)
  n <- nrow(x_scaled)
  beta <- matrix(0, M, M)
  residuals_mat <- matrix(0, n, M)
  var_matrix <- matrix(0, M, M)
  bias_sq <- matrix(0, M, M)

  for (i in seq_len(M)) {
    x_subset <- as.matrix(x_scaled[, 1:i])
    ortho <- eigen(t(x_subset) %*% x_subset)
    q_mat <- ortho$vectors
    lambda <- ortho$values
    x_tilde <- x_subset %*% q_mat %*% diag(lambda^(-0.5), i, i)
    beta_star <- t(x_tilde) %*% y
    beta_hat <- q_mat %*% diag(lambda^(-0.5), i, i) %*% beta_star
    beta[1:i, i] <- beta_hat
    residuals_mat[, i] <- y - x_tilde %*% as.matrix(beta_star)
    bias_sq[, i] <- (beta[, i] - beta_full)^2
    var_matrix_star <- diag(as.numeric((t(residuals_mat[, i]) %*% residuals_mat[, i]) / (n - i)), i, i)
    var_matrix_hat <- var_matrix_star %*% (q_mat %*% diag(lambda^(-1), i, i) %*% t(q_mat))
    var_matrix[1:i, i] <- diag(var_matrix_hat)
    var_matrix[, i] <- var_matrix[, i] + bias_sq[, i]
  }

  e_k <- residuals_mat[, M]
  sigma_hat <- as.numeric((t(e_k) %*% e_k) / (n - M))
  g_mat <- t(residuals_mat) %*% residuals_mat
  a_vec <- (sigma_hat^2) * matrix(seq_len(M), ncol = 1)

  constraint_mat <- cbind(rep(1, M), diag(M))
  constraint_vec <- c(1, rep(0, M))
  qp_solution <- quadprog::solve.QP(
    Dmat = 2 * g_mat,
    dvec = -as.numeric(a_vec),
    Amat = constraint_mat,
    bvec = constraint_vec,
    meq = 1
  )
  weights <- matrix(qp_solution$solution, ncol = 1)

  beta_scaled <- beta %*% weights
  final_beta <- beta_scaled / scale_vector
  std_scaled <- sqrt(var_matrix) %*% weights
  final_std <- std_scaled / scale_vector
  fma_results <- data.frame(
    term = output_colnames,
    coefficient = as.numeric(final_beta),
    std_error = as.numeric(final_std),
    stringsAsFactors = FALSE
  )
  fma_results$p_value <- 2 * stats::pnorm(abs(fma_results$coefficient / fma_results$std_error), lower.tail = FALSE)

  bma_out <- format_bma_out(baseline_bma_df)
  fma_out <- data.frame(
    label = c("Constant", ordered_heterogeneity_terms[-1]),
    coefficient = c(fma_results$coefficient[fma_results$term == "constant"], rep(NA_real_, length(ordered_heterogeneity_terms) - 1)),
    std_error = c(fma_results$std_error[fma_results$term == "constant"], rep(NA_real_, length(ordered_heterogeneity_terms) - 1)),
    p_value = c(fma_results$p_value[fma_results$term == "constant"], rep(NA_real_, length(ordered_heterogeneity_terms) - 1)),
    stringsAsFactors = FALSE
  )

  for (i in seq_len(nrow(fma_results))) {
    if (fma_results$term[i] == "constant") {
      next
    }
    label <- if (fma_results$term[i] %in% names(bma_labels)) bma_labels[[fma_results$term[i]]] else fma_results$term[i]
    hit <- fma_out$label == label
    fma_out$coefficient[hit] <- fma_results$coefficient[i]
    fma_out$std_error[hit] <- fma_results$std_error[i]
    fma_out$p_value[hit] <- fma_results$p_value[i]
  }

  table_4 <- merge(bma_out, fma_out, by.x = "label", by.y = "label", all = TRUE, sort = FALSE)
  names(table_4) <- c("term", "pip", "post_mean", "post_sd", "coef", "se", "p_value")
  table_4 <- table_4[match(c("Constant", ordered_heterogeneity_terms[-1]), table_4$term), ]
  write_table_csv(table_4, "table_4_heterogeneity.csv")

  # 10. Table 5: best-practice estimates --------------------------------------
  # The scenario formulas correspond to the alternative empirical settings
  # reported in the manuscript table.
  data_bpe_id <- subset(
    het_base,
    select = c(
      idstudy, elasticity, se, lcells, annual_freq, all_workers, low_exp,
      male, female, immi, dv_log_mean_wage, hourly, national, person_fe,
      skill_fe, lcit_per_year, published
    )
  )
  data_bpe_id <- data_bpe_id |>
    dplyr::group_by(idstudy) |>
    dplyr::mutate(weight = 1 / dplyr::n()) |>
    dplyr::ungroup()

  bpe_model <- stats::lm(
    elasticity ~ se + lcells + annual_freq + all_workers + low_exp + male + female +
      immi + dv_log_mean_wage + hourly + national + person_fe + skill_fe +
      lcit_per_year + published,
    data = data_bpe_id
  )

  vcov_clustered <- sandwich::vcovCL(bpe_model, cluster = ~idstudy)
  base_vars <- paste0(
    "(Intercept) + 3.945469*lcells + 1*annual_freq + ",
    "0.07267709*low_exp + 0.4093836*male + 0.09015639*female + ",
    "0.09176769*immi + 1*hourly + 1*person_fe + 1*skill_fe + ",
    "2.69443*lcit_per_year + 1*published"
  )

  scenario_formulas <- c(
    "Mean log wage (Regional)" = paste0(base_vars, " + 0.551058*all_workers = 0"),
    "Log mean wage (Regional)" = paste0(base_vars, " + 0.551058*all_workers + 1*dv_log_mean_wage = 0"),
    "Mean log wage (National)" = paste0(base_vars, " + 0.551058*all_workers + 1*national = 0"),
    "Log mean wage (National)" = paste0(base_vars, " + 0.551058*all_workers + 1*dv_log_mean_wage + 1*national = 0"),
    "Full-time only" = paste0(base_vars, " = 0"),
    "All workers (incl. part-time)" = paste0(base_vars, " + 1*all_workers = 0")
  )

  scenario_results <- lapply(names(scenario_formulas), function(name) {
    hypothesis <- multcomp::glht(bpe_model, linfct = scenario_formulas[[name]], vcov = vcov_clustered)
    estimate <- as.numeric(stats::coef(hypothesis))
    ci <- confint(hypothesis)$confint
    sigma_ci <- sort(c(-1 / ci[1, "upr"], -1 / ci[1, "lwr"]))

    data.frame(
      specification = name,
      bpe = estimate,
      ci_lower = ci[1, "lwr"],
      ci_upper = ci[1, "upr"],
      ci = paste0("[", trimws(fmt_num(ci[1, "lwr"], 2)), ", ", trimws(fmt_num(ci[1, "upr"], 2)), "]"),
      implied_sigma = -1 / estimate,
      sigma_ci = paste0("[", trimws(fmt_num(sigma_ci[1], 1)), ", ", trimws(fmt_num(sigma_ci[2], 1)), "]"),
      stringsAsFactors = FALSE
    )
  })
  write_table_csv(do.call(rbind, scenario_results), "table_5_bpe.csv")

  # 11. Completion message -----------------------------------------------------
  cat("Replication package outputs created in:\n")
  cat("  ", package_path("tables"), "\n")
  cat("  ", package_path("figures"), "\n")
}

main()
