Monte Carlo Simulation using R Code
Isaiah McCommons
9/14/2018
Example 1 - Simple DCF Calculation
In this section, the reproducible R code will be used to simple DCF calculation with sample data which will result in the display of NPV, IRR, PBP, and DPBP
Setup cf dataframe
cf_df <- data.frame(
Project1 = c(-1000, 1250, 10, 10, 20, 20),
Project2 = c(-1000, -10, 0, 10, 20, 2000)
)Calculate DCF and CF
dcf <- function(x, r, t0=FALSE){
# calculates discounted cash flows (DCF) given cash flow and discount rate
#
# x - cash flows vector
# r - vector or discount rates, in decimals. Single values will be recycled
# t0 - cash flow starts in year 0, default is FALSE, i.e. discount rate in first period is zero.
if(length(r)==1){
r <- rep(r, length(x))
if(t0==TRUE){r[1]<-0}
}
x/cumprod(1+r)
}
npv <- function(x, r, t0=FALSE){
# calculates net present value (NPV) given cash flow and discount rate
#
# x - cash flows vector
# r - discount rate, in decimals
# t0 - cash flow starts in year 0, default is FALSE
sum(dcf(x, r, t0))
}
pbp <- function(x, ...){
# calculates payback period (PBP)
#
# x - cash flows vector
# ... - ignored
i <- match(1, sign(cumsum(x)))
i-2+(-cumsum(x)[i-1]/x[i])
}
dpbp <- function(x, r, t0=FALSE){
# calculates discounted payback period (DPBP) given cash flow and discount rate
#
# x - cash flows vector
# r - discount rate, in decimals
# t0 - cash flow starts in year 0, default is FALSE
pbp(dcf(x, r, t0))
}
irr <- function(x, t0=FALSE, ...){
# calculates internal rate of return (IRR) given cash flow
#
# x - cash flows vector
# t0 - cash flow starts in year 0, default is FALSE
tryCatch(uniroot(f=function(i){sum(dcf(x, i, t0))},
interval=c(0,1))$root,
error=function(e) return(NA)
)
}Display the results of the calculation
# Combine these functions with the power of dplyr and
# you will have a toolbox capable of competing with the CFO's favorite spreadsheet:
cf_df %>%
summarise_all(funs(NPV=npv,PBP=pbp, DPBP=dpbp, IRR=irr), r=0.06, t0=TRUE) %>%
gather(key=key, value = value) %>%
separate(key, into = c("Project", "Metric")) %>%
spread(key=Project, value=value) %>%
mutate(Metric=fct_relevel(Metric, "NPV", "IRR", "PBP", "DPBP"),
Metric=fct_recode(Metric,
`Net Present Value (NPV), USD mln`="NPV",
`Internal Rate of Return (IRR), %`="IRR",
`Payback Period, years` = "PBP",
`Discounted Payback Period, years`="DPBP")) %>%
arrange(as.numeric(Metric))Example 2 - Simple Single Integer Calculation
This section will calculate a single integer or a single-point estimate. The standard DCF analysis provides financial managers single-point estimates associated with the Internal Rate of Return (IRR) or even the Present Value (PV) of the cash flow of a firm (Kim & Reinschmidt, 2012).
Calculate DCF using a single integer of NPV.
options(scipen = 999)
start.time<-Sys.time() ## start stopwatch
mydata<-lapply(1:1000, function(i) {
## Net Profit = Sales Volume * (Selling Price - Unit cost) - Fixed costs
## create random uniform integer for SalesVolume
DU = sample(x=1:3,size = 1,replace = T)
if(DU==1){
SalesVolume=100000
SellingPrice=8
}
if(DU==2){
SalesVolume=75000
SellingPrice=10
}
if(DU==3){
SalesVolume=50000
SellingPrice=11
}
## create triangular distribution for unit cost (as per example)
library(triangle)
Unitcost = rtriangle(1,5.5,7.5,c = 6.5) ## c=mode of distribution
#### hist(rtriangle(10000000,5.5,7.5))
Fixedcosts = 120000
NetProfit = SalesVolume * (SellingPrice - Unitcost) - Fixedcosts
})
stop.time<-Sys.time() ## stop stopwatch
time.taken.lapply = stop.time-start.time
NPV = unlist(mydata)Plot histogram of the results of NPV
hist(NPV)
Find the mean of the results
mean(NPV)## [1] 92883.73Summerize the results
summary(NPV)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -68673 57829 99523 92884 130777 212722Example 3 - Monte Carlo Simulation
In this example, a Monte Carlo simulation will calculate the NPV 10,000 times and find the distribution of the NPV calculations. Monte Carlo can be used to compute the distributions of project value (Balcombe & Smith, 1999).
Setting up loop for Monte Carlo simulation 10,000 times. For the purpose of brevity the R code for the calculation will not be displayed.
start.time.fl<-Sys.time() ## start stopwatch
NetProfits=vector(mode="numeric") ## mind the "s"
SalesVolumes = vector(mode="numeric")
UnitCosts = vector(mode="numeric")
SellingPrices = vector(mode="numeric")
for(i in 1:10000) {
## Net Profit = Sales Volume * (Selling Price - Unit cost) - Fixed costs
## create random uniform integer for SalesVolume
DU = sample(x=1:3,size = 1,replace = T)
if(DU==1){
SalesVolume=100000
SellingPrice=8
}
if(DU==2){
SalesVolume=75000
SellingPrice=10
}
if(DU==3){
SalesVolume=50000
SellingPrice=11
}
## create triangular distribution for unit cost (as per example)
library(triangle)
Unitcost = rtriangle(1,5.5,7.5,c = 6.5) ## c=mode of distribution
#### hist(rtriangle(10000000,5.5,7.5))
Fixedcosts = 120000
NetProfit = SalesVolume * (SellingPrice - Unitcost) - Fixedcosts
NetProfits[i]=NetProfit
SalesVolumes[i] = SalesVolume
UnitCosts[i] = Unitcost
SellingPrices[i] = SellingPrice
}
stop.time.fl<-Sys.time() ## start stopwatch
time.taken.fl <-stop.time.fl - start.time.fl
time.taken.fl## Time difference of 0.6165287 secstime.taken.lapply## Time difference of 0.07970095 secstime.taken.fl - time.taken.lapply## Time difference of 0.5368278 secsPlot the results
par(mfrow = c(2,2))
hist(SalesVolumes)
hist(SellingPrices)
hist(UnitCosts)
hist(NetProfits)
abline(v=mean(NetProfits), col=2)
Summerize the results
par(mfrow = c(1,1))
summary(SalesVolumes)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 50000 50000 75000 74450 100000 100000summary(SellingPrices)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 8.0 8.0 10.0 9.7 11.0 11.0summary(UnitCosts)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.511 6.198 6.496 6.494 6.787 7.497summary(NetProfits)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -69598 62355 102617 93669 132701 215997