# 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)

purchase_price <- ( 920199834)* stamp_duty_legal


tax_rate = 0

# Yearly debt interest - Fixed rate 
debt_interest <- 0.054

# Debt LTV ratio
debt_ltv <- 0.665


 loan = debt_ltv * purchase_price
 
# Total rent income (NOI)
noi <- c(72000000, 72000000, 72000000, 72000000, 72000000, 72000000)
# Vacancy rate + Operating Costs + letting Fees(v)
v <- c(0.06125, 0.06125, 0.06125, 0.06125, 0.06125, 0.06125)
# Growth rate (g)
g <- c( 0 ,  0.0055 ,  0.0055 ,  0.0055 ,  0.0055 ,  0.0055 )
# Discount rate (r)
r <- c(0.07853,0.07779, 0.07742,    0.07794,    0.07733,    0.07963)
                      

##########################################
# 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
  debt_balance <- loan_amount

interest_rate <-  debt_interest # Mezzanine + loan
fixed_capital_payment <- 15000000  # $15,000,000
loan_term <- inv_t

# 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: 362408711.136467
message("IRR: ", inv_irr)
## IRR: 0.117418749937534
# Create a data frame with the results
results <- data.frame(
  Year = 1:5,
  Debt_Cashflows = -repayment_schedule$Total_Payment[1:5],
  BT_Cashflows = bt_cashflows[1:inv_t],
  AT_Cashflows = at_cashflows,
  Discounted_Cashflows = discounted_cashflows
)

results_transposed <- t(results)


# Print the results table
print(results_transposed)
##                           [,1]      [,2]      [,3]      [,4]       [,5]
## Year                         1         2         3         4          5
## Debt_Cashflows       -48044376 -47234376 -46424376 -45614376 -581737266
## BT_Cashflows          67590000  67961745  68335535  68711380 1006217092
## AT_Cashflows          19545624  20727369  21911159  23097004  424479826
## Discounted_Cashflows  18122467  17831086  17494999  17108410  291851749
##########################################
# 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: 362580562.956201
message("PV SD: ", pv_simulation_results_sd)
## PV SD: 16950392.3865354
message("PV min: ", pv_simulation_results_min)
## PV min: 303624721.226985
message("PV max: ", pv_simulation_results_max)
## PV max: 431317480.861901
message("IRR mean: ", irr_simulation_results_mu)
## IRR mean: 0.117279178067217
message("IRR SD: ", irr_simulation_results_sd)
## IRR SD: 0.0109038065186467
message("IRR min: ", irr_simulation_results_min)
## IRR min: 0.0761713502118506
message("IRR max: ", irr_simulation_results_max)
## IRR max: 0.158657446924185
options(scipen = 2)
# Create a histogram of present value simulation results
histogram_pv_simulation_results <- hist(pv_simulation_results, 
                                        probability = TRUE, 
                                        breaks = 100,
                                        col = "skyblue",     # Set color of bars
                                        border = "white",    # Set color of bar borders
                                        main = "Distribution of Present Value Simulation Results",  # Add a title
                                        xlab = "Present Value",   # Label for x-axis
                                        ylab = "Density",    # Label for y-axis
                                        xlim = c(min(pv_simulation_results), max(pv_simulation_results)),  # Set 
)

# Add a smooth density line to the histogram
lines(density(pv_simulation_results), col = "red", lwd = 2)

# Add a legend
legend("topright", 
       legend = c("Histogram", "Density"),
       fill = c("skyblue", "red"),
       border = NA,
       bty = "n"
)

# Create a histogram of IRR simulation results
histogram_irr_simulation_results <- hist(irr_simulation_results, 
                                          probability = TRUE, 
                                          breaks = 100,
                                          col = "skyblue",     # Set color of bars
                                          border = "white",    # Set color of bar borders
                                          main = "Distribution of IRR Simulation Results",  # Add a title
                                          xlab = "IRR",   # Label for x-axis
                                          ylab = "Density",    # Label for y-axis
                                          xlim = c(min(irr_simulation_results), max(irr_simulation_results))
)

# Add a smooth density line to the histogram
lines(density(irr_simulation_results), col = "red", lwd = 2)

# Add a legend
legend("topright", 
       legend = c("Histogram", "Density"),
       fill = c("skyblue", "red"),
       border = NA,
       bty = "n"
)

# 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: 362408711.136467
message("IRR: ", inv_irr)
## IRR: 0.117418749937534
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: 362580562.956201
message("PV SD: ", pv_simulation_results_sd)
## PV SD: 16950392.3865354
message("PV min: ", pv_simulation_results_min)
## PV min: 303624721.226985
message("PV max: ", pv_simulation_results_max)
## PV max: 431317480.861901
message("IRR mean: ", irr_simulation_results_mu)
## IRR mean: 0.117279178067217
message("IRR SD: ", irr_simulation_results_sd)
## IRR SD: 0.0109038065186467
message("IRR min: ", irr_simulation_results_min)
## IRR min: 0.0761713502118506
message("IRR max: ", irr_simulation_results_max)
## IRR max: 0.158657446924185
message("---------------")
## ---------------

Original Research Project

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

# 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)

purchase_price <- ( 920199834)* stamp_duty_legal


tax_rate = 0

# Yearly debt interest - Fixed rate 
debt_interest <- 0.054

# Debt LTV ratio
debt_ltv <- 0.665


 loan = debt_ltv * purchase_price
 
# Total rent income (NOI)
noi <- c(72000000, 72000000, 72000000, 72000000, 72000000, 72000000)
# Vacancy rate + Operating Costs + letting Fees(v)
v <- c(0.06125, 0.06125, 0.06125, 0.06125, 0.06125, 0.06125)
# Growth rate (g)
g <- c( 0 ,  0.0055 ,  0.0055 ,  0.0055 ,  0.0055 ,  0.0055 )
# Discount rate (r)
r <- c(0.07853,0.07779, 0.07742,    0.07794,    0.07733,    0.07963)
                      



##########################################
# 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 <- debt_interest  
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: 362408711.136467
message("IRR: ", inv_irr)
## IRR: 0.117418749937534
# Set parameters
HR_mu <- 0.2964560
HR_sd <- 0.2095128
HR_min <- 0
HR_max <- 0.8
simulation_size <- 10000  # 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.255
HR_sd <- 0.2095128
HR_min <- 0
HR_max = 2
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: 604399477.702973
message("PV SD: ", pv_simulation_results_sd)
## PV SD: 175308981.87055
message("PV min: ", pv_simulation_results_min)
## PV min: 311532522.473686
message("PV max: ", pv_simulation_results_max)
## PV max: 1344844420.33327
message("IRR mean: ", irr_simulation_results_mu)
## IRR mean: 0.251371400651499
message("IRR SD: ", irr_simulation_results_sd)
## IRR SD: 0.0880766703462059
message("IRR min: ", irr_simulation_results_min)
## IRR min: 0.0815913875066876
message("IRR max: ", irr_simulation_results_max)
## IRR max: 0.553725808409047
########## 
# Graphics for simulation results

# Show present value results in a histogram
histogram_pv_simulation_results_OR <- hist(pv_simulation_results, probability = TRUE, breaks = 100)
## Warning in breaks[-1L] + breaks[-nB]: NAs producidos por enteros excedidos

# 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)
# 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()
## Warning: Removed 28 rows containing missing values (`position_stack()`).

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)
## Warning: Removed 28 rows containing missing values (`position_stack()`).

########## 
# Show all results for comparison
# Print DCF model outputs
message("---------------")
## ---------------
message("DCF model")
## DCF model
message("PV: ", pv)
## PV: 362408711.136467
message("IRR: ", inv_irr)
## IRR: 0.117418749937534
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: 604399477.702973
message("PV SD: ", pv_simulation_results_sd)
## PV SD: 175308981.87055
message("PV min: ", pv_simulation_results_min)
## PV min: 311532522.473686
message("PV max: ", pv_simulation_results_max)
## PV max: 1344844420.33327
message("IRR mean: ", irr_simulation_results_mu)
## IRR mean: 0.251371400651499
message("IRR SD: ", irr_simulation_results_sd)
## IRR SD: 0.0880766703462059
message("IRR min: ", irr_simulation_results_min)
## IRR min: 0.0815913875066876
message("IRR max: ", irr_simulation_results_max)
## IRR max: 0.553725808409047
message("---------------")
## ---------------
hist(pv_simulation_results, probability = TRUE, breaks = 100)
## Warning in breaks[-1L] + breaks[-nB]: NAs producidos por enteros excedidos