FINITE DIFFERENCE METHOD FOR THE TWO-DIMENSIONAL BLACK-SCHOLES EQUATION WITH A HYBRID BOUNDARY CONDITION

Abstract

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.

1. INTRODUCTION

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 σ₁, σ₂ 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, x_max) × (0, y_max).

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.

2. NUMERICAL SOLUTION

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

 

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)

where

\(\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, ..., N_x, j = 1, ..., N_y, and n = 0, ..., N_τ − 1,

where N_x, N_y, 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., h_i = h_j = 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(x_i , y_j), respectively, with h = 10, σ₁ = σ₂ = 0.3, r = 0.015, and ρ = 0.3. We can observe that the minimum value is obtained at the upper right domain corner.

 

FIGURE 2. (a) Mesh plot and (b) contour lines of f(x_i, y_j).

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 y_j = x_i and h_j = h_i , 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, σ₁ = σ₂ = 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 h_i 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.

 

FIGURE 3. Grid size h_i.

3. NUMERICAL EXPERIMENTS

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 = V_analytic − V,

where V_analytic 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 − K₁, y − K₂, 0},       (3.1)

where K₁ and K₂ 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, σ₁ = σ₂ = 0.3, r = 0.015, ρ = 0.3, T = 1, K₁ = K₂ = 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].

 

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).

∆_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},\)

V_x = \(\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}},\)

V_y = \(\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 K₁ and K₂ 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: σ₁ = σ₂ = 0.3, r = 0.015, ρ = 0.3, T = 1, Cash = 100, K₁ = K₂ = 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.

 

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

 

 

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.

 

 

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).

4. CONCLUSION

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.

APPENDIX