run_counterfactual_comparison.R

library(readr)
library(dplyr)
library(magrittr)
library(tidyr)
library(parallel)

setwd("~/Desktop/replicate_MSC_sims")

#####################################################################
# SHARED SETUP
#####################################################################
death_prob <<- as.matrix(read.csv("death_probabilities_from_SSA_for_simulations.csv", header=TRUE))
ret_prob <<- as.matrix(read.csv("basic_and_strategic_retirement_probabilities_by_age.csv", header=TRUE))

# Read party control data for probability calibration
party_data <- read.csv("party_control_1861_2025.csv")

# Compute transition probabilities using data through 2020 (matching baseline)
prob.cal.end <- 2021  # inauguration-year coding: 2021 captures through 2020 election

pres_regime <- party_data %>%
  filter((year - 1) %% 4 == 0) %>%
  mutate(dem_president_lag1 = lag(dem_president, 1),
         dem_president_lag2 = lag(dem_president, 2),
         pres_switch = as.numeric(dem_president != dem_president_lag1))

cong_regime <- party_data %>%
  arrange(congress) %>%
  mutate(election_type = ifelse(congress %% 2 == 0, "presidential", "midterm"),
         dem_senate_lag1 = lag(dem_senate, 1),
         election_type_lag1 = lag(election_type, 1),
         div_gov_senate_lag1 = lag(div_gov_senate, 1),
         senate_switch = as.numeric(dem_senate != dem_senate_lag1))

# Presidential probs (filter to <= prob.cal.end)
pr <- pres_regime %>% filter(year >= 1948 & year <= prob.cal.end)
switch_after1term <- mean(pr$pres_switch[pr$dem_president_lag1 != pr$dem_president_lag2], na.rm=TRUE)
switch_after2term <- mean(pr$pres_switch[pr$dem_president_lag1 == pr$dem_president_lag2], na.rm=TRUE)

Prob_demP_DD_val <- 1 - switch_after2term
Prob_demP_RR_val <- switch_after2term
Prob_demP_RD_val <- 1 - switch_after1term
Prob_demP_DR_val <- switch_after1term

# Senate probs
cr <- cong_regime %>% filter(year >= 1948 & year <= prob.cal.end)
switch.uni.pres <- mean(cr$senate_switch[cr$div_gov_senate_lag1 == 0 & cr$election_type_lag1 == "presidential"], na.rm=TRUE)
switch.div.pres <- mean(cr$senate_switch[cr$div_gov_senate_lag1 == 1 & cr$election_type_lag1 == "presidential"], na.rm=TRUE)
switch.uni.mid <- mean(cr$senate_switch[cr$div_gov_senate_lag1 == 0 & cr$election_type_lag1 == "midterm"], na.rm=TRUE)
switch.div.mid <- mean(cr$senate_switch[cr$div_gov_senate_lag1 == 1 & cr$election_type_lag1 == "midterm"], na.rm=TRUE)

cat("=== Transition probabilities ===\n")
cat("P(win 2nd term):", round(1 - switch_after1term, 4), "\n")
cat("P(win 3rd term):", round(1 - switch_after2term, 4), "\n")
cat("Switch uni pres:", round(switch.uni.pres, 4), "\n")
cat("Switch div pres:", round(switch.div.pres, 4), "\n")
cat("Switch uni mid:", round(switch.uni.mid, 4), "\n")
cat("Switch div mid:", round(switch.div.mid, 4), "\n\n")

# Shape parameters
Dem_Shape_One_val <- 3
Dem_Shape_Two_val <- 11
Rep_Shape_One_val <- 3
Rep_Shape_Two_val <- 11

n.sims <- 1000
end.year <- 2100

#####################################################################
# HELPER: set all globals and run baseline experiment
#####################################################################
run_experiment <- function(label, start_year, pres_history, court_csv, seatid_vec) {
  cat("--- Running:", label, "---\n")
  
  # Load court data
  MyData <- read.csv(court_csv, header=TRUE)
  if("justice_ideology_nemacheck" %in% names(MyData)) MyData$justice_ideology_nemacheck <- NULL
  if("justice_ideology_nsp" %in% names(MyData)) MyData <- rename(MyData, justice_ideology = justice_ideology_nsp)
  MyData$seatid <- seatid_vec
  dataset <<- subset(MyData, select=-c(justice_reliability, justice))
  
  # Set globals
  current_year <<- start_year
  president_history_init <<- pres_history
  
  Prob_demP_DD <<- Prob_demP_DD_val
  Prob_demP_RR <<- Prob_demP_RR_val
  Prob_demP_RD <<- Prob_demP_RD_val
  Prob_demP_DR <<- Prob_demP_DR_val
  
  Dem_uni_pres <<- switch.uni.pres; Rep_uni_pres <<- switch.uni.pres
  Dem_div_pres <<- switch.div.pres; Rep_div_pres <<- switch.div.pres
  Dem_uni_mid <<- switch.uni.mid; Rep_uni_mid <<- switch.uni.mid
  Dem_div_mid <<- switch.div.mid; Rep_div_mid <<- switch.div.mid
  
  Dem_Shape_One <<- Dem_Shape_One_val; Dem_Shape_Two <<- Dem_Shape_Two_val
  Rep_Shape_One <<- Rep_Shape_One_val; Rep_Shape_Two <<- Rep_Shape_Two_val
  
  Term_Limit <<- 100
  Weighting_Strategic_Retirement <<- 1
  gop_senate_advantage_parameter <<- 0
  Courtpackbyn <<- NA
  court.packing.start.year <<- start_year
  
  source("FSC_Simulator_2024.R")
  
  N <- end.year - current_year
  n_cores <- max(1, detectCores() - 1)
  local_dataset <- dataset
  local_N <- N
  
  listOfMatrices <- mclapply(1:n.sims, function(i) {
    results <- make_new_dataset(local_dataset, local_N)
    results <- core(results)
    cbind(predict_next(results), sim_number = rep(i, local_N + 1))
  }, mc.cores = n_cores, mc.set.seed = TRUE)
  
  results <- do.call(rbind, listOfMatrices)
  cat("  Done:", nrow(results), "rows\n")
  return(results)
}

#####################################################################
# Experiment 1: TRUE counterfactual — Clinton wins 2016
#####################################################################
# president_history_init: c(Bush, Obama, Obama, Obama, Obama, Clinton)
# = c(0, 1, 1, 1, 1, 1)
# Starting in 2017 with Clinton as president, Dem senate=0, divided govt
true_cf <- run_experiment(
  label = "TRUE Counterfactual (Clinton wins 2016)",
  start_year = 2017,
  pres_history = c(0, 1, 1, 1, 1, 1),  # Bush → Obama×4 → Clinton
  court_csv = "base_data_from_characteristics_for_simulations_2016_counterfactual.csv",
  seatid_vec = c(7, 9, 2, 8, 6, 1, 3, 4, 5)
)

#####################################################################
# Experiment 2: BUGGY counterfactual — real presidency, Garland court
#####################################################################
# president_history_init: c(Obama, Trump, Trump, Trump, Trump, Biden)
# = c(1, 0, 0, 0, 0, 1)
# But starting in 2017, so: history is Obama → Trump (year 1)
# Actually for 2017 start: positions map to:
# pos1 = ~2009(ish), pos2-5 = 2013-2016, pos6 = 2017
# Real world: pos1=Obama(1), pos2-5=Obama then Trump... 
# The old hardcoded was: Obama(1), Trump(0,0,0,0), Biden(1) but start was 2021
# For 2017 start: need Bush end, Obama 4 yrs, Trump yr1 = c(0,1,1,1,1,0)
buggy_cf <- run_experiment(
  label = "BUGGY Counterfactual (real presidents, Garland court)",
  start_year = 2017,
  pres_history = c(0, 1, 1, 1, 1, 0),  # Bush → Obama×4 → Trump
  court_csv = "base_data_from_characteristics_for_simulations_2016_counterfactual.csv",
  seatid_vec = c(7, 9, 2, 8, 6, 1, 3, 4, 5)
)

#####################################################################
# Load baseline for comparison
#####################################################################
load("sim_results_baseline_2024.Rdata")
bl <- baseline

#####################################################################
# ANALYSIS
#####################################################################
analyze <- function(d, label) {
  med <- d %>%
    group_by(sim_number, year) %>%
    summarise(median_ideo = median(justice_ideology), .groups="drop") %>%
    mutate(lib_med = median_ideo < 0,
           decade = paste0(floor(year/10)*10, "s"))
  
  by_decade <- med %>%
    group_by(decade) %>%
    summarise(pct_lib = round(mean(lib_med)*100, 1), .groups="drop") %>%
    filter(decade %in% c("2020s","2030s","2040s","2050s","2060s","2070s","2080s","2090s"))
  
  by_decade$label <- label
  return(by_decade)
}

true_cf_stats <- analyze(true_cf, "TRUE CF (Clinton)")
buggy_cf_stats <- analyze(buggy_cf, "BUGGY CF (Trump+Garland)")
bl_stats <- analyze(bl, "Baseline (real world)")

# Combine and pivot
all_stats <- bind_rows(true_cf_stats, buggy_cf_stats, bl_stats) %>%
  tidyr::pivot_wider(names_from = label, values_from = pct_lib)

cat("\n\n")
cat("=====================================================================\n")
cat("P(LIBERAL median) by decade\n")
cat("=====================================================================\n")
print(all_stats, n=20)

# Mirror image test
cat("\n\n=== Mirror Image Test ===\n")
cat("For 'mirror image': P(lib median | CF) ≈ P(con median | baseline) = 100 - P(lib median | baseline)\n\n")
comparison <- true_cf_stats %>%
  select(decade, pct_lib_true_cf = pct_lib) %>%
  left_join(bl_stats %>% select(decade, pct_lib_bl = pct_lib), by="decade") %>%
  mutate(
    pct_con_bl = 100 - pct_lib_bl,
    mirror_gap = abs(pct_lib_true_cf - pct_con_bl)
  )
cat("Decade | P(lib|TRUE CF) | P(con|baseline) | Gap\n")
cat("-------|----------------|-----------------|-----\n")
for(i in 1:nrow(comparison)) {
  cat(sprintf("%s  |     %5.1f%%      |     %5.1f%%       | %4.1f pp\n",
    comparison$decade[i], comparison$pct_lib_true_cf[i], 
    comparison$pct_con_bl[i], comparison$mirror_gap[i]))
}

# Also compare buggy vs true
cat("\n\n=== TRUE vs BUGGY counterfactual difference ===\n")
buggy_comparison <- true_cf_stats %>%
  select(decade, true_pct = pct_lib) %>%
  left_join(buggy_cf_stats %>% select(decade, buggy_pct = pct_lib), by="decade") %>%
  mutate(diff = true_pct - buggy_pct)
cat("Decade | TRUE CF | BUGGY CF | Diff (pp)\n")
cat("-------|---------|----------|----------\n")
for(i in 1:nrow(buggy_comparison)) {
  cat(sprintf("%s  | %5.1f%%  |  %5.1f%%  |  %+5.1f\n",
    buggy_comparison$decade[i], buggy_comparison$true_pct[i], 
    buggy_comparison$buggy_pct[i], buggy_comparison$diff[i]))
}

# Save results
save(true_cf, file="sim_results_true_counterfactual_2016.Rdata")
save(buggy_cf, file="sim_results_buggy_counterfactual_2016.Rdata")
cat("\n\nResults saved.\n")