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. 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 e ^0.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.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 x ₀ = 0 and x _i+1 = x_i + h_i for i = 0, · · · , N_x − 1, where N_x is the number of grid intervals and h_i is the grid spacing, see Fig. 1 . And we assume that x_Nx = S _max , x _Nx+1 = S _max + h _Nx−1 , and x _Nx+2 = eS _max . 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 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 = S _max [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. To accelerate the computation, we save the linear interpolation factors at each point x_i before the main time-step iterations. The factor matrix A is defined as where m _i,j = ( x_ie^ξj − x_ij ) / h_ij , x_ij ≤ x_ie^ξj < x_ij+1 , for some 0 ≤ i_j ≤ N_x + 1. Then, we have which is schematically illustrated in Fig. 2(b) . 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 u _exact ( 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. 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. 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: S _max = 300, and T = 1. The numerical solutions are computed with the uniform grids h = 2 ²⁻ⁿ and the uniform time-steps ∆ τ = 250 h ² for n = 1, 2, and 3. We denote the error as e = u − u _exact , where u and u _exact are the numerical and analytic solutions, respectively. And we calculate its discrete l ₂ -norm error as Table 1 shows the l ₂ -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. 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, S _max = 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. 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 S _max = 600, N_x = 4800 and N_τ = 1440 on uniform grid. The following Table 2 represents the l ₂ -norm error and relative CPU time on uniform and non-uniform grids, respectively Here, l ₂ -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. 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.