Simulated NPV with Time-Series
Load Environment
Let's start by making sure all the tools we nned are loaded up. You will almost always need to load these packages up anytime you are build a Monte-Carlo model.
using Distributions
using Dates
using Gadfly
using StatsBase
using MCHammer
using DataFrames
Setup Inputs and Outputs
The next critical step is to setup key inputs, arrays and other important model parameters.
#Setup the Date Range for the analysis
dr = collect(Date(2019,1,01):Dates.Year(1):Date(2023,01,01))
#Setup Global Inputs
ForecastYrs = 5
Trials = 10000
Units = [5000, 10000, 17000, 18000, 20000]
InitialInvestment = 250000
Investment = [100000, 0, 0, 25000,0] #fill(0,ForecastYrs)
#Setup Outputs
Sensitivity_Tbl = []
ProjectNPV = []
USP = []
USC =[]
DR = []
OP =[]
Annual_CashFlows =[]
Any[]
Build Simulation Model
Monte-Carlo simulation needs to generate a table of scenarios which are known as ***Trials***. A trial documents, in the form of a row, all of the inputs and calculated outputs for a particular scenario.
Using this results table allows you to runs all sorts of analysis, including sensitivity analysis and assigning probabilities to outcomes. To generate this table, you need to loop your equation/function as many times as you need and vary the inputs using probability distributions.
Another challenge to account for is that our example is a 5yr NPV model which requires building and analyzing the results over multiple periods. To extend the model, we are using MCHammer's GBMM function that allows to project a random walk forecast over how ever many periods you need, which extends automatically the model in Julia.
for i = 1:Trials
UnitSellPrice = GBMM(80, 0.2, 0.1, ForecastYrs)
UnitCost = GBMM(40, 0.1, 0.05, ForecastYrs)
#Each period the discount rate is independent. If you use an additive method instead of multiplicative, you can end up with differences. These may or may not impact the decision. For simulation it is best to use the risk free rate.
#Multiplicative Method
DiscountRate = cumprod(rand(Normal(0.02,0.0075),ForecastYrs)+fill(1,ForecastYrs))#accumulate(+,rand(Normal(0.02,0.0075),ForecastYrs))+fill(1,ForecastYrs)
#Additive Method
#DiscountRate = accumulate(+,rand(Normal(0.02,0.0075),ForecastYrs))+fill(1,ForecastYrs)
#a static DR
#DiscountRate = accumulate(+, fill(0.02,ForecastYrs))+fill(1,ForecastYrs)
#print(DiscountRate)
#DCF Elements
Annual_Sales = UnitSellPrice .* Units
Annual_COGS = UnitCost .* Units
OPEX = rand(TriangularDist(.2,0.5,0.35),ForecastYrs) .* Annual_Sales
#Constant Dollar Cashflow
#CashFlow_C = (Annual_Sales - Annual_COGS - OPEX - Investment)
#Discounted CashFLow over multpile periods. This function uses arrays and DOT functions.
CashFlow = (Annual_Sales - Annual_COGS - OPEX - Investment) ./ DiscountRate
#Calculated Output
Trial_NPV = sum(CashFlow)-InitialInvestment
#Convert Arrays to Scalars for sensitivity analysis
push!(ProjectNPV, Trial_NPV)
push!(USC, mean(UnitCost))
push!(USP, mean(UnitSellPrice))
push!(DR, mean(DiscountRate))
push!(OP, mean(OPEX))
push!(Annual_CashFlows, CashFlow)
end
Setting up data for analysis and charting
Setup inputs/outputs(above) and output tables (below) for sensitivity analysis and charting. Since correlation is based on the same math as regression, the only way to calculate sensitivity on an Array > 1 (in this case multiple years) is to condense the array into a scalar value using either mean, sum or any other transform because what ever you pick will generate a similar or identical result.
Sensitivity_Tbl = DataFrame(hcat(ProjectNPV, USC, USP, DR, OP))
names!(Sensitivity_Tbl, [:ProjectNPV, :USC, :USP, :DR, :OP])
NPV_Sensitivity = cormat(Sensitivity_Tbl,1)
5×5 Array{Float64,2}: 1.0 -0.268443 0.903869 -0.0169668 0.602936 -0.268443 1.0 -0.0100426 -0.00161571 -0.00494222 0.903869 -0.0100426 1.0 0.0116294 0.835425 -0.0169668 -0.00161571 0.0116294 1.0 0.00784858 0.602936 -0.00494222 0.835425 0.00784858 1.0
Stats
Generate model results and list all the outputs for your charting and analysis functions.
print("Project Mean: ", mean(ProjectNPV),"\n")
print("Project Std.Dev: ", std(ProjectNPV),"\n")
print("Prob. of Neg. NPV: ", GetCertainty(ProjectNPV,0,0),"\n")
print("NPV p10, p50, p90 : ", quantile(collect(Float64, ProjectNPV),[0.1,0.5,0.9]),"\n")
println("")
println("OUTPUTS: Annual_CashFlows, ProjectNPV, Sensitivity_Tbl")
println("date range = dr")
Project Mean: 2.5276551743783243e6 Project Std.Dev: 1.0689942031448558e6 Prob. of Neg. NPV: 0.0016 NPV p10, p50, p90 : [1.2033067273826795e6, 2.463327457125737e6, 3.9152107803786374e6] OUTPUTS: Annual_CashFlows, ProjectNPV, Sensitivity_Tbl date range = dr
To generate a complete list of percentiles, use the fractiles()
.
fractiles(ProjectNPV)
11×2 Array{Any,2}: "P0.0" -590221.0 "P10.0" 1.20331e6 "P20.0" 1.60518e6 "P30.0" 1.93304e6 "P40.0" 2.19838e6 "P50.0" 2.46333e6 "P60.0" 2.72835e6 "P70.0" 3.03693e6 "P80.0" 3.38385e6 "P90.0" 3.91521e6 "P100.0" 7.52754e6
##Visualizing the outputs
Looking at the Probability Distribution
histogram_chrt(ProjectNPV, "Five Year NPV")
density_chrt(ProjectNPV, "Five Year NPV")
What variables are most influential on my output distribution?
sensitivity_chrt(Sensitivity_Tbl, 1, 3)
What does my cashflow look like over time?
The trend chart is a median centered chart that establishes a 90% confidence interval for each period. Remember dr or the date range is specified at the top.
CashFlow forecast
trend_chrt(Annual_CashFlows, dr)