# '
# The file "table_multiple bandwidths.r" takes the panel data for estimation 
# from "bld/data/" and estimates the treatment effect on the hospitalization and
# death rate of different diseases, and on the potential mechanism at various bandwidths as a robustness check.
# The exact variables are specified in "src/code/specifications.json". This creates Supplementary
# Tables SI7-SI9 of the paper, respectively. 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

# 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 regressions
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"))
outcome_list <- model$outcomes
mechanisms_list <- model$mechanisms[-grep("_over_limit_relative", names(model$mechanisms))]

# Loop over all bandwidths and store result of each bandwidth in results
results <- list()
i <- 1
for (bw in c(25000, 75000, 150000, 250000, 10000000)) {
  
  # Load analysis data
  analysis <- readRDS(paste0("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 <= bw & dat$year > 2002 )

  # Treatment effects for Mechanism outcomes at multiple bandwidths
  allModels = lapply(names(mechanisms_list[c(2:length(mechanisms_list))]), function(outcome){
    if(outcome == "fires"){
      d$timeTreated <- d$time * d$treated.fires
      frm <- as.formula(paste0("asinh(", outcome, ") ~ ", paste(names(model$controls)[-1], collapse = "+"),
                               "+ time*factor(segments)+ time*mun_lon*mun_lat + treated.fires + timeTreated | mcode + year"))
      feols(frm, data = d)#, vcov = "twoway")
    }
    else{
      d$timeTreated <- d$time * d$treated
      frm <- as.formula(paste0("asinh(", outcome, ") ~ ", 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 = TRUE,
      placement = "h",
      title = "Mechanism variables",
      keep = "%timeTreated",
      se.below = TRUE,
      headers =  sapply(2:length(mechanisms_list), function(i) paste(mechanisms_list[[i]])),
      dict = c(timeTreated = paste0(bw/1000, " km")),
      label = paste0("table:", mechanisms_list$title),
      depvar = FALSE,
      notes = "Will be replaced",
      drop.section = "fixef",
      coef.just = "c",
      style.tex = style.tex("qje"),
      digits = "r3",
      digits.stats = "r3"
      
    )
  
  # Store regression results of each bandwidth
  results[[i]] <- c(star[10:11], star[13])
  i <- i + 1
}
# Append all effects from all bandwidths to one large table
brk <- "& & & & & & \\\\"
table_out <- c(star[1:9], results[[1]], brk, results[[2]], brk, results[[3]], brk, results[[4]], brk, star[10:11], star[13], star[16:length(star)]) #

# Add individual note to each table
note.latex <- "[-1.8ex]  \\end{tabular}
  \\begin{tablenotes}[para,flushleft]
    \\textit{Note:} Coefficients are given as decimal percentage 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
table_out[grepl("end{tabular}", table_out, fixed=TRUE)] <- note.latex
table_out <- sub("10000 km", "Full sample", table_out)
table_out <- sub("[h]", "", table_out, fixed = TRUE)
sink(paste0("bld/tables/multiple_bw_", mechanisms_list[[1]], ".tex"))
sink(cat(table_out[5:(length(table_out)-4)], sep = "\n"))

# Treatment effects for health outcomes at multiple bandwidths
for (outcome in c("cases_prevalence", "death_prevalence")){
  # different bandwidths
  results <- list()
  i <- 1
  for (bw in c(25000, 75000, 150000, 250000, 10000000)) {#10000

    # Load analysis data
    analysis <- readRDS(paste0("bld/data/analysis_year.RDS"))

    # Drop duplicates and keep only obs at 100km bandwidth
    dat <- subset(analysis, analysis$dist2cutoff <= bw & analysis$year > 2002)

    # Loop over hospitalization and death rate variables for estimation.
    # Loop over each variable and create model for each
    # Overall treatment effect
    allModels = lapply(seq_along(outcome_list)[-1], function(x){
      dat$timeTreated <- dat$time * dat$treated
      frm <- as.formula(paste0(outcome, " ~ ", paste(names(model$controls)[-1], collapse = "+"),
                               "+ timeTreated + factor(segments)*time + time*mun_lon*mun_lat| mcode + year"))
      feols(frm, data = subset(dat, dat$icd == names(outcome_list[x])) %>%
              distinct(mcode, year, .keep_all = TRUE))#, vcov = "twoway")
    })

    # Create regression table
    star <- 
      etable(
        allModels,
        tex = TRUE,
        placement = "h",
        title = "Outcome variables",
        keep = "%timeTreated",
        se.below = TRUE,
        headers =  sapply(1:length(compact(allModels)), function(i) paste(outcome_list[[i+1]])),
        dict = c(timeTreated = paste0(bw/1000, " km")),
        label = paste0("table:", outcome),
        depvar = FALSE,
        notes = "Will be replaced",
        drop.section = "fixef",
        coef.just = "c",
        style.tex = style.tex("qje"),
        digits = "r3",
        digits.stats = "r3"
      )

    # Store regression results of each bandwidth
    results[[i]] <- c(star[10:11], star[13])
    i <- i + 1
  }

  # Append all effects from all bandwidths to one large table
  table_out <- c(star[1:9], results[[1]], brk, results[[2]], brk, results[[3]], brk, results[[4]], brk, star[10:11], star[13], star[16:length(star)]) #
  
  # Add individual note to each table
  if (outcome == "cases_prevalence"){
    note.latex <- "[-1.8ex] \\end{tabular}
      \\begin{tablenotes}[para,flushleft]
        \\textit{Note:} The dependent variables are hospitalizations per 1,000 inhabitants associated with the named health problems
      according to the ICD-10 classification.
        All regressions use year and municipality fixed effects.
      Time variant control variables are discussed in the Methods section.
      Standard errors (in parentheses) are clustered at the municipality level. \\\\
      $^{*}$p$ < $0.1; $^{**}$p$ < $0.05; $^{***}$p$ < $0.01.\\\\
      \\end{tablenotes}"

    var.label.latex <- "\\tabularnewline\\midrule\\midrule\\\\[-1.8ex]
     & \\multicolumn{6}{c}{Hospitalizations per 1,000} \\\\"
  } else {
    note.latex <- "[-1.8ex] \\end{tabular}
      \\begin{tablenotes}[para,flushleft]
        \\textit{Note:} The dependent variables are deaths per 1,000 inhabitants associated with the named health problems
      according to the ICD-10 classification.
        All regressions use year and municipality fixed effects.
      Time variant control variables are discussed in the Methods section.
      Standard errors (in parentheses) are clustered at the municipality level. \\\\
      $^{*}$p$ < $0.1; $^{**}$p$ < $0.05; $^{***}$p$ < $0.01.\\\\
      \\end{tablenotes}"

    var.label.latex <- "\\tabularnewline\\midrule\\midrule\\\\[-1.8ex]
    & \\multicolumn{6}{c}{Deaths per 1,000} \\\\"
  }
  
  # Add note to the table and store as a tex file
  table_out[grepl("end{tabular}", table_out, fixed=TRUE)]<- note.latex
  table_out[grepl("\\tabularnewline\\midrule\\midrule", table_out, fixed=TRUE)] <- var.label.latex
  table_out <- sub("10000 km", "Full sample", table_out)
  table_out <- sub("[h]", "", table_out, fixed = TRUE)
  sink(paste0("bld/tables/multiple_bw_", outcome, ".tex"))
  cat(table_out[5:(length(table_out)-4)], sep = "\n")
  sink()
}