Financial Applications on Fractional Lévy Stochastic Processes

fractal and fractional Article Financial Applications on Fractional Lévy Stochastic Processes Reem Abdullah Aljethi 1,2 and Adem Kılıçman 1, * 1 2 * Department of Mathematics and Statistics, Institute for Mathematical Research, Universiti Putra Malaysia, UPM Serdang, Seri Kembangan 43400, Selangor, Malaysia; gs52345@student.upm.edu.my Department of Mathematics, Imam Mohammad Ibn Saud Islamic University, Riyadh 11564, Saudi Arabia Correspondence: akilic@upm.edu.my Abstract: In this present work, we perform a numerical analysis of the value of the European style options as well as a sensitivity analysis for the option price with respect to some parameters of the model when the underlying price process is driven by a fractional Lévy process. The option price is given by a deterministic representation by means of a real valued function satisfying some fractional PDE. The numerical scheme of the fractional PDE is obtained by means of a weighted and shifted Grunwald approximation. Keywords: European option pricing; Lévy process; fractional diffusion equations; sensitivity analysis 1. Introduction and Notations The pricing of financial instruments is a problem of significant interest for risk management and hedging and/or in trading derivatives expecting more income in finance industries. Financial instruments are linked to underlying assets whose dynamics are modeled by an exponential of a Lévy process of the form



 1 d(log St ) = (µ − σ2 )dt + σdLt , 2 Citation: Aljethi, R.A.; Kılıçman, A. Financial Applications on Fractional Lévy Stochastic Processes. Fractal Fract. 2022, 6, 278. https://doi.org/ 10.3390/fractalfract6050278 Academic Editor: Ricardo Almeida Received: 23 February 2022 Accepted: 13 May 2022 Published: 22 May 2022 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in (1) under the risk neutral measure Q, where µ > 0 is the average compounded growth of the stock St , dLt is the increment of Lévy and σ is the volatility of the risky asset, induced by Lévy jump. We refer to Raible (2000) [1], Schoutens [2] (2003) and Eberlein (2014) [3] and the references therein for more details on non-fractional Lévy models in finance. It is well known that the distribution of a Lévy process is characterized by its characteristic function, which is given by the Lévy–Khintchine formula. We can characterize the Lévy when a Lt dependent random variable with time is a Lévy process if and only if it has independent stationary increments with an exponential characteristic function given by the so-called Lévy–Khintchine representation. Indeed, if ( Lt )t≥0 is a Lévy process on R with characteristic triplet (m, σ, ν), then E[ei̟Lt ] = etψ(̟ ) , t ∈ R. The function ψ is called the characteristic exponent and can be expressed as published maps and institutional affiliations. Copyright: © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article ψ(̟ ) = im̟ + σ2 (i̟ )2 + 2 Z ∞ −∞ (ei̟κ − 1 − i̟I|κ |<1 )ν(dκ ), (2) where m ∈
R, σ ≥ 0 and ν is the Lévy measure that controls the jumps such that R min 1, κ 2 ν(dκ ) is finite and I A is the indicator function defined by R distributed under the terms and conditions of the Creative Commons I A ( x ) :=  1 0 if x ∈ A otherwise. Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). Fractal Fract. 2022, 6, 278. https://doi.org/10.3390/fractalfract6050278 https://www.mdpi.com/journal/fractalfract Fractal Fract. 2022, 6, 278 2 of 12 An Example of a Lévy Process: Lévy Stable Processes (LSP) α,β Let Lt be a Lévy α-stable process with skew parameter β and stability index α. We shall take β = −1 and assume that 1 < α < 2. Consequently, for each n > 0 and 0 < t < T, h  i   πα  −1 α α E exp nσLα, = exp − tn σ sec < +∞. t 2 Recall that the density function νLS ( x ) = is given by ( νLS ( x ) = ν(dx ) dx ad | x | 1+ α bd x 1+ α of the stable Lévy process with α ∈ (0, 2) for x < 0, for x > 0 where 1 = a + b and d > 0. Substituting the above density in Equation (2) yields the characteristic exponent of a LSP with β = a − b ψLS (̟ ) = i̟m − απ σα |̟ |α × [1 − iβtan( )sign(̟ )]. 2 2 An equivalent from of ψLS (̟ ) is given below ψLS (̟ ) = i̟m − σα [(i̟ )α (1 − β) + (−i̟ )α (1 + β)] 4cos( απ 2 ) where the condition β = −1, yield a = 0 and b = 1 (see Cartea et al. [4] for more details), hence 1 σα [(2)(i̟ )α ]. (3) ψLS (̟ ) = i̟m − ) 4cos( απ 2 The price of an European call option may be expressed as the risk-neutral conditional expectation V (St , t) = e−r(T −t) EQ [max(ST − K, 0)|Ft ], and the value of an European put option may be expressed as the risk-neutral conditional expectation V (St , t) = e−r(T −t) EQ [max(K − ST , 0)|Ft ], where the price of the stock follows the dynamic (1) and Ft is the completed and right continuous natural filtration associated to the Lévy stable process. Other extensions to this result have been done in many directions; see Abouagwa and Li [5] (2019) and also Wu et al. [6] (2021). Our contribution in this paper is to perform a numerical analysis of the value of the European style options as well as a sensitivity analysis for the option price with respect to the model parameters. Numerical methods have been studied using a finite difference scheme for non-fractional jump diffusion and exponential Lévy models for option pricing by Cont and Voltchkova [7] (2005). For the sensitivity analysis, we shall focus mainly on the very well-known Greeks (Delta, Gamma and Vega) and refer to a relatively recent paper by Shokrollahi et al. [8] (2015). The rest of the article is structured as follows: Section 2 is devoted to the stochastic model we are interested in as well as the deterministic representation of the option price by a real valued function satisfying some fractional PDEs. Section 3 presents the numerical scheme used to discretize the fractional PDE by means of the weighted and shifted Grunwald approximation. Section 4 deals with a numerical implementation to solve the fractional diffusion equation. Moreover, a subsection concerning the numerical computations of Greek parameters giving some information about sensitivity analysis is added. Fractal Fract. 2022, 6, 278 3 of 12 2. The Stochastic Model and Deterministic Representation −1 The dynamic of the stock price driven by the α-stable Lévy process Lα, with addit tional impact term is given for each 0 ≤ t < T by the following equation:   −1 dγ St = St (r − q)dtγ + σdLα, + λSt dβγt , S0 initial price, (4) t γ where 0 < γ < 1 and 1 < α ≤ 2, β t is the number of shares of the stock at time t and the γ term λdβ t (λ ≥ 0) stands for the impact price of the investor’s trading strategy (see Chen et al. (2014) [9] for more details and specific parameters). Throughout this paper, we only consider trading strategies that can be written in the following form γ −1 dβ t = η dtγ + ζ dLα, , (5) t
γ γ with initial number of shares β 0 and some constants η, ζ. A self-financing strategy θt , β t t≥0 represented by the wealth process (Vt )t≥0 for trader is given by γ Vt = θt St0 + β t St = V0 + Z t 0 θu dSu0 + Z t γ 0 β u dSu , where dSt0 = rSt0 dt, S00 = 1. According to the recent paper by Aljedhi and Kılıçman (2020) [10], the wealth process (Vt )t≥0 can be written as Vt = u(St , t) where u( x, t) is the solution of the fractional partial differential equation (FPDE) ∂γ u( x, t) ∂u( x, t) ∂α u( x, t) +A = ru( x, t), +B γ ∂t ∂x ∂x α (6) under the conditions u( x, 0) = h( x ), and ut ( x, 0) = f ( x ). Here, A and B are constants and α σα απ their expressions are A = r − q + λη + σ2 sec( απ 2 ) and B = 2 sec( 2 ). Notice also that based on the martingale pricing approach under a risk neural probability measure, u(St , t) can be represented as u(St , t) = EQ [φ(ST )|Ft ], where φ( x ) :=  European call option, European put option, ( x − K )+ (K − x )+ . The first goal of this paper is to find numerical values of the price of European-style option, based on implicit scheme for the (FPDE) (6) and the Euler–Maruyama scheme for the (4). The second goal is to study the sensitivity of these European options while parameters of the model change. In our numerical analysis, we combine the method of Demirci et al. (2012) [11] and the method adopted by Shen et al. (2012) [12] to solve our fractional partial differential equation. In particular, in [11], they dealt with Caputo partial differential equations with an initial condition of the form D γ y ( t ) = g ( y ( t ), t ), y (0) = y0 . It can be solved as a solution y∗ (ω ) to a non-fractional differential equation   tγ y(t) = y∗ . Γ ( γ + 1) Fractal Fract. 2022, 6, 278 4 of 12 Moreover, Jumarie (2010) [13] solved the equation, dy = g(t)dtα , t ≥ 0, 0 < α ≤ 1, y(0) = y0 , where the integral of dtα has the form Z t 0 g( x )dx α = α Z t 0 (t − x )α−1 g( x )dx (7) and have used the expression dα x = Γ(1 + α)dx. Based on our specific structure of the strategy described in (5), the dynamic of the stock price process (4) becomes   −1 −1 dγ St = St (r − q)dtγ + σdLα, + λSt (ηt dtγ + ζ t dLα, ). (8) t t It is well known that, using Itô’s formula to xt = log St , the (8) has a solution that can be written as   Z T (τ )γ α,−1 x T = xt exp , (9) (r − µλη − q) + (σ + λζ t )dLu Γ (1 + γ ) t where τ = T − t. We shall use the Euler–Maruyama scheme to find some realizations of the price process and then obtain some simulations of the underlying asset described in Equation (9) by means of the following discrete scheme for 1 < α ≤ 2: x n +1 = x n + υn +1 where x0 = log S0 , tn = nT N 1 (τn+1 )γ −1 (r − q + λη − m) + υnα+1 (σ + λζ ) Lα, n +1 , Γ (1 + γ ) for n = 0, . . . , N, τn+1 = T − tn+1 , υn+1 = −1 α,−1 α,−1 Lα, n +1 = L t n +1 − L t n . (10) σα sec( απ 2 ) and 2 3. Numerical Discretization of the Fractional PDE The price of the call or put European option associated to the underlying asset governed by the (9) is given by u(St , t) = EQ [ h(ST )|Ft ], where u( x, t) is the solution to the Fractional PDE (0 ≤ α, γ ≤ 2) ∂γ u( x, t) ∂u( x, t) ∂α u( x, t) = A + ru( x, t) + B ∂tγ ∂x ∂x α (11) with the conditions u( x, 0) = h( x ), and ut ( x, 0) = f ( x ). Here A and B are constants, and their expressions are A = r − q + λη + σα απ sec( ), 2 2 B= σα απ sec( ). 2 2 We shall use the finite difference method using the shifted Grunwald approximation of the fractional derivative. To do this, we will be working with the notation uin = u( xmin + ih, nτ ) with h as the step-size of the space variable and τ as the step-size of the time variable i = 0, . . . , Nx and n = 1, . . . , Nt . We consider now (11) at ( xi , tn ), one gets: ∂un ∂α un ∂γ uin = A i + B αi + ruin . γ ∂t ∂x ∂x (12) Fractal Fract. 2022, 6, 278 5 of 12 Before we move on to the method itself, we should first recall the definition for the Caputo time-fractional derivative [12] as 1 ∂γ u( x, tn+1 ) = γ ∂t Γ (1 − γ ) = = 1 Γ (1 − γ ) n Z t k +1 ∑ k =0 t k n u( x, tn+1 ) − u( x, tk ) ∑ τ k =0 n τ −γ Γ (2 − γ ) ∂u ∂η ( x, η ) dη (tn − η )γ Z t k +1 tk dη + O(τ 2−γ ) (tn − η )γ γ ∑ bk [u(x, tn+1−k ) − u(x, tn−k )] + O(τ2−γ ) k =0 γ where the weights bk satisfy the following formula γ γ bk = (k + 1)γ − k1−γ , b0 = 1, k = 0, 1, · · · , n. For the Riemann–Liouville space derivative with the weighted and shifted Grunwald approximation [14,15], we use 1 ∂α u( x, t) = α ∂x α h m ∑ ωkα u(x, tm−k ) + O(h), k =0 where the Grünwald–Letnikov weights ωkα = (−1)k   α , k ≥ 0 are the coefficients of the k ∞ power series of the generating function (1 − z)k = ∑ ωkα zk . These coefficients satisfy the k =0 recursive formula ωkα =  1−  α+1 ωkα−1 , k ω0α = 1, k = 1, · · · , m. Hence, the numerical scheme of (12) is given by τ −γ Γ (2 − γ ) n −1 γ ∑ bn − k − 1 k =0  
A n B ui+1 − uin−1 + α uik+1 − uik = 2h h i +1 ∑ ωkα uin−k+1 + ruin . (13) k =1 where γ γ bk = (k + 1)γ − k1−γ , b0 = 1,   α+1 ωkα−1 , w0α = 1. ωkα = 1 − k By developing the left side of the Equation (13), one gets n −1 τ −γ γ uin − ∑ (bn−k−1 − bn−k )uik − bn−1 u0i Γ (2 − γ ) k =1 ! =
A n B u − uin−1 + α 2h i+1 h i +1 ∑ ωkα uin−k+1 + ruin . (14) k =1 Rearranging and substituting uin in Equation (14) by its numerical approximation Uin , one obtains for i ∈ {1, . . . , Nx − 1} and n ∈ {1, . . . , Nt − 1} : i +1 n −1
γ Uin − µ1 Uin+1 − Uin−1 − µ2 ∑ ωkα Uin−k+1 = µ3 ∑ (bn−k−1 − bn−k )Uik + µ0 Ui0 , k =1 U0n = u( x0 , tn ), (15) k =1 n UN = u( x Nx , tn ), x (16) Fractal Fract. 2022, 6, 278 6 of 12 where µ0 = bn − 1 τ γ Γ ( 2 − γ ) τ γ AΓ(2 − γ) , µ1 = , 1 − rΓ(2 − γ) 2h(1 − rΓ(2 − γ)) and µ2 = τ γ BΓ(2 − γ) hα (1 − rΓ(2 − γ)) , µ3 = 1 . 1 − rτ γ Γ(2 − γ) In order to develop the numerical result, we must write Equation (15) in matrix form = U n−1 , which is completed with boundary conditions (16). We note that, unlike in integer-order PDEs, the coefficient matrix M is a square matrix of size Nx , and will not be tri-diagonal, which will slow down our numerical algorithm quite significantly. This is due to the non-local nature of the time-fractional and space-fractional derivatives. The structure of M is described in the following way, for i = 1, . . . , Nx − 2 and j = 1, . . . , Nx − 2 MU n Mi,j  −µ2 ωiα− j+1      −µ1 − µ2 ω2α = 1 − µ2 ω1α   α    − µ 1 − µ 2 ω0 0 if if if if if j j j j j < i−1 = i−1 =i = i+1 > i+1 To complete this matrix with the boundary conditions, we have that M0,0 = 1, M0,j = 0 for j = 1, . . . , Nx − 1, M Nx −1,Nx −1 = 1, M Nx −1,j = 0 for j = 1, . . . , Nx − 1 and the second member is given by Uin−1 = µ3 n −1 γ ∑ (bn−k−1 − bn−k )Uik + µ0 Ui0 . k =1 Now, all we have to deal with are the boundary conditions. For European-type options, the behavior of numerical scheme is similar regardless of performing call and put options. However, in the context of our Lévy stable model, this is not the case. Hence, we describe each case separately. • • For a call option, the payoff is defined by Π T = max(ST − K, 0). Since we must truncate the spatial domain in order to make it workable from a numerical point of view, we will use the following boundary conditions u( xmin , t) = 0 and u( xmax , t) = e xmax − Ke−r(T −t) where xmin = − log(4K ) and xmax = log(4K ). For a put option, the payoff is defined by Π T = max(K − ST , 0) and we will use the following conditions u( xmax , t) = 0 and u( xmin , t) = Ke−r(T −t) − e xmin where xmin = − log(4K ) and xmax = log(4K ). 4. Numerical Simulations and Discussions In this section, we present some numerical simulations in order to validate the stochastic model with associated partial differential equations and illustrate the utility of the novel model. In addition, we compute the price of call and put options in order see the difference with the classical Black and Schole’s price. Finally, we perform the sensitivity analysis, usually applied by financial engineers, that provide useful information for the investors. To this end, we use the parameter settings reported in Table 1 for the one-year option. Table 1. Parameters’ value. S0 Strike K r q λ η m σ β α 100 80 0.02 1 10 0.01 0.02 0.15 −1 1.76 Fractal Fract. 2022, 6, 278 7 of 12 Under the risk neutral equivalent probability measure and thanks to the martingale property, the value of an option at the present time can be represented as u(St , t) = EQ [(K − ST )+ |Ft )] (17) that is the conditional expected value of the discounted payoff under a chosen risk neutral measure with respect to the all the information available in the market up to time t. This representation allows us to use Monte Carlo techniques to estimate values of the option price by simulating a large number of sample values of ST . Notice that this method uses the risk-neutral valuation approach, which means that one has to take m = r. In order to compare the numerical solutions using others numerical computations, we use the Euler–Maruyama method for stochastic differential equation associated with our Lévy stable process. However, the numerical solution of the FPDE is obtained using the finite difference method cited in (15), which is reported in the Figure 1. Moreover, the numerical solution of our SDE is obtained by the Euler–Maryuma scheme. Figure 1. Comparison between numerical solutions (SDE and FPDE). Upon the comparison, using the same parameters, of the left hand side of Equation (17) obtained by a FPDE approach and the right hand side of Equation (17) computed using Monte Carlo simulation, we observe that both solutions are close to each other. In order to validate our Lévy model, we firstly solve the SDE and FPDE for call options. Secondly, we present the solutions for put option. In Tables 2 and 3, we present the call and put prices obtained by SDE and FPDE under the same parameters. Table 2. Call price obtained with Euler scheme and finite difference. Upon comparison, the two call prices are very closed to each other, with small error Emax which is defined by Emax = max |u(xi , tn ) − Uin |. Strike K Call Price with SDE Call Price with FPDE Emax 80 43.9820 43.9789 0.0040 100 40.9454 40.8985 0.0068 110 39.8865 39.8783 0.0082 Fractal Fract. 2022, 6, 278 8 of 12 Table 3. Put price obtained with Euler scheme and finite difference. Upon comparison, the two put prices are very close to each other, with small error Emax which defined by Emax = max |u(xi , tn ) − Uin |. Strike K Put Price with SDE Put Price with FPDE Emax 80 26.0659 26.1043 0.0384 100 5.5968 5.5481 0.0487 110 2.7987 2.7683 0.0304 In what follows, in Figure 2, we plot the price of the call obtained by solving the fractional PDE (12). Figure 2. This figure shows the evolution of the call option price with respect to the time to maturity using the FPDE representation. In comparison with the Black–Scholes call price when the underlying is driven by the Brownian motion, the two prices have the same shape. Secondly, we are interested in the calculus of the put price by using the FPDE and SDE. We use the same parameters setting in Table 1 completed by the following conditions 1 . u( xmax , 0) = 0 and u( xmin , 0) = Ke−rT − 4K In Figure 3, we plot the price of the put obtained by solving the fractional PDE (12). Fractal Fract. 2022, 6, 278 9 of 12 Figure 3. This figure shows the evolution of the call put option price with respect to the time to maturity using the FPDE representation. In comparison with the Black–Scholes put price when the underlying is driven by the Brownian motion, the two prices have the same shape. Approximation of Lévy Model Greeks In the Black and Scholes model, the derivation and analytic expressions for the Greeks for put and call prices can be done. We refer to De Olivera and Mordecki (2014) [16] for the computation of Greeks using the Fourier transform approach. However, due to the complexity of our model, we chose to use finite differences to approximate the derivatives. Then, we use the finite differences to find the Greeks (Delta “∆”, Vega “Vega”, Gamma “Γ”, Rho “Rho”). When this method is used to calculate Delta and Vega Greeks (first order derivative), however, the computation time is increased because the option price must be calculated more than once, multiple times in the case of the Gamma Greek (secondorder derivative). In order to approximate the first-order derivative of Greek sensitivity of the option price with respect to a parameter or to variable, we use the first-order central differences: ∆ Levy = and Vega Levy = ∂u(St , r, σ, t) u(S + dS, r, σ, t) − u(S − dS, r, σ, t) := , ∂S 2dS u(S, r, σ + dσ, t) − u(S − dS, r, σ − dσ, t) ∂u(St , r, σ, t) := , ∂σ 2dσ (18) (19) and u(S, r + dr, σ, t) − u(S, r − dr, σ, t) ∂u(St , r, σ, t) := . (20) ∂r 2dr For the second-order derivative of Greeks, such as Gamma, we use the second order central differences for a single variable: Rho Levy = Γ Levy = ∂2 u(St , r, σ, t) u(S + dS, r, σ, t) − 2u(S, r, σ, t) + u(S − dS, r, σ, t) := . ∂S2 (dS)2 (21) Fractal Fract. 2022, 6, 278 10 of 12 To illustrate the Greek sensitivity analysis for put and call option prices, we develop MATLAB codes. It is informative to present visual simulations of the Greeks (Delta, Vega, Gamma and Rho) in 3D. Indeed, the Figure 4 presents a plot of the Greeks for the considered Lévy model. In Figure 5, we plot the Greeks for the put option in 3D. We show that the Delta is constant with value equal to 1.51 for put option’s price for Lévy model, but it oscillates for Black–Scholes. We have also the same remarks for Gamma, Vega and Rho. Figure 4. Lévy model Greeks for the call option: By analyzing the Greeks plots, we observe that the call option’s price have much higher Delta values than out of the call option’s price of Black–Scholes model, and this value oscillates around 2.5, which ranges between 2.49 and 2.51. Gamma reaches its maximum when the underlying price is a little bit smaller, exactly equal to the strike of the call option, and the chart shows that Gamma is significantly constant for the Lévy model. We display also Vega as a function of the asset price and time to maturity for a call options with strike at 80 and we highlight the fact that Vega is much higher for the call option’s price than for the call option’s price for Black–Scholes. At last, Rho reaches its maximum when the underlying price is a little smaller, not exactly equal to the strike of the call option’s price, and it is significantly higher with hyper-volatility than for the call option’s price of Black–Scholes. Fractal Fract. 2022, 6, 278 11 of 12 Figure 5. Lévy model Greeks for the put option: in the figures above, we plotted the Greeks for the put option in 3D. We observe that the Delta is constant with value equal to 1.51 for put option’s price for Levy model, but it oscillates for Black–Scholes. We can see that the put–call parity is maintained for the Lévy model: Vega(Call) = Vega(Put) and Gamma(Call) = Gamma(Put). We have a negative Rho which ranges between −1000 and 0, and the figure displays its fluctuation with respect to the underlying asset. 5. Conclusions Through this work, the fractional diffusion equation in terms of the Lévy process is solved. We also discussed the Euler–Maruyama equation for solving the price of fractional financial derivatives of European options price. Lévy processes have proven to be a suitable method that strikes the right balance between the properties of the mathematical approach and the value of option evolution. Furthermore, we made a comparison between European call and put option prices based on the Lévy diffusion equation and provided numerical solutions of the fractional partial differential equation that were derived by the fractional diffusion model. Author Contributions: Both authors (R.A.A. and A.K.) contributed equally to this work. All authors have read and agreed publishing this manuscript. Funding: The second author would like to acknowledge that this research is partially supported by Ministry of Higher Education under the Fundamental Research Grants Scheme(FRGS) with project number FRGS/2/2014/SG04/UPM/01/1 and having vot number 5524674. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: Not applicable. Acknowledgments: The authors would like to thank the reviewers for their fruitful comments that have improved the presentation of the paper. Conflicts of Interest: The authors declare there is no conflict of interest. Fractal Fract. 2022, 6, 278 12 of 12 References 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. Raible, S. Lévy Processes in Finance: Theory, Numerics, and Empirical Facts. Ph.D. Thesis, University of Freiburg, Freiburg im Breisgau, Germany, 2000. Schoutens, W. Lévy Processes in Finance, 1st ed.; Wiley Series in Probability and Statistics; John Wiley & Sons: Hoboken, NJ, USA, 2003. Eberlein, E. Fourier—Based Valuation Methods in Mathematical Finance. In Quantitative Energy Finance; Benth, F.E., Kholodnyi, V.A., Laurence, P., Eds.; Springer: Berlin/Heidelberg, Germany, 2014; pp. 85–114. Cartea, A.; del-Castillo-Negrete, D. Fractional Diffusion Models of Option Prices in Markets with Jumps; Elsevier: Amsterdam, The Netherlands, 2006. Abouagwa, M.; Li, J. Stochastic fractional differential Equations driven by Lévy noise under Carathéodory conditions. J. Math. Phys. 2019, 60, 022701. [CrossRef] Wu, P.; Yang, Z.; Wang, H.; Song, R. Time fractional stochastic differential equations driven by pure jump Lévy noise. J. Math. Anal. Appl. 2021, 504, 125–412. [CrossRef] Cont, R.; Voltchkova, E. Finite difference methods for option pricing in jump diffusion and exponential Lévy models. SIAM J. Numer. Anal. 2005, 43, 1596–1626. [CrossRef] Shokrollahi, F.; Kılıçman, A.; Ibrahim, N.A. Greeks and Partial Differential Equations for some Pricing Currency Option Models. Malays. J. Math. Sci. 2015, 9, 417–442. Chen, W.; Xu, X.; Zhu, S.P. Analytical pricing European-style option under the modified Black-Scholes equation with a partialfractional derivative. Quartely Appl. Math. 2014, 72, 597–611. [CrossRef] Aljedhi, R.A.; Kılıçman, A. Fractional Partial Differential Equations Associated with Lévy Stable Process. Mathematics 2020, 8, 508. [CrossRef] Demirci, E.; Ozalp, N. A method for solving differential equations of fractional order. J. Comput. Appl. Math. 2012, 236, 2754–2762. [CrossRef] Shen, S.; Liu, F.; Chen, J.; Turner, I.; Anh, V. Numerical techniques for the variable order time fractional diffusion equation. Appl. Math. Comput. 2012, 218, 10861–10870. [CrossRef] Jumarie, G. Derivation and solutions of some farctional Black-Scholes equations in space and time. J. Comput. Math. Appl. 2010, 59, 1142–1164. [CrossRef] Yuste, S.B. Weighted average finite difference methods for fractional diffusion equations. J. Comput. Phys. 2006, 216, 264–274. [CrossRef] Lubich, C. Discretized Fractional Calculus. SIAM J. Math. Anal. 1986, 17, 704–719. [CrossRef] De Olivera, F.; Mordecki, E. Computing Greeks for Lévy Models: The Fourier Transform Approach. arXiv 2014, arXiv:1407.1343.