Advanced
A FAST AND ROBUST NUMERICAL METHOD FOR OPTION PRICES AND GREEKS IN A JUMP-DIFFUSION MODEL
A FAST AND ROBUST NUMERICAL METHOD FOR OPTION PRICES AND GREEKS IN A JUMP-DIFFUSION MODEL
The Pure and Applied Mathematics. 2015. May, 22(2): 159-168
Copyright © 2015, Korean Society of Mathematical Education
  • Received : December 05, 2014
  • Accepted : February 06, 2015
  • Published : May 31, 2015
Download
PDF
e-PUB
PubReader
PPT
Export by style
Share
Article
Author
Metrics
Cited by
TagCloud
About the Authors
DARAE JEONG
DEPARTMENT OF MATHEMATICS, KOREA UNIVERSITY, SEOUL 136-713, REPUBLIC OF KOREAEmail address:tinayoyo@korea.ac.kr
YOUNG ROCK KIM
MAJOR IN MATHEMATICS EDUCATION, HANKUK UNIVERSITY OF FOREIGN STUDIES, SEOUL, 130- 791, REPUBLIC OF KOREAEmail address:rocky777@hufs.ac.kr
SEUNGGYU LEE
DEPARTMENT OF MATHEMATICS, KOREA UNIVERSITY, SEOUL 136-713, REPUBLIC OF KOREAEmail address:sky509@korea.ac.kr
YONGHO CHOI
DEPARTMENT OF MATHEMATICS, KOREA UNIVERSITY, SEOUL 136-713, REPUBLIC OF KOREAEmail address:poohyongho@korea.ac.kr
WOONG-KI LEE
BUSINESS SCHOOL, KOREA UNIVERSITY, SEOUL 136-071, REPUBLIC OF KOREAEmail address:unggii@korea.ac.kr
JAE-MAN SHIN
DEPARTMENT OF FINANCIAL ENGINEERING, KOREA UNIVERSITY, SEOUL 136-071, REPUBLIC OF KOREAEmail address:sim00@korea.ac.kr
HYO-RIM AN
DEPARTMENT OF FINANCIAL ENGINEERING, KOREA UNIVERSITY, SEOUL 136-071, REPUBLIC OF KOREAEmail address:ahr78@korea.ac.kr
HYEONGSEOK HWANG
DEPARTMENT OF FINANCIAL ENGINEERING, KOREA UNIVERSITY, SEOUL 136-071, REPUBLIC OF KOREAEmail address:hhs288@korea.ac.kr
HJUNSEOK KIM
DEPARTMENT OF MATHEMATICS, KOREA UNIVERSITY, SEOUL 136-713, REPUBLIC OF KOREAEmail address:cfdkim@korea.ac.kr

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.
Keywords
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
PPT Slide
Lager Image
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
PPT Slide
Lager Image
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
PPT Slide
Lager Image
where
PPT Slide
Lager Image
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:
PPT Slide
Lager Image
The jump-diffusion equation is discretized on a grid defined by x 0 = 0 and x i+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 = 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
PPT Slide
Lager Image
, 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
PPT Slide
Lager Image
A non-uniform grid.
PPT Slide
Lager Image
where the first and second derivatives are defined as
PPT Slide
Lager Image
and
PPT Slide
Lager Image
is an approximate numerical quadrature of
PPT Slide
Lager Image
. Using the composite Simpson’s rule, we can rewrite Eq. (2.2) as
PPT Slide
Lager Image
where
PPT Slide
Lager Image
PPT Slide
Lager Image
PPT Slide
Lager Image
PPT Slide
Lager Image
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
PPT Slide
Lager Image
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
PPT Slide
Lager Image
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.
PPT Slide
Lager Image
where M is a positive even integer. Figure 2(a) shows a schematic illustration of the option price u ( x , τ ) and G ( xeξ ) function.
PPT Slide
Lager Image
(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
PPT Slide
Lager Image
where m i,j = ( xieξj xij ) / hij , xij xieξj < xij+1 , for some 0 ≤ ij Nx + 1. Then, we have
PPT Slide
Lager Image
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 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] :
PPT Slide
Lager Image
where
PPT Slide
Lager Image
PPT Slide
Lager Image
PPT Slide
Lager Image
PPT Slide
Lager Image
PPT Slide
Lager Image
PPT Slide
Lager Image
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.
PPT Slide
Lager Image
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.
PPT Slide
Lager Image
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: S max = 300, and T = 1. The numerical solutions are computed with the uniform grids h = 2 2−n and the uniform time-steps ∆ τ = 250 h 2 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 2 -norm error as
PPT Slide
Lager Image
Table 1 shows the l 2 -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.
l2-norm errors and convergence rates for numerical solutions on uniform grid with respect toh.
PPT Slide
Lager Image
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 .
PPT Slide
Lager Image
Gamma (Γ) is the second derivative with respect to the underlying asset x .
PPT Slide
Lager Image
Vega ( V ) is the derivative with respect to volatility σ .
PPT Slide
Lager Image
Rho ( ρ ) is the derivative with respect to the risk-free interest rate r .
PPT Slide
Lager Image
Theta (Θ) is the derivative with respect to the time t .
PPT Slide
Lager Image
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.
PPT Slide
Lager Image
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 S max = 600, Nx = 4800 and Nτ = 1440 on uniform grid. The following Table 2 represents the l 2 -norm error and relative CPU time on uniform and non-uniform grids, respectively
Comparison withl2-norm error on the interval [80, 120] and relative CPU time for numerical solutions on uniform and non-uniform grid.
PPT Slide
Lager Image
Comparison with l2-norm error on the interval [80, 120] and relative CPU time for numerical solutions on uniform and non-uniform grid.
Here, l 2 -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.
Acknowledgements
The corresponding author (J.S. Kim) was supported by The Small and Medium Business Administration.
View Fulltext  
Black F. , Scholes M. 1973 The pricing of options and corporate liabilities J. Polit. Econ. 81 637 - 654    DOI : 10.1086/260062
Briani M. , Natalini R. , Russo G. 2007 Implicit-explicit numerical schemes for jump-diffusion processes Calcolo 44 (1) 33 - 57    DOI : 10.1007/s10092-007-0128-x
Carr P. , Cousot L. 2011 A PDE approach to jump-diffusions Quant. Fin. 11 (1) 33 - 52    DOI : 10.1080/14697688.2010.531042
Cheang G. , Chiarella C. 2011 A modern view on Merton’s jump-diffusion model Quant.Fin. Res. Cent. 287
Duffy D.J. 2006 Finite Difference methods in financial engineering: a Partial Differential Equation approach John Wiley & Sons
Feng L. , Linetsky V. 2008 Pricing options in jump-diffusion models: an extrapolation approach Oper. Res. 56 (2) 304 - 325    DOI : 10.1287/opre.1070.0419
Kremer J.W. , Roenfeldt R.L. 1993 Warrant pricing: jump-diffusion vs. Black-Scholes J. Finan. Quant. Anal. 28 (2) 255 - 272    DOI : 10.2307/2331289
Kwon Y.H. , Lee Y. 2011 A second-order finite difference method for option pricing under jump-diffusion models SIAM J. Numer. Anal. 49 (6) 2598 - 2617    DOI : 10.1137/090777529
Lee S. , Li Y. , Choi Y. , Hwang H. , Kim J. 2014 Accurate and efficient computations for the Greeks of European multi-asset options J. KSIAM 18 (1) 61 - 74
Merton R.C. 1976 Option pricing when underlying stock returns are discontinuous J. Polit. Econ. 3 (1) 125 - 144
Tankov P. 2003 Financial modelling with jump processes CRC press
1998 Using MATLAB The MathWorks Inc. Natick, MA http://www.mathworks.com
Windcliff H. , Forsyth P.A. , Vetzal K.R. 2004 Analysis of the stability of the linear boundary condition for the Black-Scholes equation J. Comput. Fin. 8 65 - 92