• Received : 2019.03.07
  • Accepted : 2019.03.12
  • Published : 2019.03.25


In this paper, we develop an accurate explicit finite difference method for the two-dimensional Black-Scholes equation with a hybrid boundary condition. In general, the correlation term in multi-asset options is problematic in numerical treatments partially due to cross derivatives and numerical boundary conditions at the far field domain corners. In the proposed hybrid boundary condition, we use a linear boundary condition at the boundaries where at least one asset is zero. After updating the numerical solution by one time step, we reduce the computational domain so that we do not need boundary conditions. To demonstrate the accuracy and efficiency of the proposed algorithm, we calculate option prices and their Greeks for the two-asset European call and cash-or-nothing options. Computational results show that the proposed method is accurate and is very useful for nonlinear boundary conditions.


An option is the right to buy or sell a basic product at a specific price in the future. It is important for financial companies or customers to know the price of option. There are many mathematical models and computation methods [1, 2, 3, 4, 5, 6, 7, 8] to determine option pricing. The Black–Scholes (BS) partial differential equation (PDE) [9] is commonly used to measure the option price and its Greeks.

In [10], the authors improved the model of Mortensen [11] and showed that a rich class of recovery rate scenarios could be incorporated into the model in a computationally manageable way. Paul Glasserman [12] focused on improving efficiency and discussed applications of the Monte Carlo simulation (MCS) [13, 14] on security pricing issues. Milstein [15] calculated the option price using the MCS and then calculated the Greeks. However, the MCS takes a lot of time because it has many execution time. Because the option value is different at each time, the fluctuation of the Greeks is too high.

One of the other numerical methods is the finite difference method (FDM) [16, 17, 18, 19]. Foulon [20] used semi-discretization of the Heston PDE by using FDM on non-uniform grids and solved the equation using the alternating direction implicit (ADI [21]) on two dimensions. In [22], the authors suggested operator splitting method (OSM) for solving the linear complementarity problems arising from the pricing of American options with two assets. In [23], the authors proposed the numerical contour integral method for free-boundary PDE to solve American volatility option pricing. In [25], the authors used the Laplace Transform Method with two free boundaries for pricing American CEV strangles Option.

The problematic part in solving a two-dimensional BS equation with an implicit method is the cross differential term. Therefore, the cross differential term is typically treated as an explicit manner. In particular, at the corner point, i.e., where both the underlying assets are maximum, the solution is highly dependent on the boundary treatment. In this paper, we extend the previous work [24] for the one-dimensional no boundary method to the two-dimensional BS equation:

\(\frac{\partial u}{\partial t}=-\frac{\sigma_{1}^{2} x^{2}}{2} \frac{\partial^{2} u}{\partial x^{2}}-\frac{\sigma_{2}^{2} y^{2}}{2} \frac{\partial^{2} u}{\partial y^{2}}-\rho \sigma_{1} \sigma_{2} x y \frac{\partial^{2} u}{\partial x \partial y}-r x \frac{\partial u}{\partial x}-r y \frac{\partial u}{\partial y}+r u\)       (1.1)

for (x, y, t) ∈ (0, ∞) × (0, ∞) × [0, T)

with the final payoff condition u(x, y, T) = payoff(x, y) at expiry T, where x, y are the value of underlying assets, u(x, y, t) is the option price, r is the risk–free rate, ρ is the correlation value between x, y and σ1, σ2 are the volatility of underlying assets, respectively. After applying the change of variable τ = T − t to Eq. (1.1) and truncate the infinite domain into a finite computational domain, we have the following natural initial boundary value problem

\(\frac{\partial u}{\partial \tau}=\frac{\sigma_{1}^{2} x^{2}}{2} \frac{\partial^{2} u}{\partial x^{2}}+\frac{\sigma_{2}^{2} y^{2}}{2} \frac{\partial^{2} u}{\partial y^{2}}+\rho \sigma_{1} \sigma_{2} x y \frac{\partial^{2} u}{\partial x \partial y}+r x \frac{\partial u}{\partial x}+r y \frac{\partial u}{\partial y}-r u\)       (1.2)

for (x, y, τ ) ∈ Ω × (0, T],

where Ω = (0, xmax) × (0, ymax).

The main purpose of this paper is to develop a hybrid boundary condition for the twodimensional BS equation so that we can avoid extrapolation problem at the corner point and get stable numerical solutions.

This paper is organized as follows. In Section 2, we provide the numerical solution algorithm for the two-dimensional BS model. In Section 3, we present the computational results using the proposed hybrid boundary condition. Finally, conclusions are drawn in Section 4.


We discretize the BS equation on a non-uniform grid as shown in Fig. 1. Let hi−1 = xi −xi−1 and hj−1 = yj −yj−1 be the non-uniform space step sizes in the x- and y-directions, respectively. Here, x0 = y0 = 0. Let us denote u(xi , yj , τn ) by \(u^n_{ij}\), which is a numerical approximation of the solution, where τn = n∆τ and ∆τ is a time step. The two-dimensional BS Eq. (1.2) is discretized as

 E1TAAE_2019_v23n1_19_f0001.png 이미지

FIGURE 1. Temporal evolution of computational domain.

\(\begin{aligned} \frac{u_{i j}^{n+1}-u_{i j}^{n}}{\Delta \tau}=& \frac{\left(\sigma_{1} x_{i}\right)^{2}}{2}\left(\frac{\partial^{2} u}{\partial x^{2}}\right)_{i j}^{n}+\frac{\left(\sigma_{2} y_{j}\right)^{2}}{2}\left(\frac{\partial^{2} u}{\partial y^{2}}\right)_{i j}^{n} \\ &+\rho \sigma_{1} \sigma_{2} x_{i} y_{j}\left(\frac{\partial^{2} u}{\partial x \partial y}\right)_{i j}^{n}+r x_{i}\left(\frac{\partial u}{\partial x}\right)_{i j}^{n}+r y_{j}\left(\frac{\partial u}{\partial y}\right)_{i j}^{n}-r u_{i j}^{n} \end{aligned}\)       (2.1)


\(\begin{aligned} \left(\frac{\partial u}{\partial x}\right)_{i j}^{n} &=\frac{h_{i-1}}{h_{i}\left(h_{i-1}+h_{i}\right)} u_{i+1, j}^{n}+\frac{h_{i}-h_{i-1}}{h_{i-1} h_{i}} u_{i j}^{n}-\frac{h_{i}}{h_{i-1}\left(h_{i-1}+h_{i}\right)} u_{i-1, j,}^{n} \\ \left(\frac{\partial u}{\partial y}\right)_{i j}^{n} &=\frac{h_{j-1}}{h_{j}\left(h_{j-1}+h_{j}\right)} u_{i, j+1}^{n}+\frac{h_{j}-h_{j-1}}{h_{j-1} h_{j}} u_{i j}^{n}-\frac{h_{j}}{h_{j-1}\left(h_{j-1}+h_{j}\right)} u_{i j-1,}^{n} \end{aligned}\)

\(\begin{aligned} \left(\frac{\partial^{2} u}{\partial x^{2}}\right)_{i j}^{n} &=\frac{2}{h_{i}\left(h_{i-1}+h_{i}\right)} u_{i+1, j}^{n}-\frac{2}{h_{i-1} h_{i}} u_{i j}^{n}+\frac{2}{h_{i-1}\left(h_{i-1}+h_{i}\right)} u_{i-1, j,}^{n} \\ \left(\frac{\partial^{2} u}{\partial y^{2}}\right)_{i j}^{n} &=\frac{2}{h_{j}\left(h_{j-1}+h_{j}\right)} u_{i, j+1}^{n}-\frac{2}{h_{j-1} h_{j}} u_{i j}^{n}+\frac{2}{h_{j-1}\left(h_{j-1}+h_{j}\right)} u_{i, j-1,}^{n} \end{aligned}\)

Rewriting Eq. (2.1) gives

\(u_{i j}^{n+1}=\Delta \tau\left[\left(\frac{\left(\sigma_{1} x_{i}\right)^{2}}{h_{i-1}\left(h_{i-1}+h_{i}\right)}-\frac{r x_{i} h_{i}}{h_{i-1}\left(h_{i-1}+h_{i}\right)}-\frac{\rho \sigma_{1} \sigma_{2} x_{i} y_{j}}{h_{i-1} h_{j-1}}\right) u_{i-1, j}^{n}\right.\)

\(+\left(\frac{\left(\sigma_{2} y_{j}\right)^{2}}{h_{j-1}\left(h_{j-1}+h_{j}\right)}-\frac{r y_{j} h_{j}}{h_{j-1}\left(h_{j-1}+h_{j}\right)}-\frac{\rho \sigma_{1} \sigma_{2} x_{i} y_{j}}{h_{i-1} h_{j-1}}\right) u_{i, j-1}^{n}\)

\(+\left(\frac{1}{\Delta \tau}-\frac{\left(\sigma_{1} x_{i}\right)^{2}}{h_{i-1} h_{i}}-\frac{\left(\sigma_{2} y_{j}\right)^{2}}{h_{j-1} h_{j}}+\frac{\rho \sigma_{1} \sigma_{2} x_{i} y_{j}}{h_{i-1} h_{j-1}}+\frac{h_{i}-h_{i-1}}{h_{i-1} h_{i}}+\frac{h_{j}-h_{j-1}}{h_{j-1} h_{j}}-r\right) u_{i j}^{n}\)

\(+\left(\frac{\left(\sigma_{1} x_{i}\right)^{2}}{h_{i}\left(h_{i-1}+h_{i}\right)}+\frac{r x_{i} h_{i-1}}{h_{i}\left(h_{i-1}+h_{i}\right)}\right) u_{i+1, j}^{n}\)

\(\left.+\left(\frac{\left(\sigma_{2} y_{j}\right)^{2}}{h_{j}\left(h_{j-1}+h_{j}\right)}+\frac{r y_{j} h_{j-1}}{h_{j}\left(h_{j-1}+h_{j}\right)}\right) u_{i, j+1}^{n}+\frac{\rho \sigma_{1} \sigma_{2} x_{i} y_{j}}{h_{i-1} h_{j-1}} u_{i-1, j-1}^{n}\right],\)       (2.2)

for i = 1, ..., Nx, j = 1, ..., Ny, and n = 0, ..., Nτ − 1,

where Nx, Ny, and Nτ = T /∆τ are the number of grid points in the x-, y-, and τ - direction, respectively. To make the coefficient of \(u^n_{ij}\) be positive, the following should be satisfied:

\(\frac{1}{\Delta \tau}>\frac{\left(\sigma_{1} x_{i}\right)^{2}}{h_{i-1} h_{i}}+\frac{\left(\sigma_{2} y_{j}\right)^{2}}{h_{j-1} h_{j}}-\frac{\rho \sigma_{1} \sigma_{2} x_{i} y_{j}}{h_{i-1} h_{j-1}}-\frac{h_{i}-h_{i-1}}{h_{i-1} h_{i}}-\frac{h_{j}-h_{j-1}}{h_{j-1} h_{j}}+r.\)       (2.3)

Suppose the step size is uniform, i.e., hi = hj = h for some constant h. Then, we have

\(\Delta \tau<f\left(x_{i}, y_{j}\right)=\frac{h^{2}}{\left(\sigma_{1} x_{i}\right)^{2}+\left(\sigma_{2} y_{j}\right)^{2}+r h^{2}-\rho \sigma_{1} \sigma_{2} x_{i} y_{j}}.\)

In Fig. 2, (a) and (b) show a mesh plot and contour lines of f(xi , yj), respectively, with h = 10, σ1 = σ2 = 0.3, r = 0.015, and ρ = 0.3. We can observe that the minimum value is obtained at the upper right domain corner.

 E1TAAE_2019_v23n1_19_f0003.png 이미지

FIGURE 2. (a) Mesh plot and (b) contour lines of f(xi, yj).

Let a temporary time step be

\(\Delta \tau_{t m p}=\frac{s h^{2}}{\left(\sigma_{1} x_{I}\right)^{2}+\left(\sigma_{2} y_{J}\right)^{2}+r h^{2}-\rho \sigma_{1} \sigma_{2} x_{I} y_{J}},\)

where s (0 < s < 1) is a safety factor and I, J are maximum indices of an interested region. Let Nτ = [T /∆τtmp] + 1 be the total number of time steps, where [k] means the largest integer not greater than k. We redefine the time step [24] as

\(\Delta \tau=T / N_{\tau}.\)       (2.4)

Then, from Eq. (2.3) with yj = xi and hj = hi , we have

\(h_{i}=\frac{1}{s} \times \frac{\left(\sigma_{1}^{2}+\sigma_{2}^{2}\right) x_{i}^{2} / h_{i-1}+2}{s / \Delta \tau+\rho \sigma_{1} \sigma_{2} x_{i}^{2} / h_{i-1}^{2}+2 / h_{i-1}-r}, \text { for } i>I,\)       (2.5)

where 1/s is a safety factor. Uniform grid (h = 4) is used up to an index I and then nonuniform grid (2.5) is used. Fig. 3 shows hi against index i with s = 0.99, σ1 = σ2 = 0.3, r = 0.015, ρ = 0.3, and ∆τ = 0.0076. Now, we describe the main process of the proposed algorithm. We set the grid size hi using uniform grid and Eq. (2.5). Then, we solve the Eq. (2.2) by using a hybrid boundary condition: a linear boundary condition is applied at x = 0 or y = 0, and the computation domain is reduced by one grid in each direction at the other boundaries. Therefore, as we proceed the computation, the domain size decreases as shown in Fig. 1.

 E1TAAE_2019_v23n1_19_f0002.png 이미지

FIGURE 3. Grid size hi.


We numerically compute the European call option and the cash-or-nothing option and then compare with each analytic solution to show the accuracy of our proposed method. We use the relative error to compare the results. The relative error is defined as follows:

error = Vanalytic − V,

where Vanalytic is the analytic solution and V is the numerical solution. All tests were run on a 2.7 GHZ Intel PC with 4 GB of RAM loaded with MATLAB 2016a [26].

3.1. European call option of two-asset. First, let us consider a European call option [27] whose payoff is given as

u(x, y, 0) = max{x − K1, y − K2, 0},       (3.1)

where K1 and K2 are the strike prices of x and y, respectively. Fig. 4(a) shows the payoff function (3.1). Fig. 4(b)–(d) show the snapshots of the option price at τ = 0.7634, 0.8397, and 1, respectively. Here, we used the following parameters: s = 0.99, σ1 = σ2 = 0.3, r = 0.015, ρ = 0.3, T = 1, K1 = K2 = 100, and h = 4. Using the above-mentioned time step (2.4) and grid (2.5), we solve the European call option as shown in Fig. 4. The analytic solution is given in Appendix. We compared with operator splitting method(OSM) with conventional linear boundary condition, the method we present in the paper can get a more good result, which is shown in the Fig. 5. As shown in the Fig. 5(a), the price graph calculated by the OSM and conventional linear boundary condition shows that the value decreases toward the end of the domain. The conventional linear boundary condition has difficulty that generates the bad results at boundary when the high correlation in multi-dimensional problem or derivative with nonlinear payoff. However, in Fig. 5(b), the price graph that calculated by proposed method has good result. In addition, we compute the numerical Greeks and compare the difference with the analytic Greeks. The variables used here are the same as those used to calculate the option price. The characters in Table 1 mean that: Delta (∆) is the derivative of the option price relative to the underlying asset, \(∂u\over ∂x\) and \(∂u\over ∂y\). Gamma (Γ) is the second derivative of the option price, \(∂^2u\over ∂x^2\)\(∂^2u\over ∂y^2\), and \(∂^2u\over ∂x∂y\). Theta (Θ) is the derivative of the option price over time, \(∂u\over ∂t\) Rho (R) is the derivative of the option price with respect to the risk-free interest rate,\(∂u\over ∂r\) Vega (V) is the derivative of option price with respect to volatility, \(∂u\over ∂σ_1\) and \(∂u\over ∂σ_2\). We use finite difference schemes for evaluating the Greeks [28].

 E1TAAE_2019_v23n1_19_f0004.png 이미지

FIGURE 4. (a) Payoff of a European call option on the maximum of two-asset. (b)–(d) Snapshots of European call option at three different times, τ = 0:7634, 0:8397, and 1, respectively.


 Table 1 lists the price and the Greeks of European call option and error at current price, (x, y) =(100, 100).TABLE 1. Price and the Greeks of European call option and error at current price, (x, y) = (100, 100).

E1TAAE_2019_v23n1_19_t0001.png 이미지

x\(\frac{\partial u}{\partial x} \approx \frac{u\left(x_{i+1}, y, \tau\right)-u\left(x_{i-1}, y, \tau\right)}{2 h},\)

y\(\frac{\partial u}{\partial y} \approx \frac{u\left(x, y_{j+1}, \tau\right)-u\left(x, y_{j-1}, \tau\right)}{2 h},\)

Γxx\(\frac{\partial^{2} u}{\partial x^{2}} \approx \frac{u\left(x_{i+1}, y, \tau\right)-2 u(x, y, \tau)+u\left(x_{i-1}, y, \tau\right)}{h^{2}},\)

Γyy\(\frac{\partial^{2} u}{\partial y^{2}} \approx \frac{u\left(x, y_{j+1}, \tau\right)-2 u(x, y, \tau)+u\left(x, y_{j-1}, \tau\right)}{h^{2}},\)

Γxy\(\frac{\partial^{2} u}{\partial x \partial y} \approx \frac{u\left(x_{i+1}, y_{j+1}, \tau\right)-u\left(x_{i-1}, y_{j+1}, \tau\right)-u\left(x_{i+1}, y_{j-1}, \tau\right)+u\left(x_{i-1}, y_{j-1}, \tau\right)}{4 h^{2}},\)

Θ = \(-\frac{\partial u}{\partial \tau} \approx-\frac{u(x, y, \tau+\Delta \tau)-u(x, y, \tau-\Delta \tau)}{2 \Delta \tau},\)

R = \(\frac{\partial u}{\partial r} \approx \frac{u(x, y, \tau ; r+\Delta r)-u(x, y, \tau ; r-\Delta r)}{2 \Delta r},\)

Vx\(\frac{\partial u}{\partial \sigma_{1}} \approx \frac{u\left(x, y, \tau ; \sigma_{1}+\Delta \sigma_{1}\right)-u\left(x, y, \tau ; \sigma_{1}-\Delta \sigma_{1}\right)}{2 \Delta \sigma_{1}},\)

Vy\(\frac{\partial u}{\partial \sigma_{2}} \approx \frac{u\left(x, y, \tau ; \sigma_{2}+\Delta \sigma_{2}\right)-u\left(x, y, \tau ; \sigma_{2}-\Delta \sigma_{2}\right)}{2 \Delta \sigma_{2}}.\)

Table 1 lists the price and the Greeks of European call option and error at current price, (x, y) = (100, 100).

3.2. Two-asset Cash-or-Nothing Option. Next, we consider the two-asset cash-or-nothing option [27]. The payoff of the option is given as

\(u(x, y, 0)=\left\{\begin{array}{ll} \text { Cash } & \text { if } x \geq K_{1} \text { and } y \geq K_{2} \\ 0 & \text { otherwise } \end{array}\right.\)       (3.2)

where K1 and K2 are the strike prices of x and y, respectively (see Fig. 6 (a)). The analytic solution is obtained from a closed-form solution, which is provided in Appendix. The parameters used are: σ1 = σ2 = 0.3, r = 0.015, ρ = 0.3, T = 1, Cash = 100, K1 = K2 = 100 and h = 4. Because cash-or-nothing option has a rapid evolution at the strike price point, we wrap the strick price point and solve Eq. (2.2) by using the proposed scheme. We calculate the cash-or-nothing option and its Greeks just as we did for the European call option. Fig. 6 shows the payoff function and snapshots of the solution at three different times. Fig. 7 show the prices of the options borrowed by the OSM and proposed methods, respectively. Unlike the Fig. 5, noticeable differences between Fig. 7(a) and Fig. 7(b) are not identifiable. The reason for this comes from the payoff structure of cash or nothing option. Table 2 lists the price and the Greeks of cash-or-nothing option and error at current price, (x, y) = (100, 100). Checking the results in Tables 1 and 2 confirm that the analytic solution and numerical solution values are similar in almost all cases. This is effective because you can get similar value to the analytic solution like the above results without solving linear boundary treatment problems.

 E1TAAE_2019_v23n1_19_f0005.png 이미지

FIGURE 5. (a) European call option price using implicit finite difference scheme (b) European call option price using proposed scheme.


 E1TAAE_2019_v23n1_19_f0006.png 이미지

FIGURE 6. (a) Payoff of a cash-or-nothing option of two-asset. (b)–(d) Snapshots of cash-or-nothing option at three different times, τ = 0:8451, 0:9014, and 1, respectively.


 E1TAAE_2019_v23n1_19_f0007.png 이미지

FIGURE 7. (a) European cash or nothing option price using implicit finite difference scheme (b) European cash or nothing option price using proposed scheme.


 TABLE 2. Price and the Greeks of cash-or-nothing option and error at current price, (x, y) = (100, 100).

E1TAAE_2019_v23n1_19_t0002.png 이미지


In this article, we proposed the accurate explicit hybrid finite difference method for the twodimensional BS equation. The proposed method can avoid the numerical boundary treatment problem of the cross derivatives at the domain corner. The main idea of the algorithm is reducing the computational domain by each direction so that we do not need boundary conditions. Numerical tests for the pricing and Greeks of the two-asset European options demonstrate the accuracy and efficiency of the proposed algorithm.



  1. J. Topper, Financial Engineering with Finite Elements, John Wiley & Sons., 2005.
  2. P. Wilmott, J. Dewynneand, and S. Howisson, Option Pricing: Mathematical Models and Computation, Oxford Press, Oxford, 1993.
  3. R. Seydel, Tools for Computational Finance, Springer-Verlag, Berlin, 2003.
  4. Y. Achdou and P. Olivier, Computational methods for option pricing, Society for Industrial and Applied Mathematics, 2005.
  5. M. Mrazek, J. Pospisil and T. Sobotka, On calibration of stochastic and fractional stochastic volatility models, European Journal of Operational Research, 254 (2016), 1036-1046.
  6. M. Alghalith, Pricing the American options using the Black& Scholes pricing formula, Physica A. 2018.
  7. N.Laskin, Valuing options in shot noise market, Physica A, 502 (2018), 518-533.
  8. L. Lin, Y. Li, and J. Wu, The pricing of European options on two underlying assets with delays, Physica A, 495 (2018), 143-151.
  9. F. Black and S. Myron, The pricing of options and corporate liabilities, Journal of political economy, 81 (1973), 637-654.
  10. A. Eckner, Computational techniques for basic affine models of portfolio credit risk, Journal of Computational Finance, 13 (2009) 63-102.
  11. A. Mortensen, Semi-Analytical Valuation of Basket Credit Derivatives in Intensity-Based Models, The journal of derivatives, 13 (2006), 8-26.
  12. P. Boyle, M. Broadie and P. Glasserman, Monte Carlo methods for security pricing, Journal of Economic Dynamics & Control, 21 (1997), 1267-1321.
  13. C.Z. Mooney, Monte Carlo Simulation, Sage, London, 116, 1997.
  14. C.J. Wang and M.Y. Kao, Optimal search for parameters in Monte Carlo simulation for derivative pricing, European Journal of Operational Research, 249 (2016), 683-690.
  15. G.N. Milstein, M. V. Tretyakov, Numerical analysis of Monte Carlo evaluation of Greeks by finite differences, Journal of Computational Finance, 8 (2005), 1-34.
  16. Z. Cen and L. Anbo, A robust and accurate finite difference method for a generalized Black-Scholes equation, Journal of Computational and Applied Mathematics, 2350 (2011), 3728-3733.
  17. D. Tavella and C. Randall, Pricing Financial Instruments: the Finite Difference Method, John Wiley & Sons. 2000.
  18. D.J. Duffy, Finite Difference methods in financial engineering: a Partial Differential Equation approach, John Wiley & Sons. 2006.
  19. J. Toivanen, A high-order front-tracking finite diffrence method for pricing American options under jumpdiffsion models, The Journal of Computational Finance, 13 (2010), 61-79.
  20. K.J. Hout and S. Foulon, ADI finite difference schemes for option pricing in the Heston model with correlation, International Journal of Numerical Analysis & Modeling, 7 (2010), 303-320.
  21. S.P. Zhu and W.T. Chen, A predictor-corrector scheme based on the ADI method for pricing American puts with stochastic volatility, Computers & Mathematics with Applications, 62 (2011), 1-26.
  22. S. Ikonen and J. Toivanen, Operator splitting methods for American option pricing, Applied mathematics letters, 17 (2004), 809-814.
  23. Y. chen and J.Ma, Numerical Contour Integral Methods for Free-Boundary Partial Differential Equations Arising in American Volatility Options Pricing, Discrete Dynamics in Nature and Society, 2018 (2018), 61-74.
  24. D. Jeong, M. Yoo and J. Kim, Finite Difference Method for the Black-Scholes Equation Without Boundary Conditions, Computational Economics, (2017), 1-12.
  25. Z. Zhou and H. Wu, Laplace Transform Method for Pricing American CEV Strangles Option with Two Free Boundaries, Discrete Dynamics in Nature and Society, 2018 (2018).
  26. Using MATLAB,, The MathWorks Inc.Natick, MA, 1998.
  27. E.G. Haug, The complete guide to option pricing formulas, McGraw-Hill Companies, 2007.
  28. S. Lee, Y. Li, Y. Choi, H. Hwang and J. Kim, Accurate and efficient computations for the Greeks of European multi-asset options, Journal of Korean Society for Industrial and Applied Mathmatics, 18 (2014), 61-74.