did2s Falsification Test

did2s_falsification_test.R

library(tidyverse)
library(haven)

#' From an indicator for treatment, retrieve start year.
get_min_year = function(y, t) {
	return(y[order(y)][min(which(t[order(y)] == 1))])
}

data <- haven::read_dta("nlswork.dta")

data <- data %>%
	arrange(idcode, year) %>%
	group_by(idcode) %>%
	mutate(
		union = if_else(is.na(union), lag(union, n=1L), union),
		union = if_else(lag(union, n=1L) == 1, 1, 0),
		flag = sum(!is.na(union))
	) %>%
	filter(flag > 0) %>%
	mutate(
		ever_treated = (sum(union, na.rm=T) > 0)
	) %>%
	filter(ever_treated == 0) %>%
	ungroup()

coefs <- NULL

for(i in 1:1000){
	set.seed(i)
	print(glue::glue("Seed: {i}"))

	sample_data <- data %>% mutate(
			union = 0,
			union = runif(n()) < 0.05 & year >= 70,
			union = 1 * union
		) %>%
		group_by(idcode) %>%
		mutate(
			gvar = get_min_year(year, union),
			gvar = if_else(is.na(gvar), 0, gvar),
			union = if_else(year >= gvar & gvar > 0, 1, 0)
		) %>%
		ungroup() %>%
		mutate(
			rel_year = if_else(gvar > 0, year - gvar, 0)
		)


	first_stage = fixest::feols(ln_wage ~ 0 | year + idcode, data = sample_data %>% filter(union == 0), warn=FALSE, notes=FALSE)

	sample_data$adj = sample_data$ln_wage - stats::predict(first_stage, newdata = sample_data)

	second_stage = fixest::feols(adj ~ 0 + i(rel_year, ref = 0), data = sample_data, warn=FALSE, notes=FALSE)


	coefs <- broom::tidy(second_stage) %>%
		mutate(seed = i) %>%
		bind_rows(coefs, .)

}

mean_coefs <- coefs %>%
	group_by(term) %>%
	summarize(
		sd_estimate = sd(estimate, na.rm = TRUE),
		estimate = mean(estimate, na.rm = TRUE),
	) %>%
	mutate(
		rel_year = as.numeric(stringr::str_extract(term, "(?<=::).*")),
		group = case_when(
			rel_year < 0 ~ "Pre",
			TRUE ~ "Post"
		),
		ci_lower = estimate - 1.96 * sd_estimate,
		ci_upper = estimate + 1.96 * sd_estimate,
	)


# Plot results 

source("https://raw.githubusercontent.com/kylebutts/templates/master/ggplot_theme/theme_kyle.R")

ggplot(mean_coefs, mapping = aes(x = rel_year, y = estimate, fill = group, color = group)) +
	geom_ribbon(mapping = aes(ymin = ci_lower, ymax = ci_upper), alpha = 0.2) +
	geom_point() +
	geom_line() +
	geom_vline(xintercept = 0, linetype = "dashed") + 
	theme_kyle() + 
	scale_color_manual(values = c("steelblue", "#E04836")) +
	scale_fill_manual(values = c("steelblue", "#E04836")) +
	guides(color = FALSE, fill = FALSE) +
	labs(y = "Estimate", x = "Year Relative to Treatment")

ggsave("results.jpg", width = 8, height = 4, dpi = 300)