# 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 time1 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},
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.