# 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

- 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