# 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)
library("DT")
## Warning: package 'DT' was built under R version 4.3.2
##########################################
# 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.045

# Debt LTV ratio
debt_ltv <- 0.6


loan = debt_ltv * 864037403

mezz_loan = 112795841

loan_fee =  -11018448

# Total rent income (NOI) (this includes the revenue from parking as well)
noi <- c(72000000, 72000000, 72000000, 72000000, 72000000, 72000000)
# Vacancy rate + Operating Costs + letting Fees(v)
v <- c(0.0615, 0.0615, 0.0615, 0.0615, 0.0615, 0.0585)
# 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 <- 864037403 * debt_ltv
  debt_balance <- loan_amount

interest_rate <-  debt_interest #  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
calculate_repayments <- function(loan_amount, interest_rate, periods) {
  # Convert interest rate from percentage to decimal
  interest_rate_decimal <- as.numeric(sub("%", "", interest_rate)) / 100
  
  # Calculate yearly interest payment
  yearly_interest <- loan_amount * interest_rate_decimal
  
  # Initialize vector to store repayment values
  repayments <- numeric(periods)
  
  # Calculate yearly interest payments
  for (i in 1:(periods - 1)) {
    repayments[i] <- yearly_interest
  }
  
  # Calculate final repayment (capital + interest)
  final_repayment <- loan_amount * (1 + interest_rate_decimal)
  repayments[periods] <- final_repayment
  
  return(repayments)
}

# Example usage
loan_amount <- 112795841
interest_rate <- "8%"
periods <- 5

repayments_mezz <- calculate_repayments(loan_amount, interest_rate, periods)
repayments_mezz
## [1]   9023667   9023667   9023667   9023667 121819508
########## 
# AT cashflows (CF) 

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


# 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 + mezz_loan + loan_fee)


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


# Print DCF model outputs
message("PV: ", pv)
## PV: 352737129.873966
message("IRR: ", inv_irr)
## IRR: 0.117708417785998
# Create a data frame with the results
results <- data.frame(
  Year = 1:5,
  Debt_Cashflows = debt_cashflows,
  BT_Cashflows = bt_cashflows[1:inv_t],
  AT_Cashflows = at_cashflows,
  Discounted_Cashflows = discounted_cashflows
)

results_transposed <- t(round(results))


# Render the data table
datatable(results_transposed)
##########################################
# 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 - mezz_loan - loan_fee
  
  # 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: 352909054.613529
message("PV SD: ", pv_simulation_results_sd)
## PV SD: 16964634.4270807
message("PV min: ", pv_simulation_results_min)
## PV min: 293895907.880067
message("PV max: ", pv_simulation_results_max)
## PV max: 421734645.316957
message("IRR mean: ", irr_simulation_results_mu)
## IRR mean: 0.117554896328135
message("IRR SD: ", irr_simulation_results_sd)
## IRR SD: 0.0112994029573791
message("IRR min: ", irr_simulation_results_min)
## IRR min: 0.0748272968104169
message("IRR max: ", irr_simulation_results_max)
## IRR max: 0.160322185022093
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: 352737129.873966
message("IRR: ", inv_irr)
## IRR: 0.117708417785998
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: 352909054.613529
message("PV SD: ", pv_simulation_results_sd)
## PV SD: 16964634.4270807
message("PV min: ", pv_simulation_results_min)
## PV min: 293895907.880067
message("PV max: ", pv_simulation_results_max)
## PV max: 421734645.316957
message("IRR mean: ", irr_simulation_results_mu)
## IRR mean: 0.117554896328135
message("IRR SD: ", irr_simulation_results_sd)
## IRR SD: 0.0112994029573791
message("IRR min: ", irr_simulation_results_min)
## IRR min: 0.0748272968104169
message("IRR max: ", irr_simulation_results_max)
## IRR max: 0.160322185022093
message("---------------")
## ---------------
buy_and_hold_results <- list(
  Monte_Carlo_simulations = list(
    Number_of_simulations = simulation_size,
    PV_mean = pv_simulation_results_mu,
    PV_SD = pv_simulation_results_sd,
    PV_min = pv_simulation_results_min,
    PV_max = pv_simulation_results_max,
    IRR_mean = irr_simulation_results_mu,
    IRR_SD = irr_simulation_results_sd,
    IRR_min = irr_simulation_results_min,
    IRR_max = irr_simulation_results_max
  )
)

noi <- c(72385596, 72385596, 72385596, 72385596, 72385596, 72385596)

Original Research Project

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

# Original Resarch 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 <- (1056048956)* stamp_duty_legal


tax_rate = 0

# Yearly debt interest - Fixed rate 
debt_interest <- 0.045

# Debt LTV ratio
debt_ltv <- 0.6


loan = debt_ltv * 990421555
loan_fee = 650*1000 + loan*0.02


mezz_loan = 174331082



# Total rent income (NOI) (this includes the revenue from parking as well)
noi <- c(91290000   ,91290000,  91290000,   91290000    ,91290000)
# Vacancy rate + Operating Costs + letting Fees(v)
v <- c(0.059, 0.059, 0.059, 0.059, 0.059, 0.05)
# Growth rate (g)
g <- c( 0 ,  0.0055 ,  0.0055 ,  0.0055 ,  0.0055 ,  0.0055 )
# Discount rate (r)
r <- c(0.086031817, 0.085287304,    0.084915406,    0.085435029,    0.084826113,    0.087130711)



##########################################
# 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 
## Warning in noi * (1 - v): longitud de objeto mayor no es múltiplo de la
## longitud de uno menor
# 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

  debt_balance <- loan

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

# Initialize variables
beginning_balance <- debt_balance
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
calculate_repayments <- function(loan_amount, interest_rate, periods) {
  # Convert interest rate from percentage to decimal
  interest_rate_decimal <- as.numeric(sub("%", "", interest_rate)) / 100
  
  # Calculate yearly interest payment
  yearly_interest <- loan_amount * interest_rate_decimal
  
  # Initialize vector to store repayment values
  repayments <- numeric(periods)
  
  # Calculate yearly interest payments
  for (i in 1:(periods - 1)) {
    repayments[i] <- yearly_interest
  }
  
  # Calculate final repayment (capital + interest)
  final_repayment <- loan_amount * (1 + interest_rate_decimal)
  repayments[periods] <- final_repayment
  
  return(repayments)
}

# Example usage
mezz_loan
## [1] 174331082
interest_rate_mezz <- "8%"
periods <- 5

repayments_mezz <- calculate_repayments(mezz_loan, interest_rate_mezz, periods)
repayments_mezz
## [1]  13946487  13946487  13946487  13946487 188277569
########## 
# AT cashflows (CF) 

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


# 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 + mezz_loan - loan_fee)


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


# Print DCF model outputs
message("PV: ", pv)
## PV: 391807622.882507
message("IRR: ", inv_irr)
## IRR: 0.155402909462786
# Create a data frame with the results
results <- data.frame(
  Year = 1:5,
  Debt_Cashflows = debt_cashflows,
  BT_Cashflows = bt_cashflows[1:inv_t],
  AT_Cashflows = at_cashflows,
  Discounted_Cashflows = discounted_cashflows
)

results_transposed <- t(round(results))


# Render the data table
datatable(results_transposed)
##########################################
# 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 <- HR_mu - (1.5* HR_sd)
HR_max = HR_mu + (1.5* HR_sd)
HR_epsilon_1 <- pmax(pmin(rnorm(simulation_size, HR_mu, HR_sd), HR_max), HR_min)

noi <- c(72385596, 72385596, 72385596, 72385596, 72385596, 72385596)


##########
# 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 + mezz_loan - loan_fee)
  
  
  # 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: 383839566.982287
  message("PV SD: ", pv_simulation_results_sd)
## PV SD: 157045648.148658
  message("PV min: ", pv_simulation_results_min)
## PV min: 83856219.9868268
  message("PV max: ", pv_simulation_results_max)
## PV max: 728214504.432527
  message("IRR mean: ", irr_simulation_results_mu)
## IRR mean: 0.131910953938003
  message("IRR SD: ", irr_simulation_results_sd)
## IRR SD: 0.122560214906788
  message("IRR min: ", irr_simulation_results_min)
## IRR min: -0.190156117233456
  message("IRR max: ", irr_simulation_results_max)
## IRR max: 0.338830916835164
########## 
# 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)
# 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()

options(scipen = 3)
# 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 = "Buy and Hold Strategy"), alpha = 0.7, stat = "identity") +
  geom_bar(data = df_hist_1, aes(x = PV, y = Density, fill = "Alternative Strategy"), alpha = 0.7, stat = "identity") +
  labs(title = "Comparison of Present Value Simulation Results",
       x = "Present Value",
       y = "Density") +
  scale_fill_manual(values = c("Buy and Hold Strategy" = "skyblue", "Alternative Strategy" = "orange")) +
  theme_minimal()

# Print the plot
print(plot)

options(scipen = 3)

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 = "Buy and Hold Strategy"), alpha = 0.7, stat = "identity") +
  geom_bar(data = df_hist_1, aes(x = PV, y = Density, fill = "Alternative Strategy"), alpha = 0.7, stat = "identity") +
  labs(title = "Comparison of Present Value Simulation Results",
       x = "Present Value",
       y = "Density") +
  scale_fill_manual(values = c("Buy and Hold Strategy" = "skyblue", "Alternative Strategy" = "orange")) +
  theme_minimal() +
  theme(
    text = element_text(size = 16) # Adjust the size of all text elements
  )

# Print the plot
print(plot)

########## 
# Show all results for comparison
# Print DCF model outputs
message("---------------")
## ---------------
message("DCF model")
## DCF model
message("PV: ", pv)
## PV: 391807622.882507
message("IRR: ", inv_irr)
## IRR: 0.155402909462786
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: 383839566.982287
message("PV SD: ", pv_simulation_results_sd)
## PV SD: 157045648.148658
message("PV min: ", pv_simulation_results_min)
## PV min: 83856219.9868268
message("PV max: ", pv_simulation_results_max)
## PV max: 728214504.432527
message("IRR mean: ", irr_simulation_results_mu)
## IRR mean: 0.131910953938003
message("IRR SD: ", irr_simulation_results_sd)
## IRR SD: 0.122560214906788
message("IRR min: ", irr_simulation_results_min)
## IRR min: -0.190156117233456
message("IRR max: ", irr_simulation_results_max)
## IRR max: 0.338830916835164
message("---------------")
## ---------------
alternative_strategy_results <- list(
  DCF_model = list(
    PV = pv,
    IRR = inv_irr
  ),
  Monte_Carlo_simulations = list(
    Number_of_simulations = simulation_size,
    PV_mean = pv_simulation_results_mu,
    PV_SD = pv_simulation_results_sd,
    PV_min = pv_simulation_results_min,
    PV_max = pv_simulation_results_max,
    IRR_mean = irr_simulation_results_mu,
    IRR_SD = irr_simulation_results_sd,
    IRR_min = irr_simulation_results_min,
    IRR_max = irr_simulation_results_max
  )
)
# Create a data frame for Buy and Hold results

buy_and_hold_data <- data.frame(
  Parameter = c("Number of Simulations", "PV Mean", "PV Standard Deviation", "PV Minimum", "PV Maximum", "IRR Mean", "IRR Standard Deviation", "IRR Minimum", "IRR Maximum"),
  Value = c(buy_and_hold_results$Monte_Carlo_simulations$Number_of_simulations, 
             buy_and_hold_results$Monte_Carlo_simulations$PV_mean, 
             buy_and_hold_results$Monte_Carlo_simulations$PV_SD, 
             buy_and_hold_results$Monte_Carlo_simulations$PV_min, 
             buy_and_hold_results$Monte_Carlo_simulations$PV_max, 
             buy_and_hold_results$Monte_Carlo_simulations$IRR_mean, 
             buy_and_hold_results$Monte_Carlo_simulations$IRR_SD, 
             buy_and_hold_results$Monte_Carlo_simulations$IRR_min, 
             buy_and_hold_results$Monte_Carlo_simulations$IRR_max)
)

# Create a data frame for Alternative Strategy results
alternative_strategy_data <- data.frame(
  Parameter = c("Number of Simulations", "PV Mean", "PV Standard Deviation", "PV Minimum", "PV Maximum", "IRR Mean", "IRR Standard Deviation", "IRR Minimum", "IRR Maximum"),
  Value = c(alternative_strategy_results$Monte_Carlo_simulations$Number_of_simulations, 
             alternative_strategy_results$Monte_Carlo_simulations$PV_mean, 
             alternative_strategy_results$Monte_Carlo_simulations$PV_SD, 
             alternative_strategy_results$Monte_Carlo_simulations$PV_min, 
             alternative_strategy_results$Monte_Carlo_simulations$PV_max, 
             alternative_strategy_results$Monte_Carlo_simulations$IRR_mean, 
             alternative_strategy_results$Monte_Carlo_simulations$IRR_SD, 
             alternative_strategy_results$Monte_Carlo_simulations$IRR_min, 
             alternative_strategy_results$Monte_Carlo_simulations$IRR_max)
)

# Bind the two data frames together
combined_data <- cbind(buy_and_hold_data, alternative_strategy_data[, 2])

# Rename the columns
colnames(combined_data) <- c("Parameter", "Buy and Hold", "Alternative Strategy")

# Print the table
datatable(combined_data, 
          rownames = FALSE, 
          options = list(
            dom = 't',
            paging = FALSE,
            ordering = FALSE,
            searching = FALSE
          )
)
# Round the values to two decimal points
buy_and_hold_data$Value <- round(as.numeric(buy_and_hold_data$Value),digits =  4)
alternative_strategy_data$Value <- round(as.numeric(alternative_strategy_data$Value),digits = 4)

# Bind the two data frames together
combined_data <- cbind(buy_and_hold_data, alternative_strategy_data[, 2])

# Rename the columns
colnames(combined_data) <- c("Parameter", "Buy and Hold", "Alternative Strategy")

# Round the values in the combined_data data frame
combined_data[, c("Buy and Hold", "Alternative Strategy")] <- round(combined_data[, c("Buy and Hold", "Alternative Strategy")], 2)

# Print the table
datatable(combined_data, 
          rownames = FALSE, 
          options = list(
            dom = 't',
            paging = FALSE,
            ordering = FALSE,
            searching = FALSE
          )
)
# Set the number of simulations
simulation_size <- 1000

# Define parameters for vacancy rate
v_min <- -0.001
v_max <- 0.001

# Generate data for vacancy rate distribution
v_epsilon <- runif(simulation_size, v_min, v_max)

# Define parameters for growth rate
g_mu <- 0
g_sd <- 0.001

# Generate data for growth rate distribution
g_epsilon <- rnorm(simulation_size, g_mu, g_sd)

# Define parameters for discount rate
r_mu <- 0
r_sd <- 0.001

# Generate data for discount rate distribution
r_epsilon <- rnorm(simulation_size, r_mu, r_sd)

# Plot histograms for the distributions
par(mfrow = c(1, 3))  # Set layout to 1 row and 3 columns

# Vacancy rate histogram
hist(v_epsilon, main = "Vacancy Rate Distribution", xlab = "Vacancy Rate", breaks = 30)

# Growth rate histogram
hist(g_epsilon, main = "Growth Rate Distribution", xlab = "Growth Rate", breaks = 30)

# Discount rate histogram
hist(r_epsilon, main = "Discount Rate Distribution", xlab = "Discount Rate", breaks = 30)

# Set parameters
HR_mu <- 0.255
HR_sd <- 0.2095128
HR_min <- HR_mu - (1.5* HR_sd)
HR_max = HR_mu + (1.5* HR_sd)
HR_epsilon_1 <- pmax(pmin(rnorm(simulation_size, HR_mu, HR_sd), HR_max), HR_min)

# 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.255
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.2, 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"))