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

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

HPC-QuantLib

Fokker-Planck Equation, Feller Constraint and Boundary Conditions

Leave a Reply Cancel reply

Share this:

Like this:

Related

Leave a Reply Cancel reply