Improved test.

Author: jaehyukchoi

pkg/R/bachelier_impvol.R

@@ -60,6 +60,8 @@ BachelierImpvol <- function(
 
   v <- abs( forward - strike ) / price.straddle
 
+  # when the time value is very slightly negative, we correct to give zero vol.
+  v[price>0 & 1<v & v<1+8*.Machine$double.eps] <- 1
   nu <- ifelse(v<1e-8, 1/(1+v*v*(1/3 + v*v/5)), v/atanh(v))
 
   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]

pkg/R/blackscholes_impvol.R

@@ -33,11 +33,13 @@ BlackScholesImpvol <- function(
   intr=0, divr=0, cp=1L,
   forward=spot*exp(-divr*texp)/df, df=exp(-intr*texp)
 ){
-  optval <- (price/df - pmax(cp*(forward-strike), 0))/pmin(forward, strike)
+  timeval <- (price/df - pmax(cp*(forward-strike), 0))/pmin(forward, strike)
+  # when the time value is very slightly negative, we correct to give zero vol.
+  timeval[price>0 & -8*.Machine$double.eps<timeval & timeval<0] <- 0
 
   # we use inverse CDF of inversegaussian distribution
   mu <- 2/abs(log(strike/forward))
-  x <- statmod::qinvgauss(optval, mean=mu, lower.tail=F)
+  x <- statmod::qinvgauss(timeval, mean=mu, lower.tail=F)
   sig <- 2/sqrt(x*texp)
   return( sig )
 }

pkg/R/blackscholes_price.R

@@ -29,9 +29,9 @@ BlackScholesPrice <- function(
 ){
     stdev <- sigma*sqrt(texp)
 
-    # a trick to get the intrinsic value for negative or zero vol
+    # a trick to get the intrinsic value for zero (or slightly negative) vol
     # also avoid NAN in case forward=strike
-    stdev[stdev < 1e-32] <- 1e-32
+    stdev[abs(stdev) < .Machine$double.eps] <- .Machine$double.eps
 
     d1 <- log(forward/strike)/stdev + 0.5*stdev
     d2 <- d1 - stdev

pkg/tests/testthat/test-BS.R

@@ -0,0 +1,15 @@
+library(FER)
+
+test_that("Bachelier Price and Impvol", {
+  strike <- runif(1000, 50, 150)
+  sigma <- runif(1000, 1, 20)
+  texp <- runif(1000, 1e-4, 10)
+  intr <- runif(1000, 0, 0.1)
+  divr <- runif(1000, 0, 0.1)
+  cp <- ifelse(runif(1000)>0.5, 1L, -1L)
+
+  p <- FER::BachelierPrice(strike, 100, texp, sigma, cp=cp, intr=intr, divr=divr)
+  iv <- FER::BachelierImpvol(p, strike, 100, texp, cp=cp, intr=intr, divr=divr)
+  p2 <- FER::BachelierPrice(strike, 100, texp, iv, cp=cp, intr=intr, divr=divr)
+  expect_equal(p, p2, tolerance = 1e-8)
+})

pkg/tests/testthat/test-Bachelier.R

@@ -0,0 +1,15 @@
+library(FER)
+
+test_that("Black-Scholes Price and Impvol", {
+  strike <- runif(1000, 50, 150)
+  sigma <- runif(1000, 0.001, 2)
+  texp <- runif(1000, 1e-4, 10)
+  intr <- runif(1000, 0, 0.1)
+  divr <- runif(1000, 0, 0.1)
+  cp <- ifelse(runif(100)>0.5, 1L, -1L)
+
+  p <- FER::BlackScholesPrice(strike, 100, texp, sigma, cp=cp, intr=intr, divr=divr)
+  iv <- FER::BlackScholesImpvol(p, strike, 100, texp, cp=cp, intr=intr, divr=divr)
+  p2 <- FER::BlackScholesPrice(strike, 100, texp, iv, cp=cp, intr=intr, divr=divr)
+  expect_equal(p, p2, tolerance = 1e-8)
+})

pkg/tests/testthat/test-BlackScholes.R

@@ -0,0 +1,9 @@
+library(FER)
+
+test_that("CEV", {
+  value <- round( FER::CevPrice(
+    strike=seq(0.4, 2, 0.4)*0.05, spot=0.05, texp=1, sigma=0.4, beta=0.3),
+    digits=5 )
+  value2 <- c(0.04608, 0.04229, 0.03868, 0.03525, 0.03203)
+  expect_equal(value, value2, tolerance = 1e-6)
+})

pkg/tests/testthat/test-CEV.R

@@ -0,0 +1,25 @@
+library(FER)
+
+test_that("SabrHagan2002 BS Volatility Table 1", {
+  value <- round( FER::SabrHagan2002(
+    strike=seq(0.1, 2, 0.1), spot=1, texp=10, sigma=0.25,
+    vov=0.3, rho=-0.8, beta=0.3), digits = 4 )
+  value2 <- c(
+    0.7176, 0.5725, 0.4886, 0.4293, 0.3835,
+    0.3462, 0.3148, 0.2876, 0.2638, 0.2427,
+    0.2238, 0.2068, 0.1916, 0.1781, 0.1663,
+    0.1562, 0.1478, 0.1412, 0.136, 0.1322)
+  expect_equal(value, value2, tolerance = 1e-5)
+})
+
+test_that("SabrHagan2002 BS Volatility Table 13", {
+  value <- round( FER::SabrHagan2002(
+    strike=seq(0.1, 2, 0.1), spot=1, texp=20, sigma=0.25,
+    vov=0.3, rho=-0.5, beta=0.3), digits=4 )
+  value2 <- c(
+    0.7603, 0.5973, 0.5078, 0.4463, 0.3997,
+    0.3626, 0.332, 0.3063, 0.2844, 0.2658,
+    0.2498, 0.2364, 0.2251, 0.2158, 0.2083,
+    0.2023, 0.1976, 0.1941, 0.1914, 0.1895)
+  expect_equal(value, value2, tolerance = 1e-5)
+})