Hello,I have been using the NIG-functions in the R-languagetrying to determine some european call prices using a NIGdistribution, but seemed to have hit a wall. I came acrossthe some previous threads that touched on this topic here,but didn't addresses this particular case.Basically, the process (as I know it) is to determine the parameter for the Esscher transform, which is then used inadjusting the density function for integrating in a BS=stylefashion.A test case with parameters is included. The results do notappear 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 - theresults are rather odd, as the call(NIG) value is often lessthan the intrinsic value. There are some debugging routines defined in the followingR-script.Any help is greatly appreciated.Regards,--Dantest case and parameters------------------------library(fBasics) # for nig distributionalpha <- 17.32268beta <- 1.453270delta <- 0.01144148mu <- -0.0004846318mean <- 0.0004780999sd <- sqrt(delta*(alpha^2)/(sqrt((alpha^2)-(beta^2))^3))rf <- 0.04# determine the parameter, theta, for the riskneutral process # under the esscher transformesscher <- function(x) { result <- mu + delta*(sqrt(alpha^2 - (beta + x)^2) - sqrt(alpha^2 - (beta + x + 1)^2)) - rf return(result)}lowerlimit <- -alpha - betaupperlimit <- alpha - beta - 1xmin <- 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 NIGnigpdf <- 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&Scallbs <- 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 underlyingplotcall <- 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 <- callnig(sp, strike, maturity) cbs <- callbs(sp, strike, maturity) ip <- max(sp - 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 pricesplotdiff <- 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 <- callnig(sp, strike, maturity) cbs <- callbs(sp, 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, strike, mat[idx]) - callbs(sp, strike, mat[idx]) # call[i, idx] <- callbs0(sp, strike, mat[idx], 0.3, rf, rf) # call[i, idx] <- callbs(sp, strike, mat[idx]) # call[i, idx] <- callbs2(sp, strike, mat[idx]) # call[i, idx] <- callnig(sp, 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()

A few things that may help pairwise max is done with pmax and not max in some of your loops you are doingmax(sp-strike,0) ## line 123 for exampleAlso have you looked at Axel's stuff he had uploaded both a maple and excel sheet here with testresults. It will be easier if you try and compare yours to that. There is also some indicative data in wim's book.

Hello,Here is my answer to my problem. Change: intg2 <- integrate(nigpdf, iblower, ibupper)$value .... intg2 <- integrate(nigpdf, iblower, ibupper)$valueTo: 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)$valueMy mistake was to assume that the default arguments to the function "nigpdf" were evaluatedat runtime. In fact, they are evaluated at the time of the function defnition. Hence both integratecalls were returning identical values.Regarding the max/pmax issue that was pointed out - it seems to have been a function of cutting-n-pastinginto the bullen board editor. Those particular cases were just for plotting for debugging. Thanks to quantie for some suggestions (on and offboard). I did look at the previous thread again andexperimented around with the NIG/excel model that was posted. I did find some discrepancies for findingthe Esscher parameter. The NIG/excel implementation uses the quadratic formula - I used a root findingmethod within the theoretical bounds of the parameter:lowerlimit <- -alpha - betaupperlimit <- alpha - beta - 1xmin <- uniroot(esscher, c(lowerlimit, upperlimit), tol = 1.0e-6)Up to you. There were parameter sets extracted from my data that gave the NIG/exel model problems (alpha=81.65, beta=3.69, delta=0.01034, mu=-0.000123 and rf=0.02) as it seems to return a default value of 0.001 for the call. Regards,--Dan