library(tidyverse)
── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ── ✔ ggplot2 3.3.5 ✔ purrr 0.3.4 ✔ tibble 3.1.3 ✔ dplyr 1.0.7 ✔ tidyr 1.1.3 ✔ stringr 1.4.0 ✔ readr 2.0.0 ✔ forcats 0.5.1 ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ── ✖ dplyr::filter() masks stats::filter() ✖ dplyr::lag() masks stats::lag()
The first assumption is that a bank decides to give 1,000 loans for $180,000 for a year.
It is considered that each one who defaults causes a total loss of $200,000 for the bank, with all operational costs included.
n <- 1000
loss_per_foreclosure <- -200000
p <- 0.02
B <- 10000
losses <- replicate(B, {
defaults <- sample( c(0,1), n, prob=c(1-p, p), replace = TRUE)
sum(defaults * loss_per_foreclosure)
})
data.frame(losses_in_millions = losses/10^6) %>%
ggplot(aes(losses_in_millions)) +
geom_histogram(binwidth = 0.6, col = "black")
A Monte Carlo simulation is executed, for 10,000 times, and the distribution of the losses is illustrated in the histogram.
The Central Limit Theorem indicates that the loss distribution is approximately normal and can be defined by its expected value and standard deviation.
With these values, it is possible to calculate the exact value x, which is the interest rate that, on average, the bank breaks even:
lp + x(1 - p) = 0
x = - loss_per_foreclosure*p/(1-p)
x
x/180000 # interest rate of the $180,000 loan
With an interest rate of about 2%, on average the bank will break even, but it has 50% chance to loose money.
The bank needs to calculate an interest rate which has smaller chance for this to happen.
So a new value of interest rate must be calculated, which has probability of 1% of loosing money.
Applying the formulas of expected value, standard error and standard units, it is possible to calculate the new interest rate.
l <- loss_per_foreclosure
z <- qnorm(0.01)
x <- -l*( n*p - z*sqrt(n*p*(1-p)))/ ( n*(1-p) + z*sqrt(n*p*(1-p)))
x
x/180000 # interest rate of the $180,000 loan
With the code above, it shows that with an interest rate of about 3.5%, there is 1% of probability of the bank loosing money.
loss_per_foreclosure*p + x*(1-p) # expected value of the profit per loan
n*(loss_per_foreclosure*p + x*(1-p)) # expected value of the profit over n loans
With the interest rate calculated above, the expected profit per loan is about $2,124, resulting in a total expected profit of about $2 million a year.
A Monte Carlo simulation illustrates the approximation.
B <- 100000
profit <- replicate(B, {
draws <- sample( c(x, loss_per_foreclosure), n,
prob=c(1-p, p), replace = TRUE)
sum(draws)
})
mean(profit) # expected value of the profit over n loans
mean(profit<0) # probability of losing money
With 100,000 times, the expected profit is about $2 million and the probability of losing money is about 1%, which corresponds to the theoretical approximation.
To maximize the profits, one can rely on the law of large numbers and increase n the number of loans, so even with higher probability of default, with higher interest rates, the bank still can make profit.
The new values for this case are as follow:
p <- .04 # doubles the default probability
loss_per_foreclosure <- -200000
r <- 0.05 # increases the interest rate
x <- r*180000
loss_per_foreclosure*p + x*(1-p) # expected value of the profit per loan
Following the law of large numbers, it is possible to calculate the number of loans that represents a probability of 1% of the bank losing money, on this new conditions.
z <- qnorm(0.01)
l <- loss_per_foreclosure
n <- ceiling((z^2*(x-l)^2*p*(1-p))/(l*p + x*(1-p))^2)
n # number of loans required
n*(loss_per_foreclosure*p + x * (1-p)) # expected profit over n loans
The equations show that it is required a number of 22,163 loans so that the probability of loosing money is about 1%.
In this situation, the expected profit of the bank is about $14 million.
Again a Monte Carlo simulation is run, to illustrate the approximation.
B <- 10000
p <- 0.04
x <- 0.05 * 180000
profit <- replicate(B, {
draws <- sample( c(x, loss_per_foreclosure), n,
prob=c(1-p, p), replace = TRUE)
sum(draws)
})
mean(profit) # expected profit over n loans
mean(profit < 0) # probability of losing money
The Monte Carlo simulation agrees to the theoretical approximation.
The real problem, as ocurred in the financial meltdown in 2007, is that the law of large numbers relies on independent values, in this case, loans, to minimize the standard error. In the real world, many interconnected events may lead to completely larger errors.
The next code reruns the Monte Carlo simulation, with a slightly probability of changing the global default probability between 3% and 5%.
p <- 0.04
x <- 0.05*180000
profit <- replicate(B, {
new_p <- 0.04 + sample(seq(-0.01, 0.01, length = 100), 1)
draws <- sample( c(x, loss_per_foreclosure), n,
prob=c(1-new_p, new_p), replace = TRUE)
sum(draws)
})
mean(profit) # expected profit
mean(profit < 0) # probability of losing money
mean(profit < -10000000) # probability of losing over $10 million
The case presented here is not independent, with an external factor changing the global default probability. The expected profit is still about $14 million. However, the probability of losing money is about 35% and the probability of losing over $10 million is about 24%.
The distribution of the random variable obtained shows that it no longer follows a normal distribution.
data.frame(profit = profit/10^6) %>%
ggplot(aes(profit)) +
geom_histogram(binwidth = 6, col = "black")
