This document illustrates the use of various functions within the `decisionSupport`

R package (Luedeling, Goehring, and Schiffers 2019). We use these to develop a cost-benefit model and implement a Monte Carlo simulation with the example of agroforestry interventions promoted in Northwest Vietnam (see Do et al.Â 2020). We generate a model to calculate the Net Present Value of proposed agroforestry systems as compared to the common conventional monoculture cropping systems of maize.

We use the `decisonSupport()`

function to perform probabilistic Monte Carlo simulations. The function performs calculations using distributions and ranges which is a major advantage over traditional spreadsheet calculations. In a Monte Carlo simulation, a particular calculation is repeated many times. Each time, a random draw from defined distribution of each variable are used as inputs for calculating the function. The output then consists of a set of values instead of single estimates.

The function requires two arguments:

a data table (.csv format) specifying the names and distributions for all uncertain variables and,

an R function that predicts decision outcomes based on the variables named in the data table.

The `chance_event()`

function calculates the time series occurrence of particular events (risk events in this case). This function requires two main arguments:

the probability that a particular risk event will happen in a given year and

the time line for the simulation.

The `gompertz_yield()`

function calculates the yields of tree crops using identified parameters given an assumed pattern.

The value varier function `vv()`

produces realistic time series of data that includes year-to-year variation around a specified mean. This is determined by a desired coefficient of variation. Both mean and coefficient of variation are provided to the model via the input table.

The `discount()`

function calculated NPV in the last step of the decision modeling procedure. The function implements discounting of future costs or benefits by a decision maker-specific discount rate, which is also specified in the input tables.

The cost-benefit model `AF_benefit`

calculates the NPV:

```
AF_benefit <- function(x, varnames)
{
#Define each variable as vectors of 25 values corresponding to 25 years of simulation
mono_costs <- rep(0, n_years) #the annual cost of maize monoculture
establisment_cost <- rep(0, n_years) #costs for setting up agroforestry systems (AF)
maintenance_cost <- rep(0, n_years) #the annual cost spent on AF
tree_yield <- rep(0, n_years) #the yield of longan
AF_maize <- rep(0, n_years) #the yield of maize in AF
mono_maize <- rep(0, n_years) #the yield of maize in monocultures
grass_benefit <- rep(0, n_years) #the benefit from grass
reduced_soil_loss <- rep(0, n_years) #amount of soil loss reduced in AF compared to monoculture systems
#Simulate the chance for risk events to occur during the simulation period
tree.risk <- chance_event(tree_risk, value_if = 1, n = n_years)
maize.risk <- chance_event(maize_risk, value_if = 1, n = n_years)
#Calculate system costs and convert currency from VND to USD using currency
#exchange rate (cur_change) applied for any variable with monetary value
#(i.e. costs and benefits)
maintenance_cost[1] <- 0
maintenance_cost[2:4] <- vv(main_cost_stage1, CV_cost, 3)/cur_change
maintenance_cost[5:10] <- vv(main_cost_stage2, CV_cost, 6)/cur_change
maintenance_cost[11:25] <- vv(main_cost_stage3, CV_cost, 15)/cur_change
establisment_cost[1] <- establis_cost/cur_change
establisment_cost[2:25] <- 0
##Calculate system benefits
#benefit of longan fruit
tree_yield <- gompertz_yield(max_harvest = max_tree_harvest,
time_to_first_yield_estimate = time_first_tree_est,
time_to_second_yield_estimate = time_sec_tree_est,
first_yield_estimate_percent = first_tree_est_per,
second_yield_estimate_percent = sec_tree_est_per,
n_years = n_years,
var_CV = CV_tree_yield,
no_yield_before_first_estimate = FALSE)
tree_yield <- tree_yield*(1-tree.risk*yield_tree_risk)
tot_tree_yield <- tree_yield*num_of_tree
tree_benefit <- (tot_tree_yield*tree_price)/cur_change
#benefit of maize in AF systems
time <- 1:n_years
decay_speed_AF <- -log(1-decay_rate_AF)
AF_maize <- min_AFmaize+(max_AFmaize-min_AFmaize)*exp(-decay_speed_AF*(time-1))
tot_AF_maize <- vv(AF_maize*(1-maize.risk*yield_maize_risk), CV_maize_yield, n_years)
AF_maize_benefit <- (tot_AF_maize*vv(maize_price, CV_maize_price, n_years))/cur_change
#benefit of maize in monoculture
decay_speed_mono <- -log(1-decay_rate_mono)
mono_maize <- max_monomaize*exp(-decay_speed_mono*(time-1))
tot_mono_maize <- vv(mono_maize*(1-maize.risk*yield_maize_risk), CV_maize_yield, n_years)
mono_revenue <- (tot_mono_maize*vv(maize_price, CV_maize_price, n_years))/cur_change
mono_costs <- vv(mono_cost, CV_cost, n_years)/cur_change
mono_benefit <- mono_revenue-mono_costs
#benefit of grass
grass_benefit <- vv(grass_profit, CV_grass_profit, 25)/cur_change
#benefit of soil erosion control
reduced_soil_loss[1] <- 0
reduced_soil_loss[2:10] <- reduced_percents1*mono_soil_loss
reduced_soil_loss[10:25] <- reduced_percents2*mono_soil_loss
soil_saved <- vv(reduced_soil_loss, CV_eros_control, n = n_years)
erosion_benefit <- (soil_saved*soil_saved_payment)/cur_change
# Calculate NPVs of AF and monoculture and NPV of AF vs. monocultures of maize
AF_revenue <- tree_benefit + AF_maize_benefit + grass_benefit+erosion_benefit
AF_benefit <- AF_revenue-maintenance_cost-establisment_cost
cash_flow <- discount(AF_benefit,
discount_rate = discount_rate,
calculate_NPV = FALSE)
#cummulative cash flow of agroforestry system
cum_cash_flow <- cumsum(cash_flow)
#NPV of monoculture system
NPV_mono <- discount(mono_benefit, discount_rate = discount_rate,
calculate_NPV = TRUE )
#Benefit of choosing agroforestry over maize monoculture
tradeoff_benefit <- AF_benefit-mono_benefit
#NPV of AF system
NPV_AF <- discount(AF_benefit,
discount_rate = discount_rate,
calculate_NPV = TRUE)
#final NPV of the decision to choose AF or monocultures
NPV_tradeoff <- discount(tradeoff_benefit,
discount_rate = discount_rate,
calculate_NPV = TRUE )
#in the return list, one could indicate any outcome that they wish to see from the model
return(list(trade_off = NPV_tradeoff,
NPV_maize = NPV_mono,
NPV = NPV_AF))
}
```

The probabilistic simulation was implemented using the `decisionSupport`

function to return a Monte Carlo data table with 10, 000 values for each of the input and response variables in the return list. A VIP score generated for each input variable of the PLS regression indicate its influenctial importance contriduting to model outcome.

```
decisionSupport(inputFilePath = paste("Longan+maize+grass.csv", sep = ""),
outputPath = paste("MCResults", sep = ""),
write_table = TRUE,
welfareFunction = AF_benefit,
numberOfModelRuns = 10000,
functionSyntax = "plainNames")
```

The outcome table is then used to calculate the Value of Information for uncertain variables using the `multi_EVPI`

function in the `decisionSupport`

package (Luedeling, Goehring, and Schiffers 2019):

The `multi_EVPI()`

function calculates the expected value of perfect information (EVPI) for multiple variables in the model using the simulation outcome table from the `decisionSupport()`

function.

```
MCall <- read.table(file = "MCResults/mcSimulationResults.csv",
header = TRUE, sep = ",")
multi_EVPI(MCall, "trade_off", write_table = TRUE)
```

The full NPV distributions of agroforestry and maize monoculture acquired from MC simulation look like this:

VIP scores from the sensitivity analysis (left) and Value of Information analysis (right):

Luedeling, Eike, Lutz Goehring, and Katja Schiffers. 2019. *DecisionSupport: Quantitative Support of Decision Making Under Uncertainty*. http://www.worldagroforestry.org/.