An accurate solution for the generalized Black-Scholes equations governing option pricing

1.   Introduction

Pricing of options is an amplified area of discussion among financial practitioners. An option is a bond between two parties in which the option buyer buys the right, not the obligation to buy or sell an underlying asset at a prefixed strike price from or to the option writer within a fixed period. According to the option rights, options are classified into Call and Put options. An option that brings the owner the right to buy at a specific price is known as a call; an option that brings the right of the owner to sell at a particular price is known as a put. Option styles are classified into American and European options. American options can be exercised at any time up to and including the expiry. European options can only be exercised on the day of expiration. Fischer Black and Myron Scholes ^[1,2] have given a mathematics model under the assumptions:

1. The change in stock price $dS$ of the underlying satisfies the stochastic differential equation

where, $\mu$ is the drift rate, $D$ is the dividend yield, $\sigma$ is the market volatility and $dW$ is the increment of a standard Wiener process.

2. The risk-free rate of return $r$, the drift $\mu,$ dividend yield $D$, and the market volatility $\sigma$ are constants.

3. The market is arbitrage-free and frictionless.

In effect, the market should be complete in the sense that any financial derivative or commodity can be hedged with a portfolio of another commodity. Now using Ito's Lemma ^[3], we have

and by eliminating the market randomness, one can derive the celebrated Black-Scholes partial differential equation as

with the terminal condition

here, $T$ is the expiry and $E$ is the stike price of the commodity in the option.

The analytical solution for $(1.1)$ in a closed form ^[1,2,4], can be obtained as

which is also known as the Black-Scholes formula for European options, where

and $N(x)$ is the cumulative standard normal distribution function given as

But in the present scenario of the financial market, the parameters $\sigma, \; r$ and $D$ depend highly on the asset price $S$ and the time $\tau$. The analytical solution of generalized Black-Scholes model

with the same terminal condition, is often not available. The generalised Black-Scholes model (1.3) is numerically solved in ^{[4,5,6,7,8,9,10,11,12,13,14]}. Cubic B-spline collocation method ^[5,6] have second-order accuracy in approximating the generalized Black-Scholes model. In ^[14], cubic polynomial spline method in space direction after application of implicit Euler method in the time direction, producing second-order accuracy in space. A simultaneous application of the HODIE ^[15] and backward differentiation formulas ^[16,17,18] result in a second-order converging scheme ^[13] for generalized Black-Scholes equation. Besides of these models, fractal behavior of a stochastic process inflames the fractional Ito calculus for stochastic models and financial theories like time fractional Black-Scholes model ^[19,20]. Also numerical methods to solve these models are given in Zheng, Liu, Turner ^[21], R. D Staelen, Hendy ^[22].

In this paper, we fixated on method of lines (MOL) to the generalised linear Black-Scholes model for European call option, firstly by elementary finite difference schemes in space direction, and the generalized trapezoidal formulas (GTF $[\alpha = \frac{1}{3}])$ introduced by Chawla et al. ^{[23,24,25,26]}, in the temporal discretization for the system of ordinary differential equations. Section 2 of this article portrays the terminal value problem of linear Black- Scholes model. In Section 3, semi-discretization of the parabolic partial differential equation along with initial and the artificial boundary conditions is done, and the numerical scheme is derived. Section 4 deals with convergence and stability analysis of the numerical scheme. Numerical experimentations and error comparison with existing numerical schemes are given in Section 5, and Section 6 concludes the paper.

2.   Model equation

Let $r (S, t)$, $D (S, t)$, and $\sigma(S, t)$ be sufficiently smooth and bounded functions on the domain ($(0, \infty) \times (0, T)$). Consider the generalized Black-Scholes differential equation ^[1,3] for European call option

here, $V(S, t)$ is the value of the European call option at the the stock price $S$ (spatial variable) and at time $t$. with $V (0, t) = V_0(t), V(S, t) \; \sim\; V_{\infty}(t) \; \mbox{as}\; S\; \mbox{tends to}\; \infty$ and $V (S, T) = V_T(S)$. We proceed with the often case $V_0(t) = 0, V_{\infty}(t) = S$ and $V_T(S) = \max\{S-E, 0\}$. Here $\sigma$ denotes a statistical measure of the volatility of the underlying commodity, $E$, the exercise or striking price, $T$, the expiry time, $D$, the dividend pay and $r$, the risk-free rate of return. The parameters $'r'$, $'D'$, and $' \sigma '$ are constant functions in the case of classical Black-Scholes equation (1.1). The existence and uniqueness of a classical solution of (2.1) is well known ^{[27,28,29,30,31]}. Now, it can be seen that the above model is derived in an infinite domain $(0, \infty) \times (0, T)$, which makes difficulties in composing the numerical solutions. Thus we are insisted to consider the following model defined on a truncated domain $(0, S_{max}) \times (0, T)$, where $S_{max}$ is the suitably chosen positive number.

The existence and uniqueness of analytical solution of (2.2) can be found in ^{[27,28,29,30,31]}. Here, the boundary conditions are chosen according to ^[32]. Moreover, it is proved in ^[33] that if $V$ and $\mathscr V$ are solutions of (2.1) and (2.2) respectively, then at every point $(S, t) \in (0, S_{\max}) \times [0, T]$ satisfying

we have

where $d = \inf \left\{2 D(S, t)-2 r(S, t)+\sigma^{2}(S, t) :(S, t) \in\left(0, S_{\max }\right) \times(0, T)\right\}$ and $\Lambda = \left\{0, S_{\max }\right\}$.

Since the pay-off is not differentiable at the striking price, the resulting solution is not differentiable for the convergence of numerical approximations ^[34]. We replace $\max \{ S-E, 0 \}$ in the terminal condition by a smooth function '$\phi(S) = \varphi(S-E)$ ^[35] defined as

where $\epsilon > 0$, a small constant and $c_{i}, i = 0, 1, \ldots, 9$ are the constant coefficients to be determined.

Applying the following conditions on the funtion $\varphi(x) :$

we can uniquely determine the unknown coefficients $c_{i}, i = 0, 1, \ldots, 9$ viz.:

Figure 1 demonstrates the smoothening procedure of the pay-off terminal condition (European call option) by the function $\phi$, wherein the value of $\epsilon$ is taken as $0.5$ for the better view.

Thus we obtain another prototype in which $\mathscr W$ is treated as the worth or value of option,

with final condition

and boundary conditions

The existence and uniqueness of the analytical solution of (2.3) can be seen in ^[27], which also contains the proof of the following estimate:

There exists a positive constant '$\mathcal K$ ' independent of '$\phi (S)$ ' such that

To dispose of the degeneracy and backwardness of (2.3), we set $\tau = T-t, S = e^x$ and $\mathscr W(S, t) = u(x, \tau)$ to get

with

and

3.   Numerical methods

A numerical approach to the generalized Black-Scholes equation becomes sensible in this transformed settings (2.4). Here, we execute a method of vertical lines (MOL) on (2.4) to get a system of ODEs. Thenceforth we employ the generalized trapezoidal formulas (GTF($\frac{1}{3}$)) as a numerical time integration.

3.1.   Semi-discretization by Method of Lines (MOL)

Let $M+1$ be the number of price grid points $x_i = x_{\min}+ih, \; i = 0, 1, ..., M$, where ${h = \frac{x_{\max}-x_{\min}}{M}}$. For a positive integer $N$, Define temporal grid $\tau_j = jk, \; j = 0, 1, ..., N$, where ${k = \frac{T}{N}}$. Now let $u_{i, j} = u(x_i, \tau_j)$ and discretize spacial derivatives using classical central differences,

Suppose $\vec{U}(\tau) = (u_1(\tau), u_2(\tau), ..., u_{M-1}(\tau))'$ and $\vec{U}_M(\tau) = (0, 0, ..., u_M(\tau))_{M-1 \times 1}$, then $(3.1)$ can be expressed as

where

and

Now the system $(3.2)$ with the initial condition ${\vec{U}(0) = (\phi(x_1), \phi(x_2), ..., \phi(x_{M-1}))'}$ is an IVP in $\tau$.

3.2.   Full-discretization by Generalized Trapezoidal Formulas (GTF)

Let $F(\tau, \vec{U})$ be the right hand side of $(3.2)$, we apply the generalized trapezoidal formulas ^{[23,24,25,26]} for the time integration of $(3.2)$ to obtain,

where, ${\vec{U}_{j} = \vec{U}(\tau_j)}, \; F_{j} = F(\tau_j, \vec{U}_j)$ and $\hat{F}_{j} = F(\tau_j, \vec{U}_{j+1}-kF_{j+1})$, now by rearranging we have the generalized trapezoidal formulas for the generalised Black-Scholes equation,

with $\vec{U}_{0} = (\phi_1, \phi_2, ..., \phi_{M-1}), \; \vec{U}_{j}(0) = u^j(x_{\min}) = u(x_{\min}, \tau_j)$ and $\vec{U}_{j}(M) = u^j(x_{\max}) = u(x_{\max}, \tau_j)$, where,

and

with

4.   Error analysis

Lemma 4.1. Let $\hat{\sigma}, \; \hat{D}, \; \hat{r}$ be the parameters defined in (2.4), Assume that the spacial and temporal grid sizes $h$ and $k$ satisfies the conditions,

Then

for some $\alpha > 1$, where, $\mathbf{A}_{j}$ is the matrix given by (3.3).

Proof. By the assumptions on $h$ and $k$, the matrix $\mathbf{A}_{j}$ is a pentadiagonal diagonally dominant matrix with the minimum dominance

where $\mathbf{A}_{j} = (a_{kl}^j)$ and the $l_\infty$ bound of inverse is in agreement with Varah ^[36,37].

Lemma 4.2. Assume the conditions in lemma (4.1) for $h$ and $k$, the operator $L_h^k$ defined by

where $\mathbf{A}_{j}(i), \; \mathbf{F}_{j}(i)$ are the $i$ th rows of $\mathbf{A}_{j}, \; \mathbf{F}_{j}$ respectively, satisfies the consistency estimate

for some constant $\mathcal C > 0$.

Proof. The derivative ${\frac{\partial{u}}{\partial{\tau}}\big|_{x = x_i}}$ is given as

where, it can be seen that (By Taylor's expansion)

To calculate error of the time integration scheme, we have

and

Use Taylor's expansion for $u^{j+1}, \; \frac{\partial u}{\partial \tau}\big|^{j+1}$ to give,

Now, to calculate error of the scheme, we write (4.4) in the form

Then, an application of $(4.7)$ gives

Now under the assumption that the problem 2.4 satisfies sufficient regularity and compatibility conditions ^[27,30], we have ${\big|\frac{\partial^{m+n} u}{\partial x^m \partial \tau^n}\big|\leq \mathcal C'}$, for $0\leq n \leq 3$ and $0\leq m+n \leq 5$. Further, from the continuous problem (2.4), it can be seen that

This gives

Theorem 4.3. Let $u$ be the solution of the problem (2.4) and $\mathcal U_i^j$ be the solution for the discrete problem (4.3), then

for some $C > 0$.

Proof. The lemmas 4.1 and 4.2 says that the discretisation 4.3 is stable with

for some constant $\mathcal C > 0$. Together with the uniform bound of $\mathbf{A}_{j}^{-1}$, we obtain

for some $C > 0$ ^[38,39,40].

5.   Numerical illustrations

First we illustrate an example on the transformed Black-Scholes model (2.4) for which the closed form solution is available and is given by

where

Let $\mathcal U_{i, j}^{M, N}$ be the numerical approximation with $M$ and $N$ points in space and time directions respectively, Compute the true $L_{\infty}$ norm error (maximum absolute error, $e_{\max}^{M, N}$), $L_{2}$ norm error (root mean square error, $e_{rms }^{M, N}$) and corresponding orders of convergence $p_{max }^{M, N}$ and $p_{rms }^{M, N}$ as follows:

and

Example 1. Here we consider the Black-Scholes equation (2.4) for European call option with $\hat{\sigma} = 0.4, \hat{r} = 0.06$ $\hat{D} = 0.02, E = 1$ and $T = 1.$ For computational purpose, we assume that $x_{\min } = -2, \; x_{\max } = +2 \operatorname{and} \epsilon = 10^{-6}.$ The maximum and rms errors and their numerical order of convergence is given in Table 1. Table 2 gives the the maximum norm error and their numerical order of convergence for HODIE scheme ^[13] for the Black-Scholes equation ($\hat{\sigma} = 0.4, \hat{r} = 0.04$ $\hat{D} = 0.02, E = 1$ and $T = 1$). The solution profile is given in Figure 2.

Due to the unavailability of exact solution data of the generalized Black-Scholes equation, we are supposed to use the double mesh principle for computing the root mean square error ($e_{rms}^{M, N}$), sup norm error ($e_{\max}^{M, N}$) and corresponding orders of convergence $p_{rms}^{M, N}, \; p_{\max}^{M, N}$ and it is given by

and

In this example, we have illustrated the errors on the line $t = t_N$ which is a significant subdomain for the problem (2.4) in which the 'option premium' (Initial option price) is approximated by the generalised trapezoidal formulas.

Example 2. Consider the generalized Black-Scholes equation for European call option price (2.4) with $\hat{\sigma}(x, t) = 0.4(2+$ $t \sin (\exp(x)))$, $\hat{r}(x, t) = 0.06(1+(T-t)\exp (-\exp(x))), \; \hat{D}(x, t) = 0.02 \exp (-t-\exp(x)), T = 1$ and $E = 1.$ Assume that $x_{\text {min }} = -2, \; x_{\text {max }} = 2$ and $\epsilon = 10^{-6}.$ The numerics and the GTF solutions are displayed in Table 3 and Figure 3 respectively. Table 4 gives the maximum norm error and their numerical order of convergence for cubic B-spline collocation scheme combaining $\theta-$method ^[6].

Example 3. Consider the generalized Black-Scholes equation (2.4) for Binary European call option with $\hat{\sigma}(x, t) = 0.4(2+$ $t \sin (\exp(x)))$, $\hat{r}(x, t) = 0.06(1+(T-t)\exp (-\exp(x))), \; \hat{D}(x, t) = 0.02 \exp (-t-\exp(x)), T = 1$ and $E = 1.$ Assume that $x_{\text {min }} = -2, \; x_{\text {max }} = 2$ and $\epsilon = 10^{-6}.$ The numerics and the GTF solutions are displayed in Table 5 and Figure 4 respectively.

Here the initial condition and boundary conditions are as follows

Now we replace the initial condition $u_0(x) = u(x, 0)$ by a smooth function '$\phi(x) = \varphi(e^x-E)$ ^[35] defined as

where $\epsilon > 0$ is a small constant and $c_{i}, i = 0, 1, \ldots, 9$ are the constant coefficients to be determined.

Applying the following ten conditions on the funtion $\varphi(y) :$,

we can uniquely determine the unknown coefficients $c_{i}, i = 0, 1, \ldots, 9$ viz.:

Figure 5 shows the smoothening procedure of the pay-off in binary European call, wherein the value of $\epsilon$ taken as $0.4$ for the better view.

Example 4. Consider the generalized Black-Scholes equation 2.4 for Butter fly spread option with $\hat{\sigma}(x, t) = 0.4(2+$ $t \sin (\exp(x)))$, $\hat{r}(x, t) = 0.06(1+(T-t)\exp (-\exp(x))), \; \hat{D}(x, t) = 0.02 \exp (-t-\exp(x)), T = 1$ and three singular points $E_1 = 1, \; E_2 = 2, \; E_3 = 3.$ Assume that $x_{\text {min }} = -2, \; x_{\text {max }} = 2$ and $\epsilon = 10^{-6}.$ The numerics and the GTF solutions are displayed in Table 6 and Figure 6 respectively.

Here the initial condition and boundary conditions are as follows

The payoff for butterfly option has three singularities $E_1, \; E_2$ and $E_3$, So we replace the initial condition $u_0(x) = u(x, 0)$ by a smooth function '$\phi(x)$ ^[35] defined as follows:

where $\varphi_{1}\left(x\right) = \sum_{i = 0}^{9} c_{i}x^i$ and the coefficients $c_{i}, i = 0, 1, \ldots, 9$ are computed by solving the following ten conditions:

$\varphi_{2}\left(x\right) = \sum_{1 = 0}^{9} d_{i}x^i$ and the coefficients $d_{i}, i = 0, 1, \ldots, 9$ are computed by solving the following ten conditions:

and $\varphi_{3}\left(x\right) = \sum_{1 = 0}^{9} e_{i}x^i$ and the coefficients $e_{i}, i = 0, 1, \ldots, 9$ are computed by solving the following ten condi-

The coefficients $c_{i}, d_{i}$ and $e_{i}, i = 0, 1, \ldots, 9$ take the following values:

and

Figure 7 shows the smoothening procedure of the pay-off in butter fly spread option, wherein the value of $\epsilon$ taken as $0.5$ for the better view.

6.   Conclusions

In this paper, we have applied a computationally convenient time integration scheme for the generalized Black-Scholes equations. The method is based on a central difference spatial discretization on uniform mesh and the generalized trapezoidal formulas (GTF($\frac{1}{3})$) in time-stepping. The inverse of the matrix associated with the discrete operator is uniformly bounded by the inverse of minimum diagonal dominance and thereby stable for arbitrary volatility and interest rate. The proposed scheme is second-order consistent concerning the spatial variable and third-order in time, and by accepting the uniform bound, we obtain the convergence in the same order as the consistency. Furthermore, it can be seen that GTF $(\frac{1}{3})$ rectifies the singularities of the non-smooth payoff function by approximation with a ninth-degree polynomial. Numerical experiments, including the smoothening of payoffs and comparison with existing literature, are performed to demonstrate the efficiency of the proposed scheme.

Acknowledgments

The authors gratefully acknowledge the anonymous referees for their constructive inputs, valuable comments, and suggestions. The authors also would like to thank The Council of Scientific and Industrial Research (CSIR), India, for the financial aid provided for this research.

Conflict of interest

The authors declare no conflict of interest.