# Required packages. library(tidyquant) library(tidyr) library(Sim.DiffProc) library(dplyr) library(knitr)

install.packages('bazar') # para la función almost.equal install.packages('scatterplot3d')

# Merton assumes V follows a geometric brownian motion process. V0.sim <- function(i) { # the argument i is not used here. GBM(N = 365, T = 1, t0 = 0, x0 = 12.39578, theta = 0.05, sigma = 0.2123041) } set.seed(3) # Reproducibility paths <- sapply(1:100, V0.sim) # Create 10 paths for V. # Plot results. matplot(paths, type = "l", col = "black", lwd = 1, lty = 1, ylab = expression(paste(V[t])), xlab = "Time in days (0 to T)") abline(h = 10, col = "red", lwd = 2) points(1, 12.39578, pch = 19, cex = 1.5, col = "blue")

# V_t < 10 almost_default <- which(colSums(paths < 10) > 0) # V_T < 10 default <- which((paths[366,] < 10)) # There should be an easier way to do this: matplot(paths[,almost_default[which(almost_default %in% default == FALSE)]], type = "l", col = c(1:17), lwd = 1, lty = 1, ylab = expression(paste(V[t])), xlab = "Time in days (0 to T)") abline(h = 10, col = "red", lwd = 2) points(1, 12.7, pch = 19, cex = 1.5, col = "blue")

E0 <- 3 # Value of the equity as today, I recommend market capitalization. se <- 0.8 # Stock returns volatility. rf <- 0.05 # Risk-free rate. TT <- 1 # Maturity. D <- 2 # Value of the debt. The Bloomberg function DDIS is useful. # We need equations 24.3 and 24.4 to estimate V0 and sv. eq24.3 <- function(V0, sv) { ((V0 * pnorm((log(V0 / D) + (rf + sv^2 / 2) * TT) / (sv * sqrt(TT))) - D * exp(-rf * TT) * pnorm(((log(V0 / D) + (rf + sv^2 / 2) * TT) / (sv * sqrt(TT))) - sv * sqrt(TT)) - E0)) } eq24.4 <- function(V0, sv) { ((pnorm((log(V0 / D) + (rf + sv^2 / 2) * TT) / (sv * sqrt(TT))) * sv * V0 - se * E0)) } # Footnote 10 indicates that we should minimize F(x,y)ˆ2+G(x,y)ˆ2 min.footnote10 <- function(x) (eq24.3(x[1], x[2]))^2 + (eq24.4(x[1], x[2]))^2 # The minimization leads to the values of V0 and sv. V0_sv <- optim(c(1, 1), min.footnote10) # Define the values as parameters. V0 <- V0_sv$par[1] sv <- V0_sv$par[2] # Calculate the probability of default as a function. PD <- function(V0, sv) { pnorm(-(((log(V0 / D) + (rf + sv^2 / 2) * TT) / (sv * sqrt(TT))) - sv * sqrt(TT))) } # Calculate the probability of default given the previous parameters. pd <- PD(V0, sv) # Gather the results. ans <- data.frame(pd, V0, sv) ans d1 <- (log(V0 / D) + (rf + sv^2 / 2) * TT) / (sv * sqrt(TT)) d2 <- d1 - sv * sqrt(TT) Nd1 <- pnorm(d1) Nd2 <- pnorm(d2) ans <- data.frame(d1, Nd1, d2, Nd2) ans xseq <- seq(-4, 4, 0.01) densities <- dnorm(xseq, 0, 1) # Standard normal distribution. # Graphical evaluation of N(-d_2). plot(xseq, densities, xlab = "d values (d_2 is dotted line)", ylab = "Density", type = "l", lwd = 2, cex = 2, main = "The future assets will be high enough to pay the debt?" ) abline(h = 0) polygon(c(xseq[xseq >= d2], d2), c(densities[xseq >= d2], 0), col = "red") polygon(c(xseq[xseq < d2], d2), c(densities[xseq < d2], 0), col = "blue") legend("topleft", legend = c("N(d_2)=0.9494245 is represented in blue and N(-d_2)=0.05057548 in red. N(d_2)+N(-d_2)=1"), bg = "white", bty = "n", cex = 0.9) abline(v = d2, lty = 2, lwd = 2)

xseq <- seq(-4, 4, 0.01) densities <- dnorm(xseq, 0, 1) # Standard normal distribution. # Graphical evaluation of N(-d_2). plot(xseq, densities, xlab = "d values (d_2 is dotted line)", ylab = "Density", type = "l", lwd = 2, cex = 2, main = "The future assets will be high enough to pay the debt?" ) abline(h = 0) polygon(c(xseq[xseq >= d2], d2), c(densities[xseq >= d2], 0), col = "red") polygon(c(xseq[xseq < d2], d2), c(densities[xseq < d2], 0), col = "blue") legend("topleft", legend = c("N(d_2)=0.8730604 is represented in blue and N(-d_2)=0.1269396 in red. N(d_2)+N(-d_2)=1"), bg = "white", bty = "n", cex = 0.9) abline(v = d2, lty = 2, lwd = 2)