One way to deal with our ignorance is to perform a sensitivity analysis: how sensitive is the optimal choice to our assumptions? Below we use R to graphically depict the sensitivity of CBA to both the discount rate and time horizon: we are going to assume that the future values are known with certainty. These costs and benefits are purely hypothetical, but I have chosen them to reflect some stylized facts about sewage treatment:
We are going to compare three levels of treatment: Secondary<Tertiary<Enhanced. Note that, prior to December 15, 2020, Victoria had only primary treatment (screening) of its sewage, but primary treatment was no longer considered to be an option. To perform a sensitivity analysis in R we would start by using the crossing() function to create a dataframe that contains all the possible combinations of the variables of interest:
Next we generate some artificial future values for the costs and benefits. Note that the functions to create these future values were arbitrarily chosen (with much trial and error;) to satisfy the stylized facts above, and ensure the results are sensitive to the assumptions: i.e. each plant is the most preferred for some combination of horizon and discount rate. Note that cost benefit analysis utilizes net present value, so we need to discount these future values.
discount_rate <- seq(.055,.085,.001)
year <- 1:15
plant <- c("Secondary","Tertiary","Enhanced")
thing <- c("Cost","Benefit")
mydf <- crossing(year,plant,thing,discount_rate)%>%
mutate(future_value=case_when(plant == "Secondary" & thing == "Cost" ~ 21/year,
plant == "Secondary" & thing == "Benefit" ~ log(100000*year),
plant == "Tertiary" & thing == "Cost" ~ 39.99/year,
plant == "Tertiary" & thing == "Benefit" ~ 1.5035*log(100000*year),
plant == "Enhanced" & thing == "Cost" ~ 46.5/year,
plant == "Enhanced" & thing == "Benefit" ~ 1.67*log(100000*year)))%>%#THEN
mutate(present_value = future_value/((1+discount_rate)^year))
Below are graphs of the time series of the future value of costs and benefits, with a separate facet for each of the three types of plants.
plt <- ggplot(mydf,aes(year, future_value, colour=thing))+
geom_line()+
facet_wrap(vars(fct_reorder(plant,future_value,sum)))#over-ride default alphabetic layout.
plt+
labs(x=formatter(plt$labels$x),
y=formatter(plt$labels$y),
col=formatter(plt$labels$col))
We are interested in is seeing which plant is preferable for various combinations of time horizon and discount rate. In order to use only the data for the time horizon we are interested in we need to create various subsets of the data, that only contain the data up to the horizon of interest.
results <- tibble(Horizon=7:15)
subset_data <- function(horizon){
filter(mydf,year<=horizon)
}
results <- results%>%
mutate(data_subsets=map(Horizon,subset_data))
Once we have created the subsets of the data, we calculate the net present value for each of the horizon/discount rate combos.
get_npv <- function(df){
df%>%
group_by(plant, thing, discount_rate)%>%
summarise(total = sum(present_value))%>%
pivot_wider(names_from = thing, values_from = total, names_prefix = "sum_pv_")%>%
mutate(net_present_value = sum_pv_Benefit - sum_pv_Cost)
}
results <- results%>%
mutate(net_present_value=map(data_subsets,get_npv))%>%
unnest(net_present_value)
The plot below shows how the net present values associated with each plant depends on the discount rate (x axis) and the time horizon (the facets).
plt <- ggplot(results,aes(discount_rate,net_present_value,colour=plant))+
geom_line()+
facet_wrap(vars(Horizon), labeller = label_both)
plt+
labs(x=formatter(plt$labels$x),
y=formatter(plt$labels$y),
col=formatter(plt$labels$col))
plt <- ggplot(filter(results,Horizon==10),aes(discount_rate,net_present_value,colour=plant))+
geom_line()
plt+
labs(x=formatter(plt$labels$x),
y=formatter(plt$labels$y),
col=formatter(plt$labels$col),
title="10 year horizon")
If every combination of discount rate and time horizon are equally likely to be correct, then the expected net present values are:
exp_results <- results%>%
group_by(plant)%>%
summarize(expected_net_present_value=mean(net_present_value))%>%
arrange(desc(expected_net_present_value))
kbl(exp_results,col.names=formatter(colnames(exp_results)),digits=1)%>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = F)
| Plant | Expected Net Present Value |
|---|---|
| Enhanced | 49.5 |
| Tertiary | 49.0 |
| Secondary | 46.0 |