Heston PDE Grid: How small can we make it?

The aim is to apply the techniques used in Arbitrary Number of Stencil Points to the partial differential equation (PDE) of the Heston model with x_t = \ln S_t

\displaystyle  0 = \frac{\partial u}{\partial t} +\frac{1}{2}\nu\frac{\partial^2 u}{\partial x^2}  + \left(r_t-q_t -\frac{\nu}{2}\right)\frac{\partial u}{\partial x} + \rho\sigma\nu \frac{\partial^2 u}{\partial \nu \partial x} +\frac{1}{2}\sigma^2\nu\frac{\partial^2 u}{\partial \nu^2} + \kappa(\theta - \nu) \frac{\partial u}{\partial \nu}-r_tu

to get – for a given accuracy – to the smallest possible lattice grids. The Heston model parameters at hand are

\displaystyle S_0 = 100, v_0 = 0.04, \kappa = 1, \theta = v_0, \sigma=0.2, \rho=-0.75, r=5\%, q=0\%,

the strikes of the benchmark options are given by

\displaystyle k\in \{50, 75, 90, 100, 110, 125, 150, 200\}

and the maturity of the options should be one year. First let’s use a five point stencil in x direction. As expected the order of convergence is increasing from second to fourth order.

The same experiment but now with five point stencil also in \nu direction is shown below. Finite lattice effects show-up already at relatively small lattices due to the high convergence speed.

In time direction Richardson extrapolation can be used to increase the order of convergence from two to three.

So finally for this parameter set if the aim is to get the average pricing error below \displaystyle 10^{-3} then a lattice size of

\displaystyle x\times\nu\times T = 50\times 10 \times 20

seems to be sufficient when using five point stencil operators in x and \nu direction and Richardson extrapolation in time direction. Unfortunately for a given accuracy this technique is slower than the usual suspects, e.g. Operator splitting because the matrix inversion for the Crank-Nicolson scheme has to be carried out by an iterative solver, namely by BiCGstab. The code for the numerical experiments can be found in the test case