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()
NIG Option Pricing
4 messages · dkf@specere.com, Krishna Kumar, Martin Maechler
Hello,
I am not sure of a few things you are doing
As it turns out - neither was I. I was able to narrow the problem down and identify my programing error. I apologize for the noise on the list. Anyway, to those that might be interested (and to prevent anyone else wasting their time on this). Basically, my mistake was to assume that the default arguments to the function "nigpdf" were evaluated at runtime. In fact, they are evaluated at the time of the function defnition. Hence both integrate calls were returning identical values. The fix is, change the calls to integrate from: intg2 <- integrate(nigpdf, iblower, ibupper)$value .... intg2 <- integrate(nigpdf, iblower, ibupper)$value To: intg2 <- integrate(nigpdf, iblower, ibupper, alpha=alpha, beta=beta, delta=delta, mu=0)$value .... intg2 <- integrate(nigpdf, iblower, ibupper, alpha=alpha, beta=beta, delta=delta, mu=0)$value In the callnig() function. Regards, --Dan
The BS volatility in your function below is what you compute as "sd" below I don't think that makes sense. You'd have to imply the bs volatility from NiG call prices to compare NiG to B.S So you imply b.s. volatility and then price with that implied to get the B.S price. Also the reason your implicit values are odd is because "max" by default does the maximum across the list/array/vector object. What you probably want to do is to use "pmax" instead. Finally you would find AVT and wim schoutens posting on wilmott helpful. http://www.wilmott.com/messageview.cfm?catid=8&threadid=14313 Cheers, Kris
I am not sure of a few things you are doing The BS volatility in your function below is what you compute as "sd" below I don't think that makes sense. You'd have to imply the bs volatility from NiG call prices to compare NiG to B.S So you imply b.s. volatility and then price with that implied to get the B.S price. Also the reason your implicit values are odd is because "max" by default does the maximum across the list/array/vector object. What you probably want to do is to use "pmax" instead. Finally you would find AVT and wim schoutens posting on wilmott helpful. http://www.wilmott.com/messageview.cfm?catid=8&threadid=14313 Cheers, Kris
dkf at specere.com wrote:
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()
_______________________________________________ R-sig-finance at stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-finance
"dkf" == dkf <dkf at specere.com>
on Wed, 7 Dec 2005 00:23:34 +0700 (ICT) writes:
dkf> Hello,
>> I am not sure of a few things you are doing
dkf> As it turns out - neither was I. I was able to narrow
dkf> the problem down and identify my programing error. I
dkf> apologize for the noise on the list.
dkf> Anyway, to those that might be interested (and to
dkf> prevent anyone else wasting their time on this).
dkf> Basically, my mistake was to assume that the default
dkf> arguments to the function "nigpdf" were evaluated at
dkf> runtime. In fact, they are evaluated at the time of the
dkf> function defnition.
not really -- since this is not true for any R function!
I think your problem is that you are working with global
variables a lot, something which is *very strongly* discouraged....
In your case, I think you could pass a list(alpha=.., beta=..)
to the corresponding functions.
Regards,
Martin