Function cost_tco2e()

This function emulates the calculations made by CC to estimate the cost of reducing a TCO2e through behavioral changes.

#' cost_tco2e
#'
#' This function determines the cost per TCO2e emulating the calculations made by CC.
#' 
#' @param years 
#' @param attribution 
#' @param total_expenses
#' @param total_partners
#' @param tco2e_solar_panel
#' @param tco2e_electric_vehicule
#' @param tco2e_get_around_greener
#' @param tco2e_less_beef
#' @param tco2e_more_veggies
#' @param tco2e_reduce_food_waste
#' @param tco2e_give_to_nature
#' @param pct_solar_panel
#' @param pct_electric_vehicule
#' @param pct_get_around_greener
#' @param pct_less_beef
#' @param pct_more_veggies
#' @param pct_reduce_food_waste
#' @param pct_give_to_nature
#' @param pct_adoption
#' @param pct_activation
#' @param pct_consideration
#' @param pct_awareness
#' @param audience_size
#' @return cost per TCO2e ($)
#' @author Mariano Viz


cost_tco2e <- function(years, attribution = 100, total_expenses, total_partners, 
                  tco2e_solar_panel = 4.2, tco2e_electric_vehicule = 2.26, tco2e_get_around_greener = 0, tco2e_less_beef = 0.81, tco2e_more_veggies = 0.077,  tco2e_reduce_food_waste = 0.026, tco2e_give_to_nature = 2,
                  pct_solar_panel = 10, pct_electric_vehicule = 25, pct_get_around_greener = 10, pct_less_beef = 10, pct_more_veggies = 20, pct_reduce_food_waste = 20, pct_give_to_nature = 5,
                  pct_adoption = 57.64, pct_activation = 5.32, pct_consideration = 2.06, pct_awareness = 50,
                  audience_size = 100000000){
  
 
  attribution = attribution/100
  pct_solar_panel <- pct_solar_panel/100
  pct_electric_vehicule <- pct_electric_vehicule/100
  pct_get_around_greener <- pct_get_around_greener/100
  pct_less_beef <- pct_less_beef/100
  pct_more_veggies <- pct_more_veggies/100
  pct_reduce_food_waste <- pct_reduce_food_waste/100
  pct_give_to_nature <- pct_give_to_nature/100
  pct_adoption <- pct_adoption/100
  pct_activation <- pct_activation/100
  pct_consideration <- pct_consideration/100
  pct_awareness <- pct_awareness/100

  df <- data.frame(matrix(ncol=7, nrow=length(years)))
  colnames(df)<-c('year', 'attribution', 'total_expenses', 'total_partners', 'tco2e_per_partner', 'tco2e_per_year', 'cost_per_tco2e')
  df$year <- years
  df$attribution <- attribution
  df$total_expenses <- total_expenses
  
  people_reached = audience_size * pct_awareness
  people_engaged = people_reached * pct_consideration
  beliefs_shifted = people_engaged * pct_activation 
  behaviors_changed_per_partner = beliefs_shifted * pct_adoption
  tco2e_per_partner = (tco2e_solar_panel*pct_solar_panel*behaviors_changed_per_partner)+(tco2e_electric_vehicule*pct_electric_vehicule*behaviors_changed_per_partner)+(tco2e_get_around_greener*pct_get_around_greener*behaviors_changed_per_partner)+(tco2e_less_beef*pct_less_beef*behaviors_changed_per_partner)+(tco2e_more_veggies*pct_more_veggies*behaviors_changed_per_partner)+(tco2e_reduce_food_waste*pct_reduce_food_waste*behaviors_changed_per_partner)+(tco2e_give_to_nature*pct_give_to_nature*behaviors_changed_per_partner)
  
  df$total_partners <- total_partners
  df$tco2e_per_partner <- tco2e_per_partner
  df$tco2e_per_year <- (df$total_partners * df$tco2e_per_partner)*df$attribution
  
  
  df$cost_per_tco2e <- df$total_expenses/df$tco2e_per_year
  

  return(df$cost_per_tco2e)
  
}

Testing cost_tco2e()

The results shown in the CC spreadsheet for the years 2024 to 2027 were 16.5, 18.3, 21.5, and 22.7 respectively.

round(cost_tco2e(years = 2024:2027,
         attribution = 100,
         total_expenses = c(2161776, 2740457, 3618359, 5535243),
         total_partners = c(3.5, 4, 4.5, 6.5)), 1)
## [1] 16.5 18.3 21.5 22.7

Parameter Elasticity

initial_cost_tco2e <- cost_tco2e(years = 2024:2027,
         attribution = 100,
         total_expenses = c(2161776, 2740457, 3618359, 5535243),
         total_partners = c(3.5, 4, 4.5, 6.5))

#1 Attribution 
at_cost_tco2e <- cost_tco2e(years = 2024:2027,
         attribution = 100*1.01, # 1% increase --> 101% attribution (?)
         total_expenses = c(2161776, 2740457, 3618359, 5535243),
         total_partners = c(3.5, 4, 4.5, 6.5))

initial_output <- sum(initial_cost_tco2e)
at_final_output <- sum(at_cost_tco2e)
at_initial_input <- 100
at_final_input <- 101
at_elasticity <- round(abs(((initial_output - at_final_output)/(at_initial_input-at_final_input))*(at_initial_input/initial_output)),2) 

#2 Total Expenses
te_cost_tco2e <- cost_tco2e(years = 2024:2027,
         attribution = 100,
         total_expenses = c(2161776*1.01, 2740457*1.01, 3618359*1.01, 5535243*1.01),# 1% increase
         total_partners = c(3.5, 4, 4.5, 6.5))

te_final_output <- sum(te_cost_tco2e)
te_initial_input <- sum(2161776, 2740457, 3618359, 5535243)
te_final_input <- sum(2161776*1.01, 2740457*1.01, 3618359*1.01, 5535243*1.01)
te_elasticity <- round(abs(((initial_output - te_final_output)/(te_initial_input-te_final_input))*(te_initial_input/initial_output)),2) 

#3 Total Partners
tp_cost_tco2e <- cost_tco2e(years = 2024:2027,
         attribution = 100,
         total_expenses = c(2161776, 2740457, 3618359, 5535243),
         total_partners = c(3.5*1.01, 4*1.01, 4.5*1.01, 6.5*1.01))# 1% increase

tp_final_output <- sum(tp_cost_tco2e)
tp_initial_input <- sum(3.5, 4, 4.5, 6.5)
tp_final_input <- sum(3.5*1.01, 4*1.01, 4.5*1.01, 6.5*1.01)
tp_elasticity <- round(abs(((initial_output - tp_final_output)/(tp_initial_input-tp_final_input))*(tp_initial_input/initial_output)),2) 

#4 TCO2e/behavior/year
  
  #4.1 Install Solar Panels (10%)
  sp_cost_tco2e <- cost_tco2e(years = 2024:2027,
         attribution = 100,
         total_expenses = c(2161776, 2740457, 3618359, 5535243),
         total_partners = c(3.5, 4, 4.5, 6.5),
         tco2e_solar_panel = 4.2*1.01)# 1% increase

  sp_final_output <- sum(sp_cost_tco2e)
  sp_initial_input <- 4.2
  sp_final_input <- 4.2*1.01
  sp_elasticity <- round(abs(((initial_output -     sp_final_output)/(sp_initial_input-sp_final_input))*(sp_initial_input/initial_output)),2)
  
   #4.2 Drive an Electric Vehicle (25%)
  ev_cost_tco2e <- cost_tco2e(years = 2024:2027,
         attribution = 100,
         total_expenses = c(2161776, 2740457, 3618359, 5535243),
         total_partners = c(3.5, 4, 4.5, 6.5),
         tco2e_electric_vehicule = 2.26*1.01)# 1% increase

  ev_final_output <- sum(ev_cost_tco2e)
  ev_initial_input <- 2.26
  ev_final_input <- 2.26*1.01
  ev_elasticity <- round(abs(((initial_output -     ev_final_output)/(ev_initial_input-ev_final_input))*(ev_initial_input/initial_output)),2)

  #4.3 Get Around Greener (10%)
    #Default value = 0 TCO2e/year
  
  #4.4 Eat Less Beef (10%)
  lb_cost_tco2e <- cost_tco2e(years = 2024:2027,
         attribution = 100,
         total_expenses = c(2161776, 2740457, 3618359, 5535243),
         total_partners = c(3.5, 4, 4.5, 6.5),
         tco2e_less_beef = 0.81*1.01)# 1% increase

  lb_final_output <- sum(lb_cost_tco2e)
  lb_initial_input <- 0.81
  lb_final_input <- 0.81*1.01
  lb_elasticity <- round(abs(((initial_output -     lb_final_output)/(lb_initial_input-lb_final_input))*(lb_initial_input/initial_output)),2)
  
  #4.5 Eat More Veggies (20%)
  mv_cost_tco2e <- cost_tco2e(years = 2024:2027,
         attribution = 100,
         total_expenses = c(2161776, 2740457, 3618359, 5535243),
         total_partners = c(3.5, 4, 4.5, 6.5),
         tco2e_more_veggies = 0.077*1.01)# 1% increase

  mv_final_output <- sum(mv_cost_tco2e)
  mv_initial_input <- 0.077
  mv_final_input <- 0.077*1.01
  mv_elasticity <- round(abs(((initial_output -     mv_final_output)/(mv_initial_input-lb_final_input))*(mv_initial_input/initial_output)),2)
  
  #4.6 Reduce Food Waste (20%)
  fw_cost_tco2e <- cost_tco2e(years = 2024:2027,
         attribution = 100,
         total_expenses = c(2161776, 2740457, 3618359, 5535243),
         total_partners = c(3.5, 4, 4.5, 6.5),
         tco2e_reduce_food_waste = 0.026*1.01)# 1% increase

  fw_final_output <- sum(fw_cost_tco2e)
  fw_initial_input <- 0.026
  fw_final_input <- 0.026*1.01
  fw_elasticity <- round(abs(((initial_output -     fw_final_output)/(fw_initial_input-fw_final_input))*(fw_initial_input/initial_output)),2)
  
  #4.7 Give to Nature (5%)
  gn_cost_tco2e <- cost_tco2e(years = 2024:2027,
         attribution = 100,
         total_expenses = c(2161776, 2740457, 3618359, 5535243),
         total_partners = c(3.5, 4, 4.5, 6.5),
         tco2e_give_to_nature = 2*1.01)# 1% increase

  gn_final_output <- sum(gn_cost_tco2e)
  gn_initial_input <- 2
  gn_final_input <- 2*1.01
  gn_elasticity <- round(abs(((initial_output -     gn_final_output)/(gn_initial_input-gn_final_input))*(gn_initial_input/initial_output)),2)
  
#5 Convertion Rates
  
  #5.1 Adoption
  adop_cost_tco2e <- cost_tco2e(years = 2024:2027,
         attribution = 100,
         total_expenses = c(2161776, 2740457, 3618359, 5535243),
         total_partners = c(3.5, 4, 4.5, 6.5),
         pct_adoption = 57.64*1.01)# 1% increase

  adop_final_output <- sum(adop_cost_tco2e)
  adop_initial_input <- 57.64
  adop_final_input <- 57.64*1.01
  adop_elasticity <- round(abs(((initial_output -     adop_final_output)/(adop_initial_input-adop_final_input))*(adop_initial_input/initial_output)),2)
  
  #5.2 Activation
  act_cost_tco2e <- cost_tco2e(years = 2024:2027,
         attribution = 100,
         total_expenses = c(2161776, 2740457, 3618359, 5535243),
         total_partners = c(3.5, 4, 4.5, 6.5),
         pct_activation = 5.32*1.01)# 1% increase

  act_final_output <- sum(act_cost_tco2e)
  act_initial_input <- 5.32
  act_final_input <- 5.32*1.01
  act_elasticity <- round(abs(((initial_output -     act_final_output)/(act_initial_input-act_final_input))*(act_initial_input/initial_output)),2)
  
  #5.3 Consideration
  con_cost_tco2e <- cost_tco2e(years = 2024:2027,
         attribution = 100,
         total_expenses = c(2161776, 2740457, 3618359, 5535243),
         total_partners = c(3.5, 4, 4.5, 6.5),
         pct_consideration = 2.06*1.01)# 1% increase

  con_final_output <- sum(con_cost_tco2e)
  con_initial_input <- 2.06
  con_final_input <- 2.06*1.01
  con_elasticity <- round(abs(((initial_output -     con_final_output)/(con_initial_input-con_final_input))*(con_initial_input/initial_output)),2)
  
  #5.4 Awareness 
  aw_cost_tco2e <- cost_tco2e(years = 2024:2027,
         attribution = 100,
         total_expenses = c(2161776, 2740457, 3618359, 5535243),
         total_partners = c(3.5, 4, 4.5, 6.5),
         pct_awareness = 50*1.01)# 1% increase

  aw_final_output <- sum(aw_cost_tco2e)
  aw_initial_input <- 50
  aw_final_input <- 50*1.01
  aw_elasticity <- round(abs(((initial_output -     aw_final_output)/(aw_initial_input-aw_final_input))*(aw_initial_input/initial_output)),2)
  
#6 Audience Size
az_cost_tco2e <- cost_tco2e(years = 2024:2027,
         attribution = 100,
         total_expenses = c(2161776, 2740457, 3618359, 5535243),
         total_partners = c(3.5, 4, 4.5, 6.5),
         audience_size = 100000000*1.01)# 1% increase

  az_final_output <- sum(aw_cost_tco2e)
  az_initial_input <- 100000000
  az_final_input <- 100000000*1.01
  az_elasticity <- round(abs(((initial_output -     az_final_output)/(az_initial_input-az_final_input))*(az_initial_input/initial_output)),2)

  
# Table
parameters <- c("Attribution (%)", "Total Expenses ($)", "Total Partners", "Solar Panels (TCO2e/y)r", "Electric Vehicles (TCO2e/yr)", "Less Beef (TCO2e/yr)", "More Veggies (TCO2e/yr)", "Food Waste (TCO2e/yr)", "Give to Nature (TCO2e/yr)", "Adoption (%)", "Activation (%)", "Consideration (%)", "Awarness (%)", "Audience Size")

default_value <- list(100, c(2161776, 2740457, 3618359, 5535243), c(3.5, 4, 4.5, 6.5), 4.2, 2.26, 0.81, 0.077, 0.026, 2, 57.64, 5.32, 2.06, 50, 100000000)

elasticity <- c(at_elasticity, te_elasticity, tp_elasticity, sp_elasticity, ev_elasticity, lb_elasticity, mv_elasticity, fw_elasticity, gn_elasticity, adop_elasticity, act_elasticity, con_elasticity, aw_elasticity, az_elasticity)

table <- data.frame(parameter=parameters, default=cbind(default_value), elasticity=elasticity)
table <-table[order(table$elasticity,decreasing=TRUE),] 
kbl(table, row.names = FALSE, col.names = c("Parameter", "Default", "Elasticity")) %>% 
  kable_styling()
Parameter Default Elasticity
Total Expenses ($) 2161776, 2740457, 3618359, 5535243 1.00
Attribution (%) 100 0.99
Total Partners 3.5, 4.0, 4.5, 6.5 0.99
Adoption (%) 57.64 0.99
Activation (%) 5.32 0.99
Consideration (%) 2.06 0.99
Awarness (%) 50 0.99
Audience Size 1e+08 0.99
Electric Vehicles (TCO2e/yr) 2.26 0.47
Solar Panels (TCO2e/y)r 4.2 0.35
Give to Nature (TCO2e/yr) 2 0.08
Less Beef (TCO2e/yr) 0.81 0.07
More Veggies (TCO2e/yr) 0.077 0.00
Food Waste (TCO2e/yr) 0.026 0.00

Behavior Change Model Simulation

The following is a simulation based on the Behavior Change Model proposed by Chris. Its objective is to conceptualize the determination of the impact of CC interventions on behaviors over time.

# Behavior Model with Rare intervention (qj_r = 0.4)
r_fun <- function(df, xij0, pi=0.4, fj=0.8, qj_r=0.4, dr=0.02){
  df[1,2:3] <- xij0
  
  for(i in 2:nrow(df)) {
   df$xijt_r[i] = (df$xijt_r[i-1]*fj*(1-qj_r*pi)+qj_r*pi)
  }
  for(i in 1:nrow(df)) {
   df$xijt_r[i] = df$xijt_r[i]*((1/(1+dr))^df$t[i]) #discounting (2% rate)
  }
  
  return(list(df[1:11,2]))
}

# Behavior Model without Rare intervention (qj_nr = 0.2)
nr_fun <- function(df, xij0, pi=0.4, fj=0.8, qj_r=0.4, qj_nr=0.2, dr=0.02){
  df[1,2:3] <- xij0
  
  for(i in 2:nrow(df)) {
   df$xijt_nr[i] = df$xijt_nr[i-1]*fj*(1-qj_nr*pi)+qj_nr*pi
  }
  for(i in 1:nrow(df)) {
   df$xijt_nr[i] = df$xijt_nr[i]*((1/(1+dr))^df$t[i])
  }
  
  return(list(df[1:11,3]))
  
}

# 100 simulations over 10 years with Rare intervention
set.seed(1) # initial behavior stock (xij0) comes from a random sample of a uniform distribution from 0 to 0.5
r <- replicate(
  100,
  {
behavior_model_df <- data.frame(t = 0:10, xijt_r = NA, xijt_nr = NA)
xij0 <- runif(1, min = 0, max = 0.5)
r_fun(behavior_model_df, xij0)
  }
)
r <- print(as.data.frame(do.call(cbind, r)))
r <-   data.frame(t = 0:10, r)
r_melt <- melt(r ,  id.vars = 't', variable.name = 'xijt')

# 100 simulations over 10 years without Rare intervention
set.seed(1)
nr <- replicate(
  100,
  {
behavior_model_df <- data.frame(t = 0:10, xijt_r = NA, xijt_nr = NA)
xij0 <- runif(1, min = 0, max = 0.5)
nr_fun(behavior_model_df, xij0)
  }
)
nr <- print(as.data.frame(do.call(cbind, nr)))
nr <-   data.frame(t = 0:10, nr)
nr_melt <- melt(nr ,  id.vars = 't', variable.name = 'xijt')


# Average with and without Rare intevention
r_mean <- r
r_mean$xijt_r_mean <- rowMeans(r[,2:101])
nr_mean <- nr
nr_mean$xijt_nr_mean <- rowMeans(nr[,2:101])
xijt <- data.frame(t = 0:10, xijt_r = r_mean$xijt_r_mean, xijt_nr = nr_mean$xijt_nr_mean)

Simulation Results Graphs

ggplot(r_melt, aes(t, value)) +
  geom_line(aes(colour = xijt))+
  scale_colour_grey(start = 0.7, end = 0.7)+
  stat_summary(geom="line", fun = "mean", color="red")+
  scale_x_continuous(breaks=c(1,2,3,4,5,6,7,8,9,10), limits = c(0,10))+
  theme_minimal()+
  theme(legend.position = "none")+
  labs(x = "t",
       y = "“good behavior” stock",
       title = "Behavior Change with Rare intervention") 

ggplot(nr_melt, aes(t, value)) +
  geom_line(aes(colour = xijt))+
  scale_colour_grey(start = 0.7, end = 0.7)+
  stat_summary(geom="line", fun = "mean", color="red")+
  scale_x_continuous(breaks=c(1,2,3,4,5,6,7,8,9,10), limits = c(0,10))+
  theme_minimal()+
  theme(legend.position = "none")+
  labs(x = "t",
       y = "“good behavior” stock",
       title = "Behavior Change without Rare intervention") 

ggplot(xijt)+
  geom_line(aes(t, xijt_r))+
  geom_line(aes(t, xijt_nr))+
  scale_x_continuous(breaks=c(1,2,3,4,5,6,7,8,9,10), limits = c(0,10))+
  geom_ribbon(aes(ymin=xijt_r, ymax= xijt_nr, x = t ), fill = "red", alpha = 0.3)+
  geom_text(aes(x = 2, y = 0.40, label= "with Rare", hjust = 0), size = 3.5) +
  geom_text(aes(x = 2, y = 0.26, label= "without Rare", hjust = 0), size = 3.5) +
  theme_minimal()+
  theme(legend.position = "none")+
  labs(x = "t",
       y = "“good behavior” stock",
       title = "Behavior Change with/without Rare intervention")

Function behavior_model_ctco2e()

This function unifies the Behavior Change Model with the calculations made by CC to estimate the cost of reducing carbon emissions.

# Output = Total Behaviors changed

behavior_model_tbc <- function(df, total_expenses, total_partners,                            
                            xij0, pi, qj, fj, dr,
                            attribution = 100,
                            tco2e_solar_panel = 4.2, tco2e_electric_vehicule = 2.26, tco2e_get_around_greener = 0, tco2e_less_beef = 0.81, tco2e_more_veggies = 0.077,  tco2e_reduce_food_waste = 0.026, tco2e_give_to_nature = 2,
                            pct_solar_panel = 10, pct_electric_vehicule = 25, pct_get_around_greener = 10, pct_less_beef = 10, pct_more_veggies = 20, pct_reduce_food_waste = 20, pct_give_to_nature = 5,
                            #pct_adoption = 57.64, pct_activation = 5.32, pct_consideration = 2.06, pct_awareness = 50,
                  audience_size = 100000000){

  
  attribution = attribution/100
  pct_solar_panel <- pct_solar_panel/100
  pct_electric_vehicule <- pct_electric_vehicule/100
  pct_get_around_greener <- pct_get_around_greener/100
  pct_less_beef <- pct_less_beef/100
  pct_more_veggies <- pct_more_veggies/100
  pct_reduce_food_waste <- pct_reduce_food_waste/100
  pct_give_to_nature <- pct_give_to_nature/100

  df[1,2] <- xij0
  df$total_expenses <- total_expenses
  df$total_partners <- total_partners


  for(i in 2:nrow(df)) {
   df$xijt[i] = (df$xijt[i-1]*fj*(1-qj*pi)+qj*pi)}
  
  for(i in 1:nrow(df)) {
   df$xijt[i] = df$xijt[i]*((1/(1+dr))^df$t[i])}

  for(i in 1:nrow(df)) {
   df$behaviors_changed_per_partner[i] = df$xijt[i]*audience_size}
  
  for(i in 1:nrow(df)) {
   df$total_behaviors_changed[i] = df$behaviors_changed_per_partner[i]*df$total_partners[i]}

  for(i in 1:nrow(df)) {
   df$tco2e_per_partner[i] = (tco2e_solar_panel*pct_solar_panel*df$behaviors_changed_per_partner[i])+(tco2e_electric_vehicule*pct_electric_vehicule*df$behaviors_changed_per_partner[i])+(tco2e_get_around_greener*pct_get_around_greener*df$behaviors_changed_per_partner[i])+(tco2e_less_beef*pct_less_beef*df$behaviors_changed_per_partner[i])+(tco2e_more_veggies*pct_more_veggies*df$behaviors_changed_per_partner[i])+(tco2e_reduce_food_waste*pct_reduce_food_waste*df$behaviors_changed_per_partner[i])+(tco2e_give_to_nature*pct_give_to_nature*df$behaviors_changed_per_partner[i])}

  for(i in 1:nrow(df)) {
   df$tco2e_per_year[i] = (df$total_partners[i]*df$tco2e_per_partner[i])*attribution}

  for(i in 1:nrow(df)) {
   df$cost_per_tco2e[i] = df$total_expenses[i]/df$tco2e_per_year[i]}
  
  
  return(df$total_behaviors_changed)

}
  

# Output = TCO2e per year

behavior_model_ttco2e <- function(df, total_expenses, total_partners,                            
                            xij0, pi, qj, fj, dr,
                            attribution = 100,
                            tco2e_solar_panel = 4.2, tco2e_electric_vehicule = 2.26, tco2e_get_around_greener = 0, tco2e_less_beef = 0.81, tco2e_more_veggies = 0.077,  tco2e_reduce_food_waste = 0.026, tco2e_give_to_nature = 2,
                            pct_solar_panel = 10, pct_electric_vehicule = 25, pct_get_around_greener = 10, pct_less_beef = 10, pct_more_veggies = 20, pct_reduce_food_waste = 20, pct_give_to_nature = 5,
                            #pct_adoption = 57.64, pct_activation = 5.32, pct_consideration = 2.06, pct_awareness = 50,
                  audience_size = 100000000){

  
  attribution = attribution/100
  pct_solar_panel <- pct_solar_panel/100
  pct_electric_vehicule <- pct_electric_vehicule/100
  pct_get_around_greener <- pct_get_around_greener/100
  pct_less_beef <- pct_less_beef/100
  pct_more_veggies <- pct_more_veggies/100
  pct_reduce_food_waste <- pct_reduce_food_waste/100
  pct_give_to_nature <- pct_give_to_nature/100

  df[1,2] <- xij0
  df$total_expenses <- total_expenses
  df$total_partners <- total_partners


  for(i in 2:nrow(df)) {
   df$xijt[i] = (df$xijt[i-1]*fj*(1-qj*pi)+qj*pi)}
  
  for(i in 1:nrow(df)) {
   df$xijt[i] = df$xijt[i]*((1/(1+dr))^df$t[i])}

  for(i in 1:nrow(df)) {
   df$behaviors_changed_per_partner[i] = df$xijt[i]*audience_size}
  
  for(i in 1:nrow(df)) {
   df$total_behaviors_changed[i] = df$behaviors_changed_per_partner[i]*df$total_partners[i]}

  for(i in 1:nrow(df)) {
   df$tco2e_per_partner[i] = (tco2e_solar_panel*pct_solar_panel*df$behaviors_changed_per_partner[i])+(tco2e_electric_vehicule*pct_electric_vehicule*df$behaviors_changed_per_partner[i])+(tco2e_get_around_greener*pct_get_around_greener*df$behaviors_changed_per_partner[i])+(tco2e_less_beef*pct_less_beef*df$behaviors_changed_per_partner[i])+(tco2e_more_veggies*pct_more_veggies*df$behaviors_changed_per_partner[i])+(tco2e_reduce_food_waste*pct_reduce_food_waste*df$behaviors_changed_per_partner[i])+(tco2e_give_to_nature*pct_give_to_nature*df$behaviors_changed_per_partner[i])}

  for(i in 1:nrow(df)) {
   df$tco2e_per_year[i] = (df$total_partners[i]*df$tco2e_per_partner[i])*attribution}

  for(i in 1:nrow(df)) {
   df$cost_per_tco2e[i] = df$total_expenses[i]/df$tco2e_per_year[i]}
  
  
  return(df$tco2e_per_year)

}
  
  
# Output = Cost per TCO2e

behavior_model_ctco2e <- function(df, total_expenses, total_partners,                            
                            xij0, pi, qj, fj, dr,
                            attribution = 100,
                            tco2e_solar_panel = 4.2, tco2e_electric_vehicule = 2.26, tco2e_get_around_greener = 0, tco2e_less_beef = 0.81, tco2e_more_veggies = 0.077,  tco2e_reduce_food_waste = 0.026, tco2e_give_to_nature = 2,
                            pct_solar_panel = 10, pct_electric_vehicule = 25, pct_get_around_greener = 10, pct_less_beef = 10, pct_more_veggies = 20, pct_reduce_food_waste = 20, pct_give_to_nature = 5,
                            #pct_adoption = 57.64, pct_activation = 5.32, pct_consideration = 2.06, pct_awareness = 50,
                  audience_size = 100000000){

  
  attribution = attribution/100
  pct_solar_panel <- pct_solar_panel/100
  pct_electric_vehicule <- pct_electric_vehicule/100
  pct_get_around_greener <- pct_get_around_greener/100
  pct_less_beef <- pct_less_beef/100
  pct_more_veggies <- pct_more_veggies/100
  pct_reduce_food_waste <- pct_reduce_food_waste/100
  pct_give_to_nature <- pct_give_to_nature/100

  df[1,2] <- xij0
  df$total_expenses <- total_expenses
  df$total_partners <- total_partners


  for(i in 2:nrow(df)) {
   df$xijt[i] = (df$xijt[i-1]*fj*(1-qj*pi)+qj*pi)}
  
  for(i in 1:nrow(df)) {
   df$xijt[i] = df$xijt[i]*((1/(1+dr))^df$t[i])}

  for(i in 1:nrow(df)) {
   df$behaviors_changed_per_partner[i] = df$xijt[i]*audience_size}
  
  for(i in 1:nrow(df)) {
   df$total_behaviors_changed[i] = df$behaviors_changed_per_partner[i]*df$total_partners[i]}

  for(i in 1:nrow(df)) {
   df$tco2e_per_partner[i] = (tco2e_solar_panel*pct_solar_panel*df$behaviors_changed_per_partner[i])+(tco2e_electric_vehicule*pct_electric_vehicule*df$behaviors_changed_per_partner[i])+(tco2e_get_around_greener*pct_get_around_greener*df$behaviors_changed_per_partner[i])+(tco2e_less_beef*pct_less_beef*df$behaviors_changed_per_partner[i])+(tco2e_more_veggies*pct_more_veggies*df$behaviors_changed_per_partner[i])+(tco2e_reduce_food_waste*pct_reduce_food_waste*df$behaviors_changed_per_partner[i])+(tco2e_give_to_nature*pct_give_to_nature*df$behaviors_changed_per_partner[i])}

  for(i in 1:nrow(df)) {
   df$tco2e_per_year[i] = (df$total_partners[i]*df$tco2e_per_partner[i])*attribution}

  for(i in 1:nrow(df)) {
   df$cost_per_tco2e[i] = df$total_expenses[i]/df$tco2e_per_year[i]}
  
  
  return(df$cost_per_tco2e)

}

Testing behavior_model_ctco2e()

The results shown in the CC spreadsheet for the years 2024 to 2027 were 16.5, 18.3, 21.5, and 22.8 respectively.

behavior_model_df <- data.frame(t = 0:6, xijt = NA)
round(behavior_model_ctco2e(behavior_model_df, 
pi=1028282.215/100000000, qj=31509.893/1028282.215, xij0=(1028282.215/100000000)*(31509.893/1028282.215), fj=0, dr=0.00, audience_size=100000000,
total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856), total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18)), 1)
## [1] 16.5 18.3 21.5 22.8 20.5 18.2 15.9
pi <- 1028282.215/100000000
qj <- 31509.893/1028282.215
xij0 <- qj*pi
fj <- 0
dr <- 0.00
audience_size <- 100000000
total_expenses <- c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856)
total_partners <- c(3.5, 4, 4.5, 6.5, 9, 13, 18)
tco2e_solar_panel<-4.2 
tco2e_electric_vehicule<-2.26 
tco2e_get_around_greener<-0 
tco2e_less_beef<-0.81 
tco2e_more_veggies<-0.077  
tco2e_reduce_food_waste<-0.026 
tco2e_give_to_nature<-2
pct_solar_panel<-10/100
pct_electric_vehicule<-25/100
pct_get_around_greener<-10/100
pct_less_beef<-10/100 
pct_more_veggies<-20/100 
pct_reduce_food_waste<-20/100
pct_give_to_nature<-5/100
attribution <- 100/100

df <- data.frame(t = 0:6, xijt = NA)

df[1,2] <- xij0
df$total_expenses <- total_expenses
df$total_partners <- total_partners


for(i in 2:nrow(df)) {
   df$xijt[i] = (df$xijt[i-1]*fj*(1-qj*pi)+qj*pi)}
  
for(i in 1:nrow(df)) {
   df$xijt[i] = df$xijt[i]*((1/(1+dr))^df$t[i])}

for(i in 1:nrow(df)) {
   df$behaviors_changed_per_partner[i] = df$xijt[i]*audience_size}

for(i in 1:nrow(df)) {
   df$total_behaviors_changed[i] = df$behaviors_changed_per_partner[i]*df$total_partners[i]}

for(i in 1:nrow(df)) {
   df$tco2e_per_partner[i] = (tco2e_solar_panel*pct_solar_panel*df$behaviors_changed_per_partner[i])+(tco2e_electric_vehicule*pct_electric_vehicule*df$behaviors_changed_per_partner[i])+(tco2e_get_around_greener*pct_get_around_greener*df$behaviors_changed_per_partner[i])+(tco2e_less_beef*pct_less_beef*df$behaviors_changed_per_partner[i])+(tco2e_more_veggies*pct_more_veggies*df$behaviors_changed_per_partner[i])+(tco2e_reduce_food_waste*pct_reduce_food_waste*df$behaviors_changed_per_partner[i])+(tco2e_give_to_nature*pct_give_to_nature*df$behaviors_changed_per_partner[i])}

for(i in 1:nrow(df)) {
   df$tco2e_per_year[i] = (df$total_partners[i]*df$tco2e_per_partner[i])*attribution}

for(i in 1:nrow(df)) {
   df$cost_per_tco2e[i] = df$total_expenses[i]/df$tco2e_per_year[i]}




ggplot(df) +
  geom_line(aes(t, behaviors_changed_per_partner))+
  scale_x_continuous(breaks=c(1,2,3,4,5,6), limits = c(0,6))+
  scale_y_continuous(labels = scales::comma)+
  theme_minimal()+
  labs(x = "t",
       y = "Behaviors Changed",
       title = "Behaviors Changed per Partner") 

ggplot(df) +
  geom_area(aes(t, total_behaviors_changed), color = "grey", alpha = 0.3)+
  geom_line(aes(t, total_behaviors_changed))+
  scale_x_continuous(breaks=c(1,2,3,4,5,6), limits = c(0,6))+
  scale_y_continuous(labels = scales::comma)+
  theme_minimal()+
  labs(x = "t",
       y = "Behaviors Changed",
       title = "Total Behaviors Changed") 

ggplot(df) +
  geom_area(aes(t, tco2e_per_year), color = "grey", alpha = 0.3)+
  geom_line(aes(t, tco2e_per_year))+
  scale_x_continuous(breaks=c(1,2,3,4,5,6), limits = c(0,6))+
  scale_y_continuous(labels = scales::comma)+
  theme_minimal()+
  labs(x = "t",
       y = "TCO2e",
       title = "Total TCO2e reduction") 

ggplot(df) +
  geom_line(aes(t, cost_per_tco2e))+
  geom_point(aes(t, cost_per_tco2e))+
  scale_x_continuous(breaks=c(1,2,3,4,5,6), limits = c(0,6))+
  scale_y_continuous(labels = scales::comma)+
  theme_minimal()+
  labs(x = "t",
       y = "USD",
       title = "Cost per TCO2e")

Simulations plugging in different values: xij0, fj, qj, pi

set.seed(12)
bm_tbc <- replicate(
  100,
  {
behavior_model_df <- data.frame(t = 0:6, xijt = NA)

pi <- 1028282.215/100000000
qj <- 31509.893/1028282.215
fj <- 0
xij0 <- runif(1, min = 0, max = 0.3)

behavior_model_tbc(behavior_model_df, 
pi=pi, qj=qj, xij0=xij0, fj=fj, dr=0.02, audience_size=100000000,
total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856), total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18))
  })

bm_tbc <-   data.frame(t = 0:6, bm_tbc)
bm_tbc_melt <- melt(bm_tbc ,  id.vars = 't', variable.name = 'total_behaviors_changed')

ggplot(bm_tbc_melt, aes(t, value)) +
  geom_line(aes(colour = total_behaviors_changed))+
  scale_colour_grey(start = 0.7, end = 0.7)+
  stat_summary(geom="line", fun = "mean", color="red")+
  scale_x_continuous(breaks=c(1,2,3,4,5,6), limits = c(0,6))+
  scale_y_continuous(labels = scales::comma)+
  theme_minimal()+
  theme(legend.position = "none")+
  labs(x = "t",
       y = "Behaviors Changed",
       title = "Total Behaviors Changed")  

set.seed(12)
bm_ttco2e <- replicate(
  100,
  {
behavior_model_df <- data.frame(t = 0:6, xijt = NA)

pi <- 1028282.215/100000000
qj <- 31509.893/1028282.215
fj <- 0
xij0 <- runif(1, min = 0, max = 0.3)

behavior_model_ttco2e(behavior_model_df, 
pi=pi, qj=qj, xij0=xij0, fj=fj, dr=0.02, audience_size=100000000,
total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856), total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18))
  })

bm_ttco2e <-   data.frame(t = 0:6, bm_ttco2e)
bm_ttco2e_melt <- melt(bm_ttco2e ,  id.vars = 't', variable.name = 'tco2e_per_year')

ggplot(bm_ttco2e_melt, aes(t, value)) +
  geom_line(aes(colour = tco2e_per_year))+
  scale_colour_grey(start = 0.7, end = 0.7)+
  stat_summary(geom="line", fun = "mean", color="red")+
  scale_x_continuous(breaks=c(1,2,3,4,5,6), limits = c(0,6))+
  scale_y_continuous(labels = scales::comma)+
  theme_minimal()+
  theme(legend.position = "none")+
  labs(x = "t",
       y = "TCO2e",
       title = "Total TCO2e reduction") 

set.seed(12)
bm_ctco2e <- replicate(
  100,
  {
behavior_model_df <- data.frame(t = 0:6, xijt = NA)

pi <- 1028282.215/100000000
qj <- 31509.893/1028282.215
fj <- 0
xij0 <- runif(1, min = 0, max = 0.3)

behavior_model_ctco2e(behavior_model_df, 
pi=pi, qj=qj, xij0=xij0, fj=fj, dr=0.02, audience_size=100000000,
total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856), total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18))
  })

bm_ctco2e <-   data.frame(t = 0:6, bm_ctco2e)
bm_ctco2e_melt <- melt(bm_ctco2e ,  id.vars = 't', variable.name = 'cost_per_tco2e')

ggplot(bm_ctco2e_melt, aes(t, value)) +
  geom_line(aes(colour = cost_per_tco2e))+
  scale_colour_grey(start = 0.7, end = 0.7)+
  stat_summary(geom="line", fun = "mean", color="red")+
  scale_x_continuous(breaks=c(1,2,3,4,5,6), limits = c(0,6))+
  scale_y_continuous(labels = scales::comma)+
  theme_minimal()+
  theme(legend.position = "none")+
  labs(x = "t",
       y = "USD",
       title = "Cost per TCO2e") 

set.seed(12)
bm_tbc <- replicate(
  100,
  {
behavior_model_df <- data.frame(t = 0:6, xijt = NA)

pi <- 1028282.215/100000000
qj <- 31509.893/1028282.215
fj <- 0.5
xij0 <- runif(1, min = 0, max = 0.3)

behavior_model_tbc(behavior_model_df, 
pi=pi, qj=qj, xij0=xij0, fj=fj, dr=0.02, audience_size=100000000,
total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856), total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18))
  })

bm_tbc <-   data.frame(t = 0:6, bm_tbc)
bm_tbc_melt <- melt(bm_tbc ,  id.vars = 't', variable.name = 'total_behaviors_changed')

ggplot(bm_tbc_melt, aes(t, value)) +
  geom_line(aes(colour = total_behaviors_changed))+
  scale_colour_grey(start = 0.7, end = 0.7)+
  stat_summary(geom="line", fun = "mean", color="red")+
  scale_x_continuous(breaks=c(1,2,3,4,5,6), limits = c(0,6))+
  scale_y_continuous(labels = scales::comma)+
  theme_minimal()+
  theme(legend.position = "none")+
  labs(x = "t",
       y = "Behaviors Changed",
       title = "Total Behaviors Changed")  

set.seed(12)
bm_ttco2e <- replicate(
  100,
  {
behavior_model_df <- data.frame(t = 0:6, xijt = NA)

pi <- 1028282.215/100000000
qj <- 31509.893/1028282.215
fj <- 0.5
xij0 <- runif(1, min = 0, max = 0.3)

behavior_model_ttco2e(behavior_model_df, 
pi=pi, qj=qj, xij0=xij0, fj=fj, dr=0.02, audience_size=100000000,
total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856), total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18))
  })

bm_ttco2e <-   data.frame(t = 0:6, bm_ttco2e)
bm_ttco2e_melt <- melt(bm_ttco2e ,  id.vars = 't', variable.name = 'tco2e_per_year')

ggplot(bm_ttco2e_melt, aes(t, value)) +
  geom_line(aes(colour = tco2e_per_year))+
  scale_colour_grey(start = 0.7, end = 0.7)+
  stat_summary(geom="line", fun = "mean", color="red")+
  scale_x_continuous(breaks=c(1,2,3,4,5,6), limits = c(0,6))+
  scale_y_continuous(labels = scales::comma)+
  theme_minimal()+
  theme(legend.position = "none")+
  labs(x = "t",
       y = "TCO2e",
       title = "Total TCO2e reduction") 

set.seed(12)
bm_ctco2e <- replicate(
  100,
  {
behavior_model_df <- data.frame(t = 0:6, xijt = NA)

pi <- 1028282.215/100000000
qj <- 31509.893/1028282.215
fj <- 0.5
xij0 <- runif(1, min = 0, max = 0.3)

behavior_model_ctco2e(behavior_model_df, 
pi=pi, qj=qj, xij0=xij0, fj=fj, dr=0.02, audience_size=100000000,
total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856), total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18))
  })

bm_ctco2e <-   data.frame(t = 0:6, bm_ctco2e)
bm_ctco2e_melt <- melt(bm_ctco2e ,  id.vars = 't', variable.name = 'cost_per_tco2e')

ggplot(bm_ctco2e_melt, aes(t, value)) +
  geom_line(aes(colour = cost_per_tco2e))+
  scale_colour_grey(start = 0.7, end = 0.7)+
  stat_summary(geom="line", fun = "mean", color="red")+
  scale_x_continuous(breaks=c(1,2,3,4,5,6), limits = c(0,6))+
  scale_y_continuous(labels = scales::comma)+
  theme_minimal()+
  theme(legend.position = "none")+
  labs(x = "t",
       y = "USD",
       title = "Cost per TCO2e")

set.seed(12)
bm_tbc <- replicate(
  100,
  {
behavior_model_df <- data.frame(t = 0:6, xijt = NA)

pi <- 1028282.215/100000000
qj <- (31509.893/1028282.215)*10
fj <- 0.5
xij0 <- runif(1, min = 0, max = 0.3)

behavior_model_tbc(behavior_model_df, 
pi=pi, qj=qj, xij0=xij0, fj=fj, dr=0.02, audience_size=100000000,
total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856), total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18))
  })

bm_tbc <-   data.frame(t = 0:6, bm_tbc)
bm_tbc_melt <- melt(bm_tbc ,  id.vars = 't', variable.name = 'total_behaviors_changed')

ggplot(bm_tbc_melt, aes(t, value)) +
  geom_line(aes(colour = total_behaviors_changed))+
  scale_colour_grey(start = 0.7, end = 0.7)+
  stat_summary(geom="line", fun = "mean", color="red")+
  scale_x_continuous(breaks=c(1,2,3,4,5,6), limits = c(0,6))+
  scale_y_continuous(labels = scales::comma)+
  theme_minimal()+
  theme(legend.position = "none")+
  labs(x = "t",
       y = "Behaviors Changed",
       title = "Total Behaviors Changed")  

set.seed(12)
bm_ttco2e <- replicate(
  100,
  {
behavior_model_df <- data.frame(t = 0:6, xijt = NA)

pi <- 1028282.215/100000000
qj <- (31509.893/1028282.215)*10
fj <- 0.5
xij0 <- runif(1, min = 0, max = 0.3)

behavior_model_ttco2e(behavior_model_df, 
pi=pi, qj=qj, xij0=xij0, fj=fj, dr=0.02, audience_size=100000000,
total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856), total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18))
  })

bm_ttco2e <-   data.frame(t = 0:6, bm_ttco2e)
bm_ttco2e_melt <- melt(bm_ttco2e ,  id.vars = 't', variable.name = 'tco2e_per_year')

ggplot(bm_ttco2e_melt, aes(t, value)) +
  geom_line(aes(colour = tco2e_per_year))+
  scale_colour_grey(start = 0.7, end = 0.7)+
  stat_summary(geom="line", fun = "mean", color="red")+
  scale_x_continuous(breaks=c(1,2,3,4,5,6), limits = c(0,6))+
  scale_y_continuous(labels = scales::comma)+
  theme_minimal()+
  theme(legend.position = "none")+
  labs(x = "t",
       y = "TCO2e",
       title = "Total TCO2e reduction") 

set.seed(12)
bm_ctco2e <- replicate(
  100,
  {
behavior_model_df <- data.frame(t = 0:6, xijt = NA)

pi <- 1028282.215/100000000
qj <- (31509.893/1028282.215)*10
fj <- 0.5
xij0 <- runif(1, min = 0, max = 0.3)

behavior_model_ctco2e(behavior_model_df, 
pi=pi, qj=qj, xij0=xij0, fj=fj, dr=0.02, audience_size=100000000,
total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856), total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18))
  })

bm_ctco2e <-   data.frame(t = 0:6, bm_ctco2e)
bm_ctco2e_melt <- melt(bm_ctco2e ,  id.vars = 't', variable.name = 'cost_per_tco2e')

ggplot(bm_ctco2e_melt, aes(t, value)) +
  geom_line(aes(colour = cost_per_tco2e))+
  scale_colour_grey(start = 0.7, end = 0.7)+
  stat_summary(geom="line", fun = "mean", color="red")+
  scale_x_continuous(breaks=c(1,2,3,4,5,6), limits = c(0,6))+
  scale_y_continuous(labels = scales::comma)+
  theme_minimal()+
  theme(legend.position = "none")+
  labs(x = "t",
       y = "USD",
       title = "Cost per TCO2e")

set.seed(12)
bm_tbc <- replicate(
  100,
  {
behavior_model_df <- data.frame(t = 0:6, xijt = NA)

pi <- (1028282.215/100000000)*10
qj <- 31509.893/1028282.215
fj <- 0.5
xij0 <- runif(1, min = 0, max = 0.3)

behavior_model_tbc(behavior_model_df, 
pi=pi, qj=qj, xij0=xij0, fj=fj, dr=0.02, audience_size=100000000,
total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856), total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18))
  })

bm_tbc <-   data.frame(t = 0:6, bm_tbc)
bm_tbc_melt <- melt(bm_tbc ,  id.vars = 't', variable.name = 'total_behaviors_changed')

ggplot(bm_tbc_melt, aes(t, value)) +
  geom_line(aes(colour = total_behaviors_changed))+
  scale_colour_grey(start = 0.7, end = 0.7)+
  stat_summary(geom="line", fun = "mean", color="red")+
  scale_x_continuous(breaks=c(1,2,3,4,5,6), limits = c(0,6))+
  scale_y_continuous(labels = scales::comma)+
  theme_minimal()+
  theme(legend.position = "none")+
  labs(x = "t",
       y = "Behaviors Changed",
       title = "Total Behaviors Changed")  

set.seed(12)
bm_ttco2e <- replicate(
  100,
  {
behavior_model_df <- data.frame(t = 0:6, xijt = NA)

pi <- (1028282.215/100000000)*10
qj <- 31509.893/1028282.215
fj <- 0.5
xij0 <- runif(1, min = 0, max = 0.3)

behavior_model_ttco2e(behavior_model_df, 
pi=pi, qj=qj, xij0=xij0, fj=fj, dr=0.02, audience_size=100000000,
total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856), total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18))
  })

bm_ttco2e <-   data.frame(t = 0:6, bm_ttco2e)
bm_ttco2e_melt <- melt(bm_ttco2e ,  id.vars = 't', variable.name = 'tco2e_per_year')

ggplot(bm_ttco2e_melt, aes(t, value)) +
  geom_line(aes(colour = tco2e_per_year))+
  scale_colour_grey(start = 0.7, end = 0.7)+
  stat_summary(geom="line", fun = "mean", color="red")+
  scale_x_continuous(breaks=c(1,2,3,4,5,6), limits = c(0,6))+
  scale_y_continuous(labels = scales::comma)+
  theme_minimal()+
  theme(legend.position = "none")+
  labs(x = "t",
       y = "TCO2e",
       title = "Total TCO2e reduction") 

set.seed(12)
bm_ctco2e <- replicate(
  100,
  {
behavior_model_df <- data.frame(t = 0:6, xijt = NA)

pi <- (1028282.215/100000000)*10
qj <- 31509.893/1028282.215
fj <- 0.5
xij0 <- runif(1, min = 0, max = 0.3)

behavior_model_ctco2e(behavior_model_df, 
pi=pi, qj=qj, xij0=xij0, fj=fj, dr=0.02, audience_size=100000000,
total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856), total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18))
  })

bm_ctco2e <-   data.frame(t = 0:6, bm_ctco2e)
bm_ctco2e_melt <- melt(bm_ctco2e ,  id.vars = 't', variable.name = 'cost_per_tco2e')

ggplot(bm_ctco2e_melt, aes(t, value)) +
  geom_line(aes(colour = cost_per_tco2e))+
  scale_colour_grey(start = 0.7, end = 0.7)+
  stat_summary(geom="line", fun = "mean", color="red")+
  scale_x_continuous(breaks=c(1,2,3,4,5,6), limits = c(0,6))+
  scale_y_continuous(labels = scales::comma)+
  theme_minimal()+
  theme(legend.position = "none")+
  labs(x = "t",
       y = "USD",
       title = "Cost per TCO2e")

set.seed(12)
bm_tbc <- replicate(
  100,
  {
behavior_model_df <- data.frame(t = 0:6, xijt = NA)

pi <- (1028282.215/100000000)*10
qj <- 31509.893/1028282.215
fj <- 0.5
xij0 <- runif(1, min = 0, max = 0.3)

behavior_model_tbc(behavior_model_df, 
pi=pi, qj=qj, xij0=xij0, fj=fj, dr=0.2, audience_size=100000000,
total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856), total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18))
  })

bm_tbc <-   data.frame(t = 0:6, bm_tbc)
bm_tbc_melt <- melt(bm_tbc ,  id.vars = 't', variable.name = 'total_behaviors_changed')

ggplot(bm_tbc_melt, aes(t, value)) +
  geom_line(aes(colour = total_behaviors_changed))+
  scale_colour_grey(start = 0.7, end = 0.7)+
  stat_summary(geom="line", fun = "mean", color="red")+
  scale_x_continuous(breaks=c(1,2,3,4,5,6), limits = c(0,6))+
  scale_y_continuous(labels = scales::comma)+
  theme_minimal()+
  theme(legend.position = "none")+
  labs(x = "t",
       y = "Behaviors Changed",
       title = "Total Behaviors Changed")  

set.seed(12)
bm_ttco2e <- replicate(
  100,
  {
behavior_model_df <- data.frame(t = 0:6, xijt = NA)

pi <- (1028282.215/100000000)*10
qj <- 31509.893/1028282.215
fj <- 0.5
xij0 <- runif(1, min = 0, max = 0.3)

behavior_model_ttco2e(behavior_model_df, 
pi=pi, qj=qj, xij0=xij0, fj=fj, dr=0.2, audience_size=100000000,
total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856), total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18))
  })

bm_ttco2e <-   data.frame(t = 0:6, bm_ttco2e)
bm_ttco2e_melt <- melt(bm_ttco2e ,  id.vars = 't', variable.name = 'tco2e_per_year')

ggplot(bm_ttco2e_melt, aes(t, value)) +
  geom_line(aes(colour = tco2e_per_year))+
  scale_colour_grey(start = 0.7, end = 0.7)+
  stat_summary(geom="line", fun = "mean", color="red")+
  scale_x_continuous(breaks=c(1,2,3,4,5,6), limits = c(0,6))+
  scale_y_continuous(labels = scales::comma)+
  theme_minimal()+
  theme(legend.position = "none")+
  labs(x = "t",
       y = "TCO2e",
       title = "Total TCO2e reduction") 

set.seed(12)
bm_ctco2e <- replicate(
  100,
  {
behavior_model_df <- data.frame(t = 0:6, xijt = NA)

pi <- (1028282.215/100000000)*10
qj <- 31509.893/1028282.215
fj <- 0.5
xij0 <- runif(1, min = 0, max = 0.3)

behavior_model_ctco2e(behavior_model_df, 
pi=pi, qj=qj, xij0=xij0, fj=fj, dr=0.2, audience_size=100000000,
total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856), total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18))
  })

bm_ctco2e <-   data.frame(t = 0:6, bm_ctco2e)
bm_ctco2e_melt <- melt(bm_ctco2e ,  id.vars = 't', variable.name = 'cost_per_tco2e')

ggplot(bm_ctco2e_melt, aes(t, value)) +
  geom_line(aes(colour = cost_per_tco2e))+
  scale_colour_grey(start = 0.7, end = 0.7)+
  stat_summary(geom="line", fun = "mean", color="red")+
  scale_x_continuous(breaks=c(1,2,3,4,5,6), limits = c(0,6))+
  scale_y_continuous(labels = scales::comma)+
  theme_minimal()+
  theme(legend.position = "none")+
  labs(x = "t",
       y = "USD",
       title = "Cost per TCO2e")

Parameter Elasticity

behavior_model_df <- data.frame(t = 0:6, xijt = NA)
initial_ctco2e <- behavior_model_ctco2e(behavior_model_df, 
pi=1028282.215/100000000, qj=31509.893/1028282.215, xij0=(1028282.215/100000000)*(31509.893/1028282.215), fj=0, dr=0.00, audience_size=100000000,
total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856), total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18))


#1 Attribution 
at_ctco2e <- behavior_model_ctco2e(behavior_model_df, 
                                   pi=1028282.215/100000000, qj=31509.893/1028282.215, xij0=(1028282.215/100000000)*(31509.893/1028282.215), fj=0, 
                                   dr=0.00, 
                                   audience_size=100000000,
                                   total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856),
                                   total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18),
                                   attribution = 100*1.01)  # 1% increase --> 101% attribution (?)
  
initial_o <- sum(initial_ctco2e)
at_final_o <- sum(at_ctco2e)
at_initial_i <- 100
at_final_i <- 101
at_el <- round(abs(((initial_o - at_final_o)/(at_initial_i-at_final_i))*(at_initial_i/initial_o)),2)



#2 Total Expenses
te_ctco2e <- behavior_model_ctco2e(behavior_model_df, 
                                   pi=1028282.215/100000000, qj=31509.893/1028282.215, xij0=(1028282.215/100000000)*(31509.893/1028282.215), fj=0, 
                                   dr=0.00, 
                                   audience_size=100000000,
                                   total_expenses=c(2161776*1.01, 2740457*1.01, 3618359*1.01, 5535243*1.01, 6897765*1.01, 8856389*1.01, 10729856*1.01),# 1% increase
                                   total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18)) 

te_final_o <- sum(te_ctco2e)
te_initial_i <- sum(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856)
te_final_i <- sum(2161776*1.01, 2740457*1.01, 3618359*1.01, 5535243*1.01, 6897765*1.01, 8856389*1.01, 10729856*1.01)
te_el <- round(abs(((initial_o - te_final_o)/(te_initial_i-te_final_i))*(te_initial_i/initial_o)),2) 

#3 Total Partners
tp_ctco2e <- behavior_model_ctco2e(behavior_model_df, 
                                   pi=1028282.215/100000000, qj=31509.893/1028282.215, xij0=(1028282.215/100000000)*(31509.893/1028282.215), fj=0, 
                                   dr=0.00, 
                                   audience_size=100000000,
                                   total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856),
                                   total_partners= c(3.5*1.01, 4*1.01, 4.5*1.01, 6.5*1.01, 9*1.01, 13*1.01, 18*1.01)) # 1% increase


tp_final_o <- sum(tp_ctco2e)
tp_initial_i <- sum(3.5, 4, 4.5, 6.5, 9, 13, 18)
tp_final_i <- sum(3.5*1.01, 4*1.01, 4.5*1.01, 6.5*1.01, 9*1.01, 13*1.01, 18*1.01)
tp_el <- round(abs(((initial_o - tp_final_o)/(tp_initial_i-tp_final_i))*(tp_initial_i/initial_o)),2) 

#4 TCO2e/behavior/year
  
  #4.1 Install Solar Panels (10%)
  sp_ctco2e <- behavior_model_ctco2e(behavior_model_df, 
                                   pi=1028282.215/100000000, qj=31509.893/1028282.215, xij0=(1028282.215/100000000)*(31509.893/1028282.215), fj=0, 
                                   dr=0.00, 
                                   audience_size=100000000,
                                   total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856),
                                   total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18),
                                   tco2e_solar_panel = 4.2*1.01)# 1% increase

  sp_final_o <- sum(sp_ctco2e)
  sp_initial_i <- 4.2
  sp_final_i <- 4.2*1.01
  sp_el <- round(abs(((initial_o - sp_final_o)/(sp_initial_i-sp_final_i))*(sp_initial_i/initial_o)),2)
  
   #4.2 Drive an Electric Vehicle (25%)
  ev_ctco2e <- behavior_model_ctco2e(behavior_model_df, 
                                   pi=1028282.215/100000000, qj=31509.893/1028282.215, xij0=(1028282.215/100000000)*(31509.893/1028282.215), fj=0, 
                                   dr=0.00, 
                                   audience_size=100000000,
                                   total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856),
                                   total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18),
                                   tco2e_electric_vehicule = 2.26*1.01)# 1% increase

  ev_final_o <- sum(ev_ctco2e)
  ev_initial_i <- 2.26
  ev_final_i <- 2.26*1.01
  ev_el <- round(abs(((initial_o - ev_final_o)/(ev_initial_i-ev_final_i))*(ev_initial_i/initial_o)),2)

  #4.3 Get Around Greener (10%)
    #Default value = 0 TCO2e/year
  
  #4.4 Eat Less Beef (10%)
  lb_ctco2e <- behavior_model_ctco2e(behavior_model_df, 
                                   pi=1028282.215/100000000, qj=31509.893/1028282.215, xij0=(1028282.215/100000000)*(31509.893/1028282.215), fj=0, 
                                   dr=0.00, 
                                   audience_size=100000000,
                                   total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856),
                                   total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18),
                                   tco2e_less_beef = 0.81*1.01)# 1% increase

  lb_final_o <- sum(lb_ctco2e)
  lb_initial_i <- 0.81
  lb_final_i <- 0.81*1.01
  lb_el <- round(abs(((initial_o - lb_final_o)/(lb_initial_i-lb_final_i))*(lb_initial_i/initial_o)),2)
  
  #4.5 Eat More Veggies (20%)
  mv_ctco2e <- behavior_model_ctco2e(behavior_model_df, 
                                   pi=1028282.215/100000000, qj=31509.893/1028282.215, xij0=(1028282.215/100000000)*(31509.893/1028282.215), fj=0, 
                                   dr=0.00, 
                                   audience_size=100000000,
                                   total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856),
                                   total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18), 
                                   tco2e_more_veggies = 0.077*1.01)# 1% increase

  mv_final_o <- sum(mv_ctco2e)
  mv_initial_i <- 0.077
  mv_final_i <- 0.077*1.01
  mv_el <- round(abs(((initial_o - mv_final_o)/(mv_initial_i-lb_final_i))*(mv_initial_i/initial_o)),2)
  
  #4.6 Reduce Food Waste (20%)
  fw_ctco2e <- behavior_model_ctco2e(behavior_model_df, 
                                   pi=1028282.215/100000000, qj=31509.893/1028282.215, xij0=(1028282.215/100000000)*(31509.893/1028282.215), fj=0, 
                                   dr=0.00, 
                                   audience_size=100000000,
                                   total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856),
                                   total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18),
                                   tco2e_reduce_food_waste = 0.026*1.01)# 1% increase

  fw_final_o <- sum(fw_ctco2e)
  fw_initial_i <- 0.026
  fw_final_i <- 0.026*1.01
  fw_el <- round(abs(((initial_o - fw_final_o)/(fw_initial_i-fw_final_i))*(fw_initial_i/initial_o)),2)
  
  #4.7 Give to Nature (5%)
  gn_ctco2e <- behavior_model_ctco2e(behavior_model_df, 
                                   pi=1028282.215/100000000, qj=31509.893/1028282.215, xij0=(1028282.215/100000000)*(31509.893/1028282.215), fj=0, 
                                   dr=0.00, 
                                   audience_size=100000000,
                                   total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856),
                                   total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18),
                                   tco2e_give_to_nature = 2*1.01)# 1% increase

  gn_final_o <- sum(gn_ctco2e)
  gn_initial_i <- 2
  gn_final_i <- 2*1.01
  gn_el <- round(abs(((initial_o - gn_final_o)/(gn_initial_i-gn_final_i))*(gn_initial_i/initial_o)),2)
  

  
#5 Audience Size
az_ctco2e <-behavior_model_ctco2e(behavior_model_df, 
                                   pi=1028282.215/100000000, qj=31509.893/1028282.215, xij0=(1028282.215/100000000)*(31509.893/1028282.215), fj=0, 
                                   dr=0.00, 
                                   audience_size=100000000*1.01, # 1% increase
                                   total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856),
                                   total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18))

  az_final_o <- sum(az_ctco2e)
  az_initial_i <- 100000000
  az_final_i <- 100000000*1.01
  az_el <- round(abs(((initial_o - az_final_o)/(az_initial_i-az_final_i))*(az_initial_i/initial_o)),2)
  

#6 pi
pi_ctco2e <-behavior_model_ctco2e(behavior_model_df, 
                                  pi=(1028282.215/100000000)*1.01,  # 1% increase
                                  qj=31509.893/1028282.215, xij0=(1028282.215/100000000)*(31509.893/1028282.215), fj=0, 
                                  dr=0.00, 
                                  audience_size=100000000, # 1% increase
                                  total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856),
                                  total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18))

  pi_final_o <- sum(pi_ctco2e)
  pi_initial_i <- 1028282.215/100000000
  pi_final_i <- (1028282.215/100000000)*1.011
  pi_el <- round(abs(((initial_o - pi_final_o)/(pi_initial_i-pi_final_i))*(pi_initial_i/initial_o)),2)
  
#7 qj  
qj_ctco2e <-behavior_model_ctco2e(behavior_model_df, 
                                  pi=1028282.215/100000000,
                                  qj=(31509.893/1028282.215)*1.01, # 1% increase
                                  xij0=(1028282.215/100000000)*(31509.893/1028282.215), fj=0, 
                                  dr=0.00, 
                                  audience_size=100000000, # 1% increase
                                  total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856),
                                  total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18))

  qj_final_o <- sum(qj_ctco2e)
  qj_initial_i <- 31509.893/1028282.215
  qj_final_i <- (31509.893/1028282.215)*1.01
  qj_el <- round(abs(((initial_o - qj_final_o)/(qj_initial_i-qj_final_i))*(qj_initial_i/initial_o)),2) 
  
  
#8 xij0  
xij0_ctco2e <-behavior_model_ctco2e(behavior_model_df, 
                                  pi=1028282.215/100000000,
                                  qj=31509.893/1028282.215,
                                  xij0=((1028282.215/100000000)*(31509.893/1028282.215))*1.01, # 1% increase 
                                  fj=0, 
                                  dr=0.00, 
                                  audience_size=100000000, # 1% increase
                                  total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856),
                                  total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18))

  xij0_final_o <- sum(xij0_ctco2e)
  xij0_initial_i <- (1028282.215/100000000)*(31509.893/1028282.215)
  xij0_final_i <- ((1028282.215/100000000)*(31509.893/1028282.215))*1.01
  xij0_el <- round(abs(((initial_o - xij0_final_o)/(xij0_initial_i-xij0_final_i))*(xij0_initial_i/initial_o)),2)
  
#10 fj  
fj_initial_o <-behavior_model_ctco2e(behavior_model_df, 
                                  pi=1028282.215/100000000,
                                  qj=31509.893/1028282.215,
                                  xij0=(1028282.215/100000000)*(31509.893/1028282.215), 
                                  fj=0.1, 
                                  dr=0.00, 
                                  audience_size=100000000, 
                                  total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856),
                                  total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18))

fj_tco2e <-behavior_model_ctco2e(behavior_model_df, 
                                  pi=1028282.215/100000000,
                                  qj=31509.893/1028282.215,
                                  xij0=(1028282.215/100000000)*(31509.893/1028282.215), 
                                  fj=0.1*1.01, # 1% increase
                                  dr=0.00, 
                                  audience_size=100000000, 
                                  total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856),
                                  total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18))

  fj_final_o <- sum(fj_tco2e)
  fj_initial_o <- sum(fj_initial_o)
  fj_initial_i <- 0.1
  fj_final_i <- 0.1*1.01
  fj_el <- round(abs(((fj_initial_o - fj_final_o)/(fj_initial_i-fj_final_i))*(fj_initial_i/fj_initial_o)),2)
  
#11 dr  
dr_initial_o <-behavior_model_ctco2e(behavior_model_df, 
                                  pi=1028282.215/100000000,
                                  qj=31509.893/1028282.215,
                                  xij0=(1028282.215/100000000)*(31509.893/1028282.215), 
                                  fj=0, 
                                  dr=0.05, 
                                  audience_size=100000000, 
                                  total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856),
                                  total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18))

dr_tco2e <-behavior_model_ctco2e(behavior_model_df, 
                                  pi=1028282.215/100000000,
                                  qj=31509.893/1028282.215,
                                  xij0=(1028282.215/100000000)*(31509.893/1028282.215), 
                                  fj=0, 
                                  dr=0.05*1.01, # 1% increase
                                  audience_size=100000000, 
                                  total_expenses=c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856),
                                  total_partners= c(3.5, 4, 4.5, 6.5, 9, 13, 18))

  dr_final_o <- sum(dr_tco2e)
  dr_initial_o <- sum(dr_initial_o)
  dr_initial_i <- 0.05
  dr_final_i <- 0.05*1.01
  dr_el <- round(abs(((dr_initial_o - dr_final_o)/(dr_initial_i-dr_final_i))*(dr_initial_i/dr_initial_o)),2)

  
# Table
par <- c("Attribution (%)", "Total Expenses ($)", "Total Partners", "Solar Panels (TCO2e/y)r", "Electric Vehicles (TCO2e/yr)", "Less Beef (TCO2e/yr)", "More Veggies (TCO2e/yr)", "Food Waste (TCO2e/yr)", "Give to Nature (TCO2e/yr)", "Audience Size", "pi", "qj", "xij0", "fj", "Discount Rate")

default_v <- list(100, c(2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856), c(3.5, 4, 4.5, 6.5, 9, 13, 18), 4.2, 2.26, 0.81, 0.077, 0.026, 2, 100000000, 0.01, 0.03, 0.000315, 0, 0)

el <- c(at_el, te_el, tp_el, sp_el, ev_el, lb_el, mv_el, fw_el, gn_el, az_el, pi_el, qj_el, xij0_el, fj_el, dr_el)

tab <- data.frame(par=par, def=cbind(default_v), el=el)
tab <-tab[order(tab$el,decreasing=TRUE),] 
kbl(tab, row.names = FALSE, col.names = c("Parameter", "Default", "Elasticity")) %>% 
  kable_styling()
Parameter Default Elasticity
Total Expenses ($) 2161776, 2740457, 3618359, 5535243, 6897765, 8856389, 10729856 1.00
Attribution (%) 100 0.99
Total Partners 3.5, 4.0, 4.5, 6.5, 9.0, 13.0, 18.0 0.99
Audience Size 1e+08 0.99
qj 0.03 0.87
pi 0.01 0.79
Electric Vehicles (TCO2e/yr) 2.26 0.47
Solar Panels (TCO2e/y)r 4.2 0.35
Discount Rate 0 0.15
xij0 0.000315 0.12
fj 0 0.09
Give to Nature (TCO2e/yr) 2 0.08
Less Beef (TCO2e/yr) 0.81 0.07
More Veggies (TCO2e/yr) 0.077 0.00
Food Waste (TCO2e/yr) 0.026 0.00