Synthetic Control Method

# Synthetic Control Method: California Prop 99 → cigarette sales
# Based on Abadie, Diamond & Hainmueller (2010)

library(Synth)         # Core SCM package
library(tidyverse)
library(reshape2)

set.seed(2025)

# === STEP 1: Simulate realistic panel data ===
# 39 states (California + 38 donors), 1970-2000 (31 years)

states <- c("California", paste0("State", 1:38))
years <- 1970:2000
n_states <- 39
n_years <- 31
treatment_year <- 1989  # Prop 99 implemented

# Generate baseline cigarette sales (downward trend + state effects)
data_list <- list()

for (i in 1:n_states) {
  state_effect <- rnorm(1, 120, 15)  # State-specific baseline (packs/capita)
  time_trend <- -1.2  # General decline over time
  volatility <- runif(1, 3, 6)
  
  # Pre-treatment
  sales_pre <- state_effect + time_trend * (1:(treatment_year - 1970)) + 
               rnorm(treatment_year - 1970, 0, volatility)
  
  # Post-treatment
  if (states[i] == "California") {
    # California: additional -20 pack decline due to Prop 99
    treatment_effect <- -20
    sales_post <- state_effect + time_trend * ((treatment_year - 1970 + 1):n_years) + 
                  treatment_effect + 
                  rnorm(n_years - (treatment_year - 1970), 0, volatility)
  } else {
    # Controls: continue baseline trend
    sales_post <- state_effect + time_trend * ((treatment_year - 1970 + 1):n_years) + 
                  rnorm(n_years - (treatment_year - 1970), 0, volatility)
  }
  
  sales <- c(sales_pre, sales_post)
  
  data_list[[i]] <- data.frame(
    state = states[i],
    state_num = i,
    year = years,
    sales = pmax(sales, 20),  # Floor at 20 packs
    beer = rnorm(n_years, 25, 5),  # Predictor: beer consumption
    price = seq(1.5, 3.5, length.out = n_years) + rnorm(n_years, 0, 0.2),
    income = seq(20, 40, length.out = n_years) + rnorm(n_years, 0, 2)
  )
}

data <- bind_rows(data_list)

# Add treatment indicator
data$treated <- ifelse(data$state == "California" & data$year >= treatment_year, 1, 0)

cat("=== DATA STRUCTURE ===", "\n")
cat("States:", n_states, "(1 treated, 38 donors)\n")
cat("Years:", min(years), "-", max(years), "\n")
cat("Pre-treatment periods:", treatment_year - min(years), "\n")
cat("Post-treatment periods:", max(years) - treatment_year + 1, "\n\n")

# === STEP 2: Prepare data for Synth package ===
dataprep_out <- dataprep(
  foo = as.data.frame(data),
  predictors = c("beer", "price", "income"),
  predictors.op = "mean",
  time.predictors.prior = 1970:(treatment_year - 1),
  dependent = "sales",
  unit.variable = "state_num",
  unit.names.variable = "state",
  time.variable = "year",
  treatment.identifier = 1,  # California
  controls.identifier = 2:39,
  time.optimize.ssr = 1970:(treatment_year - 1),
  time.plot = 1970:2000
)

cat("=== SYNTHETIC CONTROL OPTIMIZATION ===", "\n")

# === STEP 3: Estimate synthetic control weights ===
synth_out <- synth(dataprep_out)

cat("\nConvergence:", ifelse(synth_out$solution.w.star[1] > 0, "SUCCESS", "FAILED"), "\n")

# Extract weights
weights <- data.frame(
  state = dataprep_out$Y0names,
  weight = synth_out$solution.w
)

cat("\n=== DONOR WEIGHTS(top 5) ===", "\n")
print(head(weights %>% arrange(desc(weight)), 10))

cat("\nNumber of donors with weight > 0.01:", sum(weights$weight > 0.01), "\n")

# === STEP 4: Assess pre-treatment fit ===
cat("\n=== PRE-TREATMENT FIT ===", "\n")

# Calculate RMSPE
pre_treatment_years <- 1970:(treatment_year - 1)
treated_pre <- dataprep_out$Y1plot[as.character(pre_treatment_years)]
synthetic_pre <- dataprep_out$Y0plot %*% synth_out$solution.w

rmspe_pre <- sqrt(mean((treated_pre - synthetic_pre)^2))
cat("Pre-treatment RMSPE:", round(rmspe_pre, 3), "packs\n")

# Predictor balance
cat("\n=== PREDICTOR BALANCE ===", "\n")
synth_tables <- synth.tab(dataprep.res = dataprep_out, synth.res = synth_out)
print(synth_tables$tab.pred)

# === STEP 5: Visualize results ===
cat("\n=== GENERATING PLOTS ===", "\n")

# Main path plot
path.plot(synth.res = synth_out, dataprep.res = dataprep_out,
          Ylab = "Per-Capita Cigarette Sales(packs)",
          Xlab = "Year",
          Legend = c("California", "Synthetic California"),
          Legend.position = "topright",
          Main = "California vs Synthetic Control: Cigarette Sales")
abline(v = treatment_year, lty = 2, col = "red")

# Gap plot (treatment effect over time)
gaps <- dataprep_out$Y1plot - (dataprep_out$Y0plot %*% synth_out$solution.w)
gaps.plot(synth.res = synth_out, dataprep.res = dataprep_out,
          Ylab = "Gap in Cigarette Sales(packs)",
          Xlab = "Year",
          Main = "Gap: California - Synthetic California")
abline(v = treatment_year, lty = 2, col = "red")
abline(h = 0, lty = 2, col = "gray")

# === STEP 6: Calculate treatment effects ===
post_treatment_years <- treatment_year:2000
post_gaps <- gaps[as.character(post_treatment_years)]

cat("\n=== TREATMENT EFFECTS ===", "\n")
cat("Average post-treatment effect:", round(mean(post_gaps), 2), "packs\n")
cat("Effect in final year(2000):", round(post_gaps[length(post_gaps)], 2), "packs\n")

# === STEP 7: In-space placebo tests (permutation inference) ===
cat("\n=== IN-SPACE PLACEBO TESTS ===", "\n")
cat("Running placebo SCM for all 38 donor states(may take 30-60 seconds)...\n")

placebo_effects <- numeric(38)
placebo_rmspe_pre <- numeric(38)
placebo_rmspe_post <- numeric(38)

for (i in 1:38) {
  # Run SCM with donor i as "treated"
  tryCatch({
    dataprep_placebo <- dataprep(
      foo = as.data.frame(data),
      predictors = c("beer", "price", "income"),
      predictors.op = "mean",
      time.predictors.prior = 1970:(treatment_year - 1),
      dependent = "sales",
      unit.variable = "state_num",
      unit.names.variable = "state",
      time.variable = "year",
      treatment.identifier = i + 1,  # Placebo treated
      controls.identifier = setdiff(2:39, i + 1),
      time.optimize.ssr = 1970:(treatment_year - 1),
      time.plot = 1970:2000
    )
    
    synth_placebo <- synth(dataprep_placebo, verbose = FALSE)
    
    # Calculate placebo effect (mean post-treatment gap)
    gaps_placebo <- dataprep_placebo$Y1plot - 
                    (dataprep_placebo$Y0plot %*% synth_placebo$solution.w)
    placebo_effects[i] <- mean(gaps_placebo[as.character(post_treatment_years)])
    
    # Pre and post RMSPE
    placebo_rmspe_pre[i] <- sqrt(mean(gaps_placebo[as.character(pre_treatment_years)]^2))
    placebo_rmspe_post[i] <- sqrt(mean(gaps_placebo[as.character(post_treatment_years)]^2))
  }, error = function(e) {
    placebo_effects[i] <- NA
  })
}

# True California effect
true_effect <- mean(post_gaps)

# p-value: proportion of placebo effects as extreme as true effect
p_value_twosided <- mean(abs(placebo_effects) >= abs(true_effect), na.rm = TRUE)
p_value_onesided <- mean(placebo_effects <= true_effect, na.rm = TRUE)

cat("\nTrue California effect:", round(true_effect, 2), "packs\n")
cat("Placebo effects range:", round(range(placebo_effects, na.rm = TRUE), 2), "\n")
cat("P-value(two-sided):", round(p_value_twosided, 3), "\n")
cat("P-value(one-sided):", round(p_value_onesided, 3), "\n")
cat("Rank:", sum(placebo_effects <= true_effect, na.rm = TRUE), "out of", sum(!is.na(placebo_effects)) + 1, "\n")

# Plot placebo distribution
hist(placebo_effects, breaks = 15, col = "lightblue", 
     main = "Placebo Distribution(In-Space Tests)",
     xlab = "Average Post-Treatment Effect(packs)",
     xlim = c(min(c(placebo_effects, true_effect), na.rm = TRUE) - 5,
              max(c(placebo_effects, true_effect), na.rm = TRUE) + 5))
abline(v = true_effect, col = "red", lwd = 3, lty = 2)
text(true_effect, par("usr")[4] * 0.9, "California", col = "red", pos = 4)

# === STEP 8: Pre/Post RMSPE ratio test ===
# Exclude placebos with poor pre-fit (RMSPE ratio > 2)
rmspe_ratio_california <- sqrt(mean(post_gaps^2)) / rmspe_pre
rmspe_ratios_placebo <- placebo_rmspe_post / placebo_rmspe_pre

cat("\n=== RMSPE RATIO TEST ===", "\n")
cat("California RMSPE ratio(post/pre):", round(rmspe_ratio_california, 2), "\n")
cat("Mean placebo RMSPE ratio:", round(mean(rmspe_ratios_placebo, na.rm = TRUE), 2), "\n")

# Filter placebos with good pre-fit
good_prefit <- which(placebo_rmspe_pre < 2 * rmspe_pre)
cat("Placebos with good pre-fit:", length(good_prefit), "/", length(placebo_effects), "\n")

if (length(good_prefit) > 0) {
  p_value_filtered <- mean(rmspe_ratios_placebo[good_prefit] >= 
                           rmspe_ratio_california, na.rm = TRUE)
  cat("P-value(RMSPE ratio, filtered):", round(p_value_filtered, 3), "\n")
}

cat("\n=== APA-STYLE REPORTING ===", "\n")
cat(paste0(
  "We used the synthetic control method to estimate the causal effect of ",
  "California's Proposition 99 tobacco control program on per-capita cigarette ",
  "sales. Using 38 US states as donors, we constructed a synthetic California ",
  "that closely matched pre-intervention sales(1970-1988, RMSPE = ",
  round(rmspe_pre, 2), " packs). ",
  "The synthetic control was a weighted combination of ",
  sum(weights$weight > 0.01), " states. ",
  "Post-intervention(1989-2000), California's sales diverged substantially from ",
  "the synthetic control, with an average reduction of ",
  abs(round(true_effect, 1)), " packs per capita. ",
  "Permutation-based inference using in-space placebo tests yielded p = ",
  round(p_value_onesided, 3), ", indicating the effect was unlikely due to chance. ",
  "These findings provide strong evidence that Proposition 99 causally reduced ",
  "cigarette consumption in California."
))