#################################################################################################### 
## Lower Delaware River WQS Benefits Analysis - WQI WTP
####################################################################################################

####################################################################################################
# PURPOSE: Estimates benefits for WQI changes, at the Census block group (CBG) level, for the Lower 
#          Delaware River; WQI values for six different scenarios each for actual and design flows
# By: Jessica Balukas, ICF
# Last Edited: 12/04/2023
#   11/20/23: Replaced 1.7% discount rate (proposed Circular A-4 Revisions) with 2% (final Circular A-4 Revisions)
#   12/04/23: Updated WQI results, based on revised criteria: 10th percentile: 66%, 50th percentile: 74% 
#                                                             10th percentile = 5.4 mg/L, 50th percentile = 6.1 mg/L
####################################################################################################

pacman::p_load(dplyr, reshape2, stringr, data.table, pastecs, openxlsx)

#################################################################################################### 
## Set macros and working directory
####################################################################################################
Scen_RunDate <-"2023-12-04"
# Define the number of analysis periods
yrsc_30 = 30
# Define the CPI adjustment factor to convert from 2019$ to 2022$
cpi_2022 = 292.655/255.657
# Define the CPI adjustment factor to convert 2021$ (from 2021 ACS) back to 2019$ (income variable to match the metadata year)
cpi_2019 = 255.657/270.970
# define the 3% and 7% annualization factor to convert NPV to annualized values, function assumes payment at beginning of period
af_2_30 = (1/(((1-(1/(1+0.02)**(yrsc_30-1)))/0.02)+1))
af_3_30 = (1/(((1-(1/(1+0.03)**(yrsc_30-1)))/0.03)+1))
af_7_30 = (1/(((1-(1/(1+0.07)**(yrsc_30-1)))/0.07)+1))
# Define low and high WQI changes
WQch_l = 7 
WQch_h = 20
lnWQch_l = log(WQch_l)
lnWQch_h = log(WQch_h)
# Define discount year
disc_year = 2024
gamefish = 0.190 #Set based on desired run; options are 0, 1, and 0.190 (mean from the meta-data)

# Set file paths:
# input path:
input_path <- "C:/Users/40179/ICF/CWA Anniversary & Econ Support - General/06_Technical Support for Economic Analysis/TD3_Delaware River WQS/WQI/WTP/Inputs/"
# WQI input path:
WQI_path <- "C:/Users/40179/ICF/CWA Anniversary & Econ Support - General/06_Technical Support for Economic Analysis/TD3_Delaware River WQS/WQI/R_programs/Outputs/"
# GIS input path:
GIS_path <- "C:/Users/40179/ICF/CWA Anniversary & Econ Support - General/06_Technical Support for Economic Analysis/TD3_Delaware River WQS/GIS/Data/R_programs/Outputs/"
# output path:
output_path <- "C:/Users/40179/ICF/CWA Anniversary & Econ Support - General/06_Technical Support for Economic Analysis/TD3_Delaware River WQS/WQI/WTP/Outputs/"

########################################################################################## 
## Import all data files used throughout the code
##########################################################################################

# import CBG-level 2021 ACS data (population, median housheold income, average household size)
CBG_2021ACS <- read.csv(paste0(GIS_path,"CBG_in_FMA_2021_5Yr_ACS.csv"), header=TRUE)

# import CBG-level sub_proportion calculations
sub_proportion <- read.csv(paste0(GIS_path,"Sub_Prop_Table.csv"), header=TRUE)

# Import the csv files with WQI values
wqi_actual <- read.csv(paste0(WQI_path,"2019_ActualDERiver_DOWQS_WQI_PercChange_23-12-04.csv"), header = TRUE)
wqi_design <- read.csv(paste0(WQI_path,"2019_DesignDERiver_DOWQS_WQI_PercChange_23-12-04.csv"), header = TRUE)

# Import text files of model coefficients (model 1 and model 2 estimated in STATA, 2022 update)
MRM2 <- read.table(paste0(input_path,"MRM2_Wtd_coef_2022.txt"), header = TRUE)
MRM2S_log <- read.table(paste0(input_path,"MRM2S_log_coef_2022.txt"), header = TRUE)

# Import county-level, SEDAC population growth rates
county_pop_rel2021 <- read.csv(paste0(input_path,"LDR_county_pop_rel2021_2023.07.28.csv"), header=TRUE)
county_pop_rel2021$GEOID <- str_pad(county_pop_rel2021$GEOID, width=5, side="left", pad="0")

########################################################################################## 
## Prepare WQI datasets
##########################################################################################

# Isolate baseline WQI and add as a separate column
wqi_actual_baseline <- wqi_actual[ which(wqi_actual$Scenario=='Baseline'),]
wqi_design_baseline <- wqi_design[ which(wqi_design$Scenario=='Baseline'),]
colnames(wqi_actual_baseline)[colnames(wqi_actual_baseline) == 'NARS_WQI'] <- 'WQI_base'
colnames(wqi_design_baseline)[colnames(wqi_design_baseline) == 'NARS_WQI'] <- 'WQI_base'
wqi_actual$WQI_base <- wqi_actual_baseline$WQI_base
wqi_design$WQI_base <- wqi_design_baseline$WQI_base
rm(wqi_actual_baseline, wqi_design_baseline)

# Create list
dataset_list <- list(wqi_actual, wqi_design)

for (i in seq_along(dataset_list)) {
  # Create Scen_Num variable to retain final scenarios (Note: Scen_Num 5 is the scenario presented in the Economic Analysis)
  dataset_list[[i]]$Scen_Num <- ifelse(dataset_list[[i]]$Scenario=="HADO_AsIs",1,
                                       ifelse(dataset_list[[i]]$Scenario=="Sensitivity_Scenario_6",2,
                                              ifelse(dataset_list[[i]]$Scenario=="Sensitivity_Scenario_5",3,
                                                     ifelse(dataset_list[[i]]$Scenario=="EPA_Standard",4,
                                                            ifelse(dataset_list[[i]]$Scenario=="Sensitivity_Scenario_4",5,
                                                                   ifelse(dataset_list[[i]]$Scenario=="Sensitivity_Scenario_3",6,7))))))
  
  dataset_list[[i]] <- dataset_list[[i]][ which(dataset_list[[i]]$Scen_Num!=7),]
  # Remove columns for individual parameters, SPARROW-based WQI values, change, and percent change 
  dataset_list[[i]] <- dataset_list[[i]][c(-3:-12,-14:-17)]

  # Rename WQI variable to WQI_scen
  colnames(dataset_list[[i]])[colnames(dataset_list[[i]]) == 'NARS_WQI'] <- 'WQI_scen'
}

wqi_actual_final <- as.data.frame(dataset_list[[1]])
wqi_design_final <- as.data.frame(dataset_list[[2]])
rm(dataset_list, wqi_actual, wqi_design)

########################################################################################## 
## Prep and merge 2021ACS and pop_rel2021 datasets
##########################################################################################

# Subset sub_proportion to needed variables only
sub_proportion <- sub_proportion[c(-1, -3,-4)]

# Merge CBG_2021ACS and sub_proportion
CBG_2021ACS <- merge(CBG_2021ACS, sub_proportion, by=c("GEOID"), all.x=TRUE)

# Create county ID in CBG_2021ACS
CBG_2021ACS <- CBG_2021ACS %>%
  mutate(GEOID_COUNTY = substr(GEOID, 1, 5))
CBG_2021ACS$GEOID <- str_pad(CBG_2021ACS$GEOID, width=12, side="left", pad="0")

# Rename GEOID in county_pop_rel2021
colnames(county_pop_rel2021)[colnames(county_pop_rel2021) == 'GEOID'] <- 'GEOID_COUNTY'

# Merge CBG_2021ACS and county_pop_rel2021
poprel2021_CBG2021ACS <- merge(CBG_2021ACS, county_pop_rel2021, by=c("GEOID_COUNTY"), all.x=TRUE)

# Remove years that are not part of the 2024-2053 30-year analysis period (2021-2023, 2054-2100)
poprel2021_CBG2021ACS <- poprel2021_CBG2021ACS[ which(poprel2021_CBG2021ACS$Year>2023 & poprel2021_CBG2021ACS$Year<2054),]

########################################################################################## 
## Prep datasets for WTP calcs (actual, design)
##########################################################################################

# Merge with WQI datasets
wqi_actual_final_merge <- cross_join(poprel2021_CBG2021ACS, wqi_actual_final)
wqi_design_final_merge <- cross_join(poprel2021_CBG2021ACS, wqi_design_final)

# Change WQI_scen to equal WQI_base for years 2024-2028 (no benefits during technology implementation years, only during O&M years)
wqi_actual_final_merge$WQI_scen <- ifelse(wqi_actual_final_merge$Year<2029,wqi_actual_final_merge$WQI_base, wqi_actual_final_merge$WQI_scen)
wqi_design_final_merge$WQI_scen <- ifelse(wqi_design_final_merge$Year<2029,wqi_design_final_merge$WQI_base, wqi_design_final_merge$WQI_scen)

# Remove interim datasets that are no longer needed
rm(CBG_2021ACS, county_pop_rel2021, wqi_actual_final, wqi_design_final, poprel2021_CBG2021ACS, sub_proportion)

# Create list
dataset_list <- list(wqi_actual_final_merge, wqi_design_final_merge)

# Calculate number of households per year, 2024-2053
for (i in seq_along(dataset_list)) {
  dataset_list[[i]]$pop_yr <- dataset_list[[i]]$POPULATION*dataset_list[[i]]$pop_rel2021
  dataset_list[[i]]$hhn_yr <- round(dataset_list[[i]]$pop_yr/dataset_list[[i]]$HOUSEHOLD_SIZE) 
# Create variables needed for the WQ meta-analysis
  dataset_list[[i]]$lnincome1_val <- log(dataset_list[[i]]$MED_INC)*cpi_2019 # convert to 2019$, matching dollar year of meta-data
  dataset_list[[i]]$OneShotVal_val <- 1
  dataset_list[[i]]$lnyear_val <- 3.6109179 # Natural log of the maximum value of the year index value for studies in the metadata
  dataset_list[[i]]$volunt_val <- 0
  dataset_list[[i]]$tax_only_val <- 1
  dataset_list[[i]]$user_cost_val <- 0
  dataset_list[[i]]$RUM_val <- 1
  dataset_list[[i]]$non_reviewed_val <- 0
  dataset_list[[i]]$thesis_val <- 0
  dataset_list[[i]]$lump_sum_val <- 0
  dataset_list[[i]]$census_south_val <- 0
  dataset_list[[i]]$census_midwest_val <- 0
  dataset_list[[i]]$census_west_val <- 0
  dataset_list[[i]]$non_users_val <- 0
  dataset_list[[i]]$swim_use_val <- 0
  dataset_list[[i]]$gamefish_val <- gamefish
  dataset_list[[i]]$ln_ar_agr_val <- -1.721003
  dataset_list[[i]]$ln_ar_ratio1_val <- 2.427232
  colnames(dataset_list[[i]])[colnames(dataset_list[[i]]) == 'sub_prop'] <- 'sub_proportion_val'
  dataset_list[[i]]$IBI_val <- 0
  dataset_list[[i]]$lnWQch_l_val <- lnWQch_l
  dataset_list[[i]]$lnWQch_h_val <- lnWQch_h
  
  # Add midpoint of WQ change for each analytic scenario
  dataset_list[[i]]$Q_val <- (dataset_list[[i]]$WQI_base + dataset_list[[i]]$WQI_scen)/2
  dataset_list[[i]]$ln_Q_val <- log(dataset_list[[i]]$Q_val)
  
  # Add water quality change variable
  dataset_list[[i]]$WQI_ch <- dataset_list[[i]]$WQI_scen - dataset_list[[i]]$WQI_base
}

wqi_actual_dataset <- as.data.frame(dataset_list[[1]])
wqi_design_dataset <- as.data.frame(dataset_list[[2]])
remove(dataset_list)

########################################################################################## 
## ESTIMATE BENEFITS USING MEAN MODEL COEFFICIENTS
##########################################################################################

### Calculate per household WTP
dataset_list2 <- list(wqi_actual_dataset, wqi_design_dataset)

for (i in seq_along(dataset_list2)) {

  # MRM2
  dataset_list2[[i]]$resid_var_over_2_MRM2 <- (MRM2$residual**2)/2;
  dataset_list2[[i]]$WTP_hh_MRM2 <- exp(MRM2$OneShotVal*dataset_list2[[i]]$OneShotVal_val + MRM2$tax_only*dataset_list2[[i]]$tax_only_val + MRM2$user_cost*dataset_list2[[i]]$user_cost_val + MRM2$RUM*dataset_list2[[i]]$RUM_val +
                                                      MRM2$IBI*dataset_list2[[i]]$IBI_val + MRM2$thesis*dataset_list2[[i]]$thesis_val + MRM2$lnyear*dataset_list2[[i]]$lnyear_val + MRM2$volunt*dataset_list2[[i]]$volunt_val + 
                                                      MRM2$non_reviewed*dataset_list2[[i]]$non_reviewed_val + MRM2$lump_sum*dataset_list2[[i]]$lump_sum_val + MRM2$census_south*dataset_list2[[i]]$census_south_val + 
                                                      MRM2$census_midwest*dataset_list2[[i]]$census_midwest_val + MRM2$census_west*dataset_list2[[i]]$census_west_val + MRM2$nonusers*dataset_list2[[i]]$non_users_val + 
                                                      MRM2$lnincome*dataset_list2[[i]]$lnincome1_val + MRM2$swim_use*dataset_list2[[i]]$swim_use_val + MRM2$gamefish*dataset_list2[[i]]$gamefish_val + MRM2$ln_ar_agr*dataset_list2[[i]]$ln_ar_agr_val + 
                                                      MRM2$ln_ar_ratio*dataset_list2[[i]]$ln_ar_ratio1_val + MRM2$sub_proportion*dataset_list2[[i]]$sub_proportion_val + MRM2$ln_Q*dataset_list2[[i]]$ln_Q_val + MRM2$intercept + 
                                                      dataset_list2[[i]]$resid_var_over_2_MRM2) * dataset_list2[[i]]$WQI_ch
  
  # MRM2S - low change
  dataset_list2[[i]]$resid_var_over_2_MRM2S <- (MRM2S_log$residual**2)/2;
  dataset_list2[[i]]$WTP_hh_lnWQch_l <- exp(MRM2S_log$OneShotVal*dataset_list2[[i]]$OneShotVal_val + MRM2S_log$tax_only*dataset_list2[[i]]$tax_only_val + MRM2S_log$user_cost*dataset_list2[[i]]$user_cost_val + 
                                                          MRM2S_log$RUM*dataset_list2[[i]]$RUM_val + MRM2S_log$IBI*dataset_list2[[i]]$IBI_val + MRM2S_log$thesis*dataset_list2[[i]]$thesis_val + MRM2S_log$lnyear*dataset_list2[[i]]$lnyear_val + 
                                                          MRM2S_log$volunt*dataset_list2[[i]]$volunt_val + MRM2S_log$non_reviewed*dataset_list2[[i]]$non_reviewed_val + MRM2S_log$lump_sum*dataset_list2[[i]]$lump_sum_val + 
                                                          MRM2S_log$census_south*dataset_list2[[i]]$census_south_val + MRM2S_log$census_midwest*dataset_list2[[i]]$census_midwest_val + MRM2S_log$census_west*dataset_list2[[i]]$census_west_val + 
                                                          MRM2S_log$nonusers*dataset_list2[[i]]$non_users_val + MRM2S_log$lnincome*dataset_list2[[i]]$lnincome1_val + MRM2S_log$swim_use*dataset_list2[[i]]$swim_use_val + 
                                                          MRM2S_log$gamefish*dataset_list2[[i]]$gamefish_val + MRM2S_log$ln_ar_agr*dataset_list2[[i]]$ln_ar_agr_val + MRM2S_log$ln_ar_ratio*dataset_list2[[i]]$ln_ar_ratio1_val + 
                                                          MRM2S_log$sub_proportion*dataset_list2[[i]]$sub_proportion_val + MRM2S_log$ln_Q*dataset_list2[[i]]$ln_Q_val + MRM2S_log$lnquality_ch*dataset_list2[[i]]$lnWQch_l_val +
                                                          MRM2S_log$intercept + dataset_list2[[i]]$resid_var_over_2_MRM2S) * dataset_list2[[i]]$WQI_ch
  
  # MRM2S - high change
  dataset_list2[[i]]$WTP_hh_lnWQch_h <- exp(MRM2S_log$OneShotVal*dataset_list2[[i]]$OneShotVal_val + MRM2S_log$tax_only*dataset_list2[[i]]$tax_only_val + MRM2S_log$user_cost*dataset_list2[[i]]$user_cost_val + 
                                                          MRM2S_log$RUM*dataset_list2[[i]]$RUM_val + MRM2S_log$IBI*dataset_list2[[i]]$IBI_val + MRM2S_log$thesis*dataset_list2[[i]]$thesis_val + MRM2S_log$lnyear*dataset_list2[[i]]$lnyear_val + 
                                                          MRM2S_log$volunt*dataset_list2[[i]]$volunt_val + MRM2S_log$non_reviewed*dataset_list2[[i]]$non_reviewed_val + MRM2S_log$lump_sum*dataset_list2[[i]]$lump_sum_val + 
                                                          MRM2S_log$census_south*dataset_list2[[i]]$census_south_val + MRM2S_log$census_midwest*dataset_list2[[i]]$census_midwest_val + MRM2S_log$census_west*dataset_list2[[i]]$census_west_val + 
                                                          MRM2S_log$nonusers*dataset_list2[[i]]$non_users_val + MRM2S_log$lnincome*dataset_list2[[i]]$lnincome1_val + MRM2S_log$swim_use*dataset_list2[[i]]$swim_use_val + 
                                                          MRM2S_log$gamefish*dataset_list2[[i]]$gamefish_val + MRM2S_log$ln_ar_agr*dataset_list2[[i]]$ln_ar_agr_val + MRM2S_log$ln_ar_ratio*dataset_list2[[i]]$ln_ar_ratio1_val + 
                                                          MRM2S_log$sub_proportion*dataset_list2[[i]]$sub_proportion_val + MRM2S_log$ln_Q*dataset_list2[[i]]$ln_Q_val + MRM2S_log$lnquality_ch*dataset_list2[[i]]$lnWQch_h_val +
                                                          MRM2S_log$intercept + dataset_list2[[i]]$resid_var_over_2_MRM2S) * dataset_list2[[i]]$WQI_ch
  
### Calculate total WTP 
  dataset_list2[[i]]$cpi_2022 <- cpi_2022
  dataset_list2[[i]]$af_2_30 <- af_2_30
  dataset_list2[[i]]$af_3_30 <- af_3_30
  dataset_list2[[i]]$af_7_30 <- af_7_30
  dataset_list2[[i]]$disc_year <- disc_year
  dataset_list2[[i]]$WTP_hh_MRM2 <- dataset_list2[[i]]$WTP_hh_MRM2*dataset_list2[[i]]$cpi_2022
  dataset_list2[[i]]$WTP_tot_MRM2 <- dataset_list2[[i]]$WTP_hh_MRM2*dataset_list2[[i]]$hhn_yr
  dataset_list2[[i]]$WTP_tot2_MRM2 <- dataset_list2[[i]]$WTP_tot_MRM2*(1/(1.02**(dataset_list2[[i]]$Year-dataset_list2[[i]]$disc_year)))
  dataset_list2[[i]]$WTP_tot3_MRM2 <- dataset_list2[[i]]$WTP_tot_MRM2*(1/(1.03**(dataset_list2[[i]]$Year-dataset_list2[[i]]$disc_year)))
  dataset_list2[[i]]$WTP_tot7_MRM2 <- dataset_list2[[i]]$WTP_tot_MRM2*(1/(1.07**(dataset_list2[[i]]$Year-dataset_list2[[i]]$disc_year)))
  dataset_list2[[i]]$WTP_hh_lnWQch_l <- dataset_list2[[i]]$WTP_hh_lnWQch_l*dataset_list2[[i]]$cpi_2022
  dataset_list2[[i]]$WTP_tot_lnWQch_l <- dataset_list2[[i]]$WTP_hh_lnWQch_l*dataset_list2[[i]]$hhn_yr
  dataset_list2[[i]]$WTP_tot2_lnWQch_l <- dataset_list2[[i]]$WTP_tot_lnWQch_l*(1/(1.02**(dataset_list2[[i]]$Year-dataset_list2[[i]]$disc_year)))
  dataset_list2[[i]]$WTP_tot3_lnWQch_l <- dataset_list2[[i]]$WTP_tot_lnWQch_l*(1/(1.03**(dataset_list2[[i]]$Year-dataset_list2[[i]]$disc_year)))
  dataset_list2[[i]]$WTP_tot7_lnWQch_l <- dataset_list2[[i]]$WTP_tot_lnWQch_l*(1/(1.07**(dataset_list2[[i]]$Year-dataset_list2[[i]]$disc_year)))
  dataset_list2[[i]]$WTP_hh_lnWQch_h <- dataset_list2[[i]]$WTP_hh_lnWQch_h*dataset_list2[[i]]$cpi_2022
  dataset_list2[[i]]$WTP_tot_lnWQch_h <- dataset_list2[[i]]$WTP_hh_lnWQch_h*dataset_list2[[i]]$hhn_yr
  dataset_list2[[i]]$WTP_tot2_lnWQch_h <- dataset_list2[[i]]$WTP_tot_lnWQch_h*(1/(1.02**(dataset_list2[[i]]$Year-dataset_list2[[i]]$disc_year)))
  dataset_list2[[i]]$WTP_tot3_lnWQch_h <- dataset_list2[[i]]$WTP_tot_lnWQch_h*(1/(1.03**(dataset_list2[[i]]$Year-dataset_list2[[i]]$disc_year)))
  dataset_list2[[i]]$WTP_tot7_lnWQch_h <- dataset_list2[[i]]$WTP_tot_lnWQch_h*(1/(1.07**(dataset_list2[[i]]$Year-dataset_list2[[i]]$disc_year)))
}

actual_wtp <- as.data.frame(dataset_list2[[1]])
design_wtp <- as.data.frame(dataset_list2[[2]])
remove(dataset_list2)

########################################################################################## 
## Compile annualized wtp tables for actual and design flows
##########################################################################################

# Annualized Values
actual_wtp_annualized <- actual_wtp %>% 
  group_by(GEOID, Scen_Num) %>%
  summarise(hhn_avg = mean(hhn_yr),
            across(starts_with("WTP_hh_"), list(avg = mean)),
            across(starts_with("WTP_tot2_"), list(tpv2 = sum)),
            across(starts_with("WTP_tot3_"), list(tpv3 = sum)),
            across(starts_with("WTP_tot7_"), list(tpv7 = sum))) %>%
  mutate(ann2_MRM2 = WTP_tot2_MRM2_tpv2*af_2_30) %>%
  mutate(ann2_lnWQch_l = WTP_tot2_lnWQch_l_tpv2*af_2_30) %>%
  mutate(ann2_lnWQch_h = WTP_tot2_lnWQch_h_tpv2*af_2_30) %>%
  mutate(ann3_MRM2 = WTP_tot3_MRM2_tpv3*af_3_30) %>%
  mutate(ann3_lnWQch_l = WTP_tot3_lnWQch_l_tpv3*af_3_30) %>%
  mutate(ann3_lnWQch_h = WTP_tot3_lnWQch_h_tpv3*af_3_30) %>%
  mutate(ann7_MRM2 = WTP_tot7_MRM2_tpv7*af_7_30) %>%
  mutate(ann7_lnWQch_l = WTP_tot7_lnWQch_l_tpv7*af_7_30) %>%
  mutate(ann7_lnWQch_h = WTP_tot7_lnWQch_h_tpv7*af_7_30)

design_wtp_annualized <- design_wtp %>% 
  group_by(GEOID, Scen_Num) %>%
  summarise(hhn_avg = mean(hhn_yr),
            across(starts_with("WTP_hh_"), list(avg = mean)),
            across(starts_with("WTP_tot2_"), list(tpv2 = sum)),
            across(starts_with("WTP_tot3_"), list(tpv3 = sum)),
            across(starts_with("WTP_tot7_"), list(tpv7 = sum))) %>%
  mutate(ann2_MRM2 = WTP_tot2_MRM2_tpv2*af_2_30) %>%
  mutate(ann2_lnWQch_l = WTP_tot2_lnWQch_l_tpv2*af_2_30) %>%
  mutate(ann2_lnWQch_h = WTP_tot2_lnWQch_h_tpv2*af_2_30) %>%
  mutate(ann3_MRM2 = WTP_tot3_MRM2_tpv3*af_3_30) %>%
  mutate(ann3_lnWQch_l = WTP_tot3_lnWQch_l_tpv3*af_3_30) %>%
  mutate(ann3_lnWQch_h = WTP_tot3_lnWQch_h_tpv3*af_3_30) %>%
  mutate(ann7_MRM2 = WTP_tot7_MRM2_tpv7*af_7_30) %>%
  mutate(ann7_lnWQch_l = WTP_tot7_lnWQch_l_tpv7*af_7_30) %>%
  mutate(ann7_lnWQch_h = WTP_tot7_lnWQch_h_tpv7*af_7_30)

# Total Values
actual_wtp_total <- actual_wtp_annualized %>% 
  group_by(Scen_Num) %>%
  summarise(hhn_avg = sum(hhn_avg),
            across(starts_with("WTP_hh_"), list(total = mean)),
            across(starts_with("WTP_tot2_"), list(total = sum)),
            across(starts_with("WTP_tot3_"), list(total = sum)),
            across(starts_with("WTP_tot7_"), list(total = sum)),
            across(starts_with("ann2_"), list(total = sum)),
            across(starts_with("ann3_"), list(total = sum)),
            across(starts_with("ann7_"), list(total = sum)))

design_wtp_total <- design_wtp_annualized %>% 
  group_by(Scen_Num) %>%
  summarise(hhn_avg = sum(hhn_avg),
            across(starts_with("WTP_hh_"), list(total = mean)),
            across(starts_with("WTP_tot2_"), list(total = sum)),
            across(starts_with("WTP_tot3_"), list(total = sum)),
            across(starts_with("WTP_tot7_"), list(total = sum)),
            across(starts_with("ann2_"), list(total = sum)),
            across(starts_with("ann3_"), list(total = sum)),
            across(starts_with("ann7_"), list(total = sum)))

########################################################################################## 
## Export results for each model
##########################################################################################
  
require(openxlsx)
  list_of_datasets <- list("Actual_CBG" = actual_wtp_annualized, "Actual_Total" = actual_wtp_total, 
                           "Design_CBG" = design_wtp_annualized, "Design_Total" = design_wtp_total)
  if (gamefish == 0) {
  write.xlsx(list_of_datasets, file = paste0(output_path,"WQI_WTP_LowerDelawareRiver_gamefish0_",Scen_RunDate,".xlsx"))
  }
  if (gamefish == 1) {
    write.xlsx(list_of_datasets, file = paste0(output_path,"WQI_WTP_LowerDelawareRiver_gamefish1_",Scen_RunDate,".xlsx"))
  }
  if (gamefish == 0.190) {
    write.xlsx(list_of_datasets, file = paste0(output_path,"WQI_WTP_LowerDelawareRiver_gamefishMean_",Scen_RunDate,".xlsx"))
  }