README.md

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

OptionPricer

The goal of OptionPricer is to provide an R-wrapper (via Rcpp) to fast implementations of finite-difference scheme solvers in c++ for solving PDEs, PIDEs, and corresponding variational inequalities arising from the Feynman-Kac connection of stochastic processes.

Installation

You can install the latest version on GitHub with devtools:

devtools::install_github("shill1729/OptionPricer")

Comparing Monte-Carlo and PDE pricers

As an numerical verification of the validity of these algorithms we compare outputs of the PDE and Monte-Carlo schemes. They should be approximately equal. We do an example for each model.

Black Scholes

Using 1.5 million samples in the Monte-Carlo routine and a space-time resolution of 400, we find the prices under the Black-Scholes model are equal within 10^-6 of the true value, in MSE for 30 day options, and for 1 year options, within 10^-5 of the true analytic value.

library(OptionPricer)
spot <- 100
maturity <- 30/252
rate <- 0.05
div <- 0.02
volat <- 0.2
what <- "call"
N <- 400
M <- 400
threshold <- 0.001
samples <- 1.5*10^6
strike_range <- 5
ch <- check_gbm_convergence(spot, maturity, rate, div, volat, what, N, M, threshold, samples, strike_range)
print(ch)

Merton's jump diffusion

We now compare the Monte-Carlo price, the PIDE price and the true semi-analytic price for European options on assets following Merton's jump diffusion.

Remarks on Monte-Carlo accuracy:

The naive Monte-Carlo scheme requries a relatively massive amount of samples for accuracy when compuating expectations of jump-diffusions, at least 5*10^7 for even slightly acceptable results within the good regimes of parameters (short-term, low number of jumps, small jump size compared to volatility). This worsens for longer maturities and when jump dynamics dominate the diffusion volatility. This is a somewhat obvious consequence: longer simulations of jump diffusions have much greater terminal volatility and thus many more samples are needed for the LLN or CLT to hold. Even that lower bound of samples is no guarantee, there is non-neglible probability of vastly large outliers for the terminal variate.

Remarks on PIDE accuracy:

For low jump regimes, a space-time resolution of 200 is typically acceptable within one thousandth tolerance for the MSE. For moderate to high jump regimes, at least 400 is required. Worse, for jump size volatilities thare are too small, the PIDE solution blows up.

library(OptionPricer)
spot <- 100
maturity <- 30/252
rate <- 0.05
div <- 0.02
volat <- 0.5
lambda <- 10
jm <- 0.0
jv <- 0.05
what <- "put"
N <- 200
M <- 200
threshold <- 0.001
samples <- 10^6
strike_range <- 5
parameters <- c(rate, div, volat, lambda, jm, jv)
ch <- check_merton_convergence(spot, maturity, parameters, what, N, M, threshold, samples, strike_range)
print(ch)

Verifying Put Call parity

We can verify the put-call parity for European options for the PDE method and compare the implied rate to the specified rate. At time-space resolution of 400, the implied rate has a MSE on the order of 10^-11 for 1 year maturities for a given range of strikes. At shorter maturities the error is even smaller, about 10^-16 for 30 DTE.

library(OptionPricer)
spot <- 100
strike <- 100
maturity <- 1
rate <- 0.02
div <- 0.0
volat <- 0.25
style <- "european"
output <- "price"
N <- 400
M <- 400
K <- seq(spot-10, spot+10, length.out = 100)
put_pde <- matrix(0, nrow = length(K))
call_pde <- matrix(0, nrow = length(K))
for(i in 1:length(K))
{
  put_pde[i] <- pde_gbm(spot, K[i], maturity, rate, div, volat, what = "put", style, output = output, N, M)  
  call_pde[i] <- pde_gbm(spot, K[i], maturity, rate, div, volat, what = "call", style, output = output, N, M)  
}
r <- -(1/maturity)*(log(spot+put_pde-call_pde)-log(K))
par(mfrow = c(2, 1))
error <- call_pde-put_pde-spot+K*exp(-rate*maturity)
plot(K, error, type = "l", main = "C-P-S+K*D")
plot(K, r, type= "l", main = "PC Parity implied rate", ylim = c(rate-0.00001, rate+0.00001))
abline(h = rate, col = "blue", lty = "dashed")
print("Numerical MSE of implied rate vs true rate:")
print(pside::ci_mean((rate-r)^2))

Calibration

Here are some examples of calibrating a pricing model to market data or synthetic data.

Kou model

Calibration of Kou jump-diffusion model to synthetic data

library(OptionPricer)
spot <- 100
maturity <- 30/252
rate <- 0.05
div <- 0
volat <- 0.2
lambda <- 10
prob <- 0.5
alpha <- 0.08
beta <- 0.08
ku <- 0.01
kd <- -0.01
what <- "put"
style <- "american"
N <- 100
M <- 100
bounds <- list(lower = c(0.01, 0, 0.2, ku, abs(kd)),
               upper = c(10, 90, 1, 0.1, 0.5))
parameters <- c(rate, div, volat, lambda, prob, alpha, beta, ku, kd)
initial_guess <- c(volat, lambda, prob, alpha, beta)*0.9
K <- seq(spot-10, spot, by = 1)
fee <- matrix(0, length(K))
for(i in 1:length(K))
{
  fee[i] <- pide_kou(spot, K[i], maturity, parameters, what, style, "price", N, M, M/2+1)
}
market <- data.frame(strike = K, fee = fee)
plot(market, type = "l")
fit <- calibrate_kou(initial_guess, market, spot, maturity, rate, div, ku, kd, what, style, N, M, TRUE, bounds = bounds)

shill1729/OptionPricer documentation built on June 11, 2020, 12:18 a.m.