A weighted finite difference method for subdiffusive Black Scholes Model

TL;DR

This paper develops a weighted finite difference numerical method for the subdiffusive Black-Scholes model, providing stability, convergence analysis, and demonstrating its effectiveness through numerical results.

Contribution

It introduces a generalized weighted finite difference scheme for the subdiffusive Black-Scholes equation with proven stability and convergence properties.

Findings

The method achieves $2-eta$ order accuracy in time.

The scheme is unconditionally stable.

Numerical tests confirm the theoretical convergence rates.

Abstract

In this paper we focus on the subdiffusive Black Scholes model. The main part of our work consists of the finite difference method as a numerical approach to the option pricing in the considered model. We derive the governing fractional differential equation and the related weighted numerical scheme being a generalization of the classical Crank-Nicolson scheme. The proposed method has $2-\alpha$ order of accuracy with respect to time where $\alpha\in(0,1)$ is the subdiffusion parameter, and $2$ with respect to space. Further, we provide the stability and convergence analysis. Finally, we present some numerical results.

Tables

Table 1. Table 1 : Order of convergence with respect to Δ ​ t Δ 𝑡 \Delta t for θ = 0 𝜃 0 \theta=0 , Δ ​ x = 0 , 2 Δ 𝑥 0 2 \Delta x=0,2 , x m ​ a ​ x = 1 subscript 𝑥 𝑚 𝑎 𝑥 1 x_{max}=1 , x m ​ i ​ n = − 1 subscript 𝑥 𝑚 𝑖 𝑛 1 x_{min}=-1 , and different α 𝛼 \alpha . We approximate u 𝑢 u by u ~ ~ 𝑢 \tilde{u} calculated for Δ ​ t = 3 , 85 × 10 − 4 Δ 𝑡 3 85 superscript 10 4 \Delta t=3,85\times 10^{-4} . Functions are evaluated at point x = − 0.01 𝑥 0.01 x=-0.01 .

$\alpha$	$h_{\Delta t}$	empirical order	theoretical value	relative error
$0.99$	$0,01$	$1,02$	$1,01$	$1,39\%$
$0.7$	$2,22\times {10}^{-3}$	$1,32$	$1,3$	$1,27\%$
$0.5$	$1,67\times {10}^{-3}$	$1,51$	$1,5$	$0,67\%$
$0.3$	$1,39\times {10}^{-3}$	$1,7$	$1,7$	$0,21\%$
$0.1$	$1,25\times {10}^{-3}$	$1,85$	$1,9$	$2,7\%$

Table 2. Table 2 : Order of convergence with respect to Δ ​ x Δ 𝑥 \Delta x for θ = 0 𝜃 0 \theta=0 , Δ ​ t = 0 , 05 Δ 𝑡 0 05 \Delta t=0,05 , x m ​ a ​ x = 10 subscript 𝑥 𝑚 𝑎 𝑥 10 x_{max}=10 , x m ​ i ​ n = − 1 subscript 𝑥 𝑚 𝑖 𝑛 1 x_{min}=-1 , and different α 𝛼 \alpha . To show that the order should not depend on α 𝛼 \alpha , in each case we take h Δ ​ x = 5 , 5 × 10 − 2 subscript ℎ Δ 𝑥 5 5 superscript 10 2 h_{\Delta x}=5,5\times 10^{-2} . We approximate u 𝑢 u by u ~ ~ 𝑢 \tilde{u} calculated for Δ ​ t = 3 , 33 × 10 − 3 Δ 𝑡 3 33 superscript 10 3 \Delta t=3,33\times 10^{-3} and Δ ​ x = 1 , 83 × 10 − 2 Δ 𝑥 1 83 superscript 10 2 \Delta x=1,83\times 10^{-2} . Functions are evaluated at point x = x m ​ a ​ x 𝑥 subscript 𝑥 𝑚 𝑎 𝑥 x=x_{max} .

$\alpha$	empirical order	theoretical value	relative error
$0.99$	$1,97$	$2$	$1,64\%$
$0.7$	$1,96$	$2$	$2,05\%$
$0.5$	$1,95$	$2$	$2,3\%$
$0.3$	$1,95$	$2$	$2,52\%$
$0.1$	$1,95$	$2$	$2,71\%$

Table 3. Table 3 : Finite difference method for different θ 𝜃 \theta , n 𝑛 n and N 𝑁 N .

(n,N)	(5000,140)	(3000,100)	(500,50)	(100,20)	(200,200)	(50,1300)
$\theta =0$	$0.63\%$	$0.87\%$	$1.74\%$	$6.31\%$	$0.58\%$	$3.17\%$
$\theta =0.25$	$0.43\%$	$0.59\%$	$1.12\%$	$4.87\%$	$0.44\%$	$3.15\%$
$\theta =0.5$	$0.23\%$	$0.31\%$	$0.61\%$	$2.04\%$	$0.3\%$	$3.13\%$
$\theta =0.6$	$6.45\times {10}^6\%$	$5.23\times {10}^4\%$	$68.34\%$	$2.58\%$	$0.25\%$	$3.12\%$
$\theta =0.9$	$1.06\times {10}^{55}\%$	$6.13\times {10}^{38}\%$	$2.09\times {10}^{17}\%$	$85\%$	$1.03\times {10}^9\%$	$3.09\%$

Table 4. Table 4 : Time related to run Finite difference method for Table 3.

(n,N)	(5000,140)	(3000,100)	(500,50)	(100,20)	(200,200)	(50,1300)
$\theta =0$	$9.7s$	$3.9s$	$0.7s$	$0.3s$	$3.6s$	$112.4s$
$\theta =0.25$	$10.3s$	$4.3s$	$0.7s$	$0.3s$	$3.6s$	$110.6s$
$\theta =0.5$	$9.7s$	$3.9s$	$0.6s$	$0.3s$	$3.6s$	$110.4s$
$\theta =0.6$	$10.5s$	$4.1s$	$0.8s$	$0.3s$	$3.6s$	$110s$
$\theta =0.9$	$10.2s$	$4.1s$	$0.7s$	$0.4s$	$3.5s$	$111s$

Equations

$$\left\{\begin{array}[]{ll}Z_{\alpha}\left(t\right)=Z\left(S_{\alpha}(t)\right),\\ Z\left(0\right)=Z_{0},\end{array}\right.$$

$$C_{BS}^{sub}\left(Z_0,K,T,\sigma ,r,\alpha \right)-P_{BS}^{sub}\left(Z_0,K,T,\sigma ,r,\alpha \right)=Z_0-K{e}^{-rT}.$$

$$\mathbb{Q}(A)={\int }_A\exp \left(-\gamma B(S_{\alpha }(T))-\frac{{\gamma }^2}{2}{S}_{\alpha }(T)\right)dP,$$

$$v\left(t\right)=C_{BS}^{sub}\left(Z_0,K,t,r,\sigma ,\alpha \right),$$

$$h\left(t\right)=C_{BS}\left(Z_0,K,t,r,\sigma \right).$$

$$h\left(t\right)=Z_0\Phi (d_{+})-K{e}^{-rt}\Phi (d_{-}),$$

$$v\left(t\right)=\mathbf{E}h(S_{\alpha }(t))={\int }_0^{\infty }h\left(z\right){\varrho }_{\alpha }\left(z,t\right)dz,$$

$$v\left(t\right)={\int }_0^{\infty }h(x)t^{-\alpha }{g}_{\alpha }(x/t^{\alpha })dx,$$

$$v\left(0\right)=h\left(0\right).$$

$$\mathcal{L}\{v\}(k)={\int }_0^{\infty }{e}^{-kt}v\left(t\right)dt={\int }_0^{\infty }{e}^{-kt}{\int }_0^{\infty }h\left(z\right){\varrho }_{\alpha }\left(z,t\right)dzdt={\int }_0^{\infty }{\int }_0^{\infty }h\left(z\right){\varrho }_{\alpha }\left(z,t\right)e^{-kt}dtdz={\int }_0^{\infty }h\left(z\right)\mathcal{L}\{{\varrho }_{\alpha }\}\left(z,k\right)dz.$$

$${\int }_0^{\infty }h\left(z\right)\mathcal{L}\{{\varrho }_{\alpha }\}\left(z,k\right)dz={\int }_0^{\infty }h\left(z\right)k^{\alpha -1}{e}^{-z{k}^{\alpha }}dz=k^{\alpha -1}{\int }_0^{\infty }h\left(z\right)e^{-z{k}^{\alpha }}dz=k^{\alpha -1}\mathcal{L}\{h\}\left(k^{\alpha }\right).$$

$$\mathcal{L}\{h\}\left(k\right)=k^{(1-\alpha )/\alpha }\mathcal{L}\{v\}\left(k^{1/\alpha }\right).$$

$$\left\{\begin{array}[]{ll}\displaystyle\frac{\partial h\left(z,t\right)}{\partial t}+\frac{1}{2}\sigma^{2}z^{2}\frac{\partial^{2}h\left(z,t\right)}{\partial z^{2}}+rz\frac{\partial h\left(z,t\right)}{\partial z}-rh\left(z,t\right)=0,\\ h\left(z,T\right)=\max(z-K,0).\end{array}\right.$$

$$k\mathcal{L}\{h\}\left(k\right)-h\left(0\right)+\frac{1}{2}{\sigma }^2{z}^2\frac{{\partial }^2\mathcal{L}\{h\}\left(k\right)}{\partial z^2}+rz\frac{\partial \mathcal{L}\{h\}\left(k\right)}{\partial z}-r\mathcal{L}\{h\}\left(k\right)=0,$$

$$\mathcal{L}\{v\}\left(k^{\frac{1}{\alpha }}\right)k^{\frac{1}{\alpha }}-v\left(0\right)+\frac{1}{2}{\sigma }^2{z}^2{k}^{\frac{1-\alpha }{\alpha }}\frac{{\partial }^2\mathcal{L}\{v\}\left(k^{\frac{1}{\alpha }}\right)}{\partial z^2}+rz{k}^{\frac{1-\alpha }{\alpha }}\frac{\partial \mathcal{L}\{v\}\left(k^{\frac{1}{\alpha }}\right)}{\partial z}-k^{\frac{1-\alpha }{\alpha }}r\mathcal{L}\{v\}\left(k^{\frac{1}{\alpha }}\right)=0.$$

$$k\mathcal{L}\{v\}\left(k\right)-v\left(0\right)=k^{1-\alpha }\left(-\frac{1}{2}{\sigma }^2{z}^2\frac{{\partial }^2\mathcal{L}\{v\}\left(k\right)}{\partial z^2}-rz\frac{\partial \mathcal{L}\{v\}\left(k\right)}{\partial z}+r\mathcal{L}\{v\}\left(k\right)\right).$$

$$\frac{\partial v\left(z,t\right)}{\partial t}={}_0^{RL}{D}_t^{1-\alpha }\left(-\frac{1}{2}{\sigma }^2{z}^2\frac{{\partial }^{2}v\left(z,t\right)}{\partial z^2}-rz\frac{\partial v\left(z,t\right)}{\partial z}+rv\left(z,t\right)\right),$$

$${}_0^{RL}{D}_t^{\alpha }g\left(t\right)=\frac{1}{\Gamma \left(1-\alpha \right)}\frac{d}{dt}{\int }_0^t{\left(t-s\right)}^{-\alpha }g\left(s\right)ds.$$

$$\mathcal{L}\{{}_0^{RL}{D}_t^{\alpha }f\left(t\right)\}=k^{\alpha }\mathcal{L}\{f\}(k)-{}_0^{RL}{D}_t^{\alpha -1}f{\left(t\right)}_{t=0},$$

$$\mathcal{L}\{{}_0^c{D}_t^{\alpha }f\left(t\right)\}=k^{\alpha }\mathcal{L}\{f\}(k)-\frac{f(0)}{k^{1-\alpha }},$$

$${}_0^c{D}_t^{\alpha }g\left(t\right)=\frac{1}{\Gamma \left(1-\alpha \right)}{\int }_0^t\frac{dg\left(s\right)}{ds}{\left(t-s\right)}^{-\alpha }ds.$$

$${}_0^c{D}_t^{\alpha }v\left(z,t\right)=-\frac{1}{2}{\sigma }^2{z}^2\frac{{\partial }^{2}v\left(z,t\right)}{\partial z^2}-rz\frac{\partial v\left(z,t\right)}{\partial z}+rv\left(z,t\right).$$

$$\{\begin{array}{l}{}_0^c{D}_t^{\alpha }v\left(z,t\right)=-\frac{1}{2}{\sigma }^2{z}^2\frac{{\partial }^{2}v\left(z,t\right)}{\partial z^2}-rz\frac{\partial v\left(z,t\right)}{\partial z}+rv\left(z,t\right),\\ v\left(z,T\right)=\max \left(z-K,0\right),\\ v\left(0,t\right)=0,\\ \underset{z\to \infty }{\lim }v\left(z,t\right)\sim z,\end{array}$$

$$x=\ln z$$

$$u\left(x,t\right)=v\left(e^x,T-t\right).$$

$$\{\begin{array}{l}{}_0^c{D}_t^{\alpha }u\left(x,t\right)=\frac{1}{2}{\sigma }^2\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}+\left(r-\frac{1}{2}{\sigma }^2\right)\frac{\partial u\left(x,t\right)}{\partial x}-ru\left(x,t\right),\\ u\left(x,0\right)=\max \left(e^x-K,0\right),\\ \underset{x\to -\infty }{\lim }u\left(x,t\right)=0,\\ \underset{x\to \infty }{\lim }u\left(x,t\right)\sim e^x,\end{array}$$

$$\{\begin{array}{l}{}_0^c{D}_t^{\alpha }u\left(x,t\right)=a\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}+b\frac{\partial u\left(x,t\right)}{\partial x}-cu\left(x,t\right),\\ u\left(x,0\right)=f\left(x\right),\\ u\left(x_{min},t\right)=p\left(t\right),\\ u\left(x_{max},t\right)=q\left(t\right),\end{array}$$

$$b_j={\left(j+1\right)}^{1-\alpha }-j^{1-\alpha },$$

$$d=\Gamma \left(2-\alpha \right)\Delta t^{\alpha },$$

$$u^k={\left(u_1^k,u_2^k,\dots ,u_{n-1}^k\right)}^T,$$

Full text

1 Introduction

Options are one of the most popular and important financial derivatives, therefore the question about their valuation has an essential meaning for financial institutions and global economy. The celebrated B-S formula for European options [3, 25] was of such great importance that the authors were awarded the Nobel Prize in Economics in 1997. After recent investigations [30] it seems that the B-S-Merton formula, despite its simplicity and clarity, cannot be used in many cases. One of the most vivid examples is the case when the dynamics of the underlying instrument has a tendency to have constant periods or sudden jumps [6]. The feature of the market stagnation can be observed e.g. in emerging markets [4] and in interest rate behavior [16]. The classical B-S model was proposed under some strict assumptions, however some improved models have been established to weaken these assumptions, such as stochastic volatility model [12], stochastic interest model [23], models with transactions costs [1, 8], jump - diffusion model [24].

In recent years the theory of fractional differential equations found important applications in econometrics and finance [22]. Also many researchers have investigated the generalization of the B-S equation into fractional case. The reason of this generalization is the fractal structure for financial market and increasing interest of fractional calculus. Usually the procedure is to replace the standard Brownian motion used in the classical model by the fractional Brownian motion [2, 17, 33]. More precisely, the order of time variation $dt$ is replaced by the Hurst exponent $H$ $(0<H<1)$. Changing the parameter of self-similarity out of the case $H=1/2$ leads to the lack of martingale property of the process describing dynamics of financial asset. It is equivalent that such generalized model allows arbitrage. We proceed with a completely different approach. We replace the Geometric Brownian Motion $Z(t)$ by the subdiffusive Geometric Brownian Motion $Z(S_{\alpha}(t))$ with the inverse subordinator $S_{\alpha}(t)$ [29] in the definition of the process describing the underlying asset. In this way we find the corresponding fractional differential equation, which is different than most of considered in the literature but the same as is given in [32] and [35]. For other results related to subdiffusive B-S model see [10, 11, 13, 34].

Many efficient numerical methods have been proposed for solving fractional differential equations, which include finite difference methods, finite element methods, finite volume methods, spectral methods and meshless methods (see [35] and references therein).

Developing numerical discretization methods for fractional integrals and derivatives is one of the important topics in fractional calculus due to its wide applications. In this work we propose the quadrature method for approximation the Caputo derivative which implies the order of accuracy equal to $2-\alpha$ . In [9, 19] authors studied approximations of order $3-\alpha$, while in [5] of $4-\alpha$. In [19] the Caputo derivative was approximated using the $r+1$-th Lagrange interpolation, and obtained a series of high-order numerical schemes with accuracy of $r+1-\alpha$ at $n$-th steps ($n\geq r$). With higher order of accuracy the level of complexity increases and the stability of the numerical scheme can be lost.

2 Subdiffusive B-S Model

2.1 Assumptions of the subdiffusive B-S Model

Let us consider a market, whose evolution is taking place up to time horizon $T$ and is contained in the probability space $(\Omega,\mathcal{F},\mathbb{P})$. Here, $\Omega$ is the sample space, $\mathcal{F}$ is filtration interpreted as the information about history of asset price which completely is available for the investor and $\mathbb{P}$ is the “objective” probability measure. The assumptions are the same as in the classical case [15] with the exception that we do not have to assume the market liquidity and that the underlying instrument instead of Geometric Brownian Motion (GBM) has to follow subdiffusive GBM [20]:

[TABLE]

where $Z_{\alpha}\left(t\right)$ - the price of the underlying asset, $Z(t)=Z\left(0\right)\exp\left(\mu t+\sigma B_{t}\right)$, $\sigma$ - volatility (constant), $\mu$ - drift (constant), $B_{t}$ - Brownian motion, $S_{\alpha}(t)$ - the inverse $\alpha$-stable subordinator defined as $S_{\alpha}(t)=\inf(\tau>0:U_{\alpha}(\tau)>t)$ [20], where $U_{\alpha}(t)$ is the $\alpha$-stable subordinator [29], $0<\alpha<1$. We assume $S_{\alpha}(t)$ is independent of $B_{t}$ for each $t\in[0,T]$.

The subdiffusive B-S is the generalization of the standard model to the cases, where the markets can be non-liquid. Then the subdiffusion parameter $\alpha\in(0,1)$ can be considered as the measure of non-liquidness. The constant periods appear more frequently with decreasing $\alpha$. With $\alpha\rightarrow 1$ the subdiffusive B-S model reduces to the classical (liquid) case. Due to its simplicity and practicality, the classical B-S Model is one of the most widely used in option pricing. Although in contrast to the subdiffusive case it does not take into account the empirical property of constant price periods. In Figure $1$ we compare sample simulation of underlying asset in classical and subdiffusive market model. Even short stagnation of a market can not be simulated by classical B-S model. As a generalization of the B-S model, its subdiffusive equivalent can be used in wide range cases - including all markets where B-S can be applied. The method of subdiffusive B-S model calibration from empirical data is described in [27].

By the classical put-call parity [26] we have following fact:

Proposition 2.1

[20]** For the fair price of European call and put options in subdiffusive B-S model we have following relationship:

[TABLE]

Here and in whole paper: $C_{BS}^{sub}$ - fair price of European call, $P_{BS}^{sub}$ - fair price of European put, $K$ - strike, $T$ - maturity, $\sigma$ - volatility, $r$ - interest rate.

One of the most expected property of the market is that there is no possibility to gain money without taking the risk. This property is called the lack of arbitrage and formally means that the self-financing strategy $\phi$ which follows to a positive profit without any probability of intermediate loss can not be constructed [7]. By the Fundamental theorem of asset pricing [7], the market model described by $(\Omega,\mathcal{F},\mathbb{P})$ and underlying instrument $Z_{\alpha}(t)$ with filtration $\mathcal{F}_{t\in[0,T]}$ is arbitrage-free if and only if there exists a probability measure $\mathbb{Q}$, (called the risk-neutral measure) equivalent to $\mathbb{P}$ such that the asset $Z_{\alpha}(t)$ is a martingale with respect to $\mathbb{Q}$. Under this measure, financial instruments have the same expected rate of return, regardless the variability of the prices. This is in contrast to the physical probability measure (the actual probability distribution of prices), under which more risky instruments have a higher expected rate of return than less risky instruments.

Let us introduce the probability measure

[TABLE]

where $\gamma=\frac{\mu+\frac{\sigma^{2}}{2}}{\sigma}$, $A\in\mathcal{F}$. As it is shown in [20] the process $Z_{\alpha}(t)$ is martingale with respect to $\mathbb{Q}$, so we have the following

Theorem 2.1

[20]** The subdiffusive B-S Model is arbitrage-free.

Another property of market model is the so-called completeness. Intuitively the market model is complete if the set of possible gambles on future states-of-the-world can be constructed with existing assets. More formally, the market model is complete if every $\mathcal{F}_{t\in[0,T]}$ -measurable random variable $X$ admits a replicating self- financing strategy $\phi$ [7]. The Second Fundamental theorem of asset pricing [7] states that a market model described by $(\Omega,\mathcal{F},\mathbb{P})$ and underlying instrument $Z_{\alpha}(t)$ with filtration $\mathcal{F}_{t\in[0,T]}$ is complete if and only if there is a unique martingale measure equivalent to $\mathbb{P}$.

Theorem 2.2

[20]** The market model in which the price of underlying instrument follows the subdiffusive GBM $Z_{\alpha}(t)$ is incomplete.

Market incompleteness means that there is no unique fair price of financial derivatives, because for different martingale measures different prices could be obtained. Despite $\mathbb{Q}$ defined in $(1)$ is not unique, in the sense of criterion of minimal relative entropy it is the “best” martingale measure. It means that the measure $\mathbb{Q}$ minimizes the distance to the measure $\mathbb{P}$ [21] . Other essential fact is that for $\alpha\to 1$, $\mathbb{Q}$ reduces to the measure of the classical B-S model which is arbitrage-free and complete. It is consistent with our intuition if we consider subdiffusive B-S model as the generalization of the standard B-S model. Thus, in whole paper we will use the martingale measure $\mathbb{Q}$ defined in $(1)$ as a reference measure.

2.2 The fair price of a call option in the subdiffusive B-S model

Let us define the fair price of a call option for subdiffusive and classical B-S model:

[TABLE]

[TABLE]

By the B-S formula [3, 25], we have

[TABLE]

with $\displaystyle d_{\pm}=\frac{\log\frac{Z_{0}}{K}+\left(r\pm\frac{1}{2}\sigma^{2}\right)t}{\sigma\sqrt{t}}.$ Here $\Phi$ denotes the cumulative distribution function of standard normal distribution.

By [20], the relation between functions $v$ and $h$ is as follow

[TABLE]

where, $\varrho_{\alpha}(x,t)$ is the PDF of $S_{\alpha}(t)$. Moreover, $v(t)$ follows

[TABLE]

where $g_{\alpha}(z)$ is given in terms of Fox function $H^{1,0}_{1,1}(z^{(1-\alpha,\alpha)}_{(0,1)})$ (see [20] and references therein).

Note that

[TABLE]

Let us consider the Laplace transform of the function $v$:

[TABLE]

In the same way as in [20] we find that $\mathcal{L}\{\varrho_{\alpha}\}=k^{\alpha-1}e^{-zk^{\alpha}}$. Hence, we have:

[TABLE]

So as a conclusion we have the following result:

[TABLE]

Let us write the B-S equation describing $h\left(z,t\right)$ [15]:

[TABLE]

Now let us take the Laplace transform with respect to $t$:

[TABLE]

Then we use formula $(2)$ and the fact that $v(0)=h(0)$, obtaining:

[TABLE]

Now let us change variable - we replace $k$ by $k^{\alpha}$:

[TABLE]

Inverting the Laplace transform, we get:

[TABLE]

where $\alpha\in\left(0,1\right)$, ${}_{0}^{RL}D_{t}^{\alpha}$ is Riemann-Louville fractional derivative defined as

[TABLE]

Above we used the fact that for $\alpha\in\left(0,1\right)$, the Laplace transform for the Riemann-Louville fractional derivative follows [28]:

[TABLE]

where $f\in C^{1}$. Note that the Laplace transform for the Caputo fractional derivative follows [28]:

[TABLE]

where $f\in C^{1}$, $\alpha\in\left(0,1\right)$ and ${}_{0}^{c}D_{t}^{\alpha}$ is defined as [28]

[TABLE]

Using basic properties of fractional derivatives, we transform $(3)$ into

[TABLE]

By this way we have found the following system:

[TABLE]

for $\left(z,t\right)\in\left(0,\infty\right)\times\left(0,T\right)$.

Let us introduce the following variable:

[TABLE]

and the function:

[TABLE]

Hence, we have:

Theorem 2.3

The fair price of a call option in the subdiffusive B-S model is equal to $v(z,t)$, where $v(z,t)$ satisfies $(4)$ and $(5)$, and $u(x,t)$ is the solution of the system:

[TABLE]

for $\left(x,t\right)\in\left(-\infty,\infty\right)\times\left(0,T\right)$.

The solution of $(6)$ exists and it is unique (see [14] and references therein).

3 Finite difference method

To solve the above model numerically we will approximate limits by finite numbers and derivatives by finite differences. After obtaining the discrete analogue of $(6)$ we will solve the problem recursively using initial-boundary conditions. We will proceed similarly, as it was done for implicit method in [35].

3.1 Weighted scheme for the subdiffusive B-S model

The system $(6)$ has the following form:

[TABLE]

where $a=\displaystyle\frac{1}{2}\sigma^{2}$, $b=\left(r-\displaystyle\frac{1}{2}\sigma^{2}\right)$, $c=r$, $\displaystyle f\left(x\right)=\max\left(\exp(x)-K,0\right)$, $\displaystyle p\left(t\right)\rightarrow 0$ if $\displaystyle x_{min}\rightarrow-\infty$, $\displaystyle q\left(t\right)\rightarrow\exp(x_{max})$ if $\displaystyle x_{max}\rightarrow\infty$. The put-call parity implies that $\displaystyle p\left(t\right)=0$, $\displaystyle q\left(t\right)=\exp(x_{max})-K\exp(-r\left(T-t\right))$.

Let us denote

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $u^{k}_{i}=u\left(x_{i},t_{k}\right)$, $f_{i}=f\left(x_{i}\right)$, $i=2,\dots,n-2$, $q^{k}=q\left(t_{k}\right)$, $k=0,1,\dots,N$, moreover $\displaystyle\Delta t=T/N$, $\displaystyle\Delta x=(x_{max}-x_{min})/n$ are time and space steps respectively. We will use the following approximations for space derivatives:

[TABLE]

We approximate the fractional-time derivative by:

[TABLE]

After omitting the truncation errors, the implicit discrete scheme can be expressed in following form:

[TABLE]

where $k\geq 1$, $A=\left(a_{ij}\right)_{\left(n-1\right)\times\left(n-1\right)}$, such that:

$a_{ij}=$ $\begin{cases}1+2\displaystyle\frac{ad}{\Delta x^{2}}+cd,&\text{for }j=i,\text{ }i=1,2\dots,n-1,\\ -\left(\displaystyle\frac{ad}{\Delta x^{2}}-\displaystyle\frac{bd}{2\Delta x}\right),&\text{for }j=i-1,\text{ }i=2\dots,n-1,\\ -\left(\displaystyle\frac{ad}{\Delta x^{2}}+\displaystyle\frac{bd}{2\Delta x}\right),&\text{for }j=i+1,\text{ }i=1\dots,n-2,\\ 0,&\text{in other cases}\end{cases}$

The corresponding initial-boundary conditions are

[TABLE]

where $k\geq 1$. In similar way let us write the explicit discrete scheme. We use approximations for space derivatives as follows:

[TABLE]

In the matrix form the explicit discrete scheme can be expressed in the following form:

[TABLE]

where $k\geq 1$, $B=\left(b_{ij}\right)_{\left(n-1\right)\times\left(n-1\right)}$, such that:

$b_{ij}=$ $\begin{cases}-\left(2\displaystyle\frac{ad}{\Delta x^{2}}+cd\right),&\text{for }j=i,\text{ }i=1\dots,n-1,\\ \displaystyle\frac{ad}{\Delta x^{2}}+\displaystyle\frac{bd}{2\Delta x},&\text{for }j=i+1,\text{ }i=1\dots,n-2,\\ \displaystyle\frac{ad}{\Delta x^{2}}-\displaystyle\frac{bd}{2\Delta x},&\text{for }j=i-1,\text{ }i=2\dots,n-1,\\ 0,&\text{in other cases}\end{cases}$

Taking the linear combination of $(7)$ and $(10)$ we obtain a weighted scheme:

[TABLE]

where $k\geq 1$, $C=\theta I+\left(1-\theta\right)A$, $\theta\in[0,1]$ and the corresponding initial-boundary conditions are defined in $(8)$. Let us denote that in the case of the classical B-S model $\theta=1/2$ defines the C-N scheme [15]. Motivated by this fact in whole paper we assume the following:

Definition 3.1

Scheme $(11)$ with $\theta=1/2$ is called the C-N discrete scheme.

3.2 Consistency of the weighted discrete scheme

In this section we will show the following

Theorem 3.1

For $\theta\in[0,1]$ and $1\leq i\leq n,1\leq j\leq N$, the truncation error $R_{i}^{j}$ follows

[TABLE]

Proof:

As it was shown in [18], the Caputo derivative can be expressed as follows

[TABLE]

where

[TABLE]

On the other hand, we can apply $(9)$. The full formulation of the discrete weighted scheme with the truncation error has the form:

[TABLE]

Here $k\geq 1$, $i=1,2,\dots,n-1$ and the corresponding initial-boundary conditions are defined in $(8)$. By the approximation of the Caputo derivative and $(9)$ we have

[TABLE]

where $C_{i}^{j}$ are constants ($1\leq i\leq n,1\leq j\leq N$).

Let us denote $\displaystyle C_{max}=\max_{1\leq i\leq n,1\leq j\leq N}C_{i}^{j}$. Then, for the truncation error it holds that

[TABLE]

Note that the parameter $\theta$ has no influence in the above analysis. $\blacksquare$

3.3 Stability of the weighted discrete scheme

We will proceed using von Neumann method. For $l=0,1,\dots n$, $k=0,1,..,N$ let us denote $u_{l}^{k}=u\left(x_{l},t_{k}\right)$ - the exact solution of numerical scheme, $\hat{u}_{l}^{k}$ - some approximation of $u_{l}^{k}$. After omitting truncation error and introducing $e_{l}^{k}=u_{l}^{k}-\hat{u}_{l}^{k}$, $(12)$ has the form:

[TABLE]

where $k\geq 1$, $i=1,2,\dots,n-1$. We introduce the following grid function:

$\displaystyle e^{k}(x)=$

[TABLE]

Because $\displaystyle e^{k}_{0}=e^{k}_{n}$, we make a periodic expansion for $e^{k}_{l}$ with period $Y=x_{max}-x_{min}$. Then $\displaystyle e^{k}(x)$ has the following Fourier series extension:

[TABLE]

where $\displaystyle v_{j}^{k}=\frac{1}{Y}\int_{0}^{Y}e^{k}(x)\exp({2j\pi xi/Y})dx$, $i=\sqrt{-1}$, $k=0,1,\dots N$. We define the norm $\left\lVert\cdot\right\rVert_{\Delta x}$ as

[TABLE]

where

[TABLE]

Because $e_{0}^{k}=e_{n}^{k}=0,$ it follows

[TABLE]

where $\left\lVert\cdot\right\rVert$ is $L^{2}[0,Y].$

Using the Parseval identity we have:

[TABLE]

$k=0,1,\dots,N.$ Based on the above analysis and the fact that $x_{l}=x_{min}+lh$, we infer that the solution of $(13)$, has the form:

[TABLE]

where $\displaystyle\lambda=\frac{2\pi l}{Y}.$ Substituting into $(13)$ we get:

[TABLE]

To continue we have to find a relation between coefficients $b_{j}$.

Proposition 3.1

Coefficients $b_{j}=(j+1)^{1-\alpha}-j^{1-\alpha}$ satisfy:

$b_{j}>0,$* $j=0,1\dots$* 2. 2.

$1=b_{0}>b_{1}>\dots>b_{k}$** 3. 3.

$\displaystyle\lim_{k\to\infty}b_{k}=0$** 4. 4.

$\displaystyle\sum_{j=0}^{k-1}(b_{j}-b_{j+1})+b_{k}=1$**

Proof:

$b_{j}=(j+1)^{1-\alpha}-j^{1-\alpha}>0,$ for $j\geq 0$ and $\alpha\in(0,1)$. 2. 2.

For $x\geq 0$ let us take consider the function $b(x)=(x+1)^{1-\alpha}-x^{1-\alpha}$. Note that $b^{{}^{\prime}}(x)=(1-\alpha)((x+1)^{-\alpha}-x^{-\alpha})<0$, so the function is strictly decreasing for $x\geq 0$. 3. 3.

It is the consequence of $(1)$ and $(2)$ because strictly decreasing sequence of positive coefficients is converging to [math]. 4. 4.

$\displaystyle\sum_{j=0}^{k-1}(b_{j}-b_{j+1})+b_{k}=(1-b_{1})+(b_{1}-b_{2})+(b_{2}-b_{3})+\dots+(b_{k-1}-b_{k})+b_{k}=1.$

$\blacksquare$

Now we will check under which conditions $\left|v^{n}\right|\leq\left|v^{0}\right|$ for each $n=1,\dots,N$. Then $\left\lVert e^{k}\right\rVert\leq\left\lVert e^{0}\right\rVert$, in other words the weighted scheme is stable.

Theorem 3.2

Let $\theta\in[0,1)$. If

1. (i)

[TABLE]

or 2. (ii)

$\displaystyle 1-\log_{2}\left(2-\frac{\theta}{1-\theta}\right)>\alpha$* or $\theta=1$, and the inequality*

[TABLE]

holds, then the scheme $(11)$ is stable.

Proof:

We have to show that $v^{n}$ defined in $(14)$ follows $\displaystyle\left|v^{n}\right|\leq\left|v^{0}\right|$ for $n=1,2,\dots,k.$ Let us denote

[TABLE]

Let us observe that $\displaystyle Re{\text{ }\zeta}=\sin^{2}\left(\frac{\lambda\Delta x}{2}\right)\frac{4ad}{\Delta x^{2}}+cd>0$. The proof of this fact is immediate because $a,d,c,\Delta x>0$.

At the beginning we will show that both statements imply

[TABLE]

Let us assume the first statement. Then $\displaystyle 1-\log_{2}\left(2-\frac{\theta}{1-\theta}\right)\leq\alpha$ is equivalent to $\displaystyle\theta-\left(1-\theta\right)\left(b_{0}-b_{1}\right)\leq 0.$

[TABLE]

So

[TABLE]

holds for each $\Delta t$, $\Delta x>0.$ Let us observe that

[TABLE]

It follows that:

[TABLE]

Note that the right-hand side expression higher than [math] is equivalent to $(17)$. Let us assume the second statement. The $(16)$ is equivalent to

[TABLE]

Let us observe that if $\displaystyle 1-\log_{2}\left(2-\frac{\theta}{1-\theta}\right)>\alpha$ or $\theta=1$, then

[TABLE]

So

[TABLE]

Similarly as it was done for the first statement, it can be shown that the right-hand side expression higher than [math] implies $(17)$.

We will follow the mathematical induction method to show that for each $n=1,2\dots,N$ there holds $\left|v^{n}\right|\leq\left|v^{0}\right|$.

$n=1$ By the identity

[TABLE]

the first equation of $(15)$ has the form

[TABLE]

It is equivalent to

[TABLE]

[TABLE]

So

[TABLE]

It is easy to check that $(17)$ implies $\displaystyle\left|\frac{1-\zeta\theta}{\zeta\left(1-\theta\right)+1}\right|\leq 1$, so

[TABLE] 2. 2.

Let us suppose that

[TABLE]

for $n=1,2,\dots,k,$ $k<N.$

To complete the proof we have to show that

[TABLE]

By the second equation of $(15)$, for $k\geq 1$ we have:

[TABLE]

it is equivalent to

[TABLE]

So

[TABLE]

Dividing by $\left|\left(1-\theta\right)\zeta+1\right|$ we get

[TABLE]

where the second inequality holds because $\displaystyle Re{\text{ }\zeta}>0$ and the latest by $(17)$. As a result we have

[TABLE]

By the mathematical induction method the proof is completed. $\blacksquare$

In particular, the implicit scheme is unconditionally stable for each $\alpha\in\left(0,1\right)$. Similarly, the explicit and the C-N schemes are conditionally stable for each $\alpha\in\left(0,1\right)$.

3.4 Convergence of the weighted discrete scheme

Let us denote $u_{l}^{k}=u(x_{l},t_{k})$ - the exact solution of $(6)$ evaluated at the grid point, $\hat{u}_{l}^{k}$ - the solution of the numerical scheme $(11)$. Let us define the error at the point $(x_{l},t_{k})$ by $E_{l}^{k}=u_{l}^{k}-\hat{u}_{l}^{k}$, $l=0,1,\dots n$, $k=0,1,..,N$.

Similarly, as $(13)$ we get the following system:

[TABLE]

where $k=1,\ldots,N$, $i=1,2,\dots,n-1$ and $\displaystyle\theta\in\left[0,1\right]$. Similarly as in case of the stability we will proceed with the von Neumann method.

We introduce the following grid functions:

$E^{k}(x)=$

[TABLE]

$R^{k}(x)=$

[TABLE]

Because $E^{k}_{0}=E^{k}_{n}$, we make a periodic expansion for $E^{k}_{l}$ with the period $Y=x_{max}-x_{min}$. Then $E^{k}(x)$ has the following Fourier series extension:

[TABLE]

where $\displaystyle w_{j}^{k}=\frac{1}{Y}\int_{0}^{Y}E^{k}(x)\exp({2j\pi xi/Y})dx$, $i=\sqrt{-1}$, $k=0,1,\dots N$.

By analogy, because $R^{k}_{0}=R^{k}_{n}$, we make a periodic expansion for $R^{k}_{l}$ with the period $Y$. Then $R^{k}(x)$ has the following Fourier series extension:

[TABLE]

where $\displaystyle r_{j}^{k}=\frac{1}{Y}\int_{0}^{Y}R^{k}(x)\exp({2j\pi xi/Y})dx$, $i=\sqrt{-1}$, $k=0,1,\dots N$.

We define the norm $\left\lVert\cdot\right\rVert_{\Delta x}$ as

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

Because $E_{0}^{k}=E_{n}^{k}=0,$ and $R_{0}^{k}=R_{n}^{k}=0,$ there holds

[TABLE]

[TABLE]

Using the Parseval identity we have:

[TABLE]

where $k=0,1,\dots,N.$ Based on the above analysis and the fact that $x_{l}=x_{min}+lh$, we suppose that the solution of $\left(18\right)$, has the form:

[TABLE]

[TABLE]

where $\displaystyle\lambda=\frac{2\pi l}{Y}.$ Substituting into $\left(18\right)$ we get:

[TABLE]

where $k=1,\ldots,N-1$. Let us denote

[TABLE]

Then, taking into account that $r^{0}=0$ and $w^{0}=0$, $(20)$ has the form

[TABLE]

where $k=1,\ldots,N-1$, $\theta\in[0,1]$ and $\zeta$ is previously defined.

Lemma 3.1

If condition $(i)$ or $(ii)$ in Theorem $3.2$ is satisfied, then $w^{k}$ follows

[TABLE]

where $k=0,1,\dots,N-1$ and the constant $C_{1}$ is independent of $\theta$, $\Delta t$ and $\Delta x$.

Proof:

Because $\displaystyle R_{l}^{k}=O\left(\Delta t^{2-\alpha}+\Delta x^{2}\right),$ so there exists a positive constant $C_{1}$, such that

[TABLE]

Then

[TABLE]

Convergence of the right series in second line of $\left(19\right)$ implies that

[TABLE]

Then by $(21)$ and Proposition 3.1 we have

[TABLE]

The first inequality is true because $Re{\text{ }\zeta}>0$. Now let us suppose that

[TABLE]

where $n=2,3,\dots,k$ and $C_{1}$ is a constant independent of $\theta$, $\Delta t$ and $\Delta x$.

By $(21)$ we have

[TABLE]

Dividing by the coefficient $\left|\left(1-\theta\right)\zeta+1\right|$ , we get

[TABLE]

The last inequality is true by $(17)$ and because $Re{\text{ }\zeta}>0$.

By the mathematical induction the proof is completed. $\blacksquare$

Theorem 3.3

If condition $(i)$ or $(ii)$ in Theorem $3.2$ is satisfied, then the discrete scheme $\left(11\right)$ is convergent and it follows that

[TABLE]

for $k=1,2,\dots,N,$ where $C_{1}$ and $C_{2}$ are positive constants independent of $\theta$, $\Delta t$ and $\Delta x$.

Proof:

Let us observe that $b_{k-1}^{-1}k^{-\alpha}\leq 1/(1-\alpha)$, $k=1,2,\dots,N$. By the Lemma $3.1$

[TABLE]

Similarly, by $\left(19\right)$ and $\left(22\right)$ we have the following:

[TABLE]

After taking $\displaystyle C_{2}=\frac{T^{\alpha}\sqrt{Y}}{1-\alpha}$ the proof is completed. $\blacksquare$

Let us observe that as a direct conclusion of Theorem 3.3 we get, that the optimal choice of $\theta$ for given $\alpha$ is such that $\displaystyle\log_{2}\left(2-\frac{\check{\theta_{\alpha}}}{1-\check{\theta_{\alpha}}}\right)=1-\alpha,$ equivalently $\displaystyle\check{\theta_{\alpha}}=\frac{2-2^{1-\alpha}}{3-2^{1-\alpha}}$. Then the lowest boundary for an error is achieved without loosing the unconditional stability/convergence. In Figure $2$ the jump into negative values is the result of the increasing error, what the consequence of lack of the stability is. Note that in the case of classical B-S model ($\alpha=1$) for implicit scheme, the method has $(\Delta x^{2}+\Delta t)$, but for $\check{\theta_{1}}=1/2$ the method has $(\Delta x^{2}+\Delta t^{2})$ order of convergence [15]. Similarly the C-N scheme is unconditionally convergent only for $\alpha=1$. Example 2 confirms that for $\alpha$ close to $1$ the C-N has the lowest numerical error without accelerating time of computation.

As it is shown in [31], in the similar class of problems the solution has a weak singularity near the initial time $t=0$. Due to this fact, the order of convergence with respect to time can fall from the value of $2-\alpha$ to $\alpha$. In our case the solution $u$ has no singularities in the considered domain. It is clear from the direct interpretation of $u$ as a price of an option, which for initial time $t=0$ is equal to a payoff function.

The case $\theta=0$ is the implicit numerical scheme investigated in [35]. Our work confirms the previous results and generalize them for other types of numerical schemes. We also find the optimal value of parameter $\theta$ for the following class of fractional-type problems

[TABLE]

for $\left(x,t\right)\in\mathbb{X}\times\left(0,T\right)$. Here, $\mathbb{X}=[x_{min},x_{max}]$ is a fixed interval, $T$ is a time horizon, $a$, $b$, $c\in\mathbb{R}$, $\alpha\in(0,1)$. We assume that $u$, $q$, $p$, $f$ are smooth enough (see [14] and references therein).

3.5 Numerical examples

Example 1

Let us take parameters $T=1$, $Z_{0}=1$, $\sigma=1$, $r=0.04$, $K=2$. Using the fact that

[TABLE]

where $p$ is the order of convergence and $\tilde{u}(h)$ is the solution of the numerical scheme $(11)$ (evaluated at the fixed point, similarly as $u$) for the length of the grid equal to $h$, we can numerically check the order of convergence (with respect to each variable) of the numerical scheme. The comparison prepared for both variables represent Table 1 and Table 2. The empirical order related to $\Delta t$ and $\Delta x$ should be close to $2-\alpha$ and $2$ respectively.

Example 2

Let us consider the case $\sigma=1$, $r=0.04$, $x_{max}=10$, $x_{min}=-20$, $T=4$, $K=2$, $Z_{0}=1$. The Tables $3$ and $4$ present the comparison of error and time of computation for European call option, in dependence of parameters $\theta$, $n$, $N$. The simulations are made for $\alpha=0.999$ and they are compared to the B-S formula’s result equal $0.593$. The most optimal choice of $\theta$ is close to $\check{\theta_{1}}=1/2$.

Example 3

Let us take the parameters $Z_{0}=\sigma=1$, $r=0.04$, $n=1000$, $N=140$, $x_{max}=10$, $x_{min}=-20$, $\theta=\check{\theta_{\alpha}}$. Since $ES_{\alpha}\left(T\right)=\frac{T^{\alpha}}{\Gamma\left(1-\alpha\right)}$, for $T>1$ the price of European call option is increasing function of $\alpha$ (see Figures $4$ and $5$), similarly for $T<1$ the price of European call option is decreasing function of $\alpha$. For $T>1$ it follows our intuition, because with decreasing value of $\alpha$, the constant periods in dynamics of underlying instrument appear more frequently. Such asset can be considered as more predictable, so value of its European call should be lower than the same options on instruments driven by higher values of $\alpha$. The dependence of the European call option on $T$ and $\alpha$ is presented in Figure $3$. In Figure $5$ we compare the finite difference method (FD) to the Monte Carlo method (MC) introduced in [20]. MC is oscillating around the FD output. Both methods are efficient and can be used to compute the prices of European options.

4 Summary

In this paper:

• –

We have shown that the solution of fractional B-S equation is equal to the fair price of European option with respect to $\mathbb{Q}$ in subdiffusive B-S model.

• –

We have introduced weighted numerical scheme for this equation. It allows us to approximate the fair price of European call option in subdiffusive B-S model.

• –

We have given condition under which the discrete scheme is stable and convergent.

• –

We have found the optimal choice of discretization parameter $\theta$ in dependence of subdiffusion parameter $\alpha$. Such numerical scheme is unconditionally stable, unconditionally convergent and has the lowest numerical error.

• –

We have presented some numerical examples to illustrate introduced theory.

We believe that the numerical techniques presented in this paper can successfully be repeated for other fractional diffusion-type problems.

Acknowledgments

The research of M.M. was partially supported by NCN-DFG Beethoven Grant No.2016/23/G/ST1/04083.