R Econ Visual Library

# Install and load packages ---------------
packages <- c(
  "tidyverse",
  "tidymodels",
  "haven",
  "lfe"
)

# Change to install = TRUE to install the required packages
pacman::p_load(packages, character.only = TRUE, install = FALSE)

# Load an example dataset ---------------
# Autor (2003): "Outsourcing at Will: The Contribution of Unjust Dismissal Doctrine to the Growth of Employment Outsourcing"
data <- read_dta("https://github.com/worldbank/r-econ-visual-library/raw/master/Library/Data/autor-jole-2003.dta")

admico_list <- c("admico_2", "admico_1", "admico0", "admico1", "admico2", "admico3", "mico4")

analysis_data <- data %>%
  filter(year >= 79, year <= 95, state != 98) %>%
  select(state, year, annemp, lnths, all_of(admico_list)) %>%
  group_by(state) %>%
  mutate(
    adoption_first_year = min(ifelse(admico0 == 1, year, Inf)),
    adoption = ifelse(year >= min(ifelse(admico0 == 1, year, Inf)), 1, 0),
    year_since_adoption = year - adoption_first_year,
    after = ifelse(year_since_adoption > 0, TRUE, FALSE)
    ) %>%
  ungroup() %>%
  mutate(
    year_since_adoption = relevel(as.factor(year_since_adoption), ref = "-1"),
    state = as.factor(state),
    year = as.factor(year)
    )

# Regression results
res_es <- felm(lnths ~ factor(year_since_adoption) | state + year | 0 | state, data = analysis_data)
res_dind <- felm(lnths ~ after | state + year | 0 | state, data = analysis_data)

# Variable label
labels <- c()
var_list <- c()
for (num in seq(-9, 15)){
  if (num != -1){
    labels <- c(labels, num)
    var_list <- c(var_list, paste0("factor(year_since_adoption)", as.character(num)))
  }
}

# Figure data
fig_data <- tibble(
  label = labels,
  es_coef = summary(res_es)$coef[var_list, "Estimate"] * 100,
  es_se = summary(res_es)$coef[var_list, "Cluster s.e."] * 100,
  dind_coef = summary(res_dind)$coef["afterTRUE", "Estimate"] * 100,
  dind_se = summary(res_dind)$coef["afterTRUE", "Cluster s.e."] * 100
  ) %>%
  add_row(label = -1, es_coef = 0, es_se = 0) %>%
  mutate(
    dind_coef = ifelse(label >= 0, dind_coef, 0),
    dind_se = ifelse(label >= 0, dind_se, 0)
    )

# Figure
ggplot(fig_data, aes(x = label, y = es_coef)) +
  geom_pointrange(aes(ymin = es_coef - 1.96 * es_se, ymax = es_coef + 1.96 * es_se), alpha = 0.7) +
  geom_vline(xintercept = -0.5, alpha = 0.3, linetype = "dashed", size = 0.3) +
  geom_line(data = fig_data, aes(x = label, y = dind_coef, group = (label >= 0)), color = "red", alpha = 0.7) +
  geom_line(data = fig_data, aes(x = label, y = dind_coef + 1.96 * dind_se, group = (label >= 0)), linetype = "dashed", color = "red", alpha = 0.7) +
  geom_line(data = fig_data, aes(x = label, y = dind_coef - 1.96 * dind_se, group = (label >= 0)), linetype = "dashed", color = "red", alpha = 0.7) +
  theme_classic() +
  geom_hline(yintercept = 0, alpha = 0.5, size = 0.5) +
  ylab("Coefficient estimates & 95% CI") +
  xlab("Year relative to adoption of implied contract exception") +
  theme(
    axis.text = element_text(size = 12),
    axis.title = element_text(size = 14)
  )