Probability Distribution of the Heston Model, Part II

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}$