1. Introduction
An option is the right, but not the obligation, to buy or sell an underlying asset at a specific price on or before a certain date at a fixed strike price. European call (put) option and American call (put) option are two types of known options. European option can be exercised only on expiration date
T. For American option, exercise is permitted at any time
. It was shown in [
1] that option pricing can be modeled using a partial differential equation of second order. The resulting equation is called the Black–Scholes equation. If the coefficients are constant or independent of the spatial variable, it can be transformed into the standard heat equation and can be solved exactly. However, in the financial market, the coefficients can depend on the time and the asset price. In such situation the analytical solution of the generalized Black–Scholes model is not available. Therefore, numerical simulations are essential to obtain information about the behavior of the solutions. It is important to develop numerical methods that reproduce the qualitative properties of the solution.
Here, we consider the generalized Black–Scholes model as given in [
2]:
with the following conditions:
Equation (
1) is a partial differential equation used for pricing options. The function
is the European call option price with asset price
x and time
t.
T is the maturity,
K is the price of strike,
is the risk-free interest rate, and
is the volatility function of the underlying asset and time. When the coefficients are constant functions, we get the classical Black–Scholes model. Following [
2] we consider the above model on a truncated domain
, where
is a suitably chosen positive number so that
The result about existence and uniqueness of a solution of (
4)–(
7) when
and
are constant functions can be found in [
3,
4]. It is proved in [
5] that if
and
are solutions of (
1)–(
3) and (
4)–(
7), respectively, then
The terminal condition (
5) is not a smooth one. Therefore, to guarantee the convergence of the numerical solution, we substitute the terminal condition by a four order smooth function (see in [
2]) and obtain the following model:
where the function
is defined as follows:
for a sufficiently small
,
satisfying
By using these conditions one can easily find that
The proof of the existence and uniqueness of a solution of (
9)–(
12) can be found in [
3]. A proof of the following result can also be found there:
Proposition 1. There exists a positive constant P independent of such thatfor . From (
8) and (
16) we have that we can take the solution of the modified models (
9)–(
12) closer to that of the original models (
1)–(
3) by choosing a sufficiently large value for
and a sufficiently small value for
.
Traditionally, important requirements for constructing a numerical method are the investigation of the consistency of the discrete scheme with the original differential equation and linear stability analysis with smooth solutions [
6,
7,
8,
9,
10]. These requirements are important, because they are necessary to ensure the convergence of the discrete solution to the exact one, but for the explicit numerical methods based on standard finite difference approaches. However, most of the time, the qualitative properties of the true solution are not reproduced by the numerical solution. Therefore, this might result in a calamitous erroneous outcome. One way to overcome this disadvantage is to use nonstandard finite difference methods (NSFDs). The origin of the NSFD methods began with Mickens in 1989 [
11]. A summary of the known results up to 1994 are given in [
12]. The simplicity in construction and the possiblity of their extensions allowed researchers to apply NSFD schemes to solve several complex differential equation models that arise at different disciplines. A scheme is called nonstandard if at least one of the following conditions is satisfied:
- 1.
The usual denominator
h or
of the discrete first or second order derivatives is replaced by a non-negative function
such that
Some examples of
satisfying these conditions are
- 2.
The nonlinear terms that appear in the differential equation are approximated in a non-local way, for instance, the nonlinear terms
and
can be replaced by (see [
13]):
In addition to the usual properties of consistency, stability and convergence, the NSDF methods give rise to numerical solutions that also maintain essential qualitative properties of the solutions [
14,
15,
16,
17,
18,
19,
20,
21].
In this paper, we introduce a new scheme for the numerical solution of the generalized Black–Scholes equation. We propose a non-standard finite difference method, where we use a non-local approximation of the reaction term of the generalized Black–Scholes equation, combined with an implicit time step technique. The method is conditionally stable, positivity preserving and of second-order of convergence with respect to the space variable. Our proposed method has some advantages: A uniform mesh is used; thus, the analysis is simple and it has better accuracy and convergence rate.
The rest of the paper is organized as follows. In
Section 2, we review the works done on the generalized Black–Scholes equation and the real options. In
Section 3, we propose a NSFD scheme for such a model and we present the analysis of the method concerning the positivity preserving property, Stability, and truncation error. Furthermore, the numerical results obtained from the new method are presented and compared with the results in [
2]. The last section describes the features of the new method, the results obtained, and ends with the presentation of ideas for future works.
2. Literature Review
In this section, we review the works done on the generalized Black–Scholes equation. In addition, at the end of the section, we will mention some of the work done related to real options. There are several methods for solving the generalized Black–Scholes equation in the literature. Cen and Le [
2] presented an accurate finite difference method, based on a central difference discretization in the space variable using a piecewise uniform mesh and an implicit method in the temporal variable. Cen et al. [
22] proposed an exponential scheme in the temporal variable combined with a central difference approach on a spatial piecewise uniform mesh. Kadalbajoo et al. [
23] suggested a cubic B-spline collocation method for the generalized model with second order accuracy in both space and time. A fitted finite volume method was used in [
24] for the generalized model, showing a second-order convergence. A Cubic spline method was proposed in [
25] for a generalized Black–Scholes equation. In that work, the implicit Euler method was used in the temporal discretization, together with a cubic polynomial spline approximation for the spatial discretization. In [
26], Mohammadi achieved third order of accuracy in space considering a quintic B-spline collocation method and first- and second-order of accuracy in time using the backward Euler method and the Crank–Nicolson method, respectively. In [
27], a high-order numerical method was obtained. In this method a two-step backward differentiation formula is used in the temporal discretization while a high-order difference approximation was used in the spatial discretization. This gives a second order of accuracy in time and third order of accuracy in space.
In [
28], Clevenhaus et al. added a penalty term to the American option problem, by doing this the problem was reduced to a PDE on a fixed domain. Their presented penalty term is bounded by an initial guess of the free boundary value and is thus independent of the solution at the last time step. The novelty and advantage of their approach consist in introducing a time-dependent penalty function which led to the construction of an efficient adaptive numerical scheme. In [
29] a new method for pricing American options by Boire et al. is presented. The numerical results show that their method successfully reduces the bias plaguing the standard importance sampling method across a wide range of moneyness and maturities, with negligible change to estimator variance.
In [
30], considering several uncertainty factors: a renewable energy production amount, an R and D investment amount, a unit price, a risk-free interest rate, and the R and D investment in the six types of renewable energy sources: biomass energy, waste energy, photovoltaic energy, solar thermal energy, marine energy, and wind power energy are evaluated. They utilize the system dynamics approach using the Black–Scholes model and conducts sensitivity analyses of the uncertainty factors. Larissa et al. [
31] investigate the contribution of adjusted net savings to sustainable economic growth for 10 Central and Eastern European and Baltic nations using data corresponding to the period 2005–2016. Their results indicated that adjusted net savings impacted on the GDP across the 10 countries analyzed. A bibliometric analysis to examine the features of Carbon capture and storage (CCS) literature including the research focus and trends as well the real options uncertainty and models, types of options, and valuation techniques in [
32] are presented. Their results provide energy and environmental policy-makers and CCS project planners with valuable insights into various aspects of CCS policy and project design. For real options in the field of agriculture, we can refer to the article [
33] in which the authors investigated a hydroponic farm considering uncertainty, sustainability, and the system’s utility in the dominant desert geography.
4. Conclusions
In this paper, the generalized Black–Scholes equation for European option pricing is studied. To numerically solve this equation, a non-standard finite difference method has been used. In the proposed method, a non-local expression to approximate the reaction sentence of the generalized Black–Scholes and an implicit time step technique is used.
The theory results show that the proposed scheme is positivity preserving, stable and the order of the scheme with respect to the space variable is two. To confirm the theory results and good performance of the new method the option price for different values of volatility and interest rates is calculated. The numerical results confirm the theoretical behavior regarding the order of convergence for Theorem 4 (see
Table 2,
Table 3,
Table 4,
Table 5,
Table 6,
Table 7 and
Table 8). For the last two examples, the option price values are plotted for all time levels, which shows that the new method is positive and non-oscillation (see
Figure 1 and
Figure 2). Furthermore, the obtained results are more accurate than other existing results in literature.
The results obtained indicate that non-standard difference schemes may be promising for solving problems which may have an impact on the stock price such as nonlinear Black–Scholes equations. In our subsequent investigation, we intend to address these problems using non-standard methods.