# [R-sig-finance] NIG Option Pricing

dkf@specere.com dkf at specere.com
Tue Dec 6 04:31:19 CET 2005

```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

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

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()

```