R/BSM.R

In czxa/FMFE: Options and Futures in Financial Mathematics and Financial Engineering

#' Black-Scholes-Merton Option pricing model and calculation of Greek value
#' @description
#' The BSM function can calculate the theoretical price and Greek value of the option according to the relevant parameters of the input option, including Delta, Gamma, Theta, Vega and Rho.
#' @param S Current stock price, default value is 41
#' @param K Execution price of option, default value is 40
#' @param sigma Annual volatility of the underlying stock price with a default value of 0.3
#' @param r Risk-free interest rate, default value 0.08
#' @param T Expiration time, default value is 1
#' @param t Start time, default value is 0, that is, default duration of options is 1
#' @param dividend Compound dividend rate, default value is 0, that is, no dividend.
#' @param type Specifies the type of option, defaulting to call
#' @examples
#' BSM(S = 41, K = 40, sigma = 0.3, r = 0.08, T = 1, t = 0, dividend = 0)
#' BSM(S = 41, K = 40, sigma = 0.3, r = 0.08, T = 1, t = 0, dividend = 0,
#'     type = "put")
#' @references
#' John C.Hull. Options, Futures, and other Derivatives 9ed
#' @export
BSM <- function (S = 41, K = 40, sigma = 0.3, r = 0.08, T = 1,t = 0, dividend = 0, type = "call"){
  d1 <- (log(S/K) + (r - dividend + sigma^2/2) * (T-t))/(sigma * sqrt(T-t))
  d2 <- (log(S/K) + (r - dividend - sigma^2/2) * (T-t))/(sigma * sqrt(T-t))
  pd1t <- - (r + sigma^2/2 - dividend)/(2*sigma*sqrt(T-t)) + log(S/K)/(2*sigma*(T-t)^1.5)
  pd2t <- - (r - sigma^2/2 - dividend)/(2*sigma*sqrt(T-t)) + log(S/K)/(2*sigma*(T-t)^1.5)
  pd1sigma <- -(log(S/K) + (r-dividend)*(T-t))/(sigma^2*sqrt(T-t)) + 0.5*sqrt(T-t)
  pd2sigma <- -(log(S/K) + (r-dividend)*(T-t))/(sigma^2*sqrt(T-t)) - 0.5*sqrt(T-t)
  pd1r <- sqrt(T-t)/sigma
  pd2r <- sqrt(T-t)/sigma
  if (type == "call"){
    phid1 <- pnorm(d1)
    price <- S * exp(-dividend * (T-t)) * phid1 - K * exp(-r * (T-t)) * pnorm(d2)
    delta <- phid1
    gamma <- delta/(S * sigma * sqrt(T-t))
    theta <- S*(exp(-d1^2/2)/sqrt(2*pi))*pd1t - K*(r*exp(-r*(T-t))*pnorm(d2) + exp(-r*(T-t))*exp(-d2^2/2)/sqrt(2*pi)*pd2t)
    vega <- S*(exp(-d1^2/2)/sqrt(2*pi))*pd1sigma - K*exp(-r*(T-t))*(exp(-d2^2/2)/sqrt(2*pi))*pd2sigma
    rho <- S*(exp(-d1^2/2)/sqrt(2*pi))*pd1r - K*(exp(-r*(T-t))*(t-T)*pnorm(d2) + exp(-d2^2/2)/sqrt(2*pi)*exp(-r*(T-t))*pd2r)
    C <- c(price, delta, gamma, theta, vega, rho)
    names(C) <- c("price", "delta", "gamma", "theta", "vega", "rho")
    return(C)
  }
  else if (type == "put"){
    phimd1 <- pnorm(-d1)
    price <- -S * exp(-dividend * (T-t)) * phimd1 + K * exp(-r * (T-t)) * pnorm(-d2)
    delta <- -phimd1
    phid1 <- 1 - phimd1
    gamma <- phid1/(S * sigma * sqrt(T-t))
    theta <- S*(exp(-d1^2/2)/sqrt(2*pi))*pd1t + K*(r*exp(-r*(T-t))*pnorm(-d2) - exp(-r*(T-t))*exp(-d2^2/2)/sqrt(2*pi)*pd2t)
    vega <- S*(exp(-d1^2/2)/sqrt(2*pi))*pd1sigma - K*exp(-r*(T-t))*(exp(-d2^2/2)/sqrt(2*pi))*pd2sigma
    rho <- S*(exp(-d1^2/2)/sqrt(2*pi))*pd1r + K*(exp(-r*(T-t))*(t-T)*pnorm(-d2) - exp(-d2^2/2)/sqrt(2*pi)*exp(-r*(T-t))*pd2r)
    P <- c(price, delta, gamma, theta, vega, rho)
    names(P) <- c("price", "delta", "gamma", "theta", "vega", "rho")
    return(P)
  }
}