An Implicit-Explicit Computational Method Based on Time Semi-Discretization for Pricing Financial Derivatives with Jumps

Yang Wang

School of Economic, Jinan University, Guangzhou, China.

DOI: 10.4236/ojs.2018.82022                

Abstract

This paper considers pricing European options under the well-known of SVJ model of Bates and related computational methods. According to the no-arbitrage principle, we first derive a partial differential equation that the value of any European contingent claim should satisfy, where the asset price obeys the SVJ model. This equation is numerically solved by using the implicit- explicit backward difference method and time semi-discretization. In order to explain the validity of our method, the stability of time semi-discretization scheme is also proved. Finally, we use a simulation example to illustrate the efficiency of the method.

Keywords

SVJ Model of Bates, Time Semi-Discretization, Stability, No-Arbitrage Principle, Implicit-Explicit Backward Difference Method

1. Introduction

It is well known in the standard Black-Scholes [1] model that the log-normal stock diffusion with constant volatility is not consistent with the market price movement. More importantly, there are evidences to indicate that the B-S model cannot describe real stock price behavior. Sometimes sudden changes of prices could happen at a random time and these changes cannot be captured by the log-normal distribution characteristic of the stock price in the Black Scholes model. In order to overcome these shortcomings of the Black-Scholes model, a variety of alternative models are proposed in the financial literature. Among these, the jump diffusion models proposed by Merton [2] and Kou [3] are widely used models. Furthermore, the SVJ model of Bates can not only make up for the shortcomings of the B-S model but also can describe the financial market in a more suitable way.

In order to reflect the effect of stochastic volatility on the market volatility, we have two ways, one assumes that the volatility is determined, and the volatility function is then determined by calibration to market price. The other assume that the stochastic volatility approach [4] [5] . In this model the volatility of the stock price is considered to be a mean reverting diffusion process, which is usually related to the stock process itself. A general approach, originally suggested by Bates, can be much more preferable to describe the financial market.

The valuation of option under jump diffusion process satisfies a partial integro-differential equation with boundary conditions. There are several numerical methods available to approximately solve the above equation. For example, in [6] , Almendral and Osterlee presented an implicit second order accurate time discretization with finite difference. Recently, Patidar [7] developed an efficient method for pricing Merton jump diffusion option. The scheme proposed by Halluin [8] required to use an iterative procedure to solve discrete equations. In particular, an approach based on implicit-explicit schemes in which integral term is treated explicitly was proposed by YongHoon Kwon [9] . This article aims to solve the pricing financial derivatives with jumps which based on time semi-discretization.

The paper is organized as follows. In Section 2, mathematical models for pricing option with jump diffusion process are given in terms of partial integro-differential equations and provide a brief review of both the Merton and Kou jump diffusion models, then the SVJ model of Bates was presented. Section 3 deals with the construction of time level implicit explicit scheme to discretize the jump diffusion model. The discretization method is called IMEX-BDF2 method with three time levels. Section 4 we consider the penalized nonlinear equation to approximate the LCP system. In Section 5, we will give a numerical simulation example. Finally the paper ends with some discussion and conclusive remarks in Section 6.

2. The Mathematical Model

In this part, we briefly discuss the mathematical model for pricing option with jump diffusion process. Consider an asset which the asset price is S, and then the movement of stock price is modeled by the following stochastic differential equation.

\ $\frac{\text{d}{S}_t}{S_{t^{-}}}=\left(\mu -\lambda \kappa \right)\text{d}t+\sigma \text{d}{W}_t+{\displaystyle {\int }_0^{\infty }\left(\eta -1\right)N\left(\text{d}\eta ,\text{d}t\right)}$ (2.1)

where μ is drift rate, t as the time to maturity, σ represents the constant volatility, $\text{d}{W}_t$ is generally satisfied with the Gauss process. The constant λ is the intensity of the independent Poisson process, κ used to express the relative size of the expected jump $\left(\eta -1\right)$ .

The price of risk asset can be described by the Brown movement. However, the Brown movement is subject to normal distribution, and its value can be negative, which does not conform to the nature of the price. According to the geometric Brown motion obeying the lognormal distribution, which avoids the defect that Brown motion may take negative value. The geometric Brown motion has been widely applied when we establish the financial asset price model.

The distribution of the risk asset which obeys the lognormal distribution is symmetric. However, the empirical analysis shows that the price is asymmetric. The log-double-exponential distribution can make up for the shortage of the lognormal distribution. It can simulate the financial market more vividly.

If we define that $g\left(\eta \right)$ is probability density function of the jump with amplitude η, under Merton’s model $g\left(\eta \right)$ is given by the log-normal density, and under Kou’s jump-diffusion model $g\left(\eta \right)$ is the following log-double-exponential density. The function $g\left(\eta \right)$ can be written as:

$g\left(\eta \right):=\{\begin{array}{l}\frac{1}{\sqrt{2\text{π}}{\sigma }_j\eta }{\text{e}}^{-\frac{{\left(\ln \eta -{\mu }_j\right)}^2}{2{\sigma }_j{}^2}}\\ \frac{1}{\eta }\left(p{\eta }_1{\text{e}}^{-{\eta }_1\ln \eta }\mathcal{H}\left(\ln \eta \right)+q{\eta }_2{\text{e}}^{{\eta }_2\ln \eta }\mathcal{H}\left(-\ln \eta \right)\right)\end{array}$ (2.2)

In Merton’s model,

$\kappa =E\left(\eta -1\right)={\text{e}}^{\left({\mu }_J+\frac{{\sigma }_J^2}{2}\right)}-1$ (2.3)

In Kou’s model, where ${\eta }_1>1,{\eta }_2>0,p>0,q=1-p$ , and $\mathcal{H}(\cdot )$ is the Heaviside function, we can show that

$\kappa =\frac{p{\eta }_1}{{\eta }_1-1}+\frac{q{\eta }_2}{{\eta }_1+1}-1$ (2.4)

although the model is well, there are still some shortcomings, there is no consideration that the volatility is also random.

In order to fit the financial market well, it should be noted that the volatility rate σ is a random fluctuation. Then the model of Merton is popularized to get the SVJ (stochastic volatility-jump) model of Bates.

$\frac{\text{d}{S}_t}{S_{t^{-}}}=\left(\mu -\lambda \kappa \right)\text{d}t+\sqrt{V_t}\text{d}{W}_t+{\displaystyle {\int }_0^{\infty }\left(\eta -1\right)N\left(\text{d}\eta ,\text{d}t\right)}$

$D{V}_t=\kappa \left(\alpha -V_t\right)\text{d}t+\sigma \sqrt{V_t}\text{d}{W}_t^v$ (2.5)

$\text{d}{W}_t\text{d}{W}_t^v=\rho \text{d}t,\ln \left(\eta \right)~N\left({\mu }_j,{\sigma }_j^2\right)$

where the $\mu ,\lambda ,\kappa ,\eta$ are same as the above, $V_t$ represents the instantaneous variance or volatility of the rate of return on assets, which is described by a square root process. It can guarantee the nonnegativity of $V_t$ . ρ is the correlation coefficient of asset price shock $W_t$ and volatility impact on the $W_t^v$ , $\rho <0$ reflects the leverage effect of asset price. According to Ito lemma, the SDE equation of the logarithmic asset price satisfaction is:

$\text{d}\ln S_t=\left(\gamma -\frac{1}{2}{V}_t\right)\text{d}t+\sqrt{V_t}\text{d}{W}_t+Y_t\text{d}{N}_t$ (2.6)

$\gamma =\mu -\lambda E\left[\eta -1\right]$ (2.7)

Let $\phi \left(S,t\right)$ represent the value of a contingent claim that depends on the underlying asset price S with current time t, then $\phi \left(S,t\right)$ satisfy following backward partial integration differential equation.

$\begin{array}{l}\frac{\partial \phi }{\partial t}+\frac{\partial \phi }{\partial \ln S_t}\left(\gamma -\frac{1}{2}{V}_t\right)+\frac{\partial \phi }{\partial V_t}\kappa \left(\alpha -V_t\right)+\frac{1}{2}\frac{{\partial }^2\phi }{\partial {\left(\ln S_t\right)}^2}{V}_t+\frac{1}{2}\frac{{\partial }^2\phi }{\partial V_t{}^2}{\sigma }^2{V}_t\\ +\frac{{\partial }^2\phi }{\partial \ln S_t\partial V_t}\rho \sigma V_t+\lambda E_t\left[\phi \left(\ln S_t+Y\right)-\phi \left(\ln S_t\right)\right]=0\end{array}$ (2.8)

where $\mathcal{X}=\ln \frac{S}{K},\phi \left(\mathcal{X},t\right)=V\left(K{\text{e}}^{\mathcal{X}},T-t\right),V\left(S,T\right)=\mathcal{G}\left(S\right)$ and the price of a European put option which the initial condition and asymptotic behavior are described by

$\phi \left(\mathcal{X},0\right)=\max \left(K-K{\text{e}}^{\mathcal{X}},0\right)$

$\phi \left(\mathcal{X},t\right)=\{\begin{array}{l}K{\text{e}}^{-rt}-K{\text{e}}^{\mathcal{X}}\mathcal{X}\to -\infty \\ 0\mathcal{X}\to \infty \end{array}$ (2.9)

3. Discretization

In this section, we will discuss the stability of some implicit-explicit (IMEX) time semi-discrete methods. For non-linear partial differential equations, a class of such IMEX methods have already been discussed in [10] . But now it is entirely different in the sense of analysis and applications. We shall construct an implicit-explicit backward difference method of order two time semi-discretization, which can fit the original equation much more properly.

Let $\left\{0=t_0<t_1<\cdots <t_N=T;t_n-t_{n-1}=k,n=1,2,\cdots ,N\right\}$ be a partition of the interval [0, T]. We will consider the following three time semi-discretization.

IMEX-BDF2 (Implicit-explicit backward difference method of order two)

$\frac{\frac{3}{2}{u}^{n+1}-2u^n+\frac{1}{2}{u}^{n-1}}{k}=L_{n+1}{u}^{n+1}+E\left({\left(\eta -1\right)}^n{u}^n+f^n\right);1\le n_0\le n\le N-1$ (3.1)

where the $u\left(\mathcal{X},t\right)=\phi \left(\mathcal{X},T-t\right)$ , $E{f}^n=2f^n-f^{n-1}$ , $f^n\left(x\right)=f\left(x,t_n\right)$ and $E{u}^n=2u^n-u^{n-1}$ , In order to prove the stability of the time semi-discretization, we will need the following lemmas.

Lemma 1: (Discrete Gronwall’s inequality) Let $\left\{a_n\right\},\left\{b_n\right\},\left\{{\alpha }_n\right\}$ and $\left\{{\beta }_n\right\}$ be four non-negative sequences, such that for non-negative integers $n_0$ & N

$a_n+b_n\le {\alpha }_n+{\displaystyle \underset{i=n_0}{\overset{n-1}{\sum }}{\beta }_i{a}_i},n=n_0,n_0+1,\cdots ,N,$ (3.2)

then, $a_n+b_n\le \left({\max }_{n_0\le l\le n}{\alpha }_l\right)\exp \left({\displaystyle {\sum }_{i=n_0}^{n-1}{\beta }_i}\right),\forall n:n_0\le n\le N$ , where the summation the ${{\displaystyle \sum }}_{i=n_0}^{n-1}$ is assumed to be zero if $n_0\ge n$ . This lemma can be easily proved by mathematical induction.

The Discrete Gronwall’s inequality illustrates that the function which satisfies the integro-differential equation, and then there is a corresponding inequality.

Lemma 2: For IMEX-BDF2, there exists an $N_0\in \mathcal{N}$ such that $\forall N\ge N_0$ & $n_0+1\le m\le N$ , we have

$\begin{array}{l}{\Vert u^m\Vert }_1^2+k\underset{n=n_0+1}{\overset{m}{{\displaystyle \sum }}}{\Vert \delta u^n\Vert }^2\\ \le C\left(k{\Vert u^{n_0-1}\Vert }_1^2+{\Vert u^{n_0}\Vert }_1^2+k{\Vert \delta u^{n_0}\Vert }^2+k\underset{n=n_0}{\overset{m-1}{{\displaystyle \sum }}}{\Vert E{f}^n\Vert }^2\right),1\le n_0\le m-1,\end{array}$ (3.3)

and

$\begin{array}{l}{\Vert u^m\Vert }_2^2+k\underset{n=n_0+1}{\overset{m}{{\displaystyle \sum }}}{\Vert \delta u^n\Vert }_1^2\\ \le C(k{\Vert u^{n_0-2}\Vert }^2+k{\Vert u^{n_0-1}\Vert }^2+{\Vert u^{n_0}\Vert }_1^2+{\Vert \delta u^{n_0-1}\Vert }^2+{\Vert \delta u^{n_0}\Vert }^2\begin{array}{c}\\ \end{array}\\ +k\underset{n=n_0}{\overset{m-1}{{\displaystyle \sum }}}\left({\Vert E{f}^n\Vert }^2+{\Vert E\delta f^n\Vert }_{-1}^2\right)),2\le n_0\le m-1.\end{array}$ (3.4)

Proof: From relation IMEX-BDF2 we see that, $\forall \upsilon \in H_0^1$ ,

$\left(3\delta u^{n+1}-\delta u^n,\upsilon \right)-2\left(L_{n+1}{u}^{n+1},\upsilon \right)=2\left(E\left(j^n{u}^n\right),\nu \right)+2\left(E{f}^n,\upsilon \right)$ (3.5)

which implies that

$\begin{array}{l}3\left(\delta u^{n+1},\upsilon \right)+2A_{n+1}\left(u^{n+1},\upsilon \right)+2B_{n+1}\left(u^{n+1},\upsilon \right)\\ =2\left(E\left(j^n{u}^n\right),\upsilon \right)+2\left(E{f}^n,\upsilon \right)+\left(\delta u^n,\upsilon \right).\end{array}$ (3.6)

By taking $\upsilon =u^{n+1}-u^n$ and using symmetry of $A_n\left(u,\upsilon \right)$ , we have

$\begin{array}{l}2k{\Vert \delta u^{n+1}\Vert }^2+A_{n+1}\left(u^{n+1},u^{n+1}\right)-A_{n+1}\left(u^n,u^n\right)+A_{n+1}\left(u^{n+1}-u^n,u^{n+1}-u^n\right)\\ =-2k{B}_{n+1}\left(u^{n+1},\delta u^{n+1}\right)+2k\left(E\left(j^n{u}^n\right),\delta u^{n+1}\right)+2k\left(E{f}^n,\delta u^{n+1}\right)+k\left(\delta u^n,\delta u^{n+1}\right).\end{array}$ (3.7)

We can get

$\begin{array}{l}3k{\Vert \delta u^{n+1}\Vert }^2+A_{n+1}\left(u^{n+1},u^{n+1}\right)\\ \le A_n\left(u^n,u^n\right)+\left(A_{n+1}\left(u^n,u^n\right)-A_n\left(u^n,u^n\right)\right)-2k{B}_{n+1}\left(u^{n+1},\delta u^{n+1}\right)\\ +2k\left(E\left(j^n{u}^n\right),\delta u^{n+1}\right)+2k\left(E{f}^n,\delta u^{n+1}\right)+k\left(\delta u^n,\delta u^{n+1}\right)\end{array}$ (3.8)

Now by using Lemma1, we have

$\begin{array}{l}2k{\Vert \delta u^{n+1}\Vert }^2+2A_{n+1}\left(u^{n+1},u^{n+1}\right)\\ \le 2A_n\left(u^n,u^n\right)+k{\Vert \delta u^n\Vert }^2+k{C}_1\left({\Vert u^{n-1}\Vert }_1^2+{\Vert u^n\Vert }_1^2+{\Vert u^{n+1}\Vert }_1^2\right)+k{C}_2{\Vert E{f}^n\Vert }^2\end{array}$ (3.9)

for some generic constants $C_1,C_2>0$ .

After summing from $n=n_0$ to $m-1$ , $2\le n_0+1\le m\le N$ , we get

$\begin{array}{l}\underset{n=n_0+1}{\overset{m}{{\displaystyle \sum }}}k{\Vert \delta u^n\Vert }^2+2M_3{\Vert u^m\Vert }_1^2\\ \le k{C}_1{\Vert u^{n_0-1}\Vert }_1^2+\left(2M+k{C}_1\right){\Vert u^{n_0}\Vert }_1^2+k{\Vert \delta u^{n_0}\Vert }^2\\ +k{C}_1\underset{n=n_0+1}{\overset{m}{{\displaystyle \sum }}}{\Vert u^n\Vert }_1^2+k{C}_2\underset{n=n_0}{\overset{m-1}{{\displaystyle \sum }}}{\Vert E{f}^n\Vert }^2\end{array}$ (3.10)

Hence there exists an $N_0\in N$ such that $\forall N\ge N_0$ above relation with the use of Discrete Gronwall’s inequality gives

${\Vert u^m\Vert }_1^2+k\underset{n=n_0+1}{\overset{m}{{\displaystyle \sum }}}{\Vert \delta u^n\Vert }^2\le C\left(k{\Vert u^{n_0-1}\Vert }_1^2+{\Vert u^{n_0}\Vert }_1^2+k{\Vert \delta u^{n_0}\Vert }^2+k\underset{n=n_0}{\overset{m-1}{{\displaystyle \sum }}}{\Vert E{f}^n\Vert }^2\right)$ (3.11)

Further, we have

$\begin{array}{l}\left(3\delta u^{n+1}-4\delta u^n+\delta u^{n-1},\upsilon \right)-2\left(L_{n+1}{u}^{n+1}-L_n{u}^n,\upsilon \right)\\ =2\left(E\left(j^n{u}^n\right)-E\left(j^{n-1}{u}^{n-1}\right),\upsilon \right)+2\left(E{f}^n-E{f}^{n-1},\upsilon \right)\end{array}$

By taking $\upsilon =\delta u^{n+1}$ and using the relation:

$\begin{array}{l}2\left(3u-4\upsilon +\omega ,u\right)\\ ={\Vert u\Vert }^2-{\Vert \upsilon \Vert }^2+{\Vert 2u-\upsilon \Vert }^2-{\Vert 2\upsilon -\omega \Vert }^2+{\Vert u-2\upsilon +\omega \Vert }^2\end{array}$

we get

$\begin{array}{l}{\Vert \delta u^{n+1}\Vert }^2-{\Vert \delta u^n\Vert }^2+{\Vert 2\delta u^{n+1}-\delta u^n\Vert }^2-{\Vert 2\delta u^n-\delta u^{n-1}\Vert }^2-4k\left(L_{n+1}\delta u^{n+1},\delta u^{n+1}\right)\\ \le 4\left(L_{n+1}{u}^n-L_n{u}^n,\delta u^{n+1}\right)+4k\left(E\delta f^n,\delta u^{n+1}\right)+4\left(E\left(j^n{u}^n\right)-E\left(j^{n-1}{u}^{n-1}\right),\delta u^{n+1}\right)\end{array}$ (3.12)

Again by using Lemma 1, we can get

$\begin{array}{l}{\Vert \delta u^{n+1}\Vert }^2-{\Vert \delta u^n\Vert }^2+{\Vert 2\delta u^{n+1}-\delta u^n\Vert }^2-{\Vert 2\delta u^n-\delta u^{n-1}\Vert }^2+4k{M}_2{\Vert \delta u^{n+1}\Vert }_1^2\\ \le 3k{M}_2{\Vert \delta u^{n+1}\Vert }_1^2+k{C}_3\left({\Vert \delta u^{n-1}\Vert }^2+{\Vert \delta u^n\Vert }^2+{\Vert \delta u^{n+1}\Vert }^2\right)\\ +k{C}_4\left({\Vert u^{n-2}\Vert }^2+{\Vert u^{n-1}\Vert }^2+{\Vert u^n\Vert }_1^2\right)+k{C}_5{\Vert E\delta f^n\Vert }_1^2\end{array}$ (3.13)

for some generic constants $C_3,C_4,C_5>0$ .

After summing from $n=n_0$ to $m-1$ , $3\le n_0+1\le m\le N$ ,

we get

$\begin{array}{l}{\Vert \delta u^m\Vert }^2+k{M}_2\underset{n=n_0+1}{\overset{m}{{\displaystyle \sum }}}{\Vert \delta u^n\Vert }_1^2\\ \le {\Vert \delta u^{n_0}\Vert }^2+{\Vert 2\delta u^{n_0}-\delta u^{n_0-1}\Vert }^2+k{C}_3\left({\Vert \delta u^{n_0-1}\Vert }^2+2\parallel {\Vert \delta u^{n_0}\Vert }^2\right)\\ +k{C}_4\left({\Vert u^{n_0-2}\Vert }^2+2{\Vert u^{n_0-1}\Vert }^2+3{\Vert u^{n_0}\Vert }_1^2\right)+3k{C}_3\underset{n=n_0+1}{\overset{m}{{\displaystyle \sum }}}{\Vert \delta u^n\Vert }^2\\ +3k{C}_4\underset{n=n_0+1}{\overset{m-1}{{\displaystyle \sum }}}{\Vert u^n\Vert }_1^2+k{C}_5\underset{n=n_0}{\overset{m-1}{{\displaystyle \sum }}}{\Vert E\delta f^n\Vert }_1^2.\end{array}$ (3.14)

Thus there exists an $N_0\in N$ such that $\forall N\ge N_0,3\le n_0+1\le m\le N$ , above relation with the use of Lemma 1 gives

$\begin{array}{l}{\Vert \delta u^m\Vert }^2+k\underset{n=n_0+1}{\overset{m}{{\displaystyle \sum }}}{\Vert \delta u^n\Vert }_1^2\\ \le C(k\left({\Vert u^{n_0-2}\Vert }^2+{\Vert u^{n_0-1}\Vert }^2\right)+{\Vert u^{n_0}\Vert }_1^2+{\Vert \delta u^{n_0-1}\Vert }^2+{\Vert \delta u^{n_0}\Vert }^2\begin{array}{c}\\ \end{array}\\ +k\underset{n=n_0}{\overset{m-1}{{\displaystyle \sum }}}\left({\Vert E{f}^n\Vert }^2+{\Vert E\delta f^n\Vert }_1^2\right)).\end{array}$ (3.15)

Finally with the help of this above relation, we have

$\begin{array}{l}{\Vert u^m\Vert }_2^2+k\underset{n=n_0+1}{\overset{m}{{\displaystyle \sum }}}{\Vert \delta u^n\Vert }_1^2\\ \le C(k\left({\Vert u^{n_0-2}\Vert }^2+{\Vert u^{n_0-1}\Vert }^2\right)+{\Vert u^{n_0}\Vert }_1^2+{\Vert \delta u^{n_0-1}\Vert }^2+{\Vert \delta u^{n_0}\Vert }^2\begin{array}{c}\\ \end{array}\\ +k\underset{n=n_0}{\overset{m-1}{{\displaystyle \sum }}}\left({\Vert E{f}^n\Vert }^2+{\Vert E\delta f\Vert }^n{}_1^2\right))\end{array}$ (3.16)

Now Lemma1, Lemma2, imply the following result: under the assumptions made in above, all the three time semi-discrete IMEX methods are stable.

4. Unconditional Stability

The types of initial and boundary conditions can be properly defined as different types of options. It is different from that the American option can be exercised at any time up to the maturity date. The types of initial and boundary conditions can not be defined easily, but it can be formulated as the linear complementarity problem (LCP) of the form.

$\{\begin{array}{l}\frac{\partial \phi }{\partial t}-\mathcal{L}\phi \ge 0\\ \phi \left(\mathcal{X},t\right)-\mathcal{G}\left(\mathcal{X}\right)\ge 0\\ \left(\frac{\partial \phi }{\partial t}-\mathcal{L}\phi \right)\left(\phi \left(\mathcal{X},t\right)-\mathcal{G}\left(\mathcal{X}\right)\right)=0\end{array}$ (4.1)

In combination with the above equation and the linear complementarity problem, we can make a transformation to simplify the operation.

$\begin{array}{l}{\mathcal{L}}_s{\phi }_s\left(x,t\right)=\frac{\partial \phi }{\partial t}+\frac{\partial \phi }{\partial \ln S_t}\left(\gamma -\frac{1}{2}{V}_t\right)+\frac{\partial \phi }{\partial V_t}\kappa \left(\alpha -V_t\right)\\ +\frac{1}{2}\frac{{\partial }^2\phi }{\partial {\left(\ln S_t\right)}^2}{V}_t+\frac{1}{2}\frac{{\partial }^2\phi }{\partial V_t{}^2}{\sigma }^2{V}_t+\frac{{\partial }^2\phi }{\partial \ln S_t\partial V_t}\rho \sigma V_t\\ +\lambda E_t\left[\phi \left(\ln S_t+Y\right)-\phi \left(\ln S_t\right)\right]\le 0\end{array}$

${\phi }_s\left(x,t\right)\ge {\phi }_s^{*}\left(x\right)\text{and}\mathcal{L}{\phi }_s\left(x,t\right)\cdot \left({\phi }_s\left(x,t\right)-{\phi }_s^{*}\left(x\right)\right)=0$

${\phi }_s\left(x_L,t\right)={\phi }_s\left(x_L\right)=\max \left(K-{\text{e}}^{x_L},0\right),t\in \left[0,T\right)$

${\phi }_s\left(x_L,t\right)={\phi }_s\left(x_L\right)=\max \left(K-{\text{e}}^{x_L},0\right),t\in \left[0,T\right)$

${\phi }_s\left(x,T\right)={\phi }_s^{*}\left(x\right)=\max \left(K-{\text{e}}^x,0\right),x\in \left(x_L,x_R\right)$ (4.2)

where K is the strike price of the options, for all $s\in S$ , and $x_L,x_R$ are the lower and upper truncated boundary respectively.

Since fractional derivative is nonlocal operator, the truncation error cannot be ignored for the nonhomogeneous Dirichlet boundary condition. The boundary condition of FPDEs can be transformed to homogeneous Dirichlet boundary condition.

Let $F_s\left(x\right)$ is defined by

$F_s\left(x\right)=\frac{{\phi }_s\left(x_R\right)-{\phi }_s\left(x_L\right)}{{\text{e}}^{x_R}-{\text{e}}^{x_L}}\left({\text{e}}^x-{\text{e}}^{x_L}\right)+{\phi }_s\left(x_L\right)$ (4.3)

using the new variable $U_s\left(x,t\right)=F_s\left(x\right)-{\phi }_s\left(x,t\right)$ , then we have

${\mathcal{L}}_s{U}_s\left(x,t\right)\ge f_s\left(x\right),$

where $f_s\left(x\right)={\mathcal{L}}_s{F}_s\left(x\right)$ , the boundary and terminal conditions will become

$f_s\left(x\right)\le {\mathcal{L}}_s{U}_s\left(x,t\right),$

$U_s\left(x,t\right)\le U_s^{*}\left(x\right),$ (4.4)

$\left(f_s\left(x\right)-\mathcal{L}{U}_s\left(x,t\right)\right)\cdot \left(U_s^{*}\left(x\right)-U_s\left(x,t\right)\right)=0.$

We have discussed that American options are described by the linear complementarity problem. Then we consider the following penalized nonlinear equation to approximate the LCP system (4.4),

$\left(f_s\left(x\right)-\mathcal{L}{U}_s^{\rho }\right)+\rho {\left[U_s^{\rho }\left(x,t\right)-U_s^{*}\left(x\right)\right]}_{+}=0,$ (4.5)

where ${\left[U_s^{\rho }\left(x,t\right)-U_s^{*}\left(x\right)\right]}_{+}=\max \left(U_s^{\rho }\left(x,t\right)-U_s^{*}\left(x\right),0\right)$ , and (4.5) converges to (4.2) as ρ goes to infinity. A semi-implicit finite difference scheme is proposed to discretize the nonlinear FPDEs system (4.5).

Now divide the intervals $\left[x_L,x_R\right]$ and $\left[0,T\right]$ into N + 1 subintervals and M subintervals, respectively. The spatial and temporal meshes are defined as follows

$x_n=x_L+nh,\text{for}n=0,1,2,\cdots ,N+1,$

$t_m=T+m\tau ,\text{for}m=0,1,2,\cdots ,M,$

where $h=\frac{x_R-x_L}{N+1},\tau =-\frac{T}{M}<0$ , $N,M$ are positive integers and ρ is chosen such that $\tau \rho =\mathcal{O}\left(1\right)$ .

Denoting $u_{n,s}^{\left(m\right)}\approx U_s^{\rho }\left(x_n,t_m\right),u_{n,s}^{*}=U_s^{*}\left(x_n\right)$ and $f_{n,s}=f_s\left(x_n\right)$ , we obtain the following finite difference scheme for the nonlinear FPDE system

$f_{n,s}-{\mathcal{L}}_s{u}_{n,s}^{\left(m+1\right)}+\rho \cdot \max \left(u_{n,s}^{\left(m+1\right)}-u_{n,s}^{*},0\right)=0$ (4.6)

Lemma 3: The nonlinear scheme above is unconditionally stable.

Proof: Let $u_s^{\left(m\right)}={\left[u_{1,s}^{\left(m\right)},u_{2,s}^{\left(m\right)},\cdots ,u_{N,s}^{\left(m\right)}\right]}^{\text{T}}$ be the solution of above in the m-th time step. Let $u_{i_0,s}^{\left(m\right)}$ be the i₀-th entry of $u_s^{\left(m\right)}$ such that

$\left|u_{i_0,s}^{\left(m\right)}\right|=\underset{1\le i\le N}{\max }\left|u_{i,s}^{\left(m\right)}\right|={\Vert u_s^{\left(m\right)}\Vert }_{\infty }$ (4.7)

thus, $\left|u_{i,s}^{\left(m\right)}\right|\le \left|u_{i_0,s}^{\left(m\right)}\right|,i=1,2,\cdots ,N$ . Then, the proof will be completed by discussing the following two cases.

Case 1: for $u_{i_0,s}^{\left(m\right)}\ge u_{i_0,s}^{*}$ , in view of above, we have

$\begin{array}{c}{\Vert u_s^{\left(m\right)}\Vert }_{\infty }=\left|u_{i_0,s}^{\left(m\right)}\right|\le \left|m_{i_0{i}_0}{u}_{i_0,s}^{\left(m\right)}\right|-\left|\underset{j\ne i_0}{{\displaystyle \sum }}{m}_{i_{0}j}{u}_{i_0,s}^{\left(m\right)}\right|\\ \le \left|m_{i_0{i}_0}{u}_{i_0,s}^{\left(m\right)}\right|-\left|\underset{j\ne i_0}{{\displaystyle \sum }}{m}_{i_{0}j}{u}_{j,s}^{\left(m\right)}\right|+\left|\tau \rho \left(u_{i_0,s}^{\left(m\right)}-u_{i_0,s}^{*}\right)\right|\\ \le \left|m_{i_0{i}_0}{u}_{i_0,s}^{\left(m\right)}+\underset{j\ne i_0}{{\displaystyle \sum }}{m}_{i_{0}j}{u}_{j,s}^{\left(m\right)}\right|+\left|\tau \rho \left(u_{i_0,s}^{\left(m\right)}-u_{i_0,s}^{*}\right)\right|\end{array}$

$\begin{array}{l}=\left|{{\displaystyle \sum }}^{\text{}}{m}_{i_{0}j}{u}_{j,s}^{\left(m\right)}\right|+\left|\tau \rho \left(u_{i_0,s}^{\left(m\right)}-u_{i_0,s}^{*}\right)\right|\\ =\left|{{\displaystyle \sum }}^{\text{}}{m}_{i_{0}j}{u}_{j,s}^{\left(m\right)}-\tau \rho \left(u_{i_0,s}^{\left(m\right)}-u_{i_0,s}^{*}\right)\right|\\ =\left|u_{i_0,s}^{\left(m-1\right)}+\tau f_{i_0,s}-\tau q_{ss}{u}_{i_0,s}^{\left(m-1\right)}-\tau \underset{s^{*}\ne s}{{\displaystyle \sum }}{q}_{s,s^{*}}{u}_{i_0,s^{*}}^{\left(m-1\right)}\right|\\ \le \left|\tau \right|{\Vert f_s\Vert }_{\infty }+\left(1+\left|\tau \right|\left|q_{s,s}\right|\right){\Vert u_s^{\left(m-1\right)}\Vert }_{\infty }+\left|\tau \right|\underset{s^{*}\ne s}{{\displaystyle \sum }}{q}_{s,s^{*}}{\Vert u_{s^{*}}^{\left(m-1\right)}\Vert }_{\infty }\end{array}$ (4.8)

Case 2: For sufficiently small $\delta t$ such that $\delta t\le \frac{1}{2+4\alpha +2\lambda C_1}$ ，we have ${\Vert {\text{e}}^m\Vert }^2\le C\left({\Vert {\text{e}}^0\Vert }^2+{\Vert {\text{e}}^1\Vert }^2+{\max }_{2\le j\le m}{\Vert {\delta }^j\Vert }^2\right),\forall 2\le m\le \frac{T}{\delta t}$ , $\alpha =\left|\frac{{\left(r-\frac{{\sigma }^2}{2}-\lambda k\right)}^2-2\left(r+\lambda \right){\sigma }^2}{2{\sigma }^2}\right|$

and C is generic constant depending on the parameter $C_1,r,\sigma ,\lambda$ and T.

This case can be proved (see Theorem1 in [11] ).

5. A Numerical Simulation Example

In this section, we present several numerical experiments to illustrate the efficiency and accuracy of proposed method.

Example: We will simulate the HS300 to fit the SVJ model. The sample number is 10,000. For the convenience of simulation, we will eliminate half of the abnormal data. The stability and effectiveness of time-discretization are verified by the simulation of the remaining data. All numerical experiments result at strike price with parameters as provided in Table 1.

The s-d and MC-error indicates the error of the simulation. In Table 1, we can find that the data are so small that it can be ignored.

Finally, through the above data simulation, we proved that the proposed method is much more efficient and accurate.

6. Discussion and Conclusion

In this article, it has proposed and analyzed three implicit-explicit (IMEX) time

Table 1. Numerical simulation analysis.

semi-discretization namely, IMEX-BDF2, for solving partial integro-differential equations which arise in option pricing theory when the underlying asset follows a jump diffusion process. All the IMEX time semi-discretization is shown to be stable. The American option whose types of initial and boundary conditions cannot be defined easily, but it can be formulated as the linear complementarity problem (LCP). The numerical simulation with European put/call under SVJ model has been carried out.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Black, F. and Scholes, M. (1973) The Pricing of Options and Corporate Liabilities. Journal of Political Economy, 81, 637-659. https://doi.org/10.1086/260062
[2]	Merton, R.C. (1976) Option Pricing when Underlying Stock Returns Are Discontinuous. Journal of Financial Economics, 3, 125-144. https://doi.org/10.1016/0304-405X(76)90022-2
[3]	Xu, W., Wu, C., Xu, W. and Li, H. (2002) A Jump-Diffusion Model for Option Pricing under Fuzzy Environments. Management Science, 48, 1086-1101. https://doi.org/10.1287/mnsc.48.8.1086.166
[4]	Hull, J. and White, A. (1987) The Pricing of Options on Assets with Stochastic Volatilities. The Journal of Finance, 42, 281-300. https://doi.org/10.1111/j.1540-6261.1987.tb02568.x
[5]	Heston, S.L. (1993) A Closed-Form Solution for Options with Stochastic Volatility with Applications to Bond and Currency Options. Review of Financial & Studies, 6, 327-343. https://doi.org/10.1093/rfs/6.2.327
[6]	Almendral, A. and Oosterlee, C.W. (2005) Numerical Valuation of Options with Jumps in the Underlying. Applied Numerical Mathematics, 53, 1-18. https://doi.org/10.1016/j.apnum.2004.08.037
[7]	Ngounda, E., Patidar, K.C. and Pindza, E. (2013) Contour Integral Method for European Options with Jumps. Communications in Nonlinear Science & Numerical Simulation, 18, 478-492. https://doi.org/10.1016/j.cnsns.2012.08.003
[8]	D’Halluin, Y., Forsyth, P.A. and Vetzal, K.R. (2003) Robust Numerical Methods for Contingent Claims under Jump Diffusion Processes. IMA Journal of Numerical Analysis, 25, 87-112.
[9]	Kwon, Y. and Lee, Y. (2011) A Second-Order Finite Difference Method for Option Pricing under Jump-Diffusion Models. SIAM Journal on Numerical Analysis, 49, 2598-2617. https://doi.org/10.1137/090777529
[10]	Ascher, U.M., Ruuth, S.J. and Spiteri, R.J. (1997) Implicit-Explicit Runge-Kutta Methods for Time-Dependent Partial Differential Equations. Elsevier Science Publishers B. V. https://doi.org/10.1016/S0168-9274(97)00056-1
[11]	Kadalbajoo, M.K., Kumar, A. and Tripathi, L.P. (2016) A Radial Basis Function Based Implicit-Explicit Method for Option Pricing under Jump-Diffusion Models. Applied Numerical Mathematics, 110, 159-173. https://doi.org/10.1016/j.apnum.2016.08.006