[R-sig-finance] NIG Option Pricing
Krishna Kumar
kriskumar at earthlink.net
Wed Dec 7 03:13:36 CET 2005
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
>
>
>
More information about the R-sig-finance
mailing list