Intraday, high Resolution Day Counters

January 4, 2015January 10, 2015 hpcquantlibLeave a comment

At the time being QuantLib’s smallest time resolution is a single day and QuantLib does not support “intraday pricing”. Close to an option’s maturity date this inaccuracy leads to wrong greeks. Especially at maturity date QuantLib gives back zero option npv and greeks. A temporary solution would be to shift the expiry date of the option one day forward but this still leads to wrong valuation of the price and the pin risk of the option.

The root of the problem is DayCounter’s yearFraction method, which acts on dates and the Date class does not support any time resolution smaller than days. On the other hand the Date class is probably the most widely used QuantLib class and therefore any solution must be backwards compatible. One solution for this problem has recently been discussed here.

Instead of adding new derived Date classes another solution will be to change the Date class itself to cope with intraday time information. The advantage of this approach is that all overloaded operators remain consistent. The date time arithmetic has been implemented using the boost::posix_time::ptime class. If one is interested in the correct treatment of time zones and day light saving one should use boost::local_time::local_date_time but this comes with a performance penalty. In order to allow backwards compatibility all signatures and behavior of the existing Date and DayCounter classes should stay the same. Additional methods are inter alia a new high-resolution date constructor

Date::Date(Day d, Month m, Year y,
           Size hours, Size minutes, Size seconds,
           Size millisec = 0, Size microsec = 0);

and methods to get the differences in days between two points in time including the fractions of the days. The maximal resolution of the methods is either micro or nano seconds depending on the underlying boost installation.

   Time Date::fractionOfDay() const
   Time daysBetween(const Date&, const Date&);

Last bit missing now is to change the corresponding day counters, e.g. the yearFraction method of the Actual365Fixed day counter becomes

Time Actual365Fixed::Impl::yearFraction(
    const Date& d1, const Date& d2, 
    const Date&, const Date&) const {

        return daysBetween(d1, d2)/365.0;
}

The new Date and DayCounter implementation is available here. It acts as a drop-in replacement for the existing classes, meaning the QuantLib test suite runs properly without any changes. Only the day counters Actual360, Actual365Fixed, ActualActual allow a strictly monotone definition of time and only these once have been adapted. The patch also contains an example of an intraday pricing of an ATM option during the last two and a half hours of the last trading day based on the Heston model and the finite difference pricing engine FDHestonVanillaEngine. The resulting Gamma and NPV of the option are shown in the diagram below.

QuantLib User Meeting 2014

December 8, 2014December 20, 2016 hpcquantlibLeave a comment

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

QuantLib-SWIG Patch for JVM/.NET Languages

May 11, 2014November 28, 2015 hpcquantlib4 Comments

Update 23.11.2015: The latest version is now part of the official QuantLib Release 1.7.

Update 22.09.2015: Please find the latest and improved version of the patch for QuantLib 1.6.2 here.

The usage of QuantLib in JVM and .NET languages (e.g. Java/Scala and C#/F#) via the SWIG interface has a known shortcoming. The implementation of QuantLib’s observer pattern does not tolerate a parallel garbage collector running in a different thread. As a result programs are randomly crashing or producing “pure virtual function calls”. A detailed description of this problem can be found e.g. here and within the references.

Please find here a patch for QuantLib 1.4 to fix this issue. It contains

Riccardo’s thread-safe singleton patch
Thread-safe Observer pattern based on boost::signals2
Irakli’s static data initialisation via function call

Installation instructions are included in the readme.txt file.

Probability Distribution of the Heston Model, Part II

February 22, 2014 hpcquantlibLeave a comment

Starting point for a semi-analytical solution of the Fokker-Planck forward equation for the Heston model is the exact sampling algorithm of Broadie and Kaya [1] (for the notation please see [2])

$\begin{array}{rcl} x_t &=& x_0 + m(t) + \sigma(t)Z \nonumber \\ \sigma^2(t) &=& (1-\rho^2)\int_0^t \nu_s ds \nonumber \\ m(t) &=& (r_t-q_t)t - \frac{1}{2}\int_0^t \nu_s ds + \rho\int_0^t \sqrt{\nu_s}dW_s^{(1)} \nonumber \end{array}$

The probability distribution function can be described as

$p(x_t, \nu_t, t) = p(\nu_t, t) p(x_t,t\mid \nu = \nu_t)$

and $p(\nu_t, t)$ is given by a noncentral chi squared distribution. The distribution $p(x_t, t \mid \nu = \nu_t)$ can be calculated using the exact simulation algorithm. In this algorithm the variable $x_t$ is given as a function of two random variables $\int_0^t \nu_s ds$ and $Z$ .

The distribution of $x_t$ can now be derived using the general transformation theorem for random variables: Let X be a random variable with probability density function f. The transformed random variable $Y=h(X)$ has the probability density function

$p(y) = f(h^{-1}(y)) \left| \det \left( \frac{\partial h^{-1}_i(y)}{\partial y_j} \right)\right|$

First step is now to rewrite the exact simulation method in terms of the two random variables

$X_1 = \int_0^t \nu_s ds \ , X_2=Z$ .

The simulation scheme then becomes

$\begin{array}{rcl} x_t &=& x_0+a(t) -\frac{1}{2}X_1 + \frac{\rho\kappa}{\sigma}X_1+\sigma(t) X_2 \nonumber \\\sigma^2(t)&=&(1-\rho^2)X_1 \nonumber \\a(t)&=&(r_t-q_t)t+\frac{\rho}{\sigma}\left( \nu_t-\nu_0-\kappa\theta t \right)\end{array}$

or in terms of the transformed random variable

$\begin{array}{rcl} Y_1 =h_1(X_1,X_2) &=& x_0+a(t)-\frac{1}{2}X_1 + \frac{\rho\kappa}{\sigma}X_1+\sqrt{\left(1-\rho^2\right)X_1}X_2 \nonumber \\ Y_2=h_2(X_1,X_2)&=& X_1 \nonumber\end{array}.$

Let $\phi(x_1)$ be the density function of $X_1$

$\phi(x_1)=\frac{2}{\pi}\int_0^\infty \cos ux_1 \mathrm{Re}(\Phi(u))\mathrm{d}u$

and $X_2$ follows by definition a normal distribution. The joint probability density function of $(X_1, X_2$ ) is then

$f\left( \begin{matrix} x_1 \\ x_2 \end{matrix}\right)=\phi(x_1) \frac{1}{\sqrt{2\pi}}e^{-\frac{x_2^2}{2}}$

with

$\begin{array}{rcl}f\left(h^{-1}(y)\right)&=&f\left(\begin{matrix}y_2 \\ \frac{1}{\sqrt{\left(1-\rho^2\right)y_2}}\left(y_1-x_0-a(t)+\frac{1}{2}y_2-\frac{\rho\kappa}{\sigma}y_2\right)\end{matrix}\right) \nonumber \\ \left|\det \left( \frac{\partial h^{-1}_i(y)}{\partial y_j} \right)\right| &=& \frac{1}{\sqrt{\left(1-\rho^2\right)y_2}} \end{array}.$

This yields to the semi-analytical formula for the solution of the Fokker-Planck equation because by definition $p(x_t,t \mid \nu = \nu_t)$ is the distribution density function of $Y_1$ , which is given by

$p(x_t,t \mid \nu = \nu_t) = \int_0^\infty \mathrm{d}y_2 p(y_1, y_2)\mid_{y_1=x_t} = \int_0^\infty \mathrm{d}y_2 \left[f\left(h^{-1}(y)\right)\left|\det \left( \frac{\partial h^{-1}_i(y)}{\partial y_j} \right)\right| \right]_{y_1=x_t}$

The integration over $y_2$ can be carried out using e.g. the Simpson integral rule together with the Cornish-Fisher expansion, which gives an upper bound for the truncation of the upper limit of the integration.

The contour plots below show the probability density function of the Heston model for some example parametrisations.

$\begin{array}{|c|c|c|c|c|c|c|c|c|} \hline {\rm Parameters} & x_0 & \nu_0 & r & q & \kappa & \theta & \sigma & \rho \\ \hline \hline a & 4.6052 & 0.4 & 5\% & 2.5\% & 1.0 & 0.4 & 0.8 & -75\% \\ \hline b & 4.6052 & 0.4 & 5\% & 2.5\% & 1.0 & 0.4 & 0.8 & \ \ 75\% \\ \hline c & 4.6052 & 0.4 & 5\% & 2.5\% & 1.0 & 0.4 & 0.4 & -75\% \\ \hline d & 4.6052 & 0.4 & 5\% & 2.5\% & 1.0 & 0.4 & 0.4 & \ \ \ \ 0\% \\ \hline \end{array}$

The example code is available here and depends on the upcoming QuantLib version 1.4.

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

[2] K. Spanderen, Probability Distribution of the Heston Model, Part I

[3] R.U. Seydel, Tools for Computational Finance, pp 86

Probability Distribution of the Heston Model, Part I

February 4, 2014September 27, 2014 hpcquantlib1 Comment

The Heston model is defined by the following stochastic differential equation of the log spot $x_t = \ln S_t$

$\begin{array}{rcl} dx_t &=& \left(r_t - q_t - \frac{\nu_t}{2}\right)dt + \sqrt\nu_t dW^{x}_t \nonumber \\ d\nu_t&=& \kappa\left(\theta-\nu_t \right ) dt + \sigma\sqrt\nu_t dW^{\nu}_t \nonumber \\ \rho dt &=& dW^{x}_tdW^{\nu}_t \end{array}$

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 $p(x_t, v_t, t)$ is given by the corresponding Fokker-Planck equation [1]

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

with the initial condition

$p(x,\nu,t=0) = \delta(x-x_0)\delta(\nu-\nu_0)$

The reduced probability density function

$p(x_t, t \mid x_0,\nu_0) = \int_0^\infty p(x_t, \nu_t, t) d\nu$

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

$\begin{array}{rcl} \Gamma &=& \kappa+i\rho\sigma p_x \\ \Omega &=& \sqrt{\Gamma^2 + \sigma^2\left(p_x^2-ip_x \right )} \\ p(x_t, t \mid \nu_0) &=& \int_{-\infty}^\infty \frac{dp_x}{2\pi}\exp\left( ip_x (x_t-x_0 -(r-q)t) -\nu_0 \frac{p_x^2-ip_x}{\Gamma + \Omega \coth\left(\Omega t/2 \right )} \right) \\ && \ \ \ \ \ \ \ \ \times \exp\left(-\frac{2\theta\kappa}{\sigma^2}\ln\left(\cosh\frac{\Omega t}{2} + \frac{\Gamma}{\Omega} \sinh \frac{\Omega t}{2}\right )+\frac{\kappa\Gamma\theta t}{\sigma^2} \right) \\ &=& \int_{-\infty}^\infty \frac{dp_x}{2\pi} \tilde{p}(p_x,t \mid \nu_0) \end{array}$

This gives the opportunity to write a pricing engine for arbitrary European payoffs. The value of an European option with payoff function $P(S_t) \in \mathbb{L}^2$ at maturity $t$ is given by

$\begin{array}{rcl} \text{npv} &=& \int_{-\infty}^\infty dx_t\int_{0}^\infty d\nu_t P(S_0 e^{x_t+(r_t-q_t)t})e^{-r_t t} p(x_t,\nu_t,t) \\ &=& e^{-r_t t}\int_{-\infty}^\infty dx_t P(S_0 e^{x_t+(r_t-q_t)t})p(x_t,t \mid x_0,\nu_0) \end{array}$

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

$| \tilde{p}(p_x, t \mid \nu_0)| = \text{QL\_EPSILON}$

determines the upper boundary for the integration over $p_x$ . The boundaries $\left[ -x_{min}, x_{max}\right]$ for the integration over $x_t$ are chosen such that the interval covers ten times the expected variance

$-x_{min} = x_{max}=10\sqrt{\int_0^{t}E\left[ \nu_t \right ] dt} = 10\sqrt{\theta t + \frac{1}{\kappa}\left(\nu_0-\theta\right)\left(1-e^{-\kappa t}\right)}$

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 $p(x_t, \nu_t, t \mid x_0, \nu_0)$ instead of the reduced density function $p(x_t, t \mid x_0, \nu_0)$ . Starting point for this algorithm is the exact solution of the Heston stochastic differential equation

$\begin{array}{rcl} x_t &=& x_o + (r_t-q_t)t - \frac{1}{2}\int_0^t \nu_s ds + \rho\int_0^t \sqrt{\nu_s} dW_s^{(1)} + \sqrt{1-\rho^2}\int_0^t \sqrt{\nu_s} dW_s^{(2)} \nonumber \\ \nu_t &=& \nu_0 + \kappa\theta t - \kappa \int_0^t \nu_s ds + \sigma \int_0^t\sqrt{\nu_s}dW_s^{(1)} \nonumber \end{array}$

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

$\nu_t = \frac{\sigma^2\left( 1-e^{-\kappa t}\right)}{4\kappa}\chi_d^{'2}\left(\frac{4\kappa e^{-\kappa t}}{\sigma^2\left(1-e^{-\kappa t}\right)}\nu_0\right), d=\frac{4\theta\kappa}{\sigma^2}$

The distribution $\Psi(0,t)$ of the integral $\int_0^t \nu_s ds$ conditional on $\nu_t$ and $\nu_0$ can be calculated via the characteristic function

$\begin{array}{rcl} \text{Pr}(\Psi(t) \le x)&=& \frac{2}{\pi}\int_0^\infty \frac{\sin ux}{u}\text{Re}(\Phi(u)) du \\ \\ \Phi(a)&=& \frac{\gamma(a)e^{-\frac{1}{2}(\gamma(a)-\kappa)t} \left(1-e^{-\kappa t}\right)} {\kappa\left(1-e^{\gamma(a)t}\right)} \exp\left( \frac{\nu_t+\nu_0}{\sigma^2} \left[ \frac{\kappa\left(1+e^{-\kappa t}\right)}{1-e^{-\kappa t}} - \frac{\gamma(a)\left(1+e^{-\gamma(a)t}\right)}{1-e^{-\gamma(a)t}} \right] \right) \\ && \times \frac{I_{0.5d-1} \left( \sqrt{\nu_0\nu_t} \frac{4\gamma (a) e^{-0.5\gamma(a)t}}{\sigma^2\left(1-e^{-\gamma(a)t}\right)}\right)}{ I_{0.5d-1} \left( \sqrt{\nu_0\nu_t} \frac{4\kappa e^{-0.5\kappa t}}{\sigma^2\left(1-e^{-\kappa t}\right)}\right)} \\ && \times \frac{\exp\left((0.5d-1) \left[-\frac{1}{2}\gamma(a)t + \ln \frac{\gamma(a)}{1-e^{-\gamma(a)t}} \right]\right)}{\left(\frac{\gamma(a) e^{-0.5\gamma(a)t} }{ 1-e^{-\gamma(a)t}}\right)^{0.5d-1}} \\ \\ \gamma(a)&=&\sqrt{\kappa^2-2 i \sigma^2 a} \end{array}$

The modified Bessel function of first kind $I_{0.5d-1}(z)$ can be evaluated using series expansion for small and medium $|z|$ or asymptotic approximation for large $|z|$ [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 $\Phi(a)$ is already a continuous version of the characteristic function and therefore the integration does not need to track the branches of $arg(z)$ 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 $\epsilon$ . The moment-generating function $\Phi(a)$ can be used to get the first, second and third moment of the distribution via finite difference quotient.

$m_n = E(X^n) = \frac{d^n}{dy^n}\Phi(x+iy)\Big|_{x=y=0}$

The next term is now fairly easy to calculate

$\int_0^t \sqrt{\nu_s} dW_s^{(1)} = \frac{1}{\sigma}\left( \nu_t - \nu_0 - \kappa\theta t+\kappa \int_0^t \nu_s ds \right)$

The log spot process can now be sampled using a standard normal random variable $Z$ 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

December 1, 2013December 1, 2013 hpcquantlibLeave a comment

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

Multi-Threading and QuantLib

July 26, 2013November 12, 2019 hpcquantlib12 Comments

Update 23.11.2015: The latest version is now part of the official QuantLib Release 1.7.

Update 22.09.2015: Please find the latest version of the patch for QuantLib 1.6.2 here.

Update 28.02.2015: Please find the latest version of the patch for QuantLib 1.5 here.

Update 11.05.2014: Please find the latest version of the patch for QuantLib 1.4 here.

QuantLib is per se not thread-safe. The standard way to utilise more than one core is to spawn several independent processes. Riccardo’s thread-safe singleton patch allows to use QuantLib within multi-threading applications as long as objects aren’t explicitly shared between different threads. In fact this patch turns the singleton pattern into a thread local singleton pattern. One possible use case of this patch is to run the test-suite with a multi-threading test runner to speed it up, e.g. on a [email protected] with four cores plus four HT cores the multi-threading test-suite runs in approx. two minutes whereas the single threaded version takes around eight minutes.

Using QuantLib in Java/Scala/C# or F# applications via the SWIG layer violates the multi-threading requirement because the garbage collector runs in a different thread and therefore QuantLib objects are shared among different threads. This creates problems with QuantLib’s implementation of the observer pattern and is discussed in detail here.

An improved implementation of the observer pattern based on boost::signals2 together with Riccardo’s thread-safe singleton patch and the multi-threading test runner can obtained from github. Under Linux/MacOS use

./configure --enable-thread-safe-observer-pattern

to enable the thread-safe singleton and the thread-safe observer pattern. Under Windows the corresponding preprocessor directives

#define QL_ENABLE_THREAD_SAFE_OBSERVER_PATTERN

are already set in the file userconfig.hpp. In order to enable the boost shared_ptr hook change the preprocessor variable BOOST_SP_ENABLE_DEBUG_HOOKS towards BOOST_SP_ENABLE_DEBUG_HOOKS_2 in the file

boost/smart_ptr/detail/sp_counted_impl.hpp

Background: the original preprocessor variable BOOST_SP_ENABLE_DEBUG_HOOKS changes the memory layout of the class boost::shared_ptr which might lead to problems with other pre-compiled libraries which also use boost::shared_ptr.

The reward for this work is a stable SWIG interface for Java/Scala/C#/F# and a thread local singleton implementation, which allows to use QuantLib within multi-threading applications as long as SWIG/QuantLib objects aren’t explicitly shared among different threads.

R/Finance 2013: Applied Finance with R

May 22, 2013May 22, 2013 hpcquantlibLeave a comment

Please find here my talk about R and QuantLib at this year’s R/Finance 2013 conference in Chicago, IL, USA. The source code of the examples is also available.

Fokker-Planck Equation, Feller Constraint and Boundary Conditions

May 4, 2013February 19, 2014 hpcquantlibLeave a comment

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

$\sigma^2 \le 2\kappa\theta$

of the square root process for the variance $\nu$

$d\nu_t = \kappa(\theta-\nu_t)dt + \sigma\sqrt{\nu_t}dW$

is violated. The corresponding Fokker-Planck forward equation

$\frac{\partial p}{\partial t} = \frac{\sigma^2}{2}\frac{\partial^2}{\partial \nu^2}(\nu p) - \frac{\partial}{\partial \nu}\left(\kappa\left(\theta-\nu \right ) p\right)$

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

$p_\nu = \frac{\eta^\eta}{\Gamma(\eta)}\frac{\nu^{\eta-1}}{\theta^\eta}e^{-\eta\nu/\theta},\ \eta=\frac{2\kappa\theta}{\sigma^2}$

diverges with

$\lim_{\nu\to 0} p_\nu^{s} \sim \nu^{\eta-1} = \nu^{-\alpha},\ \alpha=1-\frac{2\kappa\theta}{\sigma^2}$

and the zero flow boundary condition at the origin

$\left.\left[ \frac{\sigma^2}{2}\frac{\partial}{\partial \nu} (\nu p) + \kappa(\nu-\theta)p\right]\right|_{\nu=0} = 0$

becomes hard to track numerically. Rewriting the partial differential equation in terms of $q_\nu = v^\alpha p_\nu$ helps to mitigate the problem [2].

$\frac{\partial q}{\partial t} = \nu\frac{\sigma^2}{2}\frac{\partial^2}{\partial \nu^2} q + \kappa(\nu+\theta)\frac{\partial}{\partial \nu}q + \frac{2\kappa^2\theta}{\sigma^2}q$

The corresponding zero flow boundary condition in q reads as follows

$\left.\left[ \nu\frac{\sigma^2}{2}\frac{\partial}{\partial \nu} q + \kappa\nu q\right]\right|_{\nu=0} = 0$

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 $\nu_i$

$\frac{\partial q_i}{\partial t} = \alpha_iq_{i-1}+\beta_iq_i+\gamma_iq_{i+1]}$

with

$\begin{array}{rcl} h_i &=& \nu_{i+1}-\nu_i \\ \zeta_i &=& h_{i-1}h_i \\ \zeta^p_i &=& h_i(h_{i-1}+h_i) \\ \zeta^m_i &=& h_{i-1}(h_{i-1}+h_i) \\ \alpha_i &=& \frac{1}{\zeta^m_i}\left(\sigma^2\nu_i-\kappa(\theta+\nu_i)h_i)\right) \\ \beta_i &=& \frac{1}{\zeta_i}\left( -\sigma^2\nu_i+\kappa(\theta+\nu_i)(h_i-h_{i-1}) \right ) +\frac{2\kappa^2\theta}{\sigma^2} \\ \gamma_i &=& \frac{1}{\zeta^p_i}\left( \sigma^2\nu_i + \kappa(\theta+\nu_i)h_{i-1} \right ) \end{array}$

Next step is to discretize the zero flow boundary condition using a second order forward discretization [3] for the derivative $\frac{\partial q}{\partial \nu}$ .

$\left.\left[ \nu\frac{\sigma^2}{2}\frac{\partial}{\partial \nu} q + \kappa\nu q\right]\right|_{\nu=\nu_{i-1}} \approx a_iq_{i-1} + b_iq_{i} + c_iq_{i+1} = 0$

with

$\begin{array}{rcl} \eta_i &=& \frac{1}{h_{i-1}h_i(h_{i-1}+h_i)} \\ a_i &=& - \eta\frac{\sigma^2}{2}\nu_{i-1}\left((h_{i-1}+h_i)^2-h_{i-1}^2 \right ) + \kappa\nu_{i-1} \\ b_i &=& \eta\frac{\sigma^2}{2}v_{i-1}\left(h_{i-1}+h_i \right )^2 \\ c_i &=& -\eta\frac{\sigma^2}{2}v_{i-1} h_{n-1}^2 \end{array}$

Now the value $q_{0}$ at the boundary $\nu_0$ can be expressed in terms of $q_1$ and $q_2$ .

$q_{0} = \frac{-b_1q_{1} - c_1q_{2}}{a_1}$

This equation can be used to remove the boundary value $q_0$ 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 $\nu_{max}$ and for the original Fokker-Planck equation.

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

$\kappa=2.5,\ \theta=0.2,\ \sigma\in\{0.2, 2.0\}$

The uniform grid $\nu_{min} \le \nu_i \le \nu_{max}$ has 100 or 1000 grid points and the stationary solution $p^s_{\nu_i}$ of the Fokker-Planck equation is chosen as the start configuration on the grid points. The boundaries $v_{min}, v_{max}$ are defined via the inverse cumulated distribution function such that

$\int_0^{v_{min}}p^s_\nu d\nu = \int_{v_{max}}^\infty p_\nu^s d\nu = 1\%$ .

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 $q_{i,t=1} = v_i^\alpha p_{i,t=1}$ is interpolated using cubic splines to get the function $q_\nu$ . The value

$\int_{vmin}^{v_{max}}q_\nu\nu^{-\alpha} d\nu - 98\%$

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

As can be seen in the diagram below if $\sigma$ is small then the original Fokker-Planck equation shows smaller discretization errors than the transformed Fokker-Planck equation in $q=v^\alpha p$ . But if $\sigma$ 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 $\frac{2\kappa\theta}{\sigma^2} \ge 2.5$ 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

February 10, 2013September 27, 2014 hpcquantlib3 Comments

Consider the stochastic differential equation of the Heston model for the dynamics of the log-spot $x_t = \ln S_t$

The Fokker-Planck partial differential forward equation describes the time evolution of the probability density function $p(x_t,\nu_t,t)$

with the initial condition

$p(x,\nu,t=0) = \delta(x-x_0)\delta(\nu-nu_0)$

where $\delta$ 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

$\sigma^2 \le 2\kappa\theta$

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

$\left.\left[ \frac{\sigma^2}{2}\frac{\partial}{\partial \nu} (\nu p) + \kappa(\nu-\theta)p + \rho\nu\sigma\frac{\partial p}{\partial x}\right]\right|_{\nu=0} = 0, \ \forall x\in \mathbb{R}^+$

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

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

$\small{T=1.0, S_0=1,r=0.1, q=0.05, \kappa=2.0, \theta=0.4, \rho=-0.75, v_0=1.2, \sigma=0.4}$

The Feller constraint is violated for $\sigma=2.0$ 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 $\sigma=2.0$ .