# '
# The file "table_acute_exposure.r" takes the panel data for estimation
# from "bld/data/" and estimates the reduction in exposure to short-term elevated
# air pollution levels. The exact variables
# are specified in "src/code/specifications.json". This creates Supplementary
# Table SI4 of the paper. Tables are stored as tex files to "bld/tables".
# 
# '

rm(list = ls())
options(scipen = 999) # Turn off scientific notation
options(warn = 0) # 0 to turn on again, -1 to turn off
source(paste0("src/code/ggplot_theme_publication.r")) # Load a user-written ggplot theme and style.

# Load libraries
library(tidyverse) # More intuitive data wrangling
library(stargazer) # easy way to make model output look more appealing (R-inline, html, or latex)
library(fixest) # Fast two-way fixed-effects estimation
library(rjson) # Read JSON files in R

# Load variables for which the treatment effect shall be estimated from the model specification.
model <- fromJSON(file = paste0("src/code/specifications.json"))
list <- model$mechanisms[grep("_over_limit_relative", names(model$mechanisms))]

# Load the data for estimation.
analysis <- readRDS("bld/data/analysis_year.RDS")

# Drop duplicates and keep only obs at 100km bandwidth
dat <- analysis %>%
  # distinct(mcode, year, .keep_all = TRUE) %>%
  filter(icd == "J") %>%
  arrange(mcode, year)
d <- subset(dat, dat$dist2cutoff <= model$distance_to_border & dat$year > 2002) 

# Estimate acute exposure treatment effects
allModels = lapply(names(list[c(1:length(list))]), function(outcome){
  d$timeTreated <- d$time * d$treated
  frm <- as.formula(paste0(outcome, "/100 ~ ", paste(names(model$controls)[-1], collapse = "+"),
                           "+ time*factor(segments)+ time*mun_lon*mun_lat + timeTreated | mcode + year"))
  feols(frm, data = d %>% distinct(mcode, year, .keep_all = TRUE))#, vcov = "twoway")
  }
)

# Create regression table
star <- 
  etable(
    allModels,
    tex = T,
    placement = "h",
    title = "Mechanism variables",
    keep = "%timeTreated",
    se.below = TRUE,
    headers =  sapply(1:length(list), function(i) paste(list[[i]])),
    dict = c(timeTreated = "Treatment effect", mcode = "Municipality", year = "Year"),
    label = paste0("table:", list$title),
    depvar = FALSE,
    notes = "Will be replaced",
    fixef_sizes = TRUE,
    # drop.section = "fixef",
    # coef.just = "c",
    style.tex = style.tex("qje"),
    digits = "r3",
    digits.stats = "r3"
)

# Add individual note to each table
note.latex <- "[-1.8ex] \\end{tabular}
      \\begin{tablenotes}[para,flushleft]
        \\textit{Note:} Coefficients are given as percentage point changes.
        All regressions use year and municipality fixed effects.
        Time variant control variables are discussed in the Methods section.
        PM\\textsubscript{2.5} is particulate matter with a diameter of less than 2.5 $\\mu m$; $CO$ is carbon monoxide;
        NO\\textsubscript{2} is nitrogen dioxide; O\\textsubscript{3} is ozone; and SO\\textsubscript{2} is sulfur dioxide.
        Standard errors (in parentheses) are clustered at the municipality level.\\\\
        $^{*}$p$ < $0.1; $^{**}$p$ < $0.05; $^{***}$p$ < $0.01. \\\\
      \\end{tablenotes}"

# Add note to the table and store as a tex file
star[grepl("end{tabular}", star, fixed=TRUE)] <- note.latex
star <- star[c(1:14,16:length(star))]
sink(paste0("bld/tables/acute_exposure.tex"))
sink(cat(star[5:(length(star)-4)], sep = "\n"))