#===============================================================================
# 0. INSTALLATION & PACKAGE LOADING
#===============================================================================
required_packages <- c("tidyverse", "lmtest", "writexl", "car", "janitor", 
                       "fixest")

# Install required packages if necessary
new_packages <- required_packages[!(required_packages %in% installed.packages()
                                    [,"Package"])]
if(length(new_packages)) install.packages(new_packages)

# Load libraries
library(tidyverse)
library(lmtest)
library(writexl)
library(car)
library(janitor)
library(fixest)

# Set the working directory to the folder containing this script and the
# primary data file. The script reads input files and writes output files
# to this directory.
base_path <- getwd()

# ==============================================================================
# 1. DATA IMPORT AND UPDATING VARIABLE NAMES
# ==============================================================================
# Import spreadsheet data
df <- read_csv(file.path(base_path, "20260909_TreePreservation_Data_Fry.csv")) %>%
  clean_names() %>%
  rename(
    cc_11_21    = percent_percent_cc_2011_2021,
    poc         = percent_poc,
    pov         = percent_pov,
    unemp       = percent_unemp,
    child       = percent_child,
    senior      = percent_senior,
    income      = median_income_2021,
    cc_11       = x2011_cc_percent,
    mhdev_11_21 = percent_percent_m_h_dev_2011_2021,
    pop_00_10   = percent_city_pop_2000_2010,
    pop_11_21   = percent_city_pop_2011_2021
  ) %>%
  mutate(
    state_code = as.factor(state),
    code_bin1 = as.factor(code_bin1)
  )

# Drop any observations for which data is missing in core variables
reg_vars <- c("cc_11_21", "total", "poc", "pov", "unemp", "child", 
              "senior", "income", "cc_11", "mhdev_11_21", "pop_00_10", 
              "pop_11_21", "state_code")
df_initial <- df[complete.cases((df)[, reg_vars]),]

# ==============================================================================
# 2. INITIAL FILTERING (MINIMAL FILTERING)
# ==============================================================================
# Calculate 100-tile percentiles for canopy cover and development change
cc_bounds  <- quantile(df_initial$cc_11_21, probs = c(0.01, 0.99))
mhd_bounds <- quantile(df_initial$mhdev_11_21, probs = c(0.01, 0.99))

# Drop first and last percentiles
df_minimal <- df_initial %>%
  filter(cc_11_21 >= cc_bounds[1] & cc_11_21 <= cc_bounds[2]) %>%
  filter(mhdev_11_21 >= mhd_bounds[1] & mhdev_11_21 <= mhd_bounds[2])

write_csv(df_minimal, file.path(base_path, "df_minimal.csv"))

# ==============================================================================
# 3. CREATE CITY-LEVEL DATASET FROM MINIMALLY FILTERED BLOCK GROUPS
# ==============================================================================
df_city <- df_minimal %>%
  group_by(city, state_code, pop_00_10, pop_11_21, total, ordinance_time, 
           code_bin1, requirements, goals, crz, larger, fines, single, arborist,
           multi, five, stands) %>%
  summarise(
    across(
      .cols = c(poc, pov, unemp, child, senior, income, cc_11_21, cc_11, 
                mhdev_11_21),
      .fns = ~ mean(.x, na.rm = TRUE)),
    .groups = "drop")

write_csv(df_city, file.path(base_path, "df_city.csv"))

# ==============================================================================
# 4. POPULATION GROWTH FILTERING
# ==============================================================================
df_popgrowth <- df_minimal %>%
  filter((pop_00_10 <= 0.5 | is.na(pop_00_10)) & (pop_11_21 <= 0.5 | 
                                                   is.na(pop_11_21)))

write_csv(df_popgrowth, file.path(base_path, "df_popgrowth.csv"))

# ==============================================================================
# 5. TABLE 2
# ==============================================================================
create_summary_table <- function(data_frame) {
  data_frame %>%
    ungroup() %>% 
    select(where(~ is.numeric(.x) && !is.factor(.x))) %>% 
    summarise(
      across(
        .cols = everything(), 
        .fns = list(
          N    = ~ sum(!is.na(.x)),
          Mean = ~ mean(.x, na.rm = TRUE),
          SD   = ~ sd(.x, na.rm = TRUE),
          Min  = ~ min(.x, na.rm = TRUE),
          Max  = ~ max(.x, na.rm = TRUE)),
        .names = "{.col}..{.fn}")) %>%
    pivot_longer(
      cols = everything(),
      names_to = c("Variable", "Statistic"),
      names_pattern = "(.*)\\.\\.(.*)") %>%
    pivot_wider(names_from = Statistic, values_from = value) %>%
    mutate(
      across(
        .cols = c(Mean, SD, Min, Max),
        .fns = ~ format(round(.x, 4), nsmall = 4, scientific = FALSE)))}

print(create_summary_table(df_minimal), n = Inf)
print(create_summary_table(df_popgrowth), n = Inf)
print(create_summary_table(df_city), n = Inf)

#===============================================================================
# 6. FIGURE 2, BLOCK GROUP COUNTS BY ACTIONS FULFILLED
#===============================================================================
min_count <- min(df_minimal$total, na.rm = TRUE)
max_count <- max(df_minimal$total, na.rm = TRUE)

ggplot(data = df_minimal, aes(x = total))+
  geom_histogram(binwidth = 1, fill = "#f0f0f0", color = "black", alpha = 1)+
  theme_minimal(base_size = 13, base_family = "serif")+
  theme(
    panel.grid.major.x=element_blank(),
    panel.grid.minor.x=element_blank(),
    axis.line = element_line(color = "black", linewidth = 0.6),
    axis.ticks.x = element_line(color = "black", linewidth = 0.6),
    axis.ticks.length.x = unit(5, "pt"),
    axis.text.x = element_text(color = "black", 
                               vjust = -0.2, margin = margin(t = 5, b = 5)), 
    axis.text.y = element_text(color = "black"))+
  labs(x = "Number of Codes Fulfilled", y = "Total Number of Block Groups")+
  scale_x_continuous(breaks = min_count:max_count, expand = c(0,0))+
  scale_y_continuous(expand = c(0,0))

#===============================================================================
# 7. TABLE 3, PRESERVATION ACTIONS BY CITY AND BLOCK GROUP
#===============================================================================
target_codes <- c("requirements", "goals", "crz", "larger", "fines", "single",
                  "arborist", "multi", "five", "stands")

# By city
code_summary_city <- df_minimal %>%
  distinct(city, .keep_all = TRUE) %>%
  summarise(
    across(
      all_of(target_codes), 
      list(total_cities = ~sum(.x, na.rm = TRUE), 
           proportion = ~mean(.x, na.rm = TRUE)))) %>%
  pivot_longer(cols = everything(), 
               names_to = c("code_name", ".value"), 
               names_pattern = "(.*)_(total_cities|proportion)")

print(code_summary_city)

# By block group
code_summary_block <- df_minimal %>%
  summarise(
    across(
      all_of(target_codes), 
      list(total_count = ~sum(.x, na.rm = TRUE), 
           proportion = ~mean(.x, na.rm = TRUE)))) %>%
  pivot_longer(cols = everything(), 
               names_to = c("code_name", ".value"), 
               names_pattern = "(.*)_(total_count|proportion)")

print(code_summary_block)

#===============================================================================
# 8. TABLE 5, CORRELATION MATRIX
#===============================================================================
corr_vars <- c("cc_11_21", "total", "ordinance_time", "mhdev_11_21", "cc_11", 
               "pop_11_21", "poc", "pov", 
               "unemp", "child", "senior", "income")
corr_matrix <- cor(df_minimal[, corr_vars], use = "complete.obs")
round(corr_matrix, 2)

corr_df <- as.data.frame(corr_matrix) %>% 
  mutate(Variable = rownames(corr_matrix)) %>%
  select(Variable, everything())

write_xlsx(corr_df, file.path(base_path, "corr_matrix.xlsx"))

# ==============================================================================
# 9. TABLE 4, REGRESSION MODELS
# ==============================================================================
# Model 1: Baseline with time variable, standard errors clustered by city
m1 <- feols(cc_11_21 ~ total + poc + pov + unemp + child + senior + income + 
              cc_11 + mhdev_11_21 + pop_11_21 + 
              ordinance_time, data = df_minimal)
summary(m1, vcov = ~city)

# Model 2: With interactions and fixed effects
m2 <- feols(cc_11_21 ~ total + ordinance_time + poc + pov + unemp + child + 
              senior + income + cc_11*mhdev_11_21 | state_code, 
            data = df_minimal)
summary(m2, vcov = ~city)

# Model 2.5: Pooled OLS approximation of 2 for vifs
m25 <- feols(cc_11_21 ~ total + ordinance_time + poc + pov + unemp + child + 
               senior + income + cc_11*mhdev_11_21 + state_code, 
             data = df_minimal)
vif(m25)

# Model 3: City-level
m3 <- feols(cc_11_21 ~ total + pov + unemp + pop_11_21 + ordinance_time, 
            data = df_city)
summary(m3, vcov = "hetero")

# Model 4: Same as 2 but with no high-growth cities
m4 <- feols(cc_11_21 ~ total + ordinance_time + poc + pov + unemp + child + 
              senior + income + cc_11*mhdev_11_21 | state_code, 
            data = df_popgrowth)
summary(m4, vcov = ~city)

# Model 5: Same as 2 but with code bins
m5 <- feols(cc_11_21 ~ code_bin1 + poc + pov + unemp + child + senior + 
              income + cc_11*mhdev_11_21 | state_code, data = df_minimal)
summary(m5, vcov = ~city)

# Model 5.5, 5 estimated for VIFS
m55 <- feols(cc_11_21 ~ code_bin1 + poc + pov + unemp + child + senior + 
               income + cc_11*mhdev_11_21 + state_code, data = df_minimal)
vif(m55)

# Model 6, Same as model 2 but with code types
m6 <- feols(cc_11_21 ~ construction + procedure + poc + pov + unemp + child + 
              senior + income + cc_11*mhdev_11_21 | state_code, 
            data = df_minimal)
summary(m6, vcov = ~city)

# Model 6.5, 6 estimated for VIFS
m65 <- feols(cc_11_21 ~ construction + procedure + poc + pov + unemp + child + 
               senior + income + cc_11*mhdev_11_21 + state_code, 
             data = df_minimal)
vif(m65)

#===============================================================================
# 10. COUNT OF CITIES BY BIN AND BINARY VARIABLES
#===============================================================================
bin_counts <- df_minimal %>%
  group_by(code_bin1) %>%
  summarize(cities_in_bin = n_distinct(city), .groups = "drop")
print(bin_counts)

binary_city_counts <- df_minimal %>%
  distinct(city, .keep_all = TRUE) %>%
  summarise(
    cities_with_construction_1 = sum(construction == 1, na.rm = TRUE),
    cities_with_procedure_1    = sum(procedure == 1, na.rm = TRUE))
print(binary_city_counts)