Implementation of the Heston model in QuantLib

The QuantLib derivatives pricing library provides an algorithm for "analytic" pricing of European-style options under the Heston model. The description "analytic" is conventional but not very precise as the algorithm in fact involves numerical evaluation of an integral. There is however no stepping in time in this algorithm, and so it is certainly distinct from the finite-difference methods that are used to price path-dependent options.

In this article I provide some notes on the implementation of this algorithm in QuantLib.

Strategy

The overall strategy of the algorithm is to evaluate the closed-form solution for the price by integrating, using a quadrature algorithm, the Fourier integrals appearing in the Heston formula.

The caller of the algorithm can make a choice between a number of integration algorithms, including algorithms with a fixed number of points and adaptive algorithms. The default quadrature is the Gauss-Laguerre rule when the caller does not specify a maximum number of integrations, and the Gauss-Lobatto adaptive refinement algorithm when the maximum number of evaluations is provided.

Two alternative for calculating the value of the complex logarithm appearing in the forumula are provided: a branch-counting algorithm and Gatherals' algorithm which is the default.

Code structure

The bulk of code for the pricing is in files ql/pricingengines/vanilla/analytichestonengine.hpp and ql/pricingengines/vanilla/analytichestonengine.cpp. The only substantial code outside of these files is the code for the numerical integration which is in gaussianquadratures.cpp, gaussianquadratures.hpp, gausslobattointegral.cpp and gausslobattointegral.hpp, all of which are in the directory ql/trunk/ql/math/integrals.

The main wrapper class for the pricing calculation is AnalyticHestonEngine. Besides the constructors and the function for the actual pricing, its public interface contains following significant elements:

  • A sub-class Integration which defines the algorithm to be used for the numerical integration
  • The enumeration ComplexLogFormula which defines the algorithm to be used for computation of the complex logarithm
  • A protected virtual function addOnTerm which allows easy implementation of subclasses closely related to the Heston model (e.g., see AnalyticHestonHullWhiteEngine)

The actual calculation is encapsulated in the static function doCalculation. Although the function is static, it takes as an argument a pointer to an instance of the AnalyticHestonEngine class (or its subclass of course) and therefore it is not entirely decoupled from the instances of AnalyticHestonEngine.

The doCalculation makes use of an instance of the Integration subclass and a helper functor class AnalyticHestonEngine::Fj_Helper. The single Fj_Helper class defines the values of the integrand both when calculating the integral for "P1" (effectively the delta) and P2 (effectively the risk-neutral probability that spot exceed strike at expiry, see Mikhailov+Nögel), and contains most of the algebraic complexity.

Mathematical details

  1. The upper limit to which the evaluate the integral to is estimated based on the initial volatility, volatility of volatility, the mean volatility and the time to expiration