R/analytic_black_scholes.R

In shill1729/OptionPricer: An option pricer for European and American vanilla stock options

#' Black-Scholes risk-neutral European option price
#'
#' @param spot The underlying share price, S.
#' @param strike The strike price of the option, K.
#' @param maturity The time until expiration in trading years, T.
#' @param rate The risk free rate, r.
#' @param div the dividend yield rate
#' @param volat The annualized standard deviation of log-returns, sigma
#' @param what Which type of option to price.
#' @param current_time The time you want to price at, t.
#' @description {The famous formula we all know and love. The simplest analytic formula for pricing vanilla European options. See pdf (coming) for further details.}
#' @details {Either \eqn{Ke^{-r(T-t)} \Phi(x)-S_0 \Phi(x1)} for puts or \eqn{S_0 \Phi(x1)-Ke^{-r(T-t)}\Phi(x)} for calls. See the pdf (coming) for details on the arguments of the CDFs.}
#' @return numeric
analytic_gbm_price <- function(spot, strike, maturity, rate, div, volat, what = "put", current_time = 0)
{
  # Input handling:
  # current time must be less than or equal to maturity.
  time_check <- (current_time <= maturity)
  price_check <- (spot >= 0)
  vol_check <- (volat >= 0)
  input_check <- (time_check && price_check && vol_check)
  if(input_check)
  {
    # Convert to log-returns space
    beta <- log(strike/spot)
    
    # Get time until expiration T-t, where T is the maturity (trading years until expiration) and t is the current time (usually t=0)
    input_time <- maturity-current_time
    
    # The drift rates under the risk neutral measure
    drift <- (rate-div-0.5*volat^2)*input_time
    drift1 <- (rate-div+0.5*volat^2)*input_time
    
    # The volatility under the risk-neutral measure
    vola <- volat*sqrt(input_time)
    
    # The discount factor, e^{-(r-q)(T-t)}
    discount <- exp(-(rate-div)*input_time)
    
    # The forward strike price: Ke^{-r(T-t)}
    forward_strike <- strike*discount
    
    # Arguments for the standard normal CDF.
    x<-(beta-drift)/vola
    x1<-(beta-drift1)/vola
    
    # Return the appropriate price given the option type/style
    if(what=="put")
    {
      # The BSM put option price Ke^{-r(T-t)} \Phi(x)-S_0 \Phi(x1)
      price <- forward_strike*stats::pnorm(q=x)-spot*stats::pnorm(q=x1)
    } else if(what=="call")
    {
      # The BSM call option price S_0 \Phi(x1)-Ke^{-r(T-t)}\Phi(x)
      price <- spot*(1-stats::pnorm(q=x1))-forward_strike*(1-stats::pnorm(q=x))
    }
  } else
  {
    stop("Bad input: Either t>T, S_0<0, or sigma<0")
  }
  return(price)
}


#' Black-Scholes risk-neutral European option delta, sensitivity to spot
#'
#' @param spot The underlying share price, S.
#' @param strike The strike price of the option, K.
#' @param maturity The time until expiration in trading years, T.
#' @param rate The risk free rate, r.
#' @param div the dividend yield rate
#' @param volat The annualized standard deviation of log-returns, sigma
#' @param what Which type of option to price.
#' @param current_time The time you want to price at, t.
#' 
#' @description {The delta of a vanilla European option's price, i.e. the derivative with respect to the share price.}
#' @return See pdf to come...
analytic_gbm_delta <- function(spot, strike, maturity, rate, div, volat, what = "put", current_time = 0)
{
  
  # Select the input time, T-t
  input_time<-(maturity-current_time)
  # Only drift needed is the one corresponding to alpha+1; (r+0.5*sigma^2)*(T-t)
  drift1 <- (rate-div+0.5*volat^2)*input_time
  vola <- volat*sqrt(input_time)
  beta <- log(strike/spot)
  x <- (beta-drift1)/vola
  if(what == "put")
  {
    delta <- -stats::pnorm(x)
    
  } else if(what == "call")
  {
    delta <- 1-stats::pnorm(x)
  } else
  {
    stop("Analytic Black-Scholes formula is only implemented for basic options: puts and calls")
  }
  return(delta)
}


#' Black-Scholes risk-neutral European option gamma, sensitivity to delta
#'
#' @param spot The underlying share price, S.
#' @param strike The strike price of the option, K.
#' @param maturity The time until expiration in trading years, T.
#' @param rate The risk free rate, r.
#' @param div the dividend yield rate
#' @param volat The annualized standard deviation of log-returns, sigma
#' @param current_time The time you want to price at, t.
#' 
#' @description {The gamma of a vanilla European option's price, i.e. the second derivative with respect to the share price, i.e. the first derivative of delta.}
#' @return See pdf to come...
analytic_gbm_gamma <- function(spot, maturity, strike, rate, div, volat, current_time = 0)
{
  # Select the input time, T-t
  input_time <- (maturity-current_time)
  # Only drift needed is the one corresponding to alpha+1; (r+0.5*sigma^2)*(T-t)
  drift1 <- (rate-div+0.5*volat^2)*input_time
  # Volatility is sigma*sqrt(T-t)
  vola <- volat*sqrt(input_time)
  # Beta is the log-strike space point.
  beta <- log(strike/spot)
  # Both calls and puts have the same gamma
  gamma <- stats::dnorm(x=(beta-drift1)/vola)/(spot*vola)
  return(gamma)
}

#' Black-Scholes risk-neutral European option theta, sensitivity to time until expiration
#'
#' @param spot The underlying share price, S.
#' @param strike The strike price of the option, K.
#' @param maturity The time until expiration in trading years, T.
#' @param rate The risk free rate, r.
#' @param div the dividend yield rate
#' @param volat The annualized standard deviation of log-returns, sigma
#' @param what Which type of option to price.
#' @param current_time The time you want to price at, t.
#' 
#' @description {The theta of a vanilla European option's price, i.e. the derivative with respect to the current time.}
#' @return See pdf to come...
analytic_gbm_theta <- function(spot, strike, maturity, rate, div, volat, what = "put", current_time = 0)
{
  # Select the input time, T-t
  input_time<-(maturity-current_time)
  
  # The drift parameters according the risk-free MM numerarie and forward stock numeraire
  drift1 <- (rate-div-0.5*volat^2)*input_time
  drift2 <- (rate-div+0.5*volat^2)*input_time
  vola <- volat*sqrt(input_time)
  beta <- log(strike/spot)
  x <- (beta-drift1)/vola
  x2 <- (beta-drift2)/vola
  discount <- exp(-(rate-div)*input_time)
  if(what=="put")
  {
    theta <- strike*rate*discount*stats::pnorm(x)-(spot*volat*stats::dnorm(x2))/(2*sqrt(input_time))
    
  } else if(what=="call")
  {
    theta<- -strike*rate*discount*(1-stats::pnorm(x))-(spot*volat*stats::dnorm(x2))/(2*sqrt(input_time))
    
  }
  return(theta/360)
}

#' Black-Scholes risk-neutral European option kappa (or vega), sensitivity to volatility
#'
#' @param spot The underlying share price, S.
#' @param strike The strike price of the option, K.
#' @param maturity The time until expiration in trading years, T.
#' @param rate The risk free rate, r.
#' @param div the dividend yield rate
#' @param volat The annualized standard deviation of log-returns, sigma
#' @param current_time The time you want to price at, t.
#' 
#' @description {The vega of a vanilla European option's price, i.e. the derivative with respect to volat.}
#' @return See pdf to come...
analytic_gbm_vega <- function(spot, strike, maturity, rate, div, volat, current_time = 0)
{
  # Select the input time, T-t
  input_time <- (maturity-current_time)
  # Only drift needed is the one corresponding to alpha+1; (r+0.5*sigma^2)*(T-t)
  drift1 <- (rate-div+0.5*volat^2)*input_time
  vola <- volat*sqrt(input_time)
  beta <- log(strike/spot)
  x <- (beta-drift1)/vola
  # Both calls and puts have the same vega
  vega <- spot*stats::dnorm(x)*sqrt(input_time)/100
  return(vega)
}

#' Black-Scholes risk-neutral European option rho, sensitivity to interest rate
#'
#' @param spot The underlying share price, S.
#' @param strike The strike price of the option, K.
#' @param maturity The time until expiration in trading years, T.
#' @param rate The risk free rate, r.
#' @param div the dividend yield rate
#' @param volat The annualized standard deviation of log-returns, sigma
#' @param what Which type of option to price.
#' @param current_time The time you want to price at, t.
#' 
#' @description {The rho of a vanilla European option's price, i.e. the derivative with respect to risk-free rate.}
#' @return See pdf to come...
analytic_gbm_rho <- function(spot, strike, maturity, rate, div, volat, what = "put", current_time = 0)
{
  # Select the input time, T-t
  input_time <- (maturity-current_time)
  # Only drift needed is the one corresponding to alpha+1; (r-q-0.5*sigma^2)*(T-t)
  drift1 <- (rate-div-0.5*volat^2)*input_time
  vola <- volat*sqrt(input_time)
  beta <- log(strike/spot)
  x1 <- (beta-drift1)/vola
  if(what == "call")
  {
    rho <- strike*input_time*exp(-(rate-div)*input_time)*stats::pnorm(x1)
  } else if(what == "put")
  {
    rho <- -strike*input_time*exp(-(rate-div)*input_time)*stats::pnorm(-x1)
  }
  return(rho/10000)
}

#' Black-Scholes European option pricing and greeks
#'
#' @param spot the underlying share price
#' @param strike strike price of the option
#' @param maturity the maturity of the option contract
#' @param rate the risk-neutral rate
#' @param div the dividend yield
#' @param volat the volatility i.e. annualized standard deviation of log-returns
#' @param what type of option, put or call
#' @param current_time the current time, usually 0
#' @param output either "greeks" or "price
#' 
#' @description {Function calls price and greeks analytic pricing functions and returns all simultaneously in a data.frame}
#' @return data.frame
#' @export analytic_gbm
#'
#' @examples
#' spot <- 100
#' maturity <- 1
#' strike <- 100
#' rate <- 0.05
#' div <- 0
#' volat <- 0.19
#' what <- "call"
#' print(analytic_gbm(spot, strike, maturity, rate, div, volat, what))
analytic_gbm <- function(spot, strike, maturity, rate, div, volat, what, current_time = 0, output = "greeks")
{
  price <- analytic_gbm_price(spot, strike, maturity, rate, div, volat, what, current_time)  
  if(output == "price")
  {
    return(price)
  } else if(output == "greeks")
  {
    delta <- analytic_gbm_delta(spot, strike, maturity, rate, div, volat, what, current_time)
    gamma <- analytic_gbm_gamma(spot, strike, maturity, rate, div, volat, current_time)
    theta <- analytic_gbm_theta(spot, strike, maturity, rate, div, volat, what, current_time)
    vega <- analytic_gbm_vega(spot, strike, maturity, rate, div, volat, current_time)
    rho <- analytic_gbm_rho(spot, strike, maturity, rate, div, volat, what, current_time)
    model_ouput <- data.frame(price = price, delta = delta, gamma = gamma, theta = theta, vega = vega, rho = rho)
    return(model_ouput)
  }
}