A FAST AND ROBUST NUMERICAL METHOD FOR OPTION PRICES AND GREEKS IN A JUMP-DIFFUSION MODEL

Abstract

Abstract. We propose a fast and robust finite difference method for Merton's jump diffusion model, which is a partial integro-differential equation. To speed up a computational time, we compute a matrix so that we can calculate the non-local integral term fast by a simple matrix-vector operation. Also, we use non-uniform grids to increase efficiency. We present numerical experiments such as evaluation of the option prices and Greeks to demonstrate a performance of the proposed numerical method. The computational results are in good agreements with the exact solutions of the jump-diffusion model.

1. INTRODUCTION

The derivative securities pricing abilities of the jump-diffusion option pricing model generally provide more efficient and robust valuation [7, 11]. The Black–Scholes (BS) model, proposed by Fischer Black and Myron Scholes in 1973 [1], is a mathematical method for option price and is given by

where x is the underlying asset, t is the time, u(x, t) is the value of an option at t on x, σ is the volatility, and r is the risk neutral interest rate. According to the Merton’s jump-diffusion model based on a Poisson process [10], we are able to describe as the following equation: dx/x = (µ − λk)dt + σdW + (η − 1)dU, where dW is the standard Brownian motion and dU is a Poisson process of intensity λ, i.e., dU is 0 or 1 with probability 1 − λdt or λdt, respectively. Here η is taken to be a lognormally distributed jump amplitude with probability density

where k is the expectation 𝔼(η−1) given by e0.5δ2. Since the Brownian motion and Poisson process are supposed to be uncorrelated, the asset price jumps from x to xeξ if a jump occurs. Then, we are able to get a partial integro-differential equation given by

where The authors in [2] proposed implicit-explicit Runge–Kutta methods for the time integration to solve the integral term explicitly, giving higher-order accuracy schemes under weak stability time-step restrictions. Kwon and Lee developed a second-order convergent finite difference method to solve partial integro-differential equations which describe the behavior of option prices under jump-diffusion models [8]. This paper is organized as follows. In Section 2, we describe the proposed numerical scheme in detail. In Section 3, we present numerical experiments to validate accuracy and efficiency of our proposed algorithm. In Section 4, conclusions are drawn.

 

2. NUMERICAL SOLUTION

2.1. Discretization with finite differences We will use a finite difference method We will use a finite difference method [5] to numerically solve Eq. (1.1). By changing variable τ = T − t, where T is maturity, we can rewrite Eq. (1.1) as follows:

The jump-diffusion equation is discretized on a grid defined by x0 = 0 and xi+1 = xi + hi for i = 0, · · · , Nx − 1, where Nx is the number of grid intervals and hi is the grid spacing, see Fig. 1. And we assume that xNx = Smax, xNx+1 = Smax + hNx−1, and xNx+2 = eSmax. Let us denote the numerical approximation of the solution by , where ∆τ = T /Nτ is the time-step size and Nτ is the total number of time-steps. By applying an implicit scheme to Eq. (2.1), we have

Figure 1.A non-uniform grid.

where the first and second derivatives are defined as

and is an approximate numerical quadrature of . Using the composite Simpson’s rule, we can rewrite Eq. (2.2) as

where

We impose the zero Dirichlet boundary condition at x = 0 and the linear boundary condition at x = Smax [13]. The matrix form of the linear system (2.3) can be rewritten as

We solve the linear system by using the Thomas algorithm, which inverts a tridiagonal matrix directly. Next, we consider the numerical quadrature of the following integral term

To evaluate the integral term (2.4), we truncate the infinite domain (−∞, ∞) to [a, b]. Let us consider [a, b] = [−1, 1] case to be concrete. Therefore we use the composite Simpson’s rule for the numerical integration.

where M is a positive even integer. Figure 2(a) shows a schematic illustration of the option price u(x, τ ) and G(xeξ ) function.

Figure 2.(a) Schematic illustration of the option price u(x, τ ) and G(xieξ) function. (b) Linear interpolation of u(x, τ ) at x = xieξj

To accelerate the computation, we save the linear interpolation factors at each point xi before the main time-step iterations. The factor matrix A is defined as

where mi,j = (xieξj − xij) / hij, xij ≤ xieξj < xij+1, for some 0 ≤ ij ≤ Nx + 1. Then, we have which is schematically illustrated in Fig. 2(b).

 

3. NUMERICAL EXPERIMENT

In this section, we perform various numerical tests to show efficiency and robustness of the proposed numerical scheme under the jump-diffusion model. All computations were run in MATLAB version 7.7 [12].

3.1. European call option We consider a European call option whose payoff is given as u(x, 0) = max(x − K, 0), where K is strike price. We will use K = 100, r = 0.03, σ = 0.3, λ = 0.1, δ = 0.1, and uniform grid sizes h unless otherwise specified. Closed-form solution uexact(x, t) of jump-diffusion model for European call option is formulated as an infinite series of Black, Scholes, and Merton’s model [1, 10]:

where

We sum Eq. (3.1) up to n = 200. Figure 3 shows the payoff function and two numerical solutions at τ = 1 with the BS and the BS with jump-diffusion models. Here, we take h = 1 and ∆τ = 1/360. The result shows that the option values with the jump-diffusion model are larger than those with the BS model without the jump-diffusion.

Figure 3.Payoff function and two numerical solutions at τ = 1 with the BS and the BS with jump-diffusion models.

3.2. Effect of interval size on the numerical quadrature To find a suitable interval [−b, b] for numerical integration, we test the effect of b on the results with M = 50. Figure 4 shows the numerical solutions at τ = 1 with various values of b, such as 0.2, 0.3, 0.4, 0.8, 1, and 2. The result suggests that b = 1 is accurate enough. Therefore, we will use [−1, 1] in the calculation of the integral term (2.5) from now on.

Figure 4.Numerical solutions at τ = 1 with various values of b, such as, 0.2, 0.3, 0.4, 0.8, 1, and 2.

3.3. Convergence test We perform a convergence test to verify the accuracy of the proposed numerical method in this paper. In this test, we use the following parameters: Smax = 300, and T = 1. The numerical solutions are computed with the uniform grids h = 22−n and the uniform time-steps ∆τ = 250h2 for n = 1, 2, and 3. We denote the error as e = u − uexact, where u and uexact are the numerical and analytic solutions, respectively. And we calculate its discrete l2-norm error as Table 1 shows the l2-norm errors and convergence rates for numerical solutions with respect to h. This result indicates that the scheme is first-order accurate in time and second-order accurate in space.

Table 1.l2-norm errors and convergence rates for numerical solutions on uniform grid with respect to h.

3.4. Computation of the Greeks We compute the Greeks by using the proposed scheme and compare them with the closed-form solutions. Delta (∆) is the first derivative with respect to the underlying asset x.

Gamma (Γ) is the second derivative with respect to the underlying asset x.

Vega (V) is the derivative with respect to volatility σ.

Rho (ρ) is the derivative with respect to the risk-free interest rate r.

Theta (Θ) is the derivative with respect to the time t.

The parameters used in the calculations of numerical and closed-form solutions are as follows: ∆σ = 0.1, ∆r = 0.01, ∆t = 1/360, Smax = 300, and h = 1 on uniform grid. Figure 5(a), (b), (c), (d), and (e) are the results of ∆, Γ, V, ρ, and Θ from the numerical and closed-form solutions, respectively.

Figure 5.Comparison between the numerical and closed-form solutions for the Greeks. (a), (b), (c), (d), and (e) are the results of ∆, Γ, V, ρ, and Θ, respectively.

3.5. Efficiency of non-uniform grid In this section, we show the efficiency of non-uniform grid on jump-diffusion model. We first evaluate the numerical solution with Smax = 600, Nx = 4800 and Nτ = 1440 on uniform grid. The following Table 2 represents the l2-norm error and relative CPU time on uniform and non-uniform grids, respectively

Table 2.Comparison with l2-norm error on the interval [80, 120] and relative CPU time for numerical solutions on uniform and non-uniform grid.

Here, l2-norm error is evaluated on the interval [80, 120], which is indicated a neighborhood of the exercise price 100. And we have numerical solution on non-uniform grid, which is defined as x = [0 1.11 : 1.11 : 38.89 40 : 0.25 : 160 164.89 : 4.89 : 595.11 600].

In Table 2, 1 in CPU time of uniform grid stands for the relative calculation time for uniform grid. And the CPU time of non-uniform grid is scaled with that of the uniform grid. Through this test, we can confirm that non-uniform grid is more efficient than the uniform grid.

 

4. CONCLUSIONS

In this paper, we presented an efficient, fast, and robust finite difference method for the valuation of European call options with jump-diffusion processes. The mathematical model for the process is a partial integro-differential equation (1.1), Merton’s jump diffusion model. To speed up the computational time, we computed the factor matrix A in (2.6) so that we can calculate the non-local integral term fast by a simple matrix-vector operation. Also, we used non-uniform grids to increase efficiency. We presented numerical experiments such as evaluation of the option prices and the Greeks to demonstrate the performance of the proposed numerical method. The computational results were in good agreements with the exact solutions of the jump-diffusion model.