Exemption Scenarios - Requested in October 25 Meeting

knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)


library(tidyverse)
library(DBI)
library(data.table)
library(ggspatial)
library(gstat)
library(here)
library(httr)
library(jsonlite)
library(ptaxsim)
library(sf)
library(stars)
library(glue)
library(ggpattern)

# Create the DB connection with the default name expected by PTAXSIM functions
ptaxsim_db_conn <- DBI::dbConnect(RSQLite::SQLite(), "./ptaxsim.db/ptaxsim-2021.0.4.db")


options(digits=4, scipen = 999)


library(httr)
library(NatParksPalettes)

# link to the API output as a JSON file
muni_shp <- read_sf("https://gis.cookcountyil.gov/traditional/rest/services/politicalBoundary/MapServer/2/query?outFields=*&where=1%3D1&f=geojson")

cook_shp <- read_sf("https://gis.cookcountyil.gov/traditional/rest/services/plss/MapServer/1/query?outFields=*&where=1%3D1&f=geojson")



#muni_shp <- read_json("muni_shp.json")
nicknames <- readxl::read_excel("./Necessary_Files/muni_shortnames.xlsx")

class_dict <- read_csv("./Necessary_Files/class_dict_expanded.csv") %>% 
  mutate(class_code = as.character(class_code))



# `agency_dt` has all taxing agencies (but not TIFs) that existed each year and includes their total taxable base (cty_cook_eav), their levy, taxing rate, binary variables for if a municipality is home rule or not, as well as many other variables. tax_bill() uses this table for the taxable EAV that is used to  calculate the tax rates in the tax bills. For simulations, you must alter the taxable EAV or levy or other variables and then tell tax_bill() function to use the modified agency data table for simulated tax bills.
# 



# has EAV values, extensions by agency_num
agency_dt <- DBI::dbGetQuery(
  ptaxsim_db_conn,
  "SELECT *
  FROM agency
  WHERE year = 2021
  "
)


# cook_agency_names <- DBI::dbGetQuery(
#   ptaxsim_db_conn,
#   "SELECT DISTINCT agency_num, agency_name
#   FROM agency_info
#   "
# )
# 
#  
# 
# 
# # has all tax codes and the taxing agency that taxes them. Tax code rates and agency rates. 
# cook_tax_codes <- DBI::dbGetQuery(
#   ptaxsim_db_conn,
#   glue_sql("
#   SELECT*
#   FROM tax_code
#   WHERE agency_num IN ({cook_agency_names$agency_num*})
#   AND year = 2021
#   ",
#   .con = ptaxsim_db_conn
#   )
# )


muni_agency_names <- DBI::dbGetQuery(
  ptaxsim_db_conn,
  "SELECT DISTINCT agency_num, agency_name, minor_type
  FROM agency_info
  WHERE minor_type = 'MUNI'
  OR agency_num = '020060000'  

  "
)

muni_tax_codes <- DBI::dbGetQuery(
  ptaxsim_db_conn,
  glue_sql("
  SELECT*
  FROM tax_code
  WHERE agency_num IN ({muni_agency_names$agency_num*})
  AND year = 2021
  ",
  .con = ptaxsim_db_conn
  )
)  
  
tax_codes <- DBI::dbGetQuery(
  ptaxsim_db_conn,
  glue_sql("
  SELECT DISTINCT tax_code_num, tax_code_rate
  FROM tax_code
  WHERE year = 2021  
  ",
  .con = ptaxsim_db_conn
  )
)

## All tax codes. 
## tax codes within municipalities have additional info 
tc_muninames <- tax_codes %>% 
  left_join(muni_tax_codes) %>%
  left_join(muni_agency_names) %>% 
  select(-agency_rate) %>% 
  left_join(nicknames) %>% 
  select(-c(minor_type, short_name, `Column1`, `Most recent reassessed`, agency_number))


# Agency number and agency name for all TIFs
TIF_agencies <- DBI::dbGetQuery(
  ptaxsim_db_conn,
  "SELECT DISTINCT agency_num, agency_name, major_type, minor_type
  FROM agency_info
  WHERE minor_type = 'TIF'
  "
)

unique_tif_taxcodes <- DBI::dbGetQuery(
  ptaxsim_db_conn, 
  glue_sql("
  SELECT DISTINCT tax_code_num
  FROM tax_code
  WHERE agency_num IN ({TIF_agencies$agency_num*})
  AND year = 2021
  ",
  .con = ptaxsim_db_conn
  )
)


tif_distrib <- DBI::dbGetQuery(
  ptaxsim_db_conn, 
  glue_sql("
  SELECT *
  FROM tif_distribution
  WHERE tax_code_num IN ({muni_tax_codes$tax_code_num*})
  AND year = 2021
  ",
  .con = ptaxsim_db_conn
  )
) %>% mutate(tax_code_num = as.character(tax_code_num))



cross_county_lines <- c("030440000", "030585000", "030890000", "030320000", "031280000","030080000", "030560000", "031120000", "030280000", "030340000","030150000","030050000", "030180000","030500000","031210000")


cross_county_lines <- muni_agency_names %>% 
  filter(agency_num %in% cross_county_lines) %>% 
  left_join(nicknames, by = "agency_name")

Changing General Homeowner Exemption - Scenarios

Create a PIN input with modified exemption amounts, then recalculate the base by taking the difference between the real and hypothetical exemptions.

Then recalculate the tax base for each district. If increasing the exemption, the base should decrease because there is less taxable EAV.

Quartiles and Progressivity

class_dict <- read_csv("./Necessary_Files/class_dict_singlefamcodes.csv") %>% 
  mutate(class_code = as.character(class_code)) # change variable type to character so the join works.

nicknames <- readxl::read_xlsx("./Necessary_Files/muni_shortnames.xlsx")

#pin_data2 <- read_csv("./Output/4C_joined_PINs_bills_and_exemptions.csv")

muni_taxrates <- read_csv("./Output/4C_muni_taxrates.csv")

#pin_data2 <- pin_data2 %>% left_join(class_dict)

muni_TC_fullyCook <- muni_tax_codes %>%
  filter(!agency_num %in% cross_county_lines)


joined_pins <- read_csv("./Output/4C_joined_PINs_bills_and_exemptions.csv") %>%
  mutate(tax_code_num = as.character(tax_code_num)) %>%  
  left_join(tc_muninames) %>% left_join(class_dict)

# all pins in munis fully within cook county that are some form of single-family, detached home
singfam_pins <- joined_pins %>% 
  filter(tax_code %in% muni_TC_fullyCook$tax_code_num) %>% # excludes county line crossing munis
  filter(Option2 == "Single-Family")

Cook County Quartiles

Cook County quartiles are calculated fromg single family properties assessed value in 2021.

q = c(.25, .5, .75)


cook_quartiles <- singfam_pins %>%
  filter(Option2 == "Single-Family") %>%
  filter(tax_code %in% muni_TC_fullyCook$tax_code_num) %>% # excludes county line crossing munis
  arrange(av) %>%
  summarize(count_pins = n(), 
            min = min(av),
            quant25 = round(quantile(av, probs = q[1])), 
            quant50 = round(quantile(av, probs = q[2])),
            quant75 = round(quantile(av, probs = q[3])),
            max = max(av))
cook_quartiles

Scenarios & tax rates

Similar to File 5_Exemption_Scenarios.rmd.

• Calculate Class 2 Burden –> Calculate the amount of taxable EAV in the Municipality (for each scenario) and multiply it by the new composite tax rate (for each scenario).

• Burden Share = Taxable EAV within Property Class * Composite tax rate

• Composite Tax Rate = (Municipal Levy / Taxable EAV )

## Bring in tax bills and exemption data for 2021 PINs ##
# 
# joined_pins <- read_csv("./Output/4C_joined_PINs_bills_and_exemptions.csv") %>%
#   mutate(tax_code_num = as.character(tax_code_num)) %>%  
#   left_join(tc_muninames) %>% left_join(class_dict)

MuniLevy <- joined_pins %>% 
  group_by(clean_name, agency_num) %>%
  
  summarize(MuniLevy = sum(final_tax_to_dist, na.rm = TRUE), # amount billed by munis with current exemptions in place
            current_nonTIF_EAV_post_exemps = sum(final_tax_to_dist/(tax_code_rate/100), na.rm = TRUE),
            current_TIF_increment_EAV = sum(final_tax_to_tif/(tax_code_rate/100), na.rm=TRUE),  
            current_Exempt_EAV = sum(tax_amt_exe/(tax_code_rate/100), na.rm=TRUE), 
            Total_EAV = sum((tax_amt_exe+final_tax_to_dist+final_tax_to_tif)/(tax_code_rate/100), na.rm = TRUE)) %>% 
  mutate(cur_muni_comp_rate = MuniLevy / current_nonTIF_EAV_post_exemps)

MuniLevy

joined_pins <- joined_pins %>% 
  mutate(exe_neg10 = 0,
         exe_0 = ifelse(eav < 10000 & exe_homeowner!=0, eav, 
                             ifelse(eav>10000 & exe_homeowner!=0, 10000, 0 )),  #would be if there is no change in exemptions
         exe_plus10 = ifelse(eav < 20000 & exe_homeowner!=0, eav, 
                             ifelse(eav>20000 & exe_homeowner!=0, 20000, 0 )),
         exe_plus20 = ifelse(eav < 30000 & exe_homeowner!=0, eav, 
                             ifelse(eav>30000 & exe_homeowner!=0, 30000, 0 ) ),
         exe_plus30 = ifelse(eav < 40000 & exe_homeowner!=0, eav, 
                             ifelse(eav>40000 & exe_homeowner!=0, 40000, 0) ),
         exe_plus40 = ifelse(eav < 50000 & exe_homeowner!=0, eav, 
                             ifelse(eav>50000 & exe_homeowner!=0, 50000, 0) ) )

scenario_calcs <- joined_pins %>%    
  group_by(clean_name) %>%

    summarize(MuniLevy = sum(final_tax_to_dist, na.rm = TRUE), # amount billed by munis with current exemptions in place
            current_nonTIF_EAV_post_exemps = sum(final_tax_to_dist/(tax_code_rate/100), na.rm = TRUE),
            current_TIF_increment_EAV = sum(final_tax_to_tif/(tax_code_rate/100), na.rm=TRUE),  
            current_Exempt_EAV = sum(tax_amt_exe/(tax_code_rate/100), na.rm=TRUE), 
            current_GHE = sum(exe_homeowner, na.rm=TRUE),
            Total_EAV = sum((tax_amt_exe+final_tax_to_dist+final_tax_to_tif)/(tax_code_rate/100), na.rm = TRUE),
            exe_neg10 = sum(exe_neg10),
            exe_0 = sum(exe_0), # no change, for comparison
            exe_plus10 = sum(exe_plus10),
            exe_plus20 = sum(exe_plus20),
            exe_plus30 = sum(exe_plus30),
            exe_plus40 = sum(exe_plus40)) %>%

  # remove all GHE (up to 10,000 EAV added back to base per PIN), 
  # add exe_homeowner back to taxable base
  mutate(neg10_taxable_eav = Total_EAV - current_TIF_increment_EAV - current_Exempt_EAV + current_GHE,    # adds GHE exempt EAV back to taxable base and decreases tax rates
         plus10_taxable_eav = Total_EAV - current_TIF_increment_EAV - current_Exempt_EAV + current_GHE - exe_plus10, # will increase tax rates
         plus20_taxable_eav = Total_EAV - current_TIF_increment_EAV - current_Exempt_EAV + current_GHE - exe_plus20,
         plus30_taxable_eav = Total_EAV - current_TIF_increment_EAV - current_Exempt_EAV + current_GHE - exe_plus30,
         plus40_taxable_eav = Total_EAV - current_TIF_increment_EAV - current_Exempt_EAV + current_GHE - exe_plus40,
         scenario_noexemptions_taxable_eav = Total_EAV - current_TIF_increment_EAV) %>%
  
  mutate(tr_neg10 = MuniLevy / neg10_taxable_eav,
         tr_nochange = MuniLevy / current_nonTIF_EAV_post_exemps,
         tr_plus10 = MuniLevy / plus10_taxable_eav,
         tr_plus20 = MuniLevy / plus20_taxable_eav,
         tr_plus30 = MuniLevy / plus30_taxable_eav,
         tr_plus40 = MuniLevy / plus40_taxable_eav, 
         tax_rate_current = MuniLevy/current_nonTIF_EAV_post_exemps,
         taxrate_noexemps = MuniLevy /(Total_EAV - current_TIF_increment_EAV  ),
         taxrate_noTIFs = MuniLevy / (Total_EAV - current_Exempt_EAV),
         taxrate_noTIFs_orExemps = MuniLevy / Total_EAV) %>%
  select(clean_name, MuniLevy, tr_neg10:taxrate_noTIFs_orExemps, everything())

write_csv(scenario_calcs, "5b_scenario_calcs.csv")


scenario_taxrates <- scenario_calcs %>%  select(clean_name, MuniLevy, tr_neg10:taxrate_noTIFs_orExemps) 

scenario_taxrates

Class 2 Burden Shift

C2_taxableEAV <- joined_pins %>%   
  filter(class >= 200 & class <= 300) %>% 
  #left_join(scenario_taxrates, by = c("tax_code" = "tax_code_num")) %>%
  group_by(clean_name) %>%

    summarize(
      C2_av = sum(av),
      C2_eav_original = sum(equalized_AV), 
      C2_DistrictRev = sum(final_tax_to_dist, na.rm=TRUE),
      C2_current_nonTIF_EAV_post_exemps = sum(final_tax_to_dist/(tax_code_rate/100), na.rm = TRUE),
      C2_current_TIF_increment_EAV = sum(final_tax_to_tif/(tax_code_rate/100), na.rm=TRUE),  
      C2_current_Exempt_EAV = sum(tax_amt_exe/(tax_code_rate/100), na.rm=TRUE), 
      C2_current_GHE = sum(exe_homeowner, na.rm=TRUE),
      C2_Total_EAV = sum((tax_amt_exe+final_tax_to_dist+final_tax_to_tif)/(tax_code_rate/100), na.rm = TRUE),
      C2_exe_neg10 = sum(exe_neg10),
      C2_exe_0 = sum(exe_0), # no change, for comparison
      C2_exe_plus10 = sum(exe_plus10),
      C2_exe_plus20 = sum(exe_plus20),
      C2_exe_plus30 = sum(exe_plus30),
      C2_exe_plus40 = sum(exe_plus40),
      C2_PC_permuni = n())  %>% 
  left_join(MuniLevy, by = "clean_name") %>%
  mutate(C2_EAV_pct = C2_eav_original / Total_EAV)



C2_burden_shift <- C2_taxableEAV %>%
  left_join(scenario_taxrates) %>%
  mutate(C2_neg10_taxableEAV = C2_Total_EAV - C2_current_TIF_increment_EAV - C2_current_Exempt_EAV + C2_current_GHE - C2_exe_neg10,
         C2_nochange = C2_Total_EAV - C2_current_TIF_increment_EAV - C2_current_Exempt_EAV,
         C2_plus10_taxableEAV = C2_Total_EAV - C2_current_TIF_increment_EAV - C2_current_Exempt_EAV + C2_current_GHE - C2_exe_plus10,
         C2_plus20_taxableEAV = C2_Total_EAV - C2_current_TIF_increment_EAV - C2_current_Exempt_EAV + C2_current_GHE - C2_exe_plus20,
        C2_plus30_taxableEAV = C2_Total_EAV - C2_current_TIF_increment_EAV - C2_current_Exempt_EAV + C2_current_GHE - C2_exe_plus30,
        C2_plus40_taxableEAV = C2_Total_EAV - C2_current_TIF_increment_EAV - C2_current_Exempt_EAV + C2_current_GHE - C2_exe_plus40
         ) %>%
  mutate(burden_C2_neg10 = (C2_neg10_taxableEAV * tr_neg10)/ MuniLevy,
         burden_C2_nochange = C2_nochange * tax_rate_current  / MuniLevy,
         burden_C2_plus10 = (C2_plus10_taxableEAV * tr_plus10) / MuniLevy,
         burden_C2_plus20 = C2_plus20_taxableEAV * tr_plus20/ MuniLevy,
         burden_C2_plus30 = C2_plus30_taxableEAV * tr_plus30/ MuniLevy,
         burden_C2_plus40 = C2_plus40_taxableEAV * tr_plus40/ MuniLevy,
         
        burden_C2_noexemps = ( (C2_Total_EAV - C2_current_TIF_increment_EAV)*taxrate_noexemps ) / MuniLevy) %>%
  select(clean_name, C2_EAV_pct, burden_C2_neg10:burden_C2_plus40, everything())

C2_burden_shift

write_csv(C2_burden_shift, "5b_Class2_burdenshift.csv")

Scenario Tax rate graphs

scenarios_long <- scenario_taxrates  %>% 
  select(clean_name, MuniLevy, tr_neg10:taxrate_noTIFs_orExemps) %>%
  pivot_longer(cols = c(tr_neg10:taxrate_noTIFs_orExemps), names_to = "GHE_Amount")


scenario_taxrates %>% 
  select(clean_name, MuniLevy, tr_neg10:taxrate_noTIFs_orExemps) %>%
  pivot_longer(cols = c(tr_neg10:taxrate_noTIFs_orExemps), names_to = "GHE_Amount") %>%
  ggplot() + 
  geom_col(aes(x=value, y = GHE_Amount))

scenario_taxrates %>% 
  filter(clean_name %in% c("Chicago", "Dolton", "Glencoe")) %>%
  select(clean_name, MuniLevy, tr_neg10:taxrate_noTIFs_orExemps) %>%
  pivot_longer(cols = c(tr_neg10:taxrate_noTIFs_orExemps), names_to = "GHE_Amount") %>%
  ggplot() + 
  geom_col(aes(x=value, y = GHE_Amount, fill = clean_name), position = "dodge")  + 
  labs(x = "Municipality Composite Tax Rate", y = "Exemption Scenarios")

Ranked Properties and Muni Quartiles

25 v 75 Percentile Homes

q = c(.25, .5, .75)

## ranks properties that are considered single family homes in order of AV for each Muni
muni_quartiles <- joined_pins %>%
  filter(Option2 == "Single-Family") %>% 
  filter(tax_code %in% muni_TC_fullyCook$tax_code_num) %>%
  group_by(agency_name, clean_name) %>%
  arrange(av) %>%
  summarize(count_pins = n(), 
            min = min(av),
            quant25 = round(quantile(av, probs = q[1])), 
            quant50 = round(quantile(av, probs = q[2])),
            quant75 = round(quantile(av, probs = q[3])),
            max = max(av)
           ) %>% 
  arrange( desc( quant50))
muni_quartiles

## create rank variable for properties that fall within the quartiles +/- $500 range
munis_ranked <- joined_pins  %>%
  inner_join(muni_quartiles, by = c("agency_name", "clean_name")) %>% 
  mutate(rank = case_when(
    av > (quant25-500) & (av<quant25+500) ~ "q25",
    av > (quant50-500) & (av<quant50+500) ~ "q50",
    av > (quant75-500) & (av<quant75+500) ~ "q75")
    ) %>%
  select(clean_name, rank, av, pin, class, everything()) %>%
  left_join(nicknames)




munis_billchange <-  munis_ranked %>% 
  group_by(clean_name, rank) %>%
  left_join(scenario_taxrates) %>%
  arrange(av) %>%
 # group_by(agency_name, has_HO_exemp) %>% 
  mutate(#taxable_eav = final_tax_to_dist / tax_code_rate,
    # current bill = current tax rate * portion of levy billed
    
    
   # ## Made negative tax bills!! ## #
         
         bill_neg10 = tr_neg10*(equalized_AV-all_exemptions+ exe_homeowner -exe_neg10),

         bill_current = cur_comp_TC_rate/100*(equalized_AV-all_exemptions),
         bill_plus10 =  tr_plus10*(equalized_AV-all_exemptions+ exe_homeowner -exe_plus10),
         bill_plus20 = tr_plus20*(equalized_AV-all_exemptions+ exe_homeowner -exe_plus20),
         bill_plus30 = tr_plus20*(equalized_AV-all_exemptions+ exe_homeowner -exe_plus30),
         bill_plus40 = tr_plus20*(equalized_AV-all_exemptions+ exe_homeowner -exe_plus40),
    
    # ## make bills $1 if they had been negative.    
         bill_neg10 = ifelse(bill_neg10 < 1, 1, bill_neg10),
         bill_current = ifelse(bill_current < 1, 1, bill_current),
         bill_plus10 = ifelse(bill_plus10 < 1, 1, bill_plus10),
         bill_plus20 = ifelse(bill_plus20 < 1, 1, bill_plus20),
         bill_plus30 = ifelse(bill_plus30 < 1, 1, bill_plus30),
         bill_plus40 = ifelse(bill_plus40 < 1, 1, bill_plus40),
         
## Prevent tax bills from having negative values  (if exemptions > eav of home)
         # bill_neg10 = ifelse(tr_neg10*(equalized_AV-all_exemptions+ exe_homeowner -exe_neg10) > 1,
         #                              tr_neg10*(equalized_AV-all_exemptions+ exe_homeowner -exe_neg10), 1),
         # 
         # bill_current = ifelse(cur_comp_TC_rate/100*(equalized_AV-all_exemptions) > 1,
         #                       cur_comp_TC_rate/100*(equalized_AV-all_exemptions), 1),
         # 
         # bill_plus10 =  ifelse(tr_plus10*(equalized_AV-all_exemptions+ exe_homeowner -exe_plus10) > 1,
         #                       tr_plus10*(equalized_AV-all_exemptions+ exe_homeowner -exe_plus10),1),
         # 
         # bill_plus20 = ifelse(tr_plus20*(equalized_AV-all_exemptions+ exe_homeowner -exe_plus20) > 1,
         #                      tr_plus20*(equalized_AV-all_exemptions+ exe_homeowner -exe_plus20), 1),
         # 
         # bill_plus30 = ifelse(tr_plus30*(equalized_AV-all_exemptions+ exe_homeowner -exe_plus30) >1, 
         #                      tr_plus30*(equalized_AV-all_exemptions+ exe_homeowner -exe_plus30), 1),
         #                      
         #                      
         # bill_plus40 = ifelse(tr_plus40*(equalized_AV-all_exemptions+ exe_homeowner -exe_plus40) > 1,
         #                      tr_plus40*(equalized_AV-all_exemptions+ exe_homeowner -exe_plus40), 1)
         )%>%
  mutate(
         zerodol_bills_ghe0 = ifelse(bill_neg10 < 5, 1, 0),
         zerodol_bills_current = ifelse(bill_current < 5, 1, 0),
         zerodol_bills_ghe20 = ifelse(bill_plus10 < 5, 1, 0),
         zerodol_bills_ghe30 = ifelse(bill_plus20 < 5, 1, 0),
         zerodol_bills_ghe40 = ifelse(bill_plus30 < 5, 1, 0),
         zerodol_bills_ghe50 = ifelse(bill_plus40 < 5, 1, 0),
  ) %>%
  group_by(clean_name, rank, has_HO_exemp) %>% 
  summarize(median_AV = round(median(av)),
            median_EAV = round(median(eav)),
            mean_bill_neg10 = round(mean(bill_neg10, na.rm=TRUE)),
            mean_bill_cur = round(mean(bill_current, na.rm=TRUE)),
            mean_bill_plus10 = round(mean(bill_plus10, na.rm=TRUE)),
            mean_bill_plus20 = round(mean(bill_plus20, na.rm=TRUE)),
            mean_bill_plus30 = round(mean(bill_plus30, na.rm=TRUE)),
            mean_bill_plus40 = round(mean(bill_plus40, na.rm=TRUE)),
            
            # current perceived_savings = median(tax_amt_exe),
            tr_neg10 = round(mean(tr_neg10*100), digits = 2), 
            cur_comp_TC_rate = round(mean(cur_comp_TC_rate), digits = 2),
            tr_plus10 = round(mean(tr_plus10*100), digits = 2),
            tr_plus20 = round(mean(tr_plus20*100), digits = 2),
            tr_plus30 = round(mean(tr_plus30*100), digits = 2),
            tr_plus40 = round(mean(tr_plus40*100), digits = 2),
            pincount=n(),
            zerodol_bills_ghe0 = sum(zerodol_bills_ghe0),           
            zerodol_bills_current = sum(zerodol_bills_current),
            zerodol_bills_ghe20 = sum(zerodol_bills_ghe20),
            zerodol_bills_ghe30 = sum(zerodol_bills_ghe30),
            zerodol_bills_ghe40 = sum(zerodol_bills_ghe40),
            zerodol_bills_ghe50 = sum(zerodol_bills_ghe50),


  ) %>%
  arrange(has_HO_exemp, rank)


munis_billchange <- munis_billchange %>% left_join(muni_quartiles)
munis_billchange

write_csv(munis_billchange, "5b_muni_billchange_scenarios.csv")

ratios<- munis_billchange %>% 
  filter(has_HO_exemp == 1 &  !is.na(rank)) %>% # claimed exemption in 2021
  mutate(currbill_to_AV_25 = ifelse(rank == "q25", mean_bill_neg10/median_AV, NA)) %>%
  mutate(currbill_to_AV_75 = ifelse(rank == "q75", mean_bill_neg10/median_AV, NA)) %>% 
  group_by(clean_name) %>%
  summarize(GHE_0_bill_to_AV_25 = max(currbill_to_AV_25, na.rm=TRUE),
            GHE_0_bill_to_AV_75 = max(currbill_to_AV_75, na.rm=TRUE)) %>%
  mutate(muni_ratio_25to75 = GHE_0_bill_to_AV_25/GHE_0_bill_to_AV_75)


ggplot(data = ratios, aes(y = GHE_0_bill_to_AV_25, x = GHE_0_bill_to_AV_75, label = clean_name)) + 
  geom_abline(intercept = 0, slope = 1) +
  geom_point(aes(alpha = .5)) + 
  geom_text(nudge_x = .03, nudge_y=0.01, size = 3, check_overlap = TRUE)+ 
  theme_classic() + 
    scale_y_continuous(limits = c(0, .6))+

  scale_x_continuous(limits = c(0, .6))+
  theme(legend.position = "none") +
  labs(title = "Reduced GHE Amount by 10,000 EAV (0 EAV exempt from GHE)",
       subtitle = "Other exemptions still in place")

ratios<- munis_billchange %>% 
  filter(has_HO_exemp == 1 &  !is.na(rank)) %>% # claimed exemption in 2021
  mutate(currbill_to_AV_25 = ifelse(rank == "q25", mean_bill_cur/median_AV, NA)) %>%
  mutate(currbill_to_AV_75 = ifelse(rank == "q75", mean_bill_cur/median_AV, NA)) %>% 
  group_by(clean_name) %>%
  summarize(currbill_to_AV_25 = max(currbill_to_AV_25, na.rm=TRUE),
            currbill_to_AV_75 = max(currbill_to_AV_75, na.rm=TRUE)) %>%
  mutate(muni_ratio_25to75 = currbill_to_AV_25/currbill_to_AV_75)


ggplot(data = ratios, aes(y = currbill_to_AV_25, x = currbill_to_AV_75, label = clean_name)) + 
  geom_abline(intercept = 0, slope = 1) +
  geom_point(aes(alpha = .5)) + 
  geom_text(nudge_x = .03, nudge_y=0.01, size = 3, check_overlap = TRUE)+ 
  theme_classic() + 
    scale_y_continuous(limits = c(0, .6))+

  scale_x_continuous(limits = c(0, .6))+
  theme(legend.position = "none") +
  labs(title = "Current GHE Amount (up to 10,000 EAV exempt per property)")

new_ratios<- munis_billchange %>% 
  filter(has_HO_exemp == 1 &  !is.na(rank)) %>% # claimed exemption in 2021
  mutate(newbill_to_AV_25 = ifelse(rank == "q25", mean_bill_plus10/median_AV, NA)) %>%
  mutate(newbill_to_AV_75 = ifelse(rank == "q75", mean_bill_plus10/median_AV, NA)) %>% 
  group_by(clean_name) %>%
  summarize(newbill_to_AV_25 = max(newbill_to_AV_25, na.rm=TRUE),
            newbill_to_AV_75 = max(newbill_to_AV_75, na.rm=TRUE)) %>%
  mutate(muni_ratio_25to75 = newbill_to_AV_25/newbill_to_AV_75)


ggplot(data = new_ratios, aes(y = newbill_to_AV_25, x = newbill_to_AV_75, label = clean_name)) + 
  geom_abline(intercept = 0, slope = 1) +
  geom_point(aes(alpha = .5)) + 
  geom_text(nudge_x = .03, nudge_y=0.01, size = 3, check_overlap = TRUE)+ 
  theme_classic() + 
  scale_y_continuous(limits = c(0, .6))+
  scale_x_continuous(limits = c(0, .6))+
  theme(legend.position = "none") + labs(title = "Increase GHE by 10,000 EAV (for up to 20,000 EAV exempt) ", 
                                         y = "25th percentile of homes, taxbill:AV",
                                         x= "75th percentile of homes, taxbill:AV",
                                         caption =  "Uses single-family properties. Exempt EAV up to 20K per property")

new_ratios<- munis_billchange %>% 
  filter(has_HO_exemp == 1 &  !is.na(rank)) %>% # claimed exemption in 2021
  mutate(newbill_to_AV_25 = ifelse(rank == "q25", mean_bill_plus20/median_AV, NA)) %>%
  mutate(newbill_to_AV_75 = ifelse(rank == "q75", mean_bill_plus20/median_AV, NA)) %>% 
  group_by(clean_name) %>%
  summarize(newbill_to_AV_25 = max(newbill_to_AV_25, na.rm=TRUE),
            newbill_to_AV_75 = max(newbill_to_AV_75, na.rm=TRUE)) %>%
  mutate(muni_ratio_25to75 = newbill_to_AV_25/newbill_to_AV_75)


ggplot(data = new_ratios, aes(y = newbill_to_AV_25, x = newbill_to_AV_75, label = clean_name)) + 
  geom_abline(intercept = 0, slope = 1) +
  geom_point(aes(alpha = .5)) + 
  geom_text(nudge_x = .03, nudge_y=0.01, size = 3, check_overlap = TRUE)+ 
  theme_classic() + 
    scale_y_continuous(limits = c(0, .6))+

  scale_x_continuous(limits = c(0, .6))+
  theme(legend.position = "none") + labs(title = "Increase GHE by 20,000 EAV (for up to 30,000 EAV exempt) ", 
                                         y = "25th percentile of homes, taxbill:AV",
                                         x= "75th percentile of homes, taxbill:AV",
                                         caption =  "Uses single-family properties. Exempt EAV up to 30K per property")

new_ratios<- munis_billchange %>% 
  filter(has_HO_exemp == 1 &  !is.na(rank)) %>% # claimed exemption in 2021
  mutate(newbill_to_AV_25 = ifelse(rank == "q25", mean_bill_plus30/median_AV, NA)) %>%
  mutate(newbill_to_AV_75 = ifelse(rank == "q75", mean_bill_plus30/median_AV, NA)) %>% 
  group_by(clean_name) %>%
  summarize(newbill_to_AV_25 = max(newbill_to_AV_25, na.rm=TRUE),
            newbill_to_AV_75 = max(newbill_to_AV_75, na.rm=TRUE)) %>%
  mutate(muni_ratio_25to75 = newbill_to_AV_25/newbill_to_AV_75)


ggplot(data = new_ratios, aes(y = newbill_to_AV_25, x = newbill_to_AV_75, label = clean_name)) + 
  geom_abline(intercept = 0, slope = 1) +
  geom_point(aes(alpha = .5)) + 
  geom_text(nudge_x = .03, nudge_y=0.01, size = 3, check_overlap = TRUE)+ 
  theme_classic() + 
    scale_y_continuous(limits = c(0, .6))+

  scale_x_continuous(limits = c(0, .6))+
  theme(legend.position = "none") + labs(title = "Increase GHE by 30,000 EAV (for up to 40,000 EAV exempt) ", 
                                         y = "25th percentile of homes, taxbill:AV",
                                         x= "75th percentile of homes, taxbill:AV",
                                         caption =  "Uses single-family properties. Exempt EAV up to 40K per property")

new_ratios<- munis_billchange %>% 
  filter(has_HO_exemp == 1 &  !is.na(rank)) %>% # claimed exemption in 2021
  mutate(newbill_to_AV_25 = ifelse(rank == "q25", mean_bill_plus40/median_AV, NA)) %>%
  mutate(newbill_to_AV_75 = ifelse(rank == "q75", mean_bill_plus40/median_AV, NA)) %>% 
  group_by(clean_name) %>%
  summarize(newbill_to_AV_25 = max(newbill_to_AV_25, na.rm=TRUE),
            newbill_to_AV_75 = max(newbill_to_AV_75, na.rm=TRUE)) %>%
  mutate(muni_ratio_25to75 = newbill_to_AV_25/newbill_to_AV_75)


ggplot(data = new_ratios, aes(y = newbill_to_AV_25, x = newbill_to_AV_75, label = clean_name)) + 
  geom_abline(intercept = 0, slope = 1) +
  geom_point(aes(alpha = .5)) + 
  geom_text(nudge_x = .03, nudge_y=0.01, size = 3, check_overlap = TRUE)+ 
  theme_classic() + 
    scale_y_continuous(limits = c(0, .6))+

  scale_x_continuous(limits = c(0, .6))+
  theme(legend.position = "none") + labs(title = "Increase GHE by 40,000 EAV (for up to 50,000 EAV exempt) ", 
                                         y = "25th percentile of homes, taxbill:AV",
                                         x= "75th percentile of homes, taxbill:AV",
                                         caption =  "Uses single-family properties. Exempt EAV up to 50K per property")

90/10 deciles

q = c(.1, .5, .9)

## ranks properties that are considered single family homes in order of AV for each Muni
muni_quartiles <- joined_pins %>%
  filter(Option2 == "Single-Family") %>% 
  filter(tax_code %in% muni_TC_fullyCook$tax_code_num) %>%
  group_by(agency_name, clean_name) %>%
  arrange(av) %>%
  summarize(count_pins = n(), 
            min = min(av),
            quant10 = round(quantile(av, probs = q[1])), 
            quant50 = round(quantile(av, probs = q[2])),
            quant90 = round(quantile(av, probs = q[3])),
            max = max(av)
           ) %>% 
  arrange( desc( quant50))
muni_quartiles

## create rank variable for properties that fall within the quartiles +/- $500 range
munis_ranked <- joined_pins  %>%
  inner_join(muni_quartiles, by = c("agency_name", "clean_name")) %>% 
  mutate(rank = case_when(
    av > (quant10-500) & (av<quant10+500) ~ "q10",
    av > (quant50-500) & (av<quant50+500) ~ "q50",
    av > (quant90-500) & (av<quant90+500) ~ "q90")
    ) %>%
  select(clean_name, rank, av, pin, class, everything()) %>%
  left_join(nicknames)




munis_billchange <-  munis_ranked %>% 
  group_by(clean_name, rank) %>%
  left_join(scenario_taxrates) %>%
  arrange(av) %>%
 # group_by(agency_name, has_HO_exemp) %>% 
  mutate(#taxable_eav = final_tax_to_dist / tax_code_rate,
    # current bill = current tax rate * portion of levy billed
         bill_neg10 = tr_neg10*(equalized_AV-all_exemptions+ exe_homeowner -exe_neg10),

         bill_current = cur_comp_TC_rate/100*(equalized_AV-all_exemptions),
         bill_plus10 =  tr_plus10*(equalized_AV-all_exemptions+ exe_homeowner -exe_plus10),
         bill_plus20 = tr_plus20*(equalized_AV-all_exemptions+ exe_homeowner -exe_plus20),
         bill_plus30 = tr_plus30*(equalized_AV-all_exemptions+ exe_homeowner -exe_plus30),
         bill_plus40 = tr_plus40*(equalized_AV-all_exemptions+ exe_homeowner -exe_plus40),

    # ## make bills $1 if they had been negative.    
         bill_neg10 = ifelse(bill_neg10 < 1, 1, bill_neg10),
         bill_current = ifelse(bill_current < 1, 1, bill_current),
         bill_plus10 = ifelse(bill_plus10 < 1, 1, bill_plus10),
         bill_plus20 = ifelse(bill_plus20 < 1, 1, bill_plus20),
         bill_plus30 = ifelse(bill_plus30 < 1, 1, bill_plus30),
         bill_plus40 = ifelse(bill_plus40 < 1, 1, bill_plus40)) %>%

  mutate(
         zerodol_bills_ghe0 = ifelse(bill_neg10 < 5, 1, 0),
         zerodol_bills_current = ifelse(bill_current < 5, 1, 0),
         zerodol_bills_ghe20 = ifelse(bill_plus10 < 5, 1, 0),
         zerodol_bills_ghe30 = ifelse(bill_plus20 < 5, 1, 0),
         zerodol_bills_ghe40 = ifelse(bill_plus30 < 5, 1, 0),
         zerodol_bills_ghe50 = ifelse(bill_plus40 < 5, 1, 0),
  ) %>%
  group_by(clean_name, rank, has_HO_exemp) %>% 
  summarize(median_AV = round(median(av)),
            median_EAV = round(median(eav)),
            mean_bill_neg10 = round(mean(bill_neg10, na.rm=TRUE)),
            mean_bill_cur = round(mean(bill_current, na.rm=TRUE)),
            mean_bill_plus10 = round(mean(bill_plus10, na.rm=TRUE)),
            mean_bill_plus20 = round(mean(bill_plus20, na.rm=TRUE)),
            mean_bill_plus30 = round(mean(bill_plus30, na.rm=TRUE)),
            mean_bill_plus40 = round(mean(bill_plus40, na.rm=TRUE)),
            
            # current perceived_savings = median(tax_amt_exe),
            tr_neg10 = round(mean(tr_neg10*100), digits = 2), 
            cur_comp_TC_rate = round(mean(cur_comp_TC_rate), digits = 2),
            tr_plus10 = round(mean(tr_plus10*100), digits = 2),
            tr_plus20 = round(mean(tr_plus20*100), digits = 2),
            tr_plus30 = round(mean(tr_plus30*100), digits = 2),
            tr_plus40 = round(mean(tr_plus40*100), digits = 2),
            pincount=n()
  ) %>%
  arrange(has_HO_exemp, rank)


munis_billchange <- munis_billchange %>% left_join(muni_quartiles)
munis_billchange

ratios<- munis_billchange %>% 
  filter(has_HO_exemp == 1 &  !is.na(rank)) %>% # claimed exemption in 2021
  mutate(currbill_to_AV_10 = ifelse(rank == "q10", mean_bill_neg10/median_AV, NA)) %>%
  mutate(currbill_to_AV_90 = ifelse(rank == "q90", mean_bill_neg10/median_AV, NA)) %>% 
  group_by(clean_name) %>%
  summarize(GHE_0_bill_to_AV_10 = max(currbill_to_AV_10, na.rm=TRUE),
            GHE_0_bill_to_AV_90 = max(currbill_to_AV_90, na.rm=TRUE)) %>%
  mutate(muni_ratio_10to90 = GHE_0_bill_to_AV_10/GHE_0_bill_to_AV_90)


ggplot(data = ratios, aes(y = GHE_0_bill_to_AV_10, x = GHE_0_bill_to_AV_90, label = clean_name)) + 
  geom_abline(intercept = 0, slope = 1) +
  geom_point(aes(alpha = .5)) + 
  geom_text(nudge_x = .03, nudge_y=0.01, size = 3, check_overlap = TRUE)+ 
  theme_classic() + 
    scale_y_continuous(limits = c(0, .6))+

  scale_x_continuous(limits = c(0, .6))+
  theme(legend.position = "none") +
  labs(title = "Reduced GHE Amount by 10,000 EAV (0 EAV exempt from GHE)",
       subtitle = "Other exemptions still in place")

ratios<- munis_billchange %>% 
  filter(has_HO_exemp == 1 &  !is.na(rank)) %>% # claimed exemption in 2021
  mutate(currbill_to_AV_10 = ifelse(rank == "q10", mean_bill_cur/median_AV, NA)) %>%
  mutate(currbill_to_AV_90 = ifelse(rank == "q90", mean_bill_cur/median_AV, NA)) %>% 
  group_by(clean_name) %>%
  summarize(currbill_to_AV_10 = max(currbill_to_AV_10, na.rm=TRUE),
            currbill_to_AV_90 = max(currbill_to_AV_90, na.rm=TRUE)) %>%
  mutate(muni_ratio_10to90 = currbill_to_AV_10/currbill_to_AV_90)


ggplot(data = ratios, aes(y = currbill_to_AV_10, x = currbill_to_AV_90, label = clean_name)) + 
  geom_abline(intercept = 0, slope = 1) +
  geom_point(aes(alpha = .5)) + 
  geom_text(nudge_x = .03, nudge_y=0.01, size = 3, check_overlap = TRUE)+ 
  theme_classic() + 
    scale_y_continuous(limits = c(0, .6))+

  scale_x_continuous(limits = c(0, .6))+
  theme(legend.position = "none") +
  labs(title = "Current GHE Amount (up to 10,000 EAV exempt per property)")

new_ratios<- munis_billchange %>% 
  filter(has_HO_exemp == 1 &  !is.na(rank)) %>% # claimed exemption in 2021
  mutate(newbill_to_AV_10 = ifelse(rank == "q10", mean_bill_plus10/median_AV, NA)) %>%
  mutate(newbill_to_AV_90 = ifelse(rank == "q90", mean_bill_plus10/median_AV, NA)) %>% 
  group_by(clean_name) %>%
  summarize(newbill_to_AV_10 = max(newbill_to_AV_10, na.rm=TRUE),
            newbill_to_AV_90 = max(newbill_to_AV_90, na.rm=TRUE)) %>%
  mutate(muni_ratio_10to90 = newbill_to_AV_10/newbill_to_AV_90)


ggplot(data = new_ratios, aes(y = newbill_to_AV_10, x = newbill_to_AV_90, label = clean_name)) + 
  geom_abline(intercept = 0, slope = 1) +
  geom_point(aes(alpha = .5)) + 
  geom_text(nudge_x = .03, nudge_y=0.01, size = 3, check_overlap = TRUE)+ 
  theme_classic() + 
  scale_y_continuous(limits = c(0, .6))+
  scale_x_continuous(limits = c(0, .6))+
  theme(legend.position = "none") + labs(title = "Increase GHE by 10,000 EAV (for up to 20,000 EAV exempt) ", 
                                         y = "10th percentile of homes, taxbill:AV",
                                         x= "90th percentile of homes, taxbill:AV",
                                         caption =  "Uses single-family properties. Exempt EAV up to 20K per property")

new_ratios<- munis_billchange %>% 
  filter(has_HO_exemp == 1 &  !is.na(rank)) %>% # claimed exemption in 2021
  mutate(newbill_to_AV_10 = ifelse(rank == "q10", mean_bill_plus20/median_AV, NA)) %>%
  mutate(newbill_to_AV_90 = ifelse(rank == "q90", mean_bill_plus20/median_AV, NA)) %>% 
  group_by(clean_name) %>%
  summarize(newbill_to_AV_10 = max(newbill_to_AV_10, na.rm=TRUE),
            newbill_to_AV_90 = max(newbill_to_AV_90, na.rm=TRUE)) %>%
  mutate(muni_ratio_10to90 = newbill_to_AV_10/newbill_to_AV_90)


ggplot(data = new_ratios, aes(y = newbill_to_AV_10, x = newbill_to_AV_90, label = clean_name)) + 
  geom_abline(intercept = 0, slope = 1) +
  geom_point(aes(alpha = .5)) + 
  geom_text(nudge_x = .03, nudge_y=0.01, size = 3, check_overlap = TRUE)+ 
  theme_classic() + 
    scale_y_continuous(limits = c(0, .6))+

  scale_x_continuous(limits = c(0, .6))+
  theme(legend.position = "none") + labs(title = "Increase GHE by 20,000 EAV (for up to 30,000 EAV exempt) ", 
                                         y = "10th percentile of homes, taxbill:AV",
                                         x= "90th percentile of homes, taxbill:AV",
                                         caption =  "Uses single-family properties. Exempt EAV up to 30K per property")

new_ratios<- munis_billchange %>% 
  filter(has_HO_exemp == 1 &  !is.na(rank)) %>% # claimed exemption in 2021
  mutate(newbill_to_AV_10 = ifelse(rank == "q10", mean_bill_plus30/median_AV, NA)) %>%
  mutate(newbill_to_AV_90 = ifelse(rank == "q90", mean_bill_plus30/median_AV, NA)) %>% 
  group_by(clean_name) %>%
  summarize(newbill_to_AV_10 = max(newbill_to_AV_10, na.rm=TRUE),
            newbill_to_AV_90 = max(newbill_to_AV_90, na.rm=TRUE)) %>%
  mutate(muni_ratio_10to90 = newbill_to_AV_10/newbill_to_AV_90)


ggplot(data = new_ratios, aes(y = newbill_to_AV_10, x = newbill_to_AV_90, label = clean_name)) + 
  geom_abline(intercept = 0, slope = 1) +
  geom_point(aes(alpha = .5)) + 
  geom_text(nudge_x = .03, nudge_y=0.01, size = 3, check_overlap = TRUE)+ 
  theme_classic() + 
    scale_y_continuous(limits = c(0, .6))+

  scale_x_continuous(limits = c(0, .6))+
  theme(legend.position = "none") + labs(title = "Increase GHE by 30,000 EAV (for up to 40,000 EAV exempt) ", 
                                         y = "10th percentile of homes, taxbill:AV",
                                         x= "90th percentile of homes, taxbill:AV",
                                         caption =  "Uses single-family properties. Exempt EAV up to 40K per property")

new_ratios<- munis_billchange %>% 
  filter(has_HO_exemp == 1 &  !is.na(rank)) %>% # claimed exemption in 2021
  mutate(newbill_to_AV_10 = ifelse(rank == "q10", mean_bill_plus40/median_AV, NA)) %>%
  mutate(newbill_to_AV_90 = ifelse(rank == "q90", mean_bill_plus40/median_AV, NA)) %>% 
  group_by(clean_name) %>%
  summarize(newbill_to_AV_10 = max(newbill_to_AV_10, na.rm=TRUE),
            newbill_to_AV_90 = max(newbill_to_AV_90, na.rm=TRUE)) %>%
  mutate(muni_ratio_10to90 = newbill_to_AV_10/newbill_to_AV_90)


ggplot(data = new_ratios, aes(y = newbill_to_AV_10, x = newbill_to_AV_90, label = clean_name)) + 
  geom_abline(intercept = 0, slope = 1) +
  geom_point(aes(alpha = .5)) + 
  geom_text(nudge_x = .03, nudge_y=0.01, size = 3, check_overlap = TRUE)+ 
  theme_classic() + 
    scale_y_continuous(limits = c(0, .6))+

  scale_x_continuous(limits = c(0, .6))+
  theme(legend.position = "none") + labs(title = "Increase GHE by 40,000 EAV (for up to 50,000 EAV exempt) ", 
                                         y = "10th percentile of homes, taxbill:AV",
                                         x= "90th percentile of homes, taxbill:AV",
                                         caption =  "Uses single-family properties. Exempt EAV up to 50K per property")

90/10 deciles - Bigger AV range

Increased range of homes included in quantile to avoid having municipalities drop out if there were no tax bills. Not sure if it addresses the problem.

q = c(.1, .5, .9)

## ranks properties that are considered single family homes in order of AV for each Muni
muni_quartiles <- joined_pins %>%
  filter(Option2 == "Single-Family") %>% 
  filter(tax_code %in% muni_TC_fullyCook$tax_code_num) %>%
  group_by(agency_name, clean_name) %>%
  arrange(av) %>%
  summarize(count_pins = n(), 
            min = min(av),
            quant10 = round(quantile(av, probs = q[1])), 
            quant50 = round(quantile(av, probs = q[2])),
            quant90 = round(quantile(av, probs = q[3])),
            max = max(av)
           ) %>% 
  arrange( desc( quant50))
muni_quartiles

## create rank variable for properties that fall within the quartiles +/- $500 range
munis_ranked <- joined_pins  %>%
  inner_join(muni_quartiles, by = c("agency_name", "clean_name")) %>% 
  mutate(rank = case_when(
    av > (quant10-1000) & (av<quant10+1000) ~ "q10",
    av > (quant50-1000) & (av<quant50+1000) ~ "q50",
    av > (quant90-1000) & (av<quant90+1000) ~ "q90")
    ) %>%
  select(clean_name, rank, av, pin, class, everything()) %>%
  left_join(nicknames)




munis_billchange <-  munis_ranked %>% 
  group_by(clean_name, rank) %>%
  left_join(scenario_taxrates) %>%
  arrange(av) %>%
  mutate(#taxable_eav = final_tax_to_dist / tax_code_rate,
    # current bill = current tax rate * portion of levy billed
         bill_neg10 = tr_neg10*(equalized_AV-all_exemptions+ exe_homeowner -exe_neg10),

         bill_current = cur_comp_TC_rate/100*(equalized_AV-all_exemptions),
         bill_plus10 =  tr_plus10*(equalized_AV-all_exemptions+ exe_homeowner -exe_plus10),
         bill_plus20 = tr_plus20*(equalized_AV-all_exemptions+ exe_homeowner -exe_plus20),
         bill_plus30 = tr_plus30*(equalized_AV-all_exemptions+ exe_homeowner -exe_plus30),
         bill_plus40 = tr_plus40*(equalized_AV-all_exemptions+ exe_homeowner -exe_plus40),

    # ## make bills $1 if they had been negative.    
         bill_neg10 = ifelse(bill_neg10 < 1, 1, bill_neg10),
         bill_current = ifelse(bill_current < 1, 1, bill_current),
         bill_plus10 = ifelse(bill_plus10 < 1, 1, bill_plus10),
         bill_plus20 = ifelse(bill_plus20 < 1, 1, bill_plus20),
         bill_plus30 = ifelse(bill_plus30 < 1, 1, bill_plus30),
         bill_plus40 = ifelse(bill_plus40 < 1, 1, bill_plus40)) %>%

  mutate(
         zerodol_bills_ghe0 = ifelse(bill_neg10 < 5, 1, 0),
         zerodol_bills_current = ifelse(bill_current < 5, 1, 0),
         zerodol_bills_ghe20 = ifelse(bill_plus10 < 5, 1, 0),
         zerodol_bills_ghe30 = ifelse(bill_plus20 < 5, 1, 0),
         zerodol_bills_ghe40 = ifelse(bill_plus30 < 5, 1, 0),
         zerodol_bills_ghe50 = ifelse(bill_plus40 < 5, 1, 0),
  ) %>%
  group_by(clean_name, rank, has_HO_exemp) %>% 
  summarize(median_AV = round(median(av)),
            median_EAV = round(median(eav)),
            mean_bill_neg10 = round(mean(bill_neg10, na.rm=TRUE)),
            mean_bill_cur = round(mean(bill_current, na.rm=TRUE)),
            mean_bill_plus10 = round(mean(bill_plus10, na.rm=TRUE)),
            mean_bill_plus20 = round(mean(bill_plus20, na.rm=TRUE)),
            mean_bill_plus30 = round(mean(bill_plus30, na.rm=TRUE)),
            mean_bill_plus40 = round(mean(bill_plus40, na.rm=TRUE)),
            
            # current perceived_savings = median(tax_amt_exe),
            tr_neg10 = round(mean(tr_neg10), digits = 2), 
            cur_comp_TC_rate = round(mean(cur_comp_TC_rate), digits = 2),
            tr_plus10 = round(mean(tr_plus10*100), digits = 2),
            tr_plus20 = round(mean(tr_plus20*100), digits = 2),
            tr_plus30 = round(mean(tr_plus30*100), digits = 2),
            tr_plus40 = round(mean(tr_plus40*100), digits = 2),
            pincount=n()
  ) %>%
  arrange(has_HO_exemp, rank)


munis_billchange <- munis_billchange %>% left_join(muni_quartiles)
munis_billchange

ratios<- munis_billchange %>% 
  filter(has_HO_exemp == 1 &  !is.na(rank)) %>% # claimed exemption in 2021
  mutate(currbill_to_AV_10 = ifelse(rank == "q10", mean_bill_neg10/median_AV, NA)) %>%
  mutate(currbill_to_AV_90 = ifelse(rank == "q90", mean_bill_neg10/median_AV, NA)) %>% 
  group_by(clean_name) %>%
  summarize(GHE_0_bill_to_AV_10 = max(currbill_to_AV_10, na.rm=TRUE),
            GHE_0_bill_to_AV_90 = max(currbill_to_AV_90, na.rm=TRUE)) %>%
  mutate(muni_ratio_10to90 = GHE_0_bill_to_AV_10/GHE_0_bill_to_AV_90)


ggplot(data = ratios, aes(y = GHE_0_bill_to_AV_10, x = GHE_0_bill_to_AV_90, label = clean_name)) + 
  geom_abline(intercept = 0, slope = 1) +
  geom_point(aes(alpha = .5)) + 
  geom_text(nudge_x = .03, nudge_y=0.01, size = 3, check_overlap = TRUE)+ 
  theme_classic() + 
    scale_y_continuous(limits = c(0, .6))+

  scale_x_continuous(limits = c(0, .6))+
  theme(legend.position = "none") +
  labs(title = "Reduced GHE Amount by 10,000 EAV (0 EAV exempt from GHE)",
       subtitle = "Other exemptions still in place")

ratios<- munis_billchange %>% 
  filter(has_HO_exemp == 1 &  !is.na(rank)) %>% # claimed exemption in 2021
  mutate(currbill_to_AV_10 = ifelse(rank == "q10", mean_bill_cur/median_AV, NA)) %>%
  mutate(currbill_to_AV_90 = ifelse(rank == "q90", mean_bill_cur/median_AV, NA)) %>% 
  group_by(clean_name) %>%
  summarize(currbill_to_AV_10 = max(currbill_to_AV_10, na.rm=TRUE),
            currbill_to_AV_90 = max(currbill_to_AV_90, na.rm=TRUE)) %>%
  mutate(muni_ratio_10to90 = currbill_to_AV_10/currbill_to_AV_90)


ggplot(data = ratios, aes(y = currbill_to_AV_10, x = currbill_to_AV_90, label = clean_name)) + 
  geom_abline(intercept = 0, slope = 1) +
  geom_point(aes(alpha = .5)) + 
  geom_text(nudge_x = .03, nudge_y=0.01, size = 3, check_overlap = TRUE)+ 
  theme_classic() + 
    scale_y_continuous(limits = c(0, .6))+

  scale_x_continuous(limits = c(0, .6))+
  theme(legend.position = "none") +
  labs(title = "Current GHE Amount (up to 10,000 EAV exempt per property)")

new_ratios<- munis_billchange %>% 
  filter(has_HO_exemp == 1 &  !is.na(rank)) %>% # claimed exemption in 2021
  mutate(newbill_to_AV_10 = ifelse(rank == "q10", mean_bill_plus10/median_AV, NA)) %>%
  mutate(newbill_to_AV_90 = ifelse(rank == "q90", mean_bill_plus10/median_AV, NA)) %>% 
  group_by(clean_name) %>%
  summarize(newbill_to_AV_10 = max(newbill_to_AV_10, na.rm=TRUE),
            newbill_to_AV_90 = max(newbill_to_AV_90, na.rm=TRUE)) %>%
  mutate(muni_ratio_10to90 = newbill_to_AV_10/newbill_to_AV_90)


ggplot(data = new_ratios, aes(y = newbill_to_AV_10, x = newbill_to_AV_90, label = clean_name)) + 
  geom_abline(intercept = 0, slope = 1) +
  geom_point(aes(alpha = .5)) + 
  geom_text(nudge_x = .03, nudge_y=0.01, size = 3, check_overlap = TRUE)+ 
  theme_classic() + 
  scale_y_continuous(limits = c(0, .6))+
  scale_x_continuous(limits = c(0, .6))+
  theme(legend.position = "none") + labs(title = "Increase GHE by 10,000 EAV (for up to 20,000 EAV exempt) ", 
                                         y = "10th percentile of homes, taxbill:AV",
                                         x= "90th percentile of homes, taxbill:AV",
                                         caption =  "Uses single-family properties. Exempt EAV up to 20K per property")

new_ratios<- munis_billchange %>% 
  filter(has_HO_exemp == 1 &  !is.na(rank)) %>% # claimed exemption in 2021
  mutate(newbill_to_AV_10 = ifelse(rank == "q10", mean_bill_plus20/median_AV, NA)) %>%
  mutate(newbill_to_AV_90 = ifelse(rank == "q90", mean_bill_plus20/median_AV, NA)) %>% 
  group_by(clean_name) %>%
  summarize(newbill_to_AV_10 = max(newbill_to_AV_10, na.rm=TRUE),
            newbill_to_AV_90 = max(newbill_to_AV_90, na.rm=TRUE)) %>%
  mutate(muni_ratio_10to90 = newbill_to_AV_10/newbill_to_AV_90)


ggplot(data = new_ratios, aes(y = newbill_to_AV_10, x = newbill_to_AV_90, label = clean_name)) + 
  geom_abline(intercept = 0, slope = 1) +
  geom_point(aes(alpha = .5)) + 
  geom_text(nudge_x = .03, nudge_y=0.01, size = 3, check_overlap = TRUE)+ 
  theme_classic() + 
    scale_y_continuous(limits = c(0, .6))+

  scale_x_continuous(limits = c(0, .6))+
  theme(legend.position = "none") + labs(title = "Increase GHE by 20,000 EAV (for up to 30,000 EAV exempt) ", 
                                         y = "10th percentile of homes, taxbill:AV",
                                         x= "90th percentile of homes, taxbill:AV",
                                         caption =  "Uses single-family properties. Exempt EAV up to 30K per property")

new_ratios<- munis_billchange %>% 
  filter(has_HO_exemp == 1 &  !is.na(rank)) %>% # claimed exemption in 2021
  mutate(newbill_to_AV_10 = ifelse(rank == "q10", mean_bill_plus30/median_AV, NA)) %>%
  mutate(newbill_to_AV_90 = ifelse(rank == "q90", mean_bill_plus30/median_AV, NA)) %>% 
  group_by(clean_name) %>%
  summarize(newbill_to_AV_10 = max(newbill_to_AV_10, na.rm=TRUE),
            newbill_to_AV_90 = max(newbill_to_AV_90, na.rm=TRUE)) %>%
  mutate(muni_ratio_10to90 = newbill_to_AV_10/newbill_to_AV_90)


ggplot(data = new_ratios, aes(y = newbill_to_AV_10, x = newbill_to_AV_90, label = clean_name)) + 
  geom_abline(intercept = 0, slope = 1) +
  geom_point(aes(alpha = .5)) + 
  geom_text(nudge_x = .03, nudge_y=0.01, size = 3, check_overlap = TRUE)+ 
  theme_classic() + 
    scale_y_continuous(limits = c(0, .6))+

  scale_x_continuous(limits = c(0, .6))+
  theme(legend.position = "none") + labs(title = "Increase GHE by 30,000 EAV (for up to 40,000 EAV exempt) ", 
                                         y = "10th percentile of homes, taxbill:AV",
                                         x= "90th percentile of homes, taxbill:AV",
                                         caption =  "Uses single-family properties. Exempt EAV up to 40K per property")

new_ratios<- munis_billchange %>% 
  filter(has_HO_exemp == 1 &  !is.na(rank)) %>% # claimed exemption in 2021
  mutate(newbill_to_AV_10 = ifelse(rank == "q10", mean_bill_plus40/median_AV, NA)) %>%
  mutate(newbill_to_AV_90 = ifelse(rank == "q90", mean_bill_plus40/median_AV, NA)) %>% 
  group_by(clean_name) %>%
  summarize(newbill_to_AV_10 = max(newbill_to_AV_10, na.rm=TRUE),
            newbill_to_AV_90 = max(newbill_to_AV_90, na.rm=TRUE)) %>%
  mutate(muni_ratio_10to90 = newbill_to_AV_10/newbill_to_AV_90)


ggplot(data = new_ratios, aes(y = newbill_to_AV_10, x = newbill_to_AV_90, label = clean_name)) + 
  geom_abline(intercept = 0, slope = 1) +
  geom_point(aes(alpha = .5)) + 
  geom_text(nudge_x = .03, nudge_y=0.01, size = 3, check_overlap = TRUE)+ 
  theme_classic() + 
    scale_y_continuous(limits = c(0, .6))+

  scale_x_continuous(limits = c(0, .6))+
  theme(legend.position = "none") + labs(title = "Increase GHE by 40,000 EAV (for up to 50,000 EAV exempt) ", 
                                         y = "10th percentile of homes, taxbill:AV",
                                         x= "90th percentile of homes, taxbill:AV",
                                         caption =  "Uses single-family properties. Exempt EAV up to 50K per property")