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
    ) %>%
  ungroup() %>%
  mutate(
    year_since_adoption = relevel(as.factor(year_since_adoption), ref = "-1"),
    state = as.factor(state),
    year = as.factor(year)
    )

res <- felm(lnths ~ factor(year_since_adoption) | state + year | 0 | state, data = analysis_data)

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)))
  }
}

fig_data <- tibble(
  label = labels,
  coef = summary(res)$coef[var_list, "Estimate"] * 100,
  se = summary(res)$coef[var_list, "Cluster s.e."] * 100
  ) %>%
  add_row(label = -1, coef = 0, se = 0)

ggplot(fig_data, aes(x = label, y = coef)) +
  geom_point() +
  geom_ribbon(aes(ymin = coef - 1.645 * se, ymax = coef + 1.645 * se, fill = "90%"), alpha = 0.4) +
  geom_ribbon(aes(ymin = coef - 1.96 * se, ymax = coef + 1.96 * se, fill = "95%"), alpha = 0.2) +
  geom_vline(xintercept = -0.5, alpha = 0.3, linetype = "dashed", size = 0.3) +
  theme_classic() +
  geom_hline(yintercept = 0, alpha = 0.5, size = 0.5) +
  scale_fill_manual(name = "Confidence Intervals", values = c("90%" = "grey12", "95%" = "grey12")) +
  guides(fill = guide_legend(override.aes = list(alpha = c(0.4, 0.2)))) +
  ylab("Coefficient estimates & Confidence Intervals") +
  xlab("Year relative to adoption of implied contract exception") +
  theme(
    axis.text = element_text(size = 12),
    axis.title = element_text(size = 14),
    legend.title = element_text(size = 12),
    legend.text = element_text(size = 10)
  )