A robust numerical solution to a time-fractional Black–Scholes equation

Dividend paying European stock options are modeled using a time-fractional Black–Scholes (tfBS) partial differential equation (PDE). The underlying fractional stochastic dynamics explored in this work are appropriate for capturing market fluctuations in which random fractional white noise has the potential to accurately estimate European put option premiums while providing a good numerical convergence. The aim of this paper is two fold: firstly, to construct a time-fractional (tfBS) PDE for pricing European options on continuous dividend paying stocks, and, secondly, to propose an implicit finite difference method for solving the constructed tfBS PDE. Through rigorous mathematical analysis it is established that the implicit finite difference scheme is unconditionally stable. To support these theoretical observations, two numerical examples are presented under the proposed fractional framework. Results indicate that the tfBS and its proposed numerical method are very effective mathematical tools for pricing European options.


Introduction
Since the discovery of the most celebrated Black-Scholes-Merton asset pricing formula in the early 1970s, the application of Black-Scholes (BS) partial differential equations (PDEs) in valuation of derivative instruments has become very popular. The popularity of the approach is attributed to the extensively established evidence that it provides an effective asset valuation tool.
In derivative markets, once the price evolution process for a particular asset is specified, it is possible to address the question of how to price derivative contracts on this particular asset. One of the most common derivative instruments in the markets is called an option. It is a financial contract that gives the holder the right, but not obligation, to buy or sell a specified quantity of an underlying asset, e.g., a stock, at a fixed price, called the strike price, before or on the expiration date. Since it is a right and not an obligation, the holder of the contract can decide not to exercise the option and let it expire worthless. 1 The research herein was conducted while the first author was a PhD student at the University of the Western Cape.
© The Author(s) 2021. This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
There are two main types of options, standard and nonstandard. Standard options are further categorized as European or American type. The European options can only be exercised at the expiration date, whereas the American options are more flexible and can be exercised anytime at or before the expiration date. These options can further be divided into two categories, call options (which give their holders the right to buy) and put options (which give their holders the right to sell).
The study of stock markets dynamics via the classical approach is based on the wellknown efficient market hypothesis, i.e., martingale property of price movements. One of the consequences of using classical asset pricing models is that it is almost impossible to infer any information from historical price movements to predict the future. It is therefore imperative that, for one to better predict the future asset price movements, there is a need to pay attention to the repeated patterns and trends observed in financial data. The repeated patterns and trends in financial markets are mainly attributed to the following facts: (i) the markets discount everything, (ii) price moves in trends depict investors psychology (among other things), and (iii) history would mostly always repeat itself. There is empirical evidence to this effect, see, for example [17,33,40] and references therein.
The design of option valuation models is centered around the assumption that the market consists of a risky asset, for example, a stock, and a riskless asset, say a bond or a bank deposit. Emanating from its assumptions, the classical BS approach suffers from a few drawbacks, namely, (i) constant rates of return, (ii) constant volatility, and (iii) no dividends, taxes, or transactional costs. However, in real life, interest rates are bound to market forces and as such cannot remain static over a longer period of time. Based on empirical evidence, most traded assets' returns exhibit memory structures, see, for example, [17,33,40] and references therein, regimes of uneven fluctuations [6,25,37], volatility smiles and clustering [5,29,37]. Weakening some of assumptions of the classical BS approach, researchers suggested a few improved models, such as stochastic interest rate models [12,16], jump-diffusion models [14,25,26,30,46], stochastic volatility models [6,36,37,39], models with transactional costs [13,19,44,45], as well as those with dividend payments [4,29,36,37].
Bielecki et al. [5] extended the no-arbitrage pricing theory to dividend-paying securities in discrete-time markets with transaction costs. They proved that when there are no transaction costs on the dividends paid by the asset, the no-arbitrage condition under the efficient friction assumption is proved to be equivalent to the existence of a consistent pricing system. In the general case, when there are transaction costs on the dividends, the no-arbitrage condition is open. This is so because, if an asset pays dividends, the asset price falls by the amount of the dividend. No arbitrage opportunity arises because investors are compensated with the same amount of price depreciation back in cash through dividends.
Recent pool of models with dividend payments includes, but is not limited to, [4,29,37] and some references therein. Martin-Vaquero et al. [29] derived a stabilized explicit Runge-Kutta method for computing the value of an American option on multiassets with dividends. Rana and Ahmad [37] used the Glosten-Jagannathan-Runkle Generalized Autoregressive Conditional Heteroscedasticity (GJR-GARCH) forecasted volatility in pricing European call options on dividend paying stocks. Ballerster et al. [4] derived a robust numerical method for pricing vanilla options with discrete dividend payments.
Another notable setback of the classical BS approach and many of the revised versions discussed above is that their PDEs involve integer order derivatives. According to Panas [33], integer order derivatives only capture localized information (change) around a point. However, with changing market conditions which led to the evolution of a number of unusual structures in financial systems, such as repeated patterns and trends, heavy tailedness in the distributions of asset returns, volatility smiles and clustering, present challenges to the use of models involving local derivatives. The need for better approaches is thus very imperative.
The issue of memory in financial data has long been observed in a number of financial systems, such as stock prices, gross domestic products (GDPs), interest rates, foreign exchange rates, equity prices, etc., see [17,33,[40][41][42] and references therein. Garzareli et al. [17] proved the existence of memory effects in the stock price series using a conditional probability approach. Conventionally, the memory effect was measured by autocorrelation functions, and this trend continued until recently when the fractal structures of financial systems were discovered. The Hurst parameter and fractional derivative operators in models driven by fractal processes are increasingly becoming some of the most effective tools for capturing the effects of memory in financial systems. To capture memory effects under a classical BS setup, some researchers (see, for example, [21,27,28]) suggested replacing the standard Brownian motion in the stock dynamic equation with a fractional Brownian motion with Hurst parameter H ∈ [0, 1].
There has been a widespread application of fractional Brownian motions based models in finance. Edeki et al. [15], Jena et al. [20], Jumarie [21], Liang et al. [27], Mandelbrot and Cioczek-Georges [28], and Wei-Gou et al. [40], just to mention a few, are some of the authors who applied fractional calculus based models in pricing equities, wallets, options, and to general portfolio optimizations problems. Fractional Brownian motions based models have two very important properties, namely, self-similarity and long-range dependence (heredity properties), see [33,40] and references therein. These properties allows for the capture and representation of the effects of extreme asset price movements. They also help in explaining the effects of repeated patterns and trends observed in asset price movements. Jumarie [24] pointed out that stock price volatility can be well captured by a fractional Brownian motion which presents some random-like features suitable in explaining the effects of uneven fluctuations in stock price movements.
It is important to note that memory effects come in two forms, the noise memory effects and trend memory effects. Incorporating a fractional Brownian motion as the underlying process for the pricing dynamics only captures the noise memory effects [33]. Based on collective arguments in [7,32,34], fractional derivative based models are very good mathematical tools for explaining dynamics of complex processes, irregular increments, and trend memory effects which are exhibited by a number of financial instruments. To describe the trend memory effects, one replaces the integer order derivatives with their corresponding (nonlocal) derivatives of fractional orders. The resultant PDEs are often referred to as fractional partial differential equations (FPDEs). In the Black-Scholes setup, there are three main classes of FPDEs, namely, the time-fractional Black-Scholes (tfBS) PDEs, the space-fractional Black-Scholes (sfBS) PDEs, and the time-spacefractional Black-Scholes (tsfBS) PDEs. With the tfBS-PDEs, the time derivative is replaced by a corresponding fractional derivative of order α (0 < α ≤ 1), whereas in the case of sfBS-PDEs the space derivatives are replaced by their corresponding fractional derivatives of order α and β (1 < β ≤ 2). In the case of tsfBS-PDEs, one has a combination of the two.
Assuming the stock price dynamics are driven by a non-Gaussian fractional dynamic process, Jumarie [23] derived two sets of fractional BS equations, the time-fractional (tfBS) and the time-space-fractional (tsfBS) PDE. Jumarie's approach has been adapted by a number of researchers, some designed numerical schemes for solving the FPDEs and others applied the PDEs to pricing of different kind of options, see, for example, [1,2,27,29,31].
Chen et al. [10] followed an approach almost similar to that of Jumarie [23], but instead assumed that the stock price dynamics follows a standard Brownian motion as in the classical BS setup, while analogously replacing the integer derivatives with their corresponding fractional-order derivatives. Chen et al. [10] demonstrated that, using either their approach or that of Jumarie [23], one obtains similar results.
Recently, Edeki et al. [15] and Jena [20] proposed an analytical and a semianalytical solution to fractional option pricing models. In [15] an analytic solution to a time-spacefractional Black-Scholes model is proposed using the Jumarie fractional derivative operator. The approach has been proven to be consistent with actual integer and fractional data when compared with the classical Black-Scholes model. In [20] a novel semianalytical technique for solving a Schrödinger-type option pricing time-fractional model governed by a controlled Brownian motion was proposed. Their results are well in agreement with other already established numerical, analytical, and semianalytical results.
Other recent developments in fractional calculus and its applications beyond finance have been documented. Das and Samanta in [11] proposed a delayed fractional-order prey-predator system with fear (felt by prey) effect of predator on prey population incorporating prey refuge. The authors also established theoretical results on the existence, uniqueness, nonnegativity, and boundedness of the solutions to the proposed system. By drawing an analogy to the spread dynamics of an infectious disease, Graefa et al. in [18] derived a fractional-order susceptible-infected-removed (SIR) model to examine the user adoption and abandonment of online social networks, where the adoption is analogous to infection, and abandonment is analogous to recovery. Their results indicate that the fractional model has a unique user-prevailing equilibrium that is globally asymptotically stable. Their stability results also show that the infectious abandonment dynamics do not contribute to the stability of the user-free and user-prevailing equilibria, and that it only affects the location of the user-prevailing equilibrium. Another notable application of fractional calculus can be found in [38], in which memory effects in infectious diseases are modeled using a fractional susceptible-exposed-infectious-removed (SEIR) model, with the fractional derivatives defined in the Caputo sense. Their results indicate that the fractional model perform 35% better than the classical one.
Therefore, taking into account the results established in the reviewed literature, the aim of this paper is two fold: firstly, the paper provides a detailed derivation of a case-specific one-dimensional tfBS-PDE for pricing options where the underlying is governed by a fractal stochastic process. To the best of our knowledge, this approach has not yet been fully explored and still remains a subject of active research. Secondly, a robust unconditionally stable numerical method for solving the derived tfBS-PDE is proposed.
The rest of this paper is organized as follow: Sect. 2 presents some mathematical preliminaries on fractional differentiation. Specific focus is on the three common definitions of derivatives of fractional orders. Simultaneously, this section provides a brief account of work on the derivation of the tfBS-PDE for pricing options on dividend paying stocks. In Sect. 3, the derivation of an implicit finite difference scheme for solving the tfBS-PDE of Sect. 2 is presented. The theoretical analysis and discussion of the stability and convergence properties of the numerical scheme are discussed in Sect. 4. To substantiate and validate the theoretical claims regarding the proposed numerical method, extensive numerical experiments are presented in Sect. 5. Finally, some concluding remarks and scope for further research directions are indicated in Sect. 6.

Time-fractional Black-Scholes (BS) equation
This section first presents some fundamental fractional derivative definitions, namely, the Riemann-Liouville, Caputo, and Jumarie fractional derivatives. Secondly, a derivation of the time-fractional Black-Scholes PDE for pricing options on stocks that pay continuous dividends is presented.

Mathematical preliminaries on fractional differential operators
In most recent literature, derivatives of fractional orders are defined in the Caputo [9], Riemann-Liouville [43], as well as in the Jumarie (modified Riemann-Liouville) [24] sense. For a detailed treatment on different types of fractional derivatives, advantages and disadvantages of their usage, see [3] and references therein. Some other good discussions can be found in [21][22][23][24]. The following fractional derivative definitions can also be found in [9,[21][22][23]43] among other literature. Definition 2.1 (Caputo fractional derivative) Let f : R → R be a continuous, but not necessarily differentiable function. The Caputo fractional derivative of order α is defined as follows: (2.1) One of the main advantages of the Caputo fractional derivative is that it allows classical initial and boundary conditions to be included in the formulation of the problem [8,35]. Definition 2.2 (Riemann-Liouville fractional derivative) Let f : R → R be a continuous, but not necessarily differentiable function. Then, the Riemann-Liouville fractional derivative of order α is given by (2. 2) The Riemann-Liouville derivative has certain disadvantages when trying to model real-world phenomena with fractional differential equations, for example, the Riemann-Liouville derivative of a constant is not zero. In addition, if an arbitrary function is a constant at the origin, its fractional derivative has a singularity at the origin, for example, the exponential and Mittag-Leffler functions. These disadvantages reduce the field of application of the Riemann-Liouville fractional derivative [3]. Some of these disadvantages can be circumvented by modifying the Riemann-Liouville definition, for example, as done by Jumarie in [22] where he derived some interesting modifications of the Riemann-Liouville definition with the aid of fractional differencing theory coupled with his theory of Jumariefractional (generalized) Taylor series. The definition by Jumarie takes into account the existence of a fractional derivative at t = 0, which is undefined in the definitions of Caputo and Riemann-Liouville derivatives. Definition 2.3 (Jumarie (modified Riemann-Liouville) derivative) Let f : R → R be a continuous, but not necessarily differentiable function, and suppose f (t) is (i) a constant K, then its Jumarie fractional derivative of order α is defined by where η ∈ N.
Definition 2.4 (Fractional Taylor series) Let f : R → R be a continuous function that has a derivative of order αζ , ζ being any positive integer, and 0 < α ≤ 1, then where η ∈ N.

Derivation of the time-fractional BS-PDE
In deriving the time-fractional BS-PDE, let us first assume that the stock price dynamics follows the following fractional stochastic equation: where S and σ are respectively the price and volatility of the stock, r is the risk-free interest rate, and ω(t) denotes the standard Wiener process. In the presence of continuous dividends, denoted by δ, the above equation changes to Let us also consider the following important identities, which, according to Jumarie [24], are well consistent with the Jumarie fractional (generalized) Taylor series in (2.5): and Combining (2.9) and (2.10) yields a formula which allows for the conversion of integer derivatives to fractional derivatives, and vice versa: t) represents the value of a European put option, and suppose that V (S, t) satisfies the assumption Assumption 2.5 The function V (S, t) is sufficiently smooth with respect to S and its α derivative with respect to time exists for some α (0 < α ≤ 1).
Consider the risk-free investment interest rate dynamic equation dV = rV dt. (2.12) Multiplying both sides of (2.12) with (1α) yields Now, combining (2.13) with (2.9) yields the following variational fractional increment process: (2.14) Equation (2.14), along with (2.11), yields the following fractional interest rate dynamic equation: Since V (S, t) is sufficiently smooth with respect to S and its α-derivative (0 < α ≤ 1) with respect to t exists, applying the fractional Taylor series (2.5) of order α on V (S, t) up to remaining error term yields Combining this with Itô's lemma and equation (2.7) results in Using the conversion formula (2.11) but in terms of t, dt in (2.17) is replaced with to obtain Multiplying both sides of (2.19) with (1 + α) yields Using (2.15), the left-hand side of (2.20) can be rewritten as Using (2.21), along with (2.20), yields

Equation (2.22) can further be simplified into the following tfBS-PDE:
A robust numerical scheme for solving the tfBS-PDE in (2.23) coupled with the following boundary and terminal conditions is developed in the next section: lim S→∞ V (S, t) = 0, (2.24) where K is the strike price of the European put option and T is the maturity time.

Numerical method
In this section, an implicit numerical method for solving (2.23) along with (2.24) is presented. To begin, let L and N be positive integers and define h = 1/L and k = 1/N as the space and time step-sizes, respectively. Define S l = lh; l = 0, 1, 2, . . . , L and t n = nk; n = 0, 1, 2, . . . , N , such that S l ∈ [S min , S max ] and t n ∈ [0, T]. Furthermore, define V n l = V (S l , t n ) as the solution at the grid point (S l , t n ) = (lh, nk).

Analysis of the numerical method
In this section, theoretical results on the stability and convergence properties of the numerical scheme (3.9) are presented.
Define the truncation error as follows: such that ε n 0 = ε n L = 0 for all n. Now setting n = 1 in (3.9) and simplifying further results in which can be represented as Using the error equation (4.1), along with (4.3), yields and for n ≥ 2 the following error equation is obtained: (4.5) In matrix notation, (4.4) and (4.5) can be written as With the above notations and settings, the following theorem holds.

Theorem 4.1 The implicit finite difference scheme (3.9) is unconditionally stable and its global error satisfies
Proof Suppose n = 1 and let Then, using (4.4) results in which implies that which completes the proof of the theorem.

Convergence analysis
Let U n l be the exact solution of (2.23) with condition (2.24) at the grid point (S l , t n ) and define e n l = U n l -V n l , (4.12) such that e n = (e n 1 , e n 2 , . . . , e n L-1 ) T and e 0 = 0. Now since the errors e n l satisfy (3.10) and (4.3), the following is true for n ≥ 2: a nl e n l-1 + b nl e n l + c nl e n l+1 = n-1 j=1 ϕ j e n-j l + R n l , (4.13) and for n = 1, (4.14) In the above, the remainder term R n l is obtained from (3.9) by multiplying both sides of the equation by k α (2α). This gives where where C and C 1 are constants independent of h and k. Therefore, and (4.20) Substituting (4.18)-(4.20) into (4.15) and simplifying yields where C 2 and C 3 are constants independent of h and k.
In the next section, a set of numerical experiments are presented.

Numerical results
In this section, two numerical examples on pricing of standard European put options under the time-fractional BS-PDE (2.23), along with conditions (2.24), and implemented using the implicit difference scheme (3.9) are presented. A set of varying dividend yields and order (α) of fractional derivative ranging from 0.1 to 0.9 are considered. Numerical convergence and stability results are also presented.     O(h, k). The results are presented in Tables 1 and 2 for Example 5.1, and in Tables 3 and 4 for Example 5.2. Results presented in Tables 1 and 3 indicate that the proposed method is very suitable method for solving the tfBS-PDE in (2.23) for a range of values of α between 0 and 1. Though the accuracy of the results is better in the case when 1/2 ≤ α < 1 as compared to when 0 < α < 1/2, from the numerical point of view, α can be chosen small or large without substantially affecting the convergence of the method. That is so because of the fact that the proposed method is unconditionally stable, hence the choice of α does not affect the overall convergence of the method, see Tables 2 and 4. One may further note that the restriction on α to be between 0.5 and 1 is not too large a hindrance, as regular market conditions would not require extremely small values of α. This is true because, for extremely small α values, the returns on underlying stock price S t become negatively correlated, signaling antipersistent features, which in-turn violates some key fundamental principles of asset pricing theory. Similar observations are made for all considered dividend yields. To this end, it is worth noting that changing the dividend yield does have a significant effect on the option premium results. By carefully contrasting the premium profiles obtained for different dividend yields, it is observed that, regardless of the value of α, higher dividends yields result in higher put option premiums. This so because the underlying stock price is expected to drop by the amount of dividend payout, and hence a higher dividend would

Concluding remarks and directions for future research
In this paper, a time-fractional Black-Scholes PDE for pricing standard European put options written on a continuous dividend paying stock is derived. An implicit finite difference scheme for solving the derived tfBS-PDE was designed and analyzed. Theoretical results indicate that the proposed numerical method is unconditionally stable and converges with order O(h, k). Two numerical examples supporting the theoretical results were presented. Numerical results indicate that the proposed numerical method is very efficient for all values of α considered in the simulation experiments.
Overall, the results indicate that the tfBS-PDE model and the proposed numerical method can produce put option premiums which are very different from what classical theory suggests. These results are in agreement with existing results in literature that fractional calculus based models outperform their classical counterparts, see, for example, [15,18,20,38] among others. Another important feature to note is that the classical BS models are known to produce option premium curves which are similar in shape and hence may not fully reflect the actual market dynamics. However, the model formulated under the fractional calculus framework produces option premium curves which are quite sensitive to changes in associated market parameters, such as volatility, dividends, interest rates, etc. In summary, if sufficient market data is available, the proposed tfBS model can be calibrated to produce option price curves which take into account a variety of actual market conditions.
Since tractable analytical solutions to fractional Black-Scholes PDEs seldom exist, for further research directions, high-order numerical methods for the current proposed model, as well as high-order numerical solutions to models proposed in recent literature such as [15,20], can be explored.

Availability of data and materials
The datasets generated during the current study are available from the corresponding author on reasonable request.