# GY472
# Montecarlo simulation for risk assessment
# January 2024

#install.packages("jrvFinance")

##########################################
# 0. Setup environment
##########################################

#log_file <- paste0(work_dir, "log.txt")
#sink(log_file, split=TRUE)

# Import libraries
library("jrvFinance")
## Warning: package 'jrvFinance' was built under R version 4.3.2
library(ggplot2)
##########################################
# 1. Enter inputs (based on informed assumptions)
##########################################

##########
# Investment horizon
inv_t <- 5

# Time used for terminal value (perpetuity)
tv_t <- inv_t + 1

##########
# Inputs per year + terminal value

# Purchase price
stamp_duty_legal = (1+0.055)

purchase_price <- (950*1000*1000)* stamp_duty_legal


tax_rate = 0.25

# Yearly debt interest - Fixed rate
debt_interest <- 0.045 # 4.5%

# Debt LTV ratio
debt_ltv <- 0.7 # 60%

# Total rent income (NOI)
noi <- c(70293600, 70293600, 70293600, 70293600, 70293600, 70293600)
# Vacancy rate (v)
v <- c(0.0237, 0.0237, 0.0237, 0.0237, 0.0237, 0.0237)
# Growth rate (g)
g <- c(0.005, 0.005, 0.005, 0.005, 0.005, 0.005)
# Discount rate (r)
r <- c(0.0646,0.0643, 0.0643,   0.0641, 0.0643, 0.0641)


##########################################
# 2. DCF model
##########################################


##########
# Per year compounded calculations
compounded_g <- c()
compounded_r <- c()

for (t in 1:tv_t) {
  # Compounded growth factor (1+g_t)_t
  if (t == 1) {
    compounded_g[t] <- (1 + g[t])
  } else {
    compounded_g[t] <- (1 + g[t])*compounded_g[t-1]
  }
  
  # Compounded discount factor (1+g_t)_t
  if (t == 1) {
    compounded_r[t] <- (1 + r[t])
  } else {
    compounded_r[t] <- (1 + r[t])*compounded_r[t-1]
  }
}

########## 
# BT cashflows (CF) 
# IMPORTANT - This uses operations over vectors 
 bt_cashflows <- noi * (1-v) * compounded_g

# Terminal value (TV)
# Simplified Gordon growth model
# We use year 6 in our vectors as the values used to take the calculations to infinite (perpetual ownership)
tv <- bt_cashflows[inv_t+1] / (r[inv_t+1] - g[inv_t+1])

# Add TV to cashflows of the lat year of investment
bt_cashflows[inv_t] <- bt_cashflows[inv_t] + tv

Debt calculations

# Loan parameters
loan_amount <- purchase_price * debt_ltv  # $540,000,000
  debt_balance <- loan_amount

interest_rate <- 0.07  # 4.5%
fixed_capital_payment <- 15000000  # $15,000,000
loan_term <- inv_t  # 5 years

# Initialize variables
beginning_balance <- loan_amount
repayment_schedule <- data.frame(Year = 1:loan_term,
                                 Beginning_Balance = numeric(loan_term),
                                 Interest = numeric(loan_term),
                                 Fixed_Capital_Payment = rep(fixed_capital_payment, loan_term),
                                 Total_Payment = numeric(loan_term),
                                 Ending_Balance = numeric(loan_term))

# Calculate repayment schedule
for (i in 1:(loan_term - 1)) {
  repayment_schedule$Beginning_Balance[i] <- beginning_balance
  interest <- beginning_balance * interest_rate
  repayment_schedule$Interest[i] <- interest
  total_payment <- fixed_capital_payment + interest
  repayment_schedule$Total_Payment[i] <- total_payment
  ending_balance <- beginning_balance - fixed_capital_payment
  beginning_balance <- ending_balance
  repayment_schedule$Ending_Balance[i] <- ending_balance
}

# Adjust the 5th row
repayment_schedule$Beginning_Balance[loan_term] <- beginning_balance
interest <- beginning_balance * interest_rate
repayment_schedule$Interest[loan_term] <- interest
total_payment <- beginning_balance + interest  # Total payment includes remaining balance and interest
repayment_schedule$Total_Payment[loan_term] <- total_payment
repayment_schedule$Fixed_Capital_Payment[loan_term] <- beginning_balance
repayment_schedule$Ending_Balance[loan_term] <- 0
########## 
# AT cashflows (CF) 

# Debt interest yearly payments
debt_cashflows <- -repayment_schedule$Total_Payment[1:5]


# Add debt payments to BT cashflows for years 1 to the last
at_cashflows <- (bt_cashflows[1:inv_t] + debt_cashflows)

##########  
# Discounted cashflows (CF)
discounted_cashflows <- at_cashflows / compounded_r[1:inv_t]

# Present value
pv <- sum(discounted_cashflows)

# Year 0 cashflows
year_0 <- purchase_price - debt_balance


# IRR 
inv_cashflows <- c(-year_0, at_cashflows[1:inv_t])
inv_irr <- irr(inv_cashflows)


# Print DCF model outputs
message("PV: ", pv)
## PV: 448434890.807701
message("IRR: ", inv_irr)
## IRR: 0.155494895013448
##########################################
# 2. Monte Carlo simulations
##########################################
# Number of simulations
simulation_size <- 10000

# For generation of consistent pseudo-random numbers on computers, the seed value does not matter in this file
set.seed(200)

########## 
# We use random shocks based on assumed distributions for for each input
# Each shock adds a random epsilon value to the assumed input

# Vacancy rate - Uniform distribution around 0 ~ U(min, max)
v_min <- -0.001
v_max <- 0.001
v_epsilon <- runif(simulation_size, v_min, v_max)

# Growth rate - Normal distribution centered at 0 ~ N(0, sd)
g_mu <- 0
g_sd <- 0.001
g_epsilon <- rnorm(simulation_size, g_mu, g_sd)

# Discount rate - Normal distribution centered at 0 ~ N(0, sd)
r_mu <- 0
r_sd <- 0.001
r_epsilon <- rnorm(simulation_size, r_mu, r_sd)


##########
# For the simulations we run the DCF model as before, but adding a different 
# combination of random shocks
##########
pv_simulation_results <- c()
irr_simulation_results <- c()

for (s in 1:simulation_size) {
  
  # Simulation values using the random shocks
  v_simulation <- v + v_epsilon[s]
  g_simulation <- g + g_epsilon[s]
  r_simulation <- r + r_epsilon[s]
  
  # Per year compounded calculations
  compounded_g_simulation <- c()
  compounded_r_simulation <- c()
  
  for (t in 1:tv_t) {
    # Compounded growth factor (1+g_t)_t
    if (t == 1) {
      compounded_g_simulation[t] <- (1 + g_simulation[t])
    } else {
      compounded_g_simulation[t] <- (1 + g_simulation[t])*compounded_g_simulation[t-1]
    }
    
    # Compounded discount factor (1+g_t)_t
    if (t == 1) {
      compounded_r_simulation[t] <- (1 + r_simulation[t])
    } else {
      compounded_r_simulation[t] <- (1 + r_simulation[t])*compounded_r_simulation[t-1]
    }
  }

  ########## 
  # Cashflows (CF) 
  # IMPORTANT - This uses operations over vectors 
  bt_cashflows_simulation <- noi * (1-v_simulation) * compounded_g_simulation
  
  # Terminal value (TV)
  # Simplified Gordon growth model
  # We use year 6 in our vectors as the values used to take the calculations to infinite (perpetual ownership)
  tv_simulation <- bt_cashflows_simulation[inv_t+1] / (r_simulation[inv_t+1] - g_simulation[inv_t+1])
  
  # Add TV to cashflows of the lat year of investment
  bt_cashflows_simulation[inv_t] <- bt_cashflows_simulation[inv_t] + tv_simulation
  
  ########## 
  # AT cashflows (CF) 
  
  # Add debt payments to BT cashflows for years 1 to the last
  at_cashflows_simulation <- bt_cashflows_simulation[1:inv_t] + debt_cashflows
  
  
  ##########  
  # Discounted cashflows (CF)
  discounted_cashflows_simulation <- at_cashflows_simulation[1:inv_t] / compounded_r_simulation[1:inv_t]
  
  # Present value
  pv_simulation <- sum(discounted_cashflows_simulation)
  
  # Year 0 cashflows
  year_0_simulation <- purchase_price - debt_balance
  
  
  # IRR 
  inv_cashflows_simulation <- c(-year_0_simulation, at_cashflows_simulation[1:inv_t])
  inv_irr_simulation <- irr(inv_cashflows_simulation)
  
  # Store the results of this simulation
  pv_simulation_results[s] <- pv_simulation
  irr_simulation_results[s] <- inv_irr_simulation
  
}

########## 
# Means and standard deviations for simulation results
pv_simulation_results_mu <- mean(pv_simulation_results)
pv_simulation_results_sd <- sd(pv_simulation_results)
pv_simulation_results_min <- min(pv_simulation_results)
pv_simulation_results_max <- max(pv_simulation_results)
irr_simulation_results_mu <- mean(irr_simulation_results)
irr_simulation_results_sd <- sd(irr_simulation_results)
irr_simulation_results_min <- min(irr_simulation_results)
irr_simulation_results_max <- max(irr_simulation_results)

# Print Monte Carlo simulation outputs
message("Monte Carlo simulations: ", simulation_size)
## Monte Carlo simulations: 10000
message("PV mean: ", pv_simulation_results_mu)
## PV mean: 448857953.320456
message("PV SD: ", pv_simulation_results_sd)
## PV SD: 27146172.0788943
message("PV min: ", pv_simulation_results_min)
## PV min: 355997043.797927
message("PV max: ", pv_simulation_results_max)
## PV max: 561088022.814214
message("IRR mean: ", irr_simulation_results_mu)
## IRR mean: 0.155349672169076
message("IRR SD: ", irr_simulation_results_sd)
## IRR SD: 0.0136264450902011
message("IRR min: ", irr_simulation_results_min)
## IRR min: 0.104145272102177
message("IRR max: ", irr_simulation_results_max)
## IRR max: 0.2071234700224
########## 
# Graphics for simulation results

# Show present value results in a histogram
histogram_pv_simulation_results <- hist(pv_simulation_results, probability = TRUE, breaks = 100)

# Show present value results in a bar plot
#barplot_pv_simulation_results <- barplot(pv_simulation_results)
#dev.print(png, file = paste0(work_dir, "pv_simulation_barplot.png"), width = 1024)

# Show IRR results in a histogram
histogram_irr_simulation_results <- hist(irr_simulation_results, probability = TRUE, breaks = 100)

# Show IRR results in a histogram
#barplot_irr_simulation_results <- barplot(irr_simulation_results)
#dev.print(png, file = paste0(work_dir, "irr_simulation_barplot.png"), width = 1024)
library(ggplot2)

# Convert histogram object to a data frame
df_hist_1 <- data.frame(PV = histogram_pv_simulation_results$mids, 
                      Density = histogram_pv_simulation_results$density)

# Create ggplot
ggplot(df_hist_1, aes(x = PV, y = Density)) +
  geom_bar(stat = "identity", fill = "skyblue", alpha = 0.7) +
  labs(title = "Histogram of Present Value Simulation Results",
       x = "Present Value",
       y = "Density") +
  theme_minimal()

########## 
# Show all results for comparison
# Print DCF model outputs
message("---------------")
## ---------------
message("DCF model")
## DCF model
message("PV: ", pv)
## PV: 448434890.807701
message("IRR: ", inv_irr)
## IRR: 0.155494895013448
message("---------------")
## ---------------
message("Monte Carlo simulations")
## Monte Carlo simulations
message("Number of simulations: ", simulation_size)
## Number of simulations: 10000
message("PV mean: ", pv_simulation_results_mu)
## PV mean: 448857953.320456
message("PV SD: ", pv_simulation_results_sd)
## PV SD: 27146172.0788943
message("PV min: ", pv_simulation_results_min)
## PV min: 355997043.797927
message("PV max: ", pv_simulation_results_max)
## PV max: 561088022.814214
message("IRR mean: ", irr_simulation_results_mu)
## IRR mean: 0.155349672169076
message("IRR SD: ", irr_simulation_results_sd)
## IRR SD: 0.0136264450902011
message("IRR min: ", irr_simulation_results_min)
## IRR min: 0.104145272102177
message("IRR max: ", irr_simulation_results_max)
## IRR max: 0.2071234700224
message("---------------")
## ---------------

Formula

calculate_dcf <- function(purchase_price, debt_interest, debt_ltv, noi, v, g, r, inv_t) {
  
  # Investment horizon
  tv_t <- inv_t + 1
  
  # Per year compounded calculations
  compounded_g <- c()
  compounded_r <- c()
  
  for (t in 1:tv_t) {
    # Compounded growth factor (1+g_t)_t
    if (t == 1) {
      compounded_g[t] <- (1 + g[t])
    } else {
      compounded_g[t] <- (1 + g[t])*compounded_g[t-1]
    }
    
    # Compounded discount factor (1+g_t)_t
    if (t == 1) {
      compounded_r[t] <- (1 + r[t])
    } else {
      compounded_r[t] <- (1 + r[t])*compounded_r[t-1]
    }
  }
  
  # BT cashflows (CF) 
  bt_cashflows <- noi * (1-v) * compounded_g
  
  # Terminal value (TV)
  tv <- bt_cashflows[inv_t+1] / (r[inv_t+1] - g[inv_t+1])
  
  # Add TV to cashflows of the last year of investment
  bt_cashflows[inv_t] <- bt_cashflows[inv_t] + tv
  
  # AT cashflows (CF) 
  # Debt balance
  debt_balance <- purchase_price * debt_ltv
  
  # Debt interest yearly payments
  debt_cashflows <- rep(-(debt_balance * debt_interest), each=inv_t)
  
  # Add debt balance to last year for full debt repayment
  debt_cashflows[inv_t] <- debt_cashflows[inv_t] - debt_balance
  
  # Add debt payments to BT cashflows for years 1 to the last
  at_cashflows <- bt_cashflows[1:inv_t] + debt_cashflows
  
  # Discounted cashflows (CF)
  discounted_cashflows <- at_cashflows / compounded_r[1:inv_t]
  
  # Present value
  pv <- sum(discounted_cashflows)
  
  # Year 0 cashflows
  year_0 <- purchase_price - debt_balance
  
  # IRR 
  inv_cashflows <- c(-year_0, at_cashflows[1:inv_t])
  inv_irr <- irr(inv_cashflows)
  
  # Return PV and IRR
  return(list(PV = pv, IRR = inv_irr))
}
result <- calculate_dcf(purchase_price = (922*1000*1000),
                        debt_interest = 0.045,
                        debt_ltv = 0.6,
                        noi = c((70*1000*1000),(70*1000*1000),(70*1000*1000), (70*1000*1000), (70*1000*1000), (70*1000*1000)),
                        v = c(0.023, 0.023, 0.023, 0.023, 0.023, 0.023),
                        g = c(0.02, 0.02, 0.02, 0.02, 0.02, 0.02),
                        r = c(0.05, 0.05, 0.05, 0.08, 0.08, 0.08),
                        inv_t = 5)

result
## $PV
## [1] 743622579
## 
## $IRR
## [1] 0.2476977

Original Research Project

We are going to start by adding the new values in the original research project.

#
inv_t <- 5

# Time used for terminal value (perpetuity)
tv_t <- inv_t + 1

##########
# Inputs per year + terminal value

# Purchase price
stamp_duty_legal = (1+0.055)

purchase_price <- (950*1000*1000)* stamp_duty_legal


tax_rate = 0.25

# Yearly debt interest - Fixed rate
debt_interest <- 0.045 # 4.5%

# Debt LTV ratio
debt_ltv <- 0.7 # 60%

# Total rent income (NOI)
noi <- c(70293600, 70293600, 70293600, 70293600, 70293600, 70293600)*1.3
# Vacancy rate (v)
v <- c(0.0237, 0.0237, 0.0237, 0.0237, 0.0237, 0.0237)
# Growth rate (g)
g <- c(0.005, 0.005, 0.005, 0.005, 0.005, 0.005)
# Discount rate (r)
r <- c(0.0646,0.0643, 0.0643,   0.0641, 0.0643, 0.0641)


##########################################
# 2. DCF model
##########################################


##########
# Per year compounded calculations
compounded_g <- c()
compounded_r <- c()

for (t in 1:tv_t) {
  # Compounded growth factor (1+g_t)_t
  if (t == 1) {
    compounded_g[t] <- (1 + g[t])
  } else {
    compounded_g[t] <- (1 + g[t])*compounded_g[t-1]
  }
  
  # Compounded discount factor (1+g_t)_t
  if (t == 1) {
    compounded_r[t] <- (1 + r[t])
  } else {
    compounded_r[t] <- (1 + r[t])*compounded_r[t-1]
  }
}

########## 
# BT cashflows (CF) 
# IMPORTANT - This uses operations over vectors 
 bt_cashflows <- noi * (1-v) * compounded_g

# Terminal value (TV)
# Simplified Gordon growth model
# We use year 6 in our vectors as the values used to take the calculations to infinite (perpetual ownership)
tv <- bt_cashflows[inv_t+1] / (r[inv_t+1] - g[inv_t+1])

# Add TV to cashflows of the lat year of investment
bt_cashflows[inv_t] <- bt_cashflows[inv_t] + tv

Debt calculations

# Loan parameters
loan_amount <- purchase_price * debt_ltv  # $540,000,000
  debt_balance <- loan_amount

interest_rate <- 0.07  # 4.5%
fixed_capital_payment <- 15000000  # $15,000,000
loan_term <- inv_t  # 5 years

# Initialize variables
beginning_balance <- loan_amount
repayment_schedule <- data.frame(Year = 1:loan_term,
                                 Beginning_Balance = numeric(loan_term),
                                 Interest = numeric(loan_term),
                                 Fixed_Capital_Payment = rep(fixed_capital_payment, loan_term),
                                 Total_Payment = numeric(loan_term),
                                 Ending_Balance = numeric(loan_term))

# Calculate repayment schedule
for (i in 1:(loan_term - 1)) {
  repayment_schedule$Beginning_Balance[i] <- beginning_balance
  interest <- beginning_balance * interest_rate
  repayment_schedule$Interest[i] <- interest
  total_payment <- fixed_capital_payment + interest
  repayment_schedule$Total_Payment[i] <- total_payment
  ending_balance <- beginning_balance - fixed_capital_payment
  beginning_balance <- ending_balance
  repayment_schedule$Ending_Balance[i] <- ending_balance
}

# Adjust the 5th row
repayment_schedule$Beginning_Balance[loan_term] <- beginning_balance
interest <- beginning_balance * interest_rate
repayment_schedule$Interest[loan_term] <- interest
total_payment <- beginning_balance + interest  # Total payment includes remaining balance and interest
repayment_schedule$Total_Payment[loan_term] <- total_payment
repayment_schedule$Fixed_Capital_Payment[loan_term] <- beginning_balance
repayment_schedule$Ending_Balance[loan_term] <- 0
########## 
# AT cashflows (CF) 

# Debt interest yearly payments
debt_cashflows <- -repayment_schedule$Total_Payment[1:5]


# Add debt payments to BT cashflows for years 1 to the last
at_cashflows <- (bt_cashflows[1:inv_t] + debt_cashflows) 
##########  
# Discounted cashflows (CF)
discounted_cashflows <- at_cashflows / compounded_r[1:inv_t]

# Present value
pv <- sum(discounted_cashflows)

# Year 0 cashflows
year_0 <- purchase_price - debt_balance


# IRR 
inv_cashflows <- c(-year_0, at_cashflows[1:inv_t])
inv_irr <- irr(inv_cashflows)


# Print DCF model outputs
message("PV: ", pv)
## PV: 798202883.644842
message("IRR: ", inv_irr)
## IRR: 0.316223966317731
# Set parameters
HR_mu <- 0.2964560
HR_sd <- 0.2095128
HR_min <- 0
HR_max <- 0.8
simulation_size <- 1000  # You can adjust this as needed

# Generate random values for HR_epsilon_1 (conservative impact)
HR_epsilon_1 <- pmax(pmin(rnorm(simulation_size, HR_mu, HR_sd), HR_max), HR_min)

# Coefficients from the model
flower_coef <- 0.2964560
flower_se <- 0.2095128

# Simulate HR_epsilon_2 using the normal distribution with mean = coefficient and standard deviation = standard error
HR_epsilon_2 <- rnorm(simulation_size, mean = flower_coef, sd = flower_se)

# Create a data frame for ggplot
df <- data.frame(HR_epsilon_1 = HR_epsilon_1, HR_epsilon_2 = HR_epsilon_2)

# Create overlaid histograms with ggplot
ggplot(df, aes(x = HR_epsilon_1, fill = "Conservative Impact")) +
  geom_histogram(alpha = 0.6, position = "identity", bins = 30) +
  geom_histogram(aes(x = HR_epsilon_2, fill = "Standard Deviation as is"), alpha = 0.3, position = "identity", bins = 30) +
  labs(title = "Histogram of HR_epsilon",
       x = "HR_epsilon",
       y = "Frequency",
       fill = "Variable Impact") +
  scale_fill_manual(values = c("Conservative Impact" = "skyblue", "Standard Deviation as is" = "red"))

##########################################
# 2. Monte Carlo simulations
##########################################
# Number of simulations
simulation_size <- 10000

# For generation of consistent pseudo-random numbers on computers, the seed value does not matter in this file
set.seed(200)

########## 
# We use random shocks based on assumed distributions for for each input
# Each shock adds a random epsilon value to the assumed input

# Vacancy rate - Uniform distribution around 0 ~ U(min, max)
v_min <- -0.001
v_max <- 0.001
v_epsilon <- runif(simulation_size, v_min, v_max)

# Growth rate - Normal distribution centered at 0 ~ N(0, sd)
g_mu <- 0
g_sd <- 0.001
g_epsilon <- rnorm(simulation_size, g_mu, g_sd)

# Discount rate - Normal distribution centered at 0 ~ N(0, sd)
r_mu <- 0
r_sd <- 0.001
r_epsilon <- rnorm(simulation_size, r_mu, r_sd)

#Hedonic regression impact rate

HR_mu <- 0.2964560
HR_sd <- 0.2095128
HR_min <- 0
HR_max <- 0.8
HR_epsilon_1 <- pmax(pmin(rnorm(simulation_size, HR_mu, HR_sd), HR_max), HR_min)



##########
# For the simulations we run the DCF model as before, but adding a different 
# combination of random shocks
##########
pv_simulation_results <- c()
irr_simulation_results <- c()

for (s in 1:simulation_size) {
  
  # Simulation values using the random shocks
  v_simulation <- v + v_epsilon[s]
  g_simulation <- g + g_epsilon[s]
  r_simulation <- r + r_epsilon[s]
  
  # Per year compounded calculations
  compounded_g_simulation <- c()
  compounded_r_simulation <- c()
  
  for (t in 1:tv_t) {
    # Compounded growth factor (1+g_t)_t
    if (t == 1) {
      compounded_g_simulation[t] <- (1 + g_simulation[t])
    } else {
      compounded_g_simulation[t] <- (1 + g_simulation[t])*compounded_g_simulation[t-1]
    }
    
    # Compounded discount factor (1+g_t)_t
    if (t == 1) {
      compounded_r_simulation[t] <- (1 + r_simulation[t])
    } else {
      compounded_r_simulation[t] <- (1 + r_simulation[t])*compounded_r_simulation[t-1]
    }
  }

  ########## 
  # Cashflows (CF) 
  # IMPORTANT - This uses operations over vectors 
  bt_cashflows_simulation <- noi * (1-v_simulation) * compounded_g_simulation * (1 + HR_epsilon_1[s])
  
  # Terminal value (TV)
  # Simplified Gordon growth model
  # We use year 6 in our vectors as the values used to take the calculations to infinite (perpetual ownership)
  tv_simulation <- bt_cashflows_simulation[inv_t+1] / (r_simulation[inv_t+1] - g_simulation[inv_t+1])
  
  # Add TV to cashflows of the lat year of investment
  bt_cashflows_simulation[inv_t] <- bt_cashflows_simulation[inv_t] + tv_simulation
  
  ########## 
  # AT cashflows (CF) 
  
  # Add debt payments to BT cashflows for years 1 to the last
  at_cashflows_simulation <- bt_cashflows_simulation[1:inv_t] + debt_cashflows 
  
  
  ##########  
  # Discounted cashflows (CF)
  discounted_cashflows_simulation <- at_cashflows_simulation[1:inv_t] / compounded_r_simulation[1:inv_t]
  
  # Present value
  pv_simulation <- sum(discounted_cashflows_simulation)
  
  # Year 0 cashflows
  year_0_simulation <- purchase_price - debt_balance
  
  
  # IRR 
  inv_cashflows_simulation <- c(-year_0_simulation, at_cashflows_simulation[1:inv_t])
  inv_irr_simulation <- irr(inv_cashflows_simulation)
  
  # Store the results of this simulation
  pv_simulation_results[s] <- pv_simulation
  irr_simulation_results[s] <- inv_irr_simulation
  
}

########## 
# Means and standard deviations for simulation results
pv_simulation_results_mu <- mean(pv_simulation_results)
pv_simulation_results_sd <- sd(pv_simulation_results)
pv_simulation_results_min <- min(pv_simulation_results)
pv_simulation_results_max <- max(pv_simulation_results)
irr_simulation_results_mu <- mean(irr_simulation_results)
irr_simulation_results_sd <- sd(irr_simulation_results)
irr_simulation_results_min <- min(irr_simulation_results)
irr_simulation_results_max <- max(irr_simulation_results)

# Print Monte Carlo simulation outputs
message("Monte Carlo simulations: ", simulation_size)
## Monte Carlo simulations: 10000
message("PV mean: ", pv_simulation_results_mu)
## PV mean: 1253716704.69632
message("PV SD: ", pv_simulation_results_sd)
## PV SD: 297381997.316952
message("PV min: ", pv_simulation_results_min)
## PV min: 694434503.872494
message("PV max: ", pv_simulation_results_max)
## PV max: 2192854372.46808
message("IRR mean: ", irr_simulation_results_mu)
## IRR mean: 0.466424373354835
message("IRR SD: ", irr_simulation_results_sd)
## IRR SD: 0.0903364336441466
message("IRR min: ", irr_simulation_results_min)
## IRR min: 0.281117206197885
message("IRR max: ", irr_simulation_results_max)
## IRR max: 0.703576791455973
########## 
# Graphics for simulation results

# Show present value results in a histogram
histogram_pv_simulation_results_OR <- hist(pv_simulation_results, probability = TRUE, breaks = 100)

# Show present value results in a bar plot
#barplot_pv_simulation_results <- barplot(pv_simulation_results)
#dev.print(png, file = paste0(work_dir, "pv_simulation_barplot.png"), width = 1024)

# Show IRR results in a histogram
histogram_irr_simulation_results_OR <- hist(irr_simulation_results, probability = TRUE, breaks = 100)

# Show IRR results in a histogram
#barplot_irr_simulation_results <- barplot(irr_simulation_results)
#dev.print(png, file = paste0(work_dir, "irr_simulation_barplot.png"), width = 1024)
library(ggplot2)

# Convert histogram object to a data frame
df_hist <- data.frame(PV = histogram_pv_simulation_results_OR$mids, 
                      Density = histogram_pv_simulation_results_OR$density)

# Create ggplot
ggplot(df_hist, aes(x = PV, y = Density)) +
  geom_bar(stat = "identity", fill = "skyblue", alpha = 0.7) +
  labs(title = "Histogram of Present Value Simulation Results",
       x = "Present Value",
       y = "Density") +
  theme_minimal()

ggplot(df_hist_1, aes(x = PV, y = Density)) +
  geom_bar(stat = "identity", fill = "skyblue", alpha = 0.7) +
  labs(title = "Histogram of Present Value Simulation Results",
       x = "Present Value",
       y = "Density") +
  theme_minimal()

library(ggplot2)

# Convert histogram objects to data frames
df_hist <- data.frame(PV = histogram_pv_simulation_results$mids, 
                      Density = histogram_pv_simulation_results$density)
df_hist_1 <- data.frame(PV = histogram_pv_simulation_results_OR$mids, 
                        Density = histogram_pv_simulation_results_OR$density)

# Create ggplot for the first histogram
plot <- ggplot() +
  geom_bar(data = df_hist, aes(x = PV, y = Density), fill = "skyblue", alpha = 0.7, stat = "identity") +
  geom_bar(data = df_hist_1, aes(x = PV, y = Density), fill = "orange", alpha = 0.7, stat = "identity") +
  labs(title = "Comparison of Present Value Simulation Results",
       x = "Present Value",
       y = "Density") +
  theme_minimal()

# Print the plot
print(plot)

########## 
# Show all results for comparison
# Print DCF model outputs
message("---------------")
## ---------------
message("DCF model")
## DCF model
message("PV: ", pv)
## PV: 798202883.644842
message("IRR: ", inv_irr)
## IRR: 0.316223966317731
message("---------------")
## ---------------
message("Monte Carlo simulations")
## Monte Carlo simulations
message("Number of simulations: ", simulation_size)
## Number of simulations: 10000
message("PV mean: ", pv_simulation_results_mu)
## PV mean: 1253716704.69632
message("PV SD: ", pv_simulation_results_sd)
## PV SD: 297381997.316952
message("PV min: ", pv_simulation_results_min)
## PV min: 694434503.872494
message("PV max: ", pv_simulation_results_max)
## PV max: 2192854372.46808
message("IRR mean: ", irr_simulation_results_mu)
## IRR mean: 0.466424373354835
message("IRR SD: ", irr_simulation_results_sd)
## IRR SD: 0.0903364336441466
message("IRR min: ", irr_simulation_results_min)
## IRR min: 0.281117206197885
message("IRR max: ", irr_simulation_results_max)
## IRR max: 0.703576791455973
message("---------------")
## ---------------
hist(pv_simulation_results, probability = TRUE, breaks = 100)