Re-implemented functions

Author: jaehyukchoi

pkg/R/bachelier_impvol.R

@@ -1,37 +1,37 @@
 #' Calculate Bachelier model implied volatility
 #'
-#' @param type option type either 'call' or 'put'
-#' @param price Price
-#' @param spot current stock price
-#' @param forward forward stock price
-#' @param strike strike price
-#' @param texp time to expiry
+#' @param price (vector of) option price
+#' @param strike (vector of) strike price
+#' @param spot (vector of) spot price
+#' @param texp (vector of) time to expiry
 #' @param intr interest rate
 #' @param divr dividend rate
-#' @return implied vol
+#' @param cpsign call/put sign. 1 for call, -1 for put.
+#' @param forward forward price. If given, forward overrides spot
+#' @param df discount factor. If given, df overrides intr
+#' @return Bachelier implied volatility
+#'
+#' @references Choi, J., Kim, K., & Kwak, M. (2009). Numerical Approximation of the Implied Volatility Under Arithmetic Brownian Motion. Applied Mathematical Finance, 16(3), 261–268. https://doi.org/10.1080/13504860802583436
+#'
+#' @export
+#'
 #' @examples
 #' spot <- 100
 #' strike <- 100
 #' texp <- 1.2
-#' sigma <- 0.2
+#' sigma <- 20
 #' intr <- 0.05
 #' price <- 20
-#' vol <- FER::BachelierImpvol(price=price, spot=spot, strike=strike, texp=texp, intr=intr)
-#' @export
+#' vol <- FER::BachelierImpvol(price, strike, spot, texp, intr=intr)
+#'
 BachelierImpvol <- function(
-  type = "call", price, spot, forward = spot*exp((intr-divr)*texp),
-  strike = forward, texp = 1, intr = 0, divr = 0
+  price, strike = forward, spot, texp = 1,
+  intr = 0, divr = 0, cpsign = 1,
+  forward = spot*exp(-divr*texp)/df,
+  df = exp(-intr*texp)
 ){
 
-  price.forward = price * exp(intr*texp)
-
-  if( type == "call" ) {
-    price.straddle <- 2*price.forward - (forward - strike)
-  } else if( type == "put" ) {
-    price.straddle <- 2*price.forward + (forward - strike)
-  } else if( type == "straddle") {
-    price.straddle <- price.forward
-  }
+  price.straddle = 2*price/df - cpsign*(forward - strike)
 
   # vectors a and b used for rational Chebyshev approximation
   a <- c(3.994961687345134e-1,
@@ -54,24 +54,14 @@ BachelierImpvol <- function(
       -2.067719486400926e2,
        1.174240599306013e1)
 
-  #implied volatility when current stock price
-  #is different from the strike price
-
-  #variable v which is bounded in the range [-1, 1],
-  #since the straddle price is always worth more
-  #than the intrinsic value |F-K|.
   v <- abs( forward - strike ) / price.straddle
-  #transformation of v used for better
-  #approximation of h
 
   nu <- ifelse(v<1e-8, 1/(1+v*v*(1/3 + v*v/5)), v/atanh(v))
 
-  poly.a <- (((((((a[8]*nu+a[7])*nu+a[6])*nu+a[5]))*nu+a[4])*nu+a[3])*nu+a[2])*nu+a[1]
-  poly.b <- (((((((((b[10]*nu+b[9])*nu+b[8])*nu+b[7]))*nu+b[6])*nu+b[5])*nu+b[4])*nu+b[3])*nu+b[2])*nu+b[1]
+  poly.nu <- (((((((a[8]*nu+a[7])*nu+a[6])*nu+a[5]))*nu+a[4])*nu+a[3])*nu+a[2])*nu+a[1]
+  poly.de <- (((((((((b[10]*nu+b[9])*nu+b[8])*nu+b[7]))*nu+b[6])*nu+b[5])*nu+b[4])*nu+b[3])*nu+b[2])*nu+b[1]
 
-  #approximation of h(n)
-  #implied volatility
-  vol <- sqrt(pi*nu/(2*texp)) * price.straddle * (poly.a/poly.b)
+  vol.norm <- sqrt(pi*nu/(2*texp)) * price.straddle * (poly.nu/poly.de)
 
-  return(vol)
+  return(vol.norm)
 }

pkg/R/bachelier_price.R

@@ -1,89 +1,37 @@
 #' Calculate Bachelier model option price
 #'
-#' @param type option type either "call" or "put"
-#' @param price Price
-#' @param spot current stock price
-#' @param forward forward stock price
-#' @param texp time to expiry
+#' @param strike (vector of) strike price
+#' @param spot (vector of) spot price
+#' @param texp (vector of) time to expiry
+#' @param sigma (vector of) volatility
 #' @param intr interest rate
 #' @param divr dividend rate
-#' @param sigma volatility
+#' @param cpsign call/put sign. 1 for call, -1 for put.
+#' @param forward forward price. If given, forward overrides spot
+#' @param df discount factor. If given, df overrides intr
 #' @return option price
+#'
+#' @export
+#'
 #' @examples
 #' spot <- 100
 #' strike <- seq(80,125,5)
 #' texp <- 1.2
-#' sigma <- 0.2
+#' sigma <- 20
 #' intr <- 0.05
-#' price <- FER::BachelierPrice(spot=spot, texp = texp, sigma=sigma, strike=strike, intr=intr)
-#' @export
+#' price <- FER::BachelierPrice(strike, spot, texp, sigma, intr=intr)
+#'
 BachelierPrice <- function(
-  type = "call", spot, forward = spot*exp((intr-divr)*texp),
-  strike = forward, texp = 1, intr = 0, divr = 0, sigma
+  strike = forward, spot, texp = 1, sigma,
+  intr = 0, divr = 0, cpsign=1,
+  forward = spot*exp(-divr*texp)/df,
+  df = exp(-intr*texp)
 ){
-
-  #------------------------------------------------
-  #------------------------------------------------
-  #Inputs:
-  #       type: "call","put","straddle","digital"
-  #       s:  spot price of the underlying asset
-  #       k: strike price
-  #       t: time to maturity
-  #       r: risk-free rate
-  #       sigma: volatility
-  #Outputs:
-  #       price: option price
-  #       delta: option delta
-  #       gamma: option gamma
-  #       vega:  option vega
-  #------------------------------------------------
-  #------------------------------------------------
-
   stdev <- sigma*sqrt(texp)
-  d1 <- (forward-strike) / stdev
-  pnorm.d1 <- pnorm(d1) # normal CDF
-  dnorm.d1 <- dnorm(d1) # normal PDF =exp(-d1*d1/2)/sqrt(2*pi)
-  disc.factor <- exp(-intr*texp)
-
-  if(type=="call"){
-    #Option Price
-    price <- (forward-strike)*pnorm.d1 + stdev*dnorm.d1
-    #Greeks
-    delta <- pnorm.d1 #Delta
-    gamma <- 1 / (sigma*texp)/sqrt(2*pi)*exp(-d1^2/2)#Gamma
-    vega <- sqrt(texp)/sqrt(2*pi)*exp(-d1^2/2)#Vega
-
-  }else if (type=="put"){
-    #Option Price
-    price <- (strike-forward)*(1-pnorm.d1) + stdev*dnorm.d1
-    #Greeks
-    delta <- pnorm.d1 - 1 #Delta
-    gamma <- 1 / (sigma*texp)/sqrt(2*pi)*exp(-d1^2/2)#Gamma
-    vega <- sqrt(texp)/sqrt(2*pi)*exp(-d1^2/2)#Vega
+  stdev[stdev < 1e-32] <- 1e-32
+  dn <- cpsign*(forward-strike) / stdev
 
-  } else if (type=="straddle"){
-    #Straddle price
-    price <- (forward-strike)*(2*pnorm.d1-1) + 2*stdev*dnorm.d1
-    delta <- 2*pnorm.d1 - 1
-    gamma <- 2 / (sigma*texp)/sqrt(2*pi)*exp(-d1^2/2)
-    vega <- 2 * sqrt(texp)/sqrt(2*pi)*exp(-d1^2/2)#Vega
-
-  } else if (type== "digital call"){
-    #Digital call price
-    price <- pnorm(d1)
-    delta <- 1 / sqrt(2*pi)*exp(-d1^2/2)/(sigma*sqrt(texp))
-    gamma <- 1 / sqrt(2*pi)*exp(-d1^2/2)/(sigma^2*texp)
-    vega <- -1 / sqrt(2*pi)*exp(-d1^2/2)*(spot-strike)/(sigma^2*sqrt(texp))
-
-  } else if (type== "digital put"){
-    #Digital put price
-    price <- pnorm(-d1)
-    delta <- -1 / sqrt(2*pi)*exp(-d1^2/2)/(sigma*sqrt(texp))
-    gamma <- 1 / sqrt(2*pi)*exp(-d1^2/2)/(sigma^2*texp)
-    vega <- -1 / sqrt(2*pi)*exp(-d1^2/2)*(spot-strike)/(sigma^2*sqrt(texp))
-  } else {
-    cat("Error! Please input: call/put/straddle/digital call/digital put")
-  }
-  return( disc.factor * price )
+  #Option Price
+  price <- df*(cpsign*(forward-strike)*pnorm(dn) + stdev*dnorm(dn))
+  return( price )
 }
-

pkg/R/blackscholes_impvol.R

@@ -1,50 +1,38 @@
 #' Calculate Black-Scholes implied volatility
 #'
-#' @param type option type either 'call' or 'put'
-#' @param price Price
-#' @param spot current stock price
-#' @param forward forward stock price
-#' @param strike strike price
-#' @param texp time to expiry
+#' @param price (vector of) option price
+#' @param strike (vector of) strike price
+#' @param spot (vector of) spot price
+#' @param texp (vector of) time to expiry
 #' @param intr interest rate
 #' @param divr dividend rate
-#' @return Implied vol
+#' @param cpsign call/put sign. 1 for call, -1 for put.
+#' @param forward forward price. If given, forward overrides spot
+#' @param df discount factor. If given, df overrides intr
+#' @return Black-Scholes implied volatility
+#'
+#' @export
+#'
 #' @examples
 #' spot <- 100
 #' strike <- 100
 #' texp <- 1.2
 #' sigma <- 0.2
 #' intr <- 0.05
 #' price <- 20
-#' vol <- FER::BlackScholesImpvol(price=price, spot=spot, strike=strike, texp=texp, intr=intr)
-#' @export
+#' vol <- FER::BlackScholesImpvol(price, strike, spot, texp, intr=intr)
+#'
 BlackScholesImpvol <- function(
-  type = "call", price, spot, forward = spot*exp((intr-div)*texp),
-  strike = forward, texp = 1, intr = 0, divr = 0
+  price, strike = forward, spot, texp = 1,
+  intr = 0, divr = 0, cpsign = 1,
+  forward = spot*exp(-divr*texp)/df,
+  df = exp(-intr*texp)
 ){
-    price.forward = price * exp( r*texp)
-
-    n.price = length(price.forward)
-    n.strike = length(strike)
-
-    vol <- rep(NA, n.price )
-
-    if( length(strike) > 1L ) {
-      stopifnot(n.price == n.strike )
-      strike.vec <- strike
-    } else {
-      strike.vec <- rep( strike, n.price )
-    }
+  optval = (price/df - pmax(cpsign*(forward-strike), 0))/pmin(forward, strike)
 
-    for(k in 1:n.price) {
-      # Be careful here.... Chekc how functional works in R.
-      sub <- function(sigma){
-        f <- CalcBsmPrice(
-          type = type, forward = forward, strike = strike.vec[k], texp = texp, sigma = sigma
-        )[1] - price.forward[k]
-        return(f)
-      }
-      vol[k] <- uniroot(f = sub,interval = c(0,10), tol = .Machine$double.eps * 1e4)$root
-    }
-    return(vol)
+  # we use inverse CDF of inversegaussian distribution
+  mu <- 2/abs(log(strike/forward))
+  x <- statmod::qinvgauss(optval, mean=mu, lower.tail=F)
+  sig <- 2/sqrt(x*texp)
+  return( sig )
 }

pkg/R/blackscholes_price.R

@@ -1,25 +1,31 @@
 #' Calculate Black-Scholes option price
 #'
-#' @param type option type either 'call' or 'put'
-#' @param spot current stock price
-#' @param forward forward stock price
-#' @param strike strike price
-#' @param texp time to expiry
+#' @param strike (vector of) strike price
+#' @param spot (vector of) spot price
+#' @param texp (vector of) time to expiry
+#' @param sigma (vector of) volatility
 #' @param intr interest rate
 #' @param divr dividend rate
-#' @param sigma volatility
+#' @param type option type either 'call' or 'put'
+#' @param forward forward price. If given, forward overrides spot
+#' @param df discount factor. If given, df overrides intr
 #' @return option price
+#'
+#' @export
+#'
 #' @examples
 #' spot <- 100
 #' strike <- seq(80,125,5)
 #' texp <- 1.2
 #' sigma <- 0.2
 #' intr <- 0.05
-#' price <- FER::BlackScholesPrice(spot=spot, texp = texp, sigma=sigma, strike=strike, intr=intr)
-#' @export
+#' price <- FER::BlackScholesPrice(strike, spot, texp, sigma, intr=intr)
+#'
 BlackScholesPrice <- function(
-  type = "call", spot, forward = spot*exp((intr-divr)*texp),
-  strike = forward, texp = 1, intr = 0, divr = 0, sigma
+  strike = forward, spot, texp = 1, sigma,
+  intr = 0, divr = 0, cpsign = 1,
+  forward = spot*exp(-divr*texp)/df,
+  df = exp(-intr*texp)
 ){
     stdev <- sigma*sqrt(texp)
 
@@ -30,22 +36,6 @@ BlackScholesPrice <- function(
     d1 <- log(forward/strike)/stdev + 0.5*stdev
     d2 <- d1 - stdev
     disc.factor <- exp(-intr*texp)
-
-    pnorm.d1 <- pnorm(d1)
-    pnorm.d2 <- pnorm(d2)
-
-    if (type == "call" ){
-        price <- forward*pnorm.d1 - strike*pnorm.d2
-        delta <- pnorm.d1
-    }else if (type == "put"){
-        price <- strike*(1-pnorm.d2) - forward*(1-pnorm.d1)
-        delta <- pnorm.d1 - 1
-    }else if (type == "straddle"){
-        price <- forward*(2*pnorm.d1 - 1) - strike*(2*pnorm.d2 - 1)
-        delta <- 2*pnorm.d1 - 1
-    }else if (type == "digit"){
-        price <- pnorm.d2
-        delta <- 1/sqrt(2*pi)*exp(-d2^2/2)/(s*stdev)
-    }
-    return( disc.factor * price )
+    price <- df * cpsign*(forward*pnorm(cpsign*d1) - strike*pnorm(cpsign*d2))
+    return( price )
 }