An approximation of small-time probability density functions in a general jump diffusion model

Abstract

We propose a method for approximating probability density functions related to multidimensional jump diffusion processes. For small-time horizons, a closed-form approximation of the characteristic function is derived based on the Itô–Taylor expansion. The probability density function is then approximated numerically by inverting the characteristic function using fast Fourier transform. As application we consider a general stochastic volatility model, which involves time-/state-dependent drift and diffusion functions as well as jump components. We test our approach under the Heston model and the Bates model and show that our method provides accurate approximations.

Introduction

In recent years, general jump diffusion volatility models that combine stochastic volatility (e.g., [1]) as well as local volatility (see [2]) or jumps (e.g., [3]) have become popular in the financial industry. For such mixed models (e.g., the universal volatility models in [4]), it remains a challenging task to find the joint probability distribution of the asset price (or return) and its volatility at a certain time point, since the corresponding joint probability density function (PDF) is generally not given in an analytical form and is also difficult to approximate in an efficient numerical way. One popular approach is to obtain the (transition) PDF of the two-dimensional Markov processes that govern the dynamics of the asset price and its volatility by solving the Fokker–Planck equation. The Fokker–Planck equation is either a partial differential equation (PDE), in case the model does not have a jump component, or it is a partial integro-differential equations (PIDE), if jumps are added to the model. Determining the (transition) PDF is useful for likelihood-based statistical inference and for pricing contingent claims. Jensen and Poulsen [5] compare the numerical solution of the Fokker–Planck equation using a finite difference algorithm with various techniques (e.g., the Hermite polynomial approximation to the transition PDF in [6]) for approximating the transition PDF of a time-homogeneous one-dimensional diffusion process. Based on finite difference methods, Andreasen and Huge [7] study the numerical solution to the Fokker–Planck PDE for the transition PDF as well as that to its dual PDE for the option price in a stochastic local volatility model with zero correlation between spot price of the underlying asset and volatility. In [8], [9] the numerical solution to PDEs or PIDEs is conducted for the pricing of contingent claims. These techniques might also be employed on solving the Fokker–Planck equation for the transition PDF.

Authors that derive analytical formulations for the transition PDF exploit the structure of the resulting Fokker–Planck equations, which can be originally attributed to the specific form of the given models. Drǎgulescu and Yakovenko [10] derive an analytical expression for the marginal probability distribution function of log-returns in the Heston model by solving the corresponding Fokker–Planck equation exactly. Lipton [11] solves the Fokker–Planck PDE for the transition PDF in the Heston model through Fourier transform techniques. His result has been extended in [12] to a more general model with price and volatility jumps.

Solving the PDE or PIDE problem numerically suffers usually from the curse of dimensionality. The first and second derivatives in the Fokker–Planck equation have to be discretized numerically and the grid size for discretization approximations must be small enough, for keeping the solution sufficiently accurate. Especially for the multidimensional case, the system of equations to be solved at each iteration tends to be very large. Therefore, it is quite hard in practice to find an efficient way for solving the Fokker–Planck equation numerically. In our study, we focus on approximating the joint PDF in a general jump diffusion model through the corresponding characteristic function, instead of trying to solve PDE or PIDE problems.

Often we do not know the distribution function of given random variables in closed form, but we can derive the characteristic function associated to these random variables explicitly. Let $Z:={(Z_1,\dots ,Z_d)}^T$ be a d-dimensional real random variable with PDF f_Z, the characteristic function ${\Phi }_Z:{\mathbb{R}}^d\to \mathbb{C}$ corresponding to f_Z is defined as Fourier transform of f_Z,

$${\Phi }_Z\left(\Theta \right):=\mathbb{E}\left[\exp \{ı{\Theta }^{T}Z\}\right]={\int }_{{\mathbb{R}}^d}\exp \{ı{\Theta }^T\tilde{z}\}{f}_Z\left(\tilde{z}\right)\mathrm{d}\tilde{z},$$

where Θ is a d-dimensional real vector and ı denotes the imaginary unit. Suppose that the characteristic function Φ_Z and the PDF f_Z are both integrable and Φ_Z is known, we are then able to recover the density function $f_Z\in \mathbb{R}$ through the inverse Fourier transform of Φ_Z:

$$f_Z\left(z\right)=\frac{1}{{\left(2\pi \right)}^d}{\int }_{{\mathbb{R}}^d}\exp \{-ı{\tilde{\Theta }}^{T}z\}{\Phi }_Z\left(\tilde{\Theta }\right)\mathrm{d}\tilde{\Theta }.$$

Accordingly, it can be interpreted that the characteristic function Φ_Z contains exactly the same information as the PDF f_Z. For more details we refer the reader to Shephard [13], who provides a unified framework for the discussion and study of inverting characteristic functions in order to compute distribution functions or density functions. Once the characteristic function is available, the PDF can be obtained by evaluating the integral on the right-hand side of (2). An analytical evaluation of this integral is usually the exception and therefore a numerical approximation is required. For the numerical integration of the integral in (2), the fast Fourier transform (FFT) (see, for instance, [14]) is known as a computationally efficient algorithm.

The question left is then how to get the characteristic function Φ_Z. There are some models that are analytically tractable in the sense that their characteristic functions can be expressed in closed-form. One class of them consists of the so-called affine jump diffusions introduced in [15]. In an affine jump diffusion model, the drifts, covariances and jump intensities are affine in the state vector. Usually, the characteristic function of the state vector in the affine jump diffusion can be derived in closed-form by solving analytically a system of ordinary differential equations of Riccati type. Lévy processes can be regarded as another special class with explicit formulation of the characteristic function given by the Lévy–Khintchine representation. Several models or their transformations fall into either of these two categories, for example, the Heston model can be regarded as a jump-free special case of the affine jump diffusion. However, a broad scope of volatility models belongs to neither of them, e.g., the SABR model introduced in [16] and the universal volatility model in [4].

In our study, we apply the Itô–Taylor expansion to a fairly general multidimensional diffusion setting and obtain a closed-form expression for the characteristic function of the resulting approximation. Subsequently, the formula for the characteristic function for models with jump components is derived. The Itô–Taylor expansion is a stochastic extension of the Taylor formula and is used for approximating functions of stochastic processes. Preston and Wood [17] approximate the moment generating function of one-dimensional diffusion processes by using the Itô–Taylor expansion. The moment generating functions are then used for computing the approximate transition densities by applying the saddlepoint approximation. Studer [18] derives a so-called stochastic delta-gamma approximation, based on the Itô–Taylor expansion, for multidimensional diffusion processes and determines the moment generating function of this approximation, which can then be used by the saddlepoint approximation for calculating the PDF. In the multidimensional case, it is usually impossible to get an analytical solution to so-called saddlepoint equations, because of the high nonlinearity of these equations. In this case, one needs to resort to some efficient solver of systems of nonlinear equations. For more details on the topic of saddlepoint approximations see [19].

In this article, we restrict ourselves to getting the characteristic function or its approximation in a fairly general multivariate jump diffusion model on small-time horizons. Observe that the techniques we are using are inappropriate for the one-dimensional case, for which various approximations are available in the literature, see, e.g., [6], [17]. Our result for the transition PDF can be extended to larger time horizons under the Markovian setting by applying the Chapman–Kolmogorov equation.

The remainder of the paper is organized as follows. Section 2 introduces the general idea of Itô–Taylor expansions as well as two important approximations, the Euler–Maruyama approximation and the Milstein approximation. For both we derive a closed-form expression for the characteristic function of the multidimensional diffusion process. After that, we generalize the results by adding jumps to the diffusion process. In Section 3 we apply our results to a general two-dimensional stochastic volatility model. Section 4 provides computational details for calculating numerically PDFs with known characteristic functions by means of FFT. Section 5 illustrates our method by performing numerical tests under the Heston model and the Bates model, for which the joint characteristic functions of the asset log-return and its variance are known in closed form. Under the Heston model, we also compare the performance of our approach with that of the closed-form transition density approximation proposed in [20].