R/merton.R

In Skumin/baringaav: Asset valuation functions from the Climate Change suite of models

# Merton model solvers
d1_calc <- function(asset_value, default_point, sigmaA, rf, horizon = 1) {
  return((log(asset_value / default_point) + horizon * (rf + sigmaA^2 / 2))/(sigmaA * sqrt(horizon)))
}

d2_calc <- function(asset_value, default_point, sigmaA, rf, horizon = 1) {
  return(d1_calc(asset_value, default_point, sigmaA, rf, horizon) - sigmaA)
}

equity_value_calc <- function(asset_value, default_point, sigmaA, rf, horizon = 1) {
  d1 <- d1_calc(asset_value, default_point, sigmaA, rf, horizon)
  d2 <- d2_calc(asset_value, default_point, sigmaA, rf, horizon)
  return(asset_value * pnorm(d1) - default_point * exp(-rf) * pnorm(d2))
}

equity_value_barrier_calc <- function(asset_value, default_point, sigmaA, rf, horizon = 1) {
  lambda <- 0.5 + rf / (sigmaA ^ 2)
  vanilla_part <- equity_value_calc(asset_value, default_point, sigmaA, rf, horizon)
  correction_term <- (default_point / asset_value) ^ (2 * lambda - 2) * equity_value_calc((default_point ^ 2) / asset_value, default_point, sigmaA, rf, horizon)
  return(vanilla_part - correction_term)
}

hedge_ratio <- function(asset_value, default_point, sigmaA, rf) {
  hr <- exp(-((rf + 0.5 * sigmaA^2 + log(asset_value / default_point))^2 / (2 * sigmaA^2))) / (sqrt(2 * pi) * sigmaA) -
    (default_point * exp(-rf - 0.5 * (sigmaA - (rf + 0.5 * sigmaA^2 + log(asset_value / default_point)) / sigmaA)^2)) /
    (asset_value * sqrt(2 * pi) * sigmaA) + 0.5 * VGAM::erfc(-((rf + 0.5 * sigmaA^2 + log(asset_value / default_point)) / (sqrt(2) * sigmaA)))
  return(hr)
}

sigmaE_calc <- function(asset_value, equity_value, default_point, sigmaA, rf) {
  return(sigmaA * asset_value / equity_value * hedge_ratio(asset_value, default_point, sigmaA, rf))
}

dd_calc <- function(asset_value, default_point, sigmaA, muA, horizon = 1) {
  return((log(asset_value / default_point) + horizon * (muA - 0.5 * sigmaA^2)) / (sigmaA * sqrt(horizon)))
}

pd_calc <- function(asset_value, default_point, sigmaA, muA, horizon = 1) {
  return(pnorm(-dd_calc(asset_value, default_point, sigmaA, muA, horizon)))
}

pd_barrier_calc <- function(asset_value, default_point, sigmaA, muA, horizon = 1) {
  dd_part <- pnorm(-dd_calc(asset_value, default_point, sigmaA, muA, horizon))
  correction_term <- (default_point / asset_value) ^ (2 * muA / sigmaA^2 - 1) * pnorm(dd_calc(default_point, asset_value, sigmaA, muA, horizon))
  return(dd_part + correction_term)
}

rm_solver <- function(equity_value, default_point, sigmaE, rf) {
  stopifnot(equity_value > 0 & default_point >= 0 & sigmaE > 0 & rf > 0)

  # Scale the inputs first
  scalar <- floor(min(log10(equity_value), log10(default_point)))
  equity_value_sc <- equity_value / (10 ^ scalar)
  default_point_sc <- default_point / (10 ^ scalar)

  # Bound the asset vols
  sigmaA_lb <- sqrt(2 * rf / 85)
  sigmaA_ub <- 0.5

  obj_function <- function(x) {
    y <- numeric(2)

    # Price the equity
    y[1] <- equity_value_calc(x[1], default_point_sc, (sigmaA_ub - sigmaA_lb) / 2 * tanh(x[2]) - ((sigmaA_ub - sigmaA_lb) / 2 - sigmaA_ub), rf) - equity_value_sc

    # Match the vols
    y[2] <- sigmaE_calc(x[1], equity_value_sc, default_point_sc, (sigmaA_ub - sigmaA_lb) / 2 * tanh(x[2]) - ((sigmaA_ub - sigmaA_lb) / 2 - sigmaA_ub), rf) - sigmaE

    return(y)
  }

  res <- nleqslv::nleqslv(c(equity_value_sc + default_point_sc, 0.15), obj_function, control = list(ftol = 1e-16, xtol = 1e-16, maxit = 5000), method = 'Broyden', global = 'dbldog')

  if((res$termcd != 1 & (abs(res$fvec[1]) > 1e-6 | abs(res$fvec[2]) > 1e-6))) {
    while((res$termcd != 1 & (abs(res$fvec[1]) > 1e-6 | abs(res$fvec[2]) > 1e-6)) & sigmaE <= 2.0) {
      sigmaE <- sigmaE + 0.01

      obj_function <- function(x) {
        y <- numeric(2)

        # Price the equity
        y[1] <- equity_value_calc(x[1], default_point_sc, (sigmaA_ub - sigmaA_lb) / 2 * tanh(x[2]) - ((sigmaA_ub - sigmaA_lb) / 2 - sigmaA_ub), rf) - equity_value_sc

        # Match the vols
        y[2] <- sigmaE_calc(x[1], equity_value_sc, default_point_sc, (sigmaA_ub - sigmaA_lb) / 2 * tanh(x[2]) - ((sigmaA_ub - sigmaA_lb) / 2 - sigmaA_ub), rf) - sigmaE

        return(y)
      }

      res <- nleqslv::nleqslv(c(equity_value_sc + default_point_sc, 0.15), obj_function, control = list(ftol = 1e-16, xtol = 1e-16, maxit = 5000), method = 'Broyden', global = 'dbldog')
    }
  }

  return(list(asset_value = res$x[1] * (10 ^ scalar), debt_value = res$x[1] * (10 ^ scalar) - equity_value, sigmaA = (sigmaA_ub - sigmaA_lb) / 2 * tanh(res$x[2]) - ((sigmaA_ub - sigmaA_lb) / 2 - sigmaA_ub), sigmaE = sigmaE, pd = pd_calc(res$x[1], default_point_sc, (sigmaA_ub - sigmaA_lb) / 2 * tanh(res$x[2]) - ((sigmaA_ub - sigmaA_lb) / 2 - sigmaA_ub), rf), conv = res$termcd, equity_value_resid = res$fvec[1], sigmaE_resid = res$fvec[2]))
}

am_solver <- function(equity_value, default_point, pd, rf, sigmaA_bounds = NULL) {
  stopifnot(equity_value > 0 & default_point >= 0 & pd > 0 & rf > 0)

  if(length(sigmaA_bounds) > 0) {
    stopifnot(length(sigmaA_bounds) == 2)
    stopifnot(sigmaA_bounds[1] < sigmaA_bounds[2])
  }

  # Scale the inputs first
  scalar <- floor(min(log10(equity_value), log10(default_point)))
  equity_value_sc <- equity_value / (10 ^ scalar)
  default_point_sc <- default_point / (10 ^ scalar)

  if(length(sigmaA_bounds) > 0) {
    sigmaA_lb <- sigmaA_bounds[1]
    sigmaA_ub <- sigmaA_bounds[2]

    obj_function <- function(x) {
      y <- numeric(2)

      # Price the equity
      y[1] <- equity_value_calc(x[1], default_point_sc, (sigmaA_ub - sigmaA_lb) / 2 * tanh(x[2]) - ((sigmaA_ub - sigmaA_lb) / 2 - sigmaA_ub), rf) - equity_value_sc

      # Match the PD
      y[2] <- pd_calc(x[1], default_point_sc, (sigmaA_ub - sigmaA_lb) / 2 * tanh(x[2]) - ((sigmaA_ub - sigmaA_lb) / 2 - sigmaA_ub), rf) - pd

      return(y)
    }
  } else {
    obj_function <- function(x) {
      y <- numeric(2)

      # Price the equity
      y[1] <- equity_value_calc(x[1], default_point_sc, x[2], rf) - equity_value_sc

      # Match the PD
      y[2] <- pd_calc(x[1], default_point_sc, x[2], rf) - pd

      return(y)
    }
  }

  res <- nleqslv::nleqslv(c(equity_value_sc + default_point_sc, 0.3), obj_function, control = list(ftol = 1e-16, xtol = 1e-16, maxit = 5000), method = 'Broyden', global = 'dbldog')
  first_pass_conv <- 1

  # Try different starting values if there's no convergence
  if(res$termcd != 1 & (abs(res$fvec[1]) > 1e-6) | abs(res$fvec[2]) > 1e-6) {
    first_pass_conv <- 0

    # pd is in the pre-made mapping table
    if(pd %in% baringaav::mapping_table[, unique(PD)]) {
      closest_pd <- pd

    # pd is not in the pre-made mapping table
    } else {
      warning('Initial PD is not in the standard Baringa asset valuation PD master scale.')
      closest_pd_id <- findInterval(pd, c(0, sort(unique(baringaav::mapping_table$PD)), 1))
      if(closest_pd_id == 1) {
        closest_pd <- min(unique(baringaav::mapping_table$PD))
      } else if(closest_pd_id == 20) {
        closest_pd <- max(unique(baringaav::mapping_table$PD))
      } else {
        closest_pd <- c(0, sort(unique(baringaav::mapping_table$PD)), 1)[closest_pd_id]
      }
    }

    # Get the closest equity value multiplier from the mapping table
    market_leverage <- equity_value_sc / default_point_sc
    closest_equity_id <- baringaav::mapping_table[PD == closest_pd, findInterval(market_leverage, equity_multiplier)]

    if(closest_equity_id == 0) {
      closest_equity <- baringaav::mapping_table[PD == closest_pd, min(equity_multiplier)]
    } else {
      closest_equity <- baringaav::mapping_table[PD == closest_pd, equity_multiplier[closest_equity_id]]
    }
    starting_sigmaA <- baringaav::mapping_table[PD == closest_pd & equity_multiplier == closest_equity, sigmaA_mapped]
    if(length(sigmaA_bounds) > 0) {
      starting_sigmaA <- pmax(sigmaA_lb, pmin(sigmaA_ub, starting_sigmaA))
    }

    # Re-run the solver with a different initial guess for the asset volatility
    res <- nleqslv::nleqslv(c(equity_value_sc + default_point_sc, starting_sigmaA), obj_function, control = list(ftol = 1e-16, xtol = 1e-16, maxit = 5000), method = 'Broyden', global = 'dbldog')
  }

  # Return based on whether the tanh transformation was used
  if(length(sigmaA_bounds) > 0) {
    return(list(asset_value = res$x[1] * (10 ^ scalar), debt_value = res$x[1] * (10 ^ scalar) - equity_value, sigmaA = (sigmaA_ub - sigmaA_lb) / 2 * tanh(res$x[2]) - ((sigmaA_ub - sigmaA_lb) / 2 - sigmaA_ub), sigmaE = sigmaE_calc(res$x[1] * (10 ^ scalar), equity_value, default_point, (sigmaA_ub - sigmaA_lb) / 2 * tanh(res$x[2]) - ((sigmaA_ub - sigmaA_lb) / 2 - sigmaA_ub), rf), pd = pd, conv = res$termcd, equity_value_resid = res$fvec[1], pd_resid = res$fvec[2], first_pass_conv = first_pass_conv))
  } else {
    return(list(asset_value = res$x[1] * (10 ^ scalar), debt_value = res$x[1] * (10 ^ scalar) - equity_value, sigmaA = res$x[2], sigmaE = sigmaE_calc(res$x[1] * (10 ^ scalar), equity_value, default_point, res$x[2], rf), pd = pd, conv = res$termcd, equity_value_resid = res$fvec[1], pd_resid = res$fvec[2], first_pass_conv = first_pass_conv))
  }
}

am_barrier_solver <- function(equity_value, default_point, pd, rf, sigmaA_bounds = NULL) {
  stopifnot(equity_value > 0 & default_point >= 0 & pd > 0 & rf > 0)

  if(length(sigmaA_bounds) > 0) {
    stopifnot(length(sigmaA_bounds) == 2)
    stopifnot(sigmaA_bounds[1] < sigmaA_bounds[2])
  }

  # Scale the inputs first
  scalar <- floor(min(log10(equity_value), log10(default_point)))
  equity_value_sc <- equity_value / (10 ^ scalar)
  default_point_sc <- default_point / (10 ^ scalar)

  if(length(sigmaA_bounds) > 0) {
    sigmaA_lb <- sigmaA_bounds[1]
    sigmaA_ub <- sigmaA_bounds[2]

    obj_function <- function(x) {
      y <- numeric(2)

      # Price the equity
      y[1] <- equity_value_barrier_calc(x[1], default_point_sc, (sigmaA_ub - sigmaA_lb) / 2 * tanh(x[2]) - ((sigmaA_ub - sigmaA_lb) / 2 - sigmaA_ub), rf) - equity_value_sc

      # Match the PD
      y[2] <- pd_barrier_calc(x[1], default_point_sc, (sigmaA_ub - sigmaA_lb) / 2 * tanh(x[2]) - ((sigmaA_ub - sigmaA_lb) / 2 - sigmaA_ub), rf) - pd

      return(y)
    }
  } else {
    obj_function <- function(x) {
      y <- numeric(2)

      # Price the equity
      y[1] <- equity_value_barrier_calc(x[1], default_point_sc, x[2], rf) - equity_value_sc

      # Match the PD
      y[2] <- pd_barrier_calc(x[1], default_point_sc, x[2], rf) - pd

      return(y)
    }
  }

  res <- nleqslv::nleqslv(c(equity_value_sc + default_point_sc, 0.3), obj_function, control = list(ftol = 1e-16, xtol = 1e-16, maxit = 5000), method = 'Broyden', global = 'dbldog')
  first_pass_conv <- 1

  # Try different starting values if there's no convergence
  if(res$termcd != 1 & (abs(res$fvec[1]) > 1e-6) | abs(res$fvec[2]) > 1e-6) {
    first_pass_conv <- 0

    # pd is in the pre-made mapping table
    if(pd %in% baringaav::mapping_table[, unique(PD)]) {
      closest_pd <- pd

      # pd is not in the pre-made mapping table
    } else {
      warning('Initial PD is not in the standard Baringa asset valuation PD master scale.')
      closest_pd_id <- findInterval(pd, c(0, sort(unique(baringaav::mapping_table$PD)), 1))
      if(closest_pd_id == 1) {
        closest_pd <- min(unique(baringaav::mapping_table$PD))
      } else if(closest_pd_id == 20) {
        closest_pd <- max(unique(baringaav::mapping_table$PD))
      } else {
        closest_pd <- c(0, sort(unique(baringaav::mapping_table$PD)), 1)[closest_pd_id]
      }
    }

    # Get the closest equity value multiplier from the mapping table
    market_leverage <- equity_value_sc / default_point_sc
    closest_equity_id <- baringaav::mapping_table[PD == closest_pd, findInterval(market_leverage, equity_multiplier)]

    if(closest_equity_id == 0) {
      closest_equity <- baringaav::mapping_table[PD == closest_pd, min(equity_multiplier)]
    } else {
      closest_equity <- baringaav::mapping_table[PD == closest_pd, equity_multiplier[closest_equity_id]]
    }
    starting_sigmaA <- baringaav::mapping_table[PD == closest_pd & equity_multiplier == closest_equity, sigmaA_mapped]
    if(length(sigmaA_bounds) > 0) {
      starting_sigmaA <- pmax(sigmaA_lb, pmin(sigmaA_ub, starting_sigmaA)) * 0.8 # The mapping table is based on the vanilla option, this is an ad-hoc multiplier
    }

    # Re-run the solver with a different initial guess for the asset volatility
    res <- nleqslv::nleqslv(c(equity_value_sc + default_point_sc, starting_sigmaA), obj_function, control = list(ftol = 1e-16, xtol = 1e-16, maxit = 5000), method = 'Broyden', global = 'dbldog')
  }

  # Return based on whether the tanh transformation was used
  if(length(sigmaA_bounds) > 0) {
    return(list(asset_value = res$x[1] * (10 ^ scalar), debt_value = res$x[1] * (10 ^ scalar) - equity_value, sigmaA = (sigmaA_ub - sigmaA_lb) / 2 * tanh(res$x[2]) - ((sigmaA_ub - sigmaA_lb) / 2 - sigmaA_ub), sigmaE = NA, pd = pd, conv = res$termcd, equity_value_resid = res$fvec[1], pd_resid = res$fvec[2], first_pass_conv = first_pass_conv))
  } else {
    return(list(asset_value = res$x[1] * (10 ^ scalar), debt_value = res$x[1] * (10 ^ scalar) - equity_value, sigmaA = res$x[2], sigmaE = NA, pd = pd, conv = res$termcd, equity_value_resid = res$fvec[1], pd_resid = res$fvec[2], first_pass_conv = first_pass_conv))
  }
}