Computing option prices for the Heston model

This note sketches how to compute option prices under the Heston model for an array of strikes and maturities. The example is based on the implementation described in Gilli and Schumann (2010/2011). (Many more explanations are in Gilli (2019), tough mostly with MATLAB code.)

The NMOF manual has for a long time¹ contained a function priceMatrix, with an example for pricing with the characteristic function under Black–Scholes–Merton .

priceMatrix <- function(cf, S, Xvec, tauvec,
                        r, q = 0, ...,
                        nodes = NULL, weights = NULL,
                        n = 200) {

    if (is.null(nodes)) {
        tmp <- xwGauss(n)
        tmp <- changeInterval(tmp$nodes, tmp$weights,
                              oldmin = -1, oldmax = 1,
                              newmin =  0, newmax = 200)
        nodes <- tmp$nodes
        weights <- tmp$weights
    }

    callprices <- array(NA_real_, dim = c(length(Xvec), length(tauvec)))
    tmpmat <- array(NA_real_, dim = c(length(Xvec), length(weights)))

    inodes <- 1i * nodes
    itau <- 0L
    for (tau in tauvec) {
        itau <- itau + 1L
        cfi <- S * exp((r - q) * tau)
        cf1 <- cf(nodes - 1i, S, tau, r, q, ...)/inodes/cfi
        cf2 <- cf(nodes,      S, tau, r, q, ...)/inodes
        iX <- 0L
        for (X in Xvec) {
            iX <- iX + 1L
            if (itau == 1L)
                tmpmat[iX, ] <- exp(-inodes * log(X))
            P1 <- 0.5 + weights %*% Re(tmpmat[iX, ] * cf1)/pi
            P2 <- 0.5 + weights %*% Re(tmpmat[iX, ] * cf2)/pi
            callprices[iX, itau] <-
                exp(-q*tau) * S * P1 - exp(-r*tau) * X * P2
        }
   }
    callprices
}

We fix settings for the Heston model.

library("NMOF")

  Xvec <- 80:120
    tauvec <- c(c(3, 6, 9)/12,  ## 3, 6, 9 months
                1, 2, 3, 4, 5)  ## 1..5 years

    S <- 100
    r <- 0
    q <- 0
    v0 <- 0.2^2  ## variance, not volatility
    vT <- 0.2^2  ## variance, not volatility
    v <- vT
    rho <- -0.3
    k <- 0.5
    sigma <- 0.9

    n <- 50  ## number of nodes

Now for pricing.

## VARIANT 1: using priceMatrix
ans1 <- priceMatrix(cf = cfHeston, S, Xvec, tauvec, r, q = q,
                    v0 = v0, vT = vT, rho = rho, k=k, sigma = sigma,
                    n = n)
colnames(ans1) <- paste("tau", tauvec)
rownames(ans1) <- Xvec

## VARIANT 2: using callHestoncf
ans2 <- array(NA_real_, dim = c(length(Xvec), length(tauvec)))
colnames(ans2) <- paste("tau", tauvec)
rownames(ans2) <- Xvec


for (x in Xvec)
    for (tau in tauvec)
        ans2[x == Xvec, tau == tauvec] <-
            callHestoncf(X = x, tau = tau,
                         S = S, r = r, q =q, v0 = v0, vT = vT,
                         rho = rho, k = k, sigma = sigma,
                         rel.tol = .Machine$double.eps^0.5)

message("Max. rel. difference ",
        round(max(abs(ans1 - ans2)/ans2), 6))

Max. rel. difference 0.000609

## SPEED COMPARISON
library("rbenchmark")
benchmark(
  with_priceMatrix = priceMatrix(cf = cfHeston, S, Xvec, tauvec, r, q = q,
                                 v0 = v0, vT = vT, rho = rho, k = k, sigma = sigma,
                                 n = n),
  with_callHestoncf = for (x in Xvec)
                        for (tau in tauvec)
                          callHestoncf(X = x, tau = tau,
                                       S = S, r = r, q =q, v0 = v0, vT = vT,
                                       rho = rho, k = k, sigma = sigma,
                                       rel.tol = .Machine$double.eps^0.5),
  replications = 100, order = "relative")[, 1:4]

               test replications elapsed relative
1  with_priceMatrix          100   0.353    1.000
2 with_callHestoncf          100  24.954   70.691

@INCOLLECTION{Gilli2010,
  author       = {Manfred Gilli and Enrico Schumann},
  title        = {Calibrating the {H}eston Model with {D}ifferential
		  {E}volution},
  booktitle    = {EvoApplications 2010, Part II},
  publisher    = {Springer},
  year         = 2010,
  editor       = {Di Chio, Cecilia and Brabazon, Anthony and Di Caro,
		  Gianni A.  and Ebner, Marc and Farooq, Muddassar and
		  Fink, Andreas and Grahl, J{\"o}rn and Greenfield,
		  Gary and Machado, Penousal and O'Neill, Michael and
		  Tarantino, Ernesto and Urquhart, Neil},
  number       = 6025,
  series       = {Lecture Notes in Computer Science},
  pages        = {242--250}
}

@INCOLLECTION{Gilli2011,
  author       = {Manfred Gilli and Enrico Schumann},
  title        = {Calibrating Option Pricing Models with Heuristics},
  booktitle    = {Natural Computing in Computational Finance},
  publisher    = {Springer},
  year         = 2011,
  editor       = {Brabazon, Anthony and O'Neill, Michael and Maringer, Dietmar},
  volume       = 4,
}

@BOOK{Gilli2019,
  author       = {Gilli, Manfred and Maringer, Dietmar and Schumann, Enrico},
  title        = {Numerical Methods and Optimization in Finance},
  publisher    = {Elsevier/Academic Press},
  year         = 2019,
  edition      = 2
}

Footnotes:

The first commit in Git in 2012 (before that, I had used Subversion) already contains the function, and git log -S informs me that I have not touched the function since.