[R-sig-finance] NIG Option Pricing

[dkf@specere.com]

Hello,

I have been using the NIG-functions in the fBasics package
trying to determine some european call prices using a NIG
distribution, but seemed to have hit a wall.

Basically, the process (as I know it) is to determine the
parameter for the Esscher transform, which is then used in
adjusting the density function for integrating in a BS=style
fashion.

A test case with parameters is included.  The results do not
appear correct as the classic "W" shape of "call(NIG) - call(BS)"
is not observed.  further, if you compare the different call
prices vs. stock/strike, as well as the implicit value - the
results are rather odd, as the call(NIG) value is often less
than the intrinsic value.

There are some debugging routines defined in the following
R-script.

Any help is greatly appreciated.

Regards,

--Dan

test case and parameters
------------------------


library(fBasics)   # for nig distribution

alpha <- 17.32268
beta <- 1.453270
delta <- 0.01144148
mu <- -0.0004846318

mean <- 0.0004780999
sd <- sqrt(delta*(alpha^2)/(sqrt((alpha^2)-(beta^2))^3))

rf <- 0.04

#  determine the parameter, theta, for the riskneutral process
#  under the esscher transform

esscher <- function(x) {
  result <- mu + delta*(sqrt(alpha^2 - (beta + x)^2) - sqrt(alpha^2 -
(beta + x + 1)^2)) - rf
  return(result)
}

lowerlimit <- -alpha - beta
upperlimit <- alpha - beta - 1
xmin <- uniroot(esscher, c(lowerlimit, upperlimit), tol = 1.0e-6)

theta <- xmin$root

# visual test of distributions
#

showpdfs <- function() {
  #
  #  global: alpha, beta, delta, mu, theta, rf
  windows()

  span.min = qnig(0.001, alpha, beta, delta, mu)
  span.max = qnig(0.999, alpha, beta, delta, mu)
  span <- seq(span.min, span.max, length=200)
  yfit <- dnig(span, alpha = alpha, beta = beta, delta = delta, mu = mu)
  yfit2 <- dnig(span, alpha = alpha, beta = beta + theta, delta = delta,
mu = 0)
  yfitnorm <- dnorm(span, mean=mean, sd=sd)

  ylim <- log(c(min(yfit), max(yfit)))  # only bounding on fitted
distribution
  xlim <- c(span[1], span[length(span)])
  plot(span, log(yfit), ylim=ylim, xlim=xlim, type="l", lty=2, col="blue",
ylab="", xlab="logreturns")
  par(new=TRUE)
  plot(span, log(yfit2), ylim=ylim, xlim=xlim, type="l", lty=3, col="red",
ylab="", xlab="logreturns")
  par(new=TRUE)
  plot(span, log(yfitnorm), ylim=ylim, xlim=xlim, type="l", lty=4,
col="green", ylab="", xlab="logreturns")
  legend(x="topleft", lty=c(2, 3, 4), col=c("blue", "red", "green"),
c("objective", "riskneutral", "normal"))
}

# define european call under NIG

nigpdf <- function(x) {
  #
  #  global: alpha, beta, delta, mu, theta, rf
  dnig(x, alpha=alpha, beta=beta, delta=delta, mu=mu)
}

callnig <- function(stock, strike, maturity) {
  #
  #  global: alpha, beta, delta, mu, theta, rf
  if (maturity <= 0) {
    return(max(stock - strike, 0))
  }
  iblower <- log(strike/stock)
  ibupper <- Inf

  # adjust distribution estimates
  muold <- mu
  mu <- 0
  delta <- delta*maturity
  beta <- beta + theta
  intg2 <- integrate(nigpdf, iblower, ibupper)$value
  beta <- beta + 1
  intg1 <- integrate(nigpdf, iblower, ibupper)$value
  result <- stock*intg1 - strike*exp(-rf*maturity)*intg2

  # unwind adjustments
  beta <- beta - theta - 1
  delta <- delta/maturity
  mu <- muold

  return(result)
}

# define european call under B&S

callbs <- function(stock, strike, maturity) {
  #
  # global:  sd
  if (maturity <= 0) {
    return(max(stock - strike, 0))
  }
  d1 <- (log(stock/strike) + ((rf +
(0.5*(sd^2)))*maturity))/(sd*sqrt(maturity))
  d2 <- d1 - sd*sqrt(maturity)
  result <- stock*pnorm(d1) - strike*exp(-rf*maturity)*pnorm(d2)

  return(result)
}

# define function for evaluating calls
#

pfmt5 <- function(x) {
  t1 <- sprintf("%5g", x)
  return(as.numeric(t1))
}

# plot call vs underlying

plotcall <- function(from=27.5, to=32.5, strike=30, maturity=1) {
  npoints <- 200
  sp <- seq(from=from, to=to, length=npoints)
  cbs <- double(npoints)
  cnig <- double(npoints)
  ip <- double(npoints)

  windows()  # open new window for plotting

  for (i in 1:npoints) {
    cnig[i] <- callnig(sp[i], strike, maturity)
    cbs[i] <- callbs(sp[i], strike, maturity)
    ip[i] <- max(sp[i] - strike, 0)   # intrinsic value of call
  }
  plot(sp/strike, ip, ylim=c(0,4), ylab="", xlab="",
    type="l", lty=2, col="blue")
  par(new=TRUE)
  plot(sp/strike, cbs, ylim=c(0,4), ylab="", xlab="",
    type="l", lty=3, col="red")
  par(new=TRUE)
  plot(sp/strike, cnig, ylim=c(0,4), ylab="call price", xlab="(stock
price)/strike",
    type="l", lty=4, col="green")
  legend(x="bottomright", lty=c(2, 3, 4), col=c("blue", "red", "green"),
c("implicit value", "black-scholes", "nig"))
  y <- 3.7
  yinc <- 0.15
  text(0.92, y, adj=0, paste("NIG Parameters:")); y <- y - yinc
  text(0.93, y, adj=0, substitute(alpha==v, list(v=pfmt5(alpha)))); y <- y
- yinc
  text(0.93, y, adj=0, substitute(beta==v, list(v=pfmt5(beta)))); y <- y -
yinc
  text(0.93, y, adj=0, substitute(delta==v, list(v=pfmt5(delta)))); y <- y
- yinc
  text(0.93, y, adj=0, substitute(mu==v, list(v=pfmt5(mu)))); y <- y - yinc
  text(0.93, y, adj=0, substitute(theta==v, list(v=pfmt5(theta)))); y <- y
- yinc
  y <- y - yinc
  text(0.92, y, adj=0, paste("Normal Parameters:")); y <- y - yinc
  text(0.93, y, adj=0, substitute(mu==v, list(v=pfmt5(mean)))); y <- y - yinc
  text(0.93, y, adj=0, substitute(sigma==v, list(v=pfmt5(sd)))); y <- y -
yinc

  title("call(NIG), call(BS), implicit value vs. stock/strike")
}

# plot the difference of the two prices

plotdiff <- function(from=27.5, to=32.5, strike=30, maturity=1) {
  npoints <- 200
  sp <- seq(from=from, to=to, length=npoints)
  cbs <- double(npoints)
  cnig <- double(npoints)
  ip <- double(npoints)

  windows()  # open new window for plotting

  for (i in 1:npoints) {
    cnig[i] <- callnig(sp[i], strike, maturity)
    cbs[i] <- callbs(sp[i], strike, maturity)
  }
  plot((sp/strike), (cnig - cbs), type="l", col="red")
  y <- -0.2
  yinc <- 0.025
  text(0.92, y, adj=0, paste("NIG Parameters:")); y <- y - yinc
  text(0.93, y, adj=0, substitute(alpha==v, list(v=pfmt5(alpha)))); y <- y
- yinc
  text(0.93, y, adj=0, substitute(beta==v, list(v=pfmt5(beta)))); y <- y -
yinc
  text(0.93, y, adj=0, substitute(delta==v, list(v=pfmt5(delta)))); y <- y
- yinc
  text(0.93, y, adj=0, substitute(mu==v, list(v=pfmt5(mu)))); y <- y - yinc
  text(0.93, y, adj=0, substitute(theta==v, list(v=pfmt5(theta)))); y <- y
- yinc
  y <- y - yinc
  text(0.92, y, adj=0, paste("Normal Parameters:")); y <- y - yinc
  text(0.93, y, adj=0, substitute(mu==v, list(v=pfmt5(mean)))); y <- y - yinc
  text(0.93, y, adj=0, substitute(sigma==v, list(v=pfmt5(sd)))); y <- y -
yinc

  title("call(NIG) - call(BS)")
}

plotsurfacediff <- function() {

  windows()

  sp <- seq(from=15, to=45, length=100)
  call <- matrix(nrow=100, ncol=25)
  mat <- seq(from=.2, to=10, length=25)

  strike <- 30

  for (idx in 1:25) {
    for (i in 1:100) {
      # print(paste("idx = ", idx))
      # print(paste("i = ", i))
      call[i, idx] <- callnig(sp[i], strike, mat[idx]) - callbs(sp[i],
strike, mat[idx])
      # call[i, idx] <- callbs0(sp[i], strike, mat[idx], 0.3, rf, rf)
      # call[i, idx] <- callbs(sp[i], strike, mat[idx])
      # call[i, idx] <- callbs2(sp[i], strike, mat[idx])
      # call[i, idx] <- callnig(sp[i], strike, mat[idx])
    }
  }

  persp(sp, mat, call, theta=30, phi=30, expand=0.8, col="lightblue",
ticktype = "detailed",
           xlab = "stock price", ylab = "maturity", zlab = "call price")
}

#  these just do the default.

plotcall()

plotdiff()

showpdfs()

plotsurfacediff()