Please find here my talk about the Collocating Local Volatility Model at this year’s QuantLib User Meeting in Düsseldorf, Germany. The source code is also available.

# Category: Finite Difference Methods

# The Collocating Local Volatility Model

In a recent paper [1] Lech Grzelak has introduced his Collocating Local Volatility Model (CLV). This model utilises the so called collocation method [2] to map the cumulative distribution function of an arbitrary kernel process onto the true cumulative distribution function (CDF) extracted from option market quotes.

Starting point for the collocating local volatility model is the market implied CDF of an underlying at time :

The prices can also be given by another calibrated pricing model, e.g. the Heston model or the SABR model. To increase numerical stability it is best to use OTM calls and puts.

The dynamics of the spot process should be given by some stochastic process and a deterministic mapping function such that

The mapping function ensures that the terminal distribution of given by the CLV model matches the market implied CDF. The model then reads

The choice of the stochastic process does not influence the model prices of vanilla European options – they are given by the market implied terminal CDF – but influences the dynamics of the forward volatility skew generated by the model and therefore the prices of more exotic options. It is also preferable to choose an analytical trackable process to reduce the computational efforts.

The collocation methods outlined in [2] defines an efficient algorithm to approximate the mapping function based on a set of collocation points for a given set of maturities and interpolation points per maturity . Let be the market implied CDF for a given expiry . Then we get

for the collocation points with .

The optimal collocation points are given by the abscissas of the Gaussian quadrature for the terminal distribution of . The simplest choice is a normally distribute kernel process like the Ornstein-Uhlenbeck process

.

The corresponding collocation points of the Normal-CLV model are then given by

in which the collocation points of the standard normal distribution can be calculated by QuantLib’s Gauss-Hermite quadrature implementation

Array abscissas = std::sqrt(2)*GaussHermiteIntegration(n).x()

Lagrange polynomials [3] are an efficient interpolation scheme to interpolate the mapping function between the collocation points

Strictly speaking Lagrange polynomials do not preserve monotonicity and one could also use monotonic interpolation schemes supported by QuantLib’s spline interpolation routines. As outlined in [2] this method can also be used to approximate the inverse CDF of an “expensive” distributions.

Calibration of the Normal-CLV model to market prices is therefore pretty fast and straight forward as it takes the calibration of .

Monte-Carlo pricing means simulating the trackable process and evaluate the Lagrange polynomial if the value of the spot process is needed. Pricing via partial differential equation involves the one dimensinal PDE

with the terminal condition at maturity time

For plain vanilla options the upper and lower boundary condition is

**Example 1: Pricing error for plain vanilla options**

- Market prices are given by the Black-Scholes-Merton model with

.

- Normal-CLV process parameters are given by

Ten collocation points are used to define the mapping function and the time to maturity is one year. The diagram below shows the deviation of the implied volatility of the Normal-CLV price based on the PDE solution from the true value of 25%

Even ten collocation points are already enough to obtain a very small pricing error. The advice in [2] to stretch the collocation grid has turned out to be very efficient if the number of collocation points gets larger.

**Example 2: Forward volatility skew**

- Market prices are given by the Heston model with

.

- Normal-CLV process parameters are given by

The diagram below shows the implied volatility of an forward starting European option with moneyness varying from 0.5 to 2 and maturity date six month after the reset date.

The shape of the forward volatility surface of the Normal-CLV model shares important similarities with the surfaces of the more complex Heston or Heston Stochastic Local Volatility (Heston-SLV) model with large mixing angles . But by the very nature of the Normal-CLV model, the forward volatility does not depend on the values of or , which limits the variety of different forward skew dynamics this model can create. CLV models with non-normal kernel processes will support a greater variety.

**Example 3: Pricing of Double-no-Touch options**

- Market prices are given by the Heston model with

.

- Normal-CLV process parameters are given by different values and

Unsurprisingly the prices of 1Y Double-no-Touch options exhibit similar patterns with the Normal-CLV model and the Heston-SLV model as shown below in the “Moustache” graph. But the computational efforts using the Normal-CLV model are much smaller than the efforts for calibrating and solving the Heston-SLV model.

The QuantLib implementation of the Normal-CLV model is available as a pull request #117, the Rcpp based package Rclv contains the R interface to the QuantLib implementation and the demo code for all three examples.

[1] A. Grzelak, 2016, The CLV Framework – A Fresh Look at Efficient Pricing with Smile

[2] L.A. Grzelak, J.A.S. Witteveen, M.Suárez-Taboada, C.W. Oosterlee,

The Stochastic Collocation Monte Carlo Sampler: Highly efficient sampling from “expensive” distributions

[3] J-P. Berrut, L.N. Trefethen, Barycentric Lagrange interpolation, SIAM Review, 46(3):501–517, 2004.

# QuantLib User Meeting 2014

Please find here Johannes and my talk about “Stochastic Local Volatility Calibration in QuantLib” given at year’s QuantLib User Meeting in Düsseldorf

# Probability Distribution of the Heston Model, Part I

The Heston model is defined by the following stochastic differential equation of the log spot

To a significant extent the popularity of the Heston model is based on the fact that semi-closed formulas for vanilla European options exist using the characteristic function of the model. The time evolution of the probability density function is given by the corresponding Fokker-Planck equation [1]

with the initial condition

The reduced probability density function

for this initial value problem can be calculated using a semi-closed integral formula [2]

This gives the opportunity to write a pricing engine for arbitrary European payoffs. The value of an European option with payoff function at maturity is given by

The calculation needs two nested integrations which can be carried out efficiently using e. g. the Gauss-Lobatto algorithm. The solution of the equation

determines the upper boundary for the integration over . The boundaries for the integration over are chosen such that the interval covers ten times the expected variance

Obviously the nested integration makes this algorithm more tricky than the standard ways to price plain vanilla European options but it is not limited to vanilla payoffs. The implementation of this algorithm can be found here within the QuantLib trunk on Github.

Broadie and Kaya [1] have outlined an algorithm to sample from the full probability density function instead of the reduced density function . Starting point for this algorithm is the exact solution of the Heston stochastic differential equation

The probability density function of the variance process is given by a noncentral chi-squared distribution

The distribution of the integral conditional on and can be calculated via the characteristic function

The modified Bessel function of first kind can be evaluated using series expansion for small and medium or asymptotic approximation for large [5]. Unfortunately Boost provides only real versions of the Bessel functions and the copyright status of older complex valued Fortran77 routine is vague. Therefore QuantLib comes with its own implementation.

Please notice that is already a continuous version of the characteristic function and therefore the integration does not need to track the branches of when calculating the complex valued Bessel function [4].

The integration over the characteristic function is best done using either Gauss-Laguerre, Gauss-Lobatto or trapezoidal rule. The two later algorithms need to truncate the integration at some upper bound. First guess for a truncation limit can be taken from the Cornish-Fisher expansion for some very small . The moment-generating function can be used to get the first, second and third moment of the distribution via finite difference quotient.

The next term is now fairly easy to calculate

The log spot process can now be sampled using a standard normal random variable and

This sampling algorithm is exact even for very large time steps and therefore gives some advantages for quasi random Monte-Carlo methods but the inversion of the integration of the characteristic function is also very slow. The algorithm is implemented within the HestonProcess class.

[1] I. Clark, Foreign Exchange Option Pricing: A Practitioners Guide, p. 113

[2] A. Dragulescu, V. Yakovenko, Probability distribution of returns in the Heston model with stochastic volatility

[3] M. Broadie, Ö. Kaya, Exact Simulation of Stochastic Volatility and other Affine Jump Diffusion Processes

[4] R. Lord, Efficient pricing algorithms for exotic derivatives, p. 40

[5] J.R. Culham, Bessel Functions of the First and Second Kind

# QuantLib User Meeting 2013

Please find here my talk “*Beyond Simple Monte-Carlo: Parallel Computing with QuantLib*” at this year’s QuantLib User Meeting in Düsseldorf.

# Fokker-Planck Equation, Feller Constraint and Boundary Conditions

The Fokker-Planck forward equation is an important tool to calibrate local volatility extensions of stochastic volatility models, e.g. the local vol extension of the Heston model. But the treatment of the boundary conditions – especially at zero instantaneous variance – is notoriously difficult if the Feller constraint

of the square root process for the variance

is violated. The corresponding Fokker-Planck forward equation

describes the time evolution of the probability density function *p *and the boundary at the origin is instantaneously attainable if the Feller constraint is violated. In this case the stationary solution [1] of the Fokker-Planck equation

diverges with

and the zero flow boundary condition at the origin

becomes hard to track numerically. Rewriting the partial differential equation in terms of helps to mitigate the problem [2].

The corresponding zero flow boundary condition in *q *reads as follows

The first step in order to apply the boundary condition is to discretize the transformed Fokker-Planck equation in space direction on a non-uniform grid

with

Next step is to discretize the zero flow boundary condition using a second order forward discretization [3] for the derivative .

with

Now the value at the boundary can be expressed in terms of and .

This equation can be used to remove the boundary value from the discretization of the transformed Fokker-Planck equation. In the same manner the zero flow boundary condition can be implemented for the upper boundary and for the original Fokker-Planck equation.

The following parameters of the square-root process will be used to test the different boundary conditions.

The uniform grid has 100 or 1000 grid points and the stationary solution of the Fokker-Planck equation is chosen as the start configuration on the grid points. The boundaries are defined via the inverse cumulated distribution function such that

.

The grid covers 98% of the distribution and this value would not change over time if the boundary conditions are satisfied without discretiszation errors. To test the magnitude of the discretization errors we let the initial solution evolve over one year using 100 time steps with the Crank-Nicolson scheme. The resulting solution is interpolated using cubic splines to get the function . The value

servers as a quality indicator for the discretization errors of the boundary condition.

As can be seen in the diagram below if is small then the original Fokker-Planck equation shows smaller discretization errors than the transformed Fokker-Planck equation in . But if becomes larger and especially if the Feller constraint is violated then the solution of the transformed Fokker-Planck equation clearly outperforms the original solution. As a rule of thumb based on this example one should use the original Fokker-Planck equation if and the transformed equation otherwise.

The same techniques can be used to solve the Fokker-Planck equation for the Heston model. The code for the numerical tests is available in the latest QuantLib version from the SVN trunk. The diagram is based on the test FDHestonTest::testSquareRootEvolveWithStationaryDensity.

[1] A. Dragulescu, V. Yakovenko, Probability distribution of returns in the Heston model with stochastic volatility

[2] V. Lucic, Boundary Conditions for Computing Densities in Hybrid Models via PDE Methods

[3] P. Lermusiaux, Numerical Fluid Mechanics, Lecture Note 14, Page 5

# Fokker-Planck Forward Equation for the Heston Model

Consider the stochastic differential equation of the Heston model for the dynamics of the log-spot

The Fokker-Planck partial differential forward equation describes the time evolution of the probability density function

with the initial condition

where denotes the Dirac delta distribution. A semi-closed form solution for this problem is presented in [1]. When solving the Fokker-Planck forward equation special care must be taken with respect to the boundary conditions, especially if the Feller constraint

is violated. In this case the boundary at the origin is instantaneously attainable. A generalisation of Feller’s zero-flux boundary condition should be applied at the origin [2].

A three point forward differentiation formula can be used to calculate a second order accurate approximation of the partial derivative for [3].

The diagram below shows the solution of the forward equation for the model

The Feller constraint is violated for and this changes the shape of the solution completely.

The code for this example is available here and it is based on the latest QuantLib version from the SVN trunk. It also depends on RInside and Rcpp to generate the plots. In addition the zip contains a short movie clip of the time evolution of the solution for .

[1] A. Dragulescu, V. Yakovenko, Probability distribution of returns in the Heston model with stochastic volatility

[2] V. Lucic, Boundary Conditions for Computing Densities in Hybrid Models via PDE Methods

[3] K. A. Kopecky, Numerical Differentiation