On a parabolic partial differential equation and system modeling a production planning problem

1.   Introduction

Assuming that the price of underlying assets satisfies the geometric Brownian motion, the Black-Scholes option pricing model was firstly proposed in 1973 which depends only on the risk-free interest rate and the volatility ^[1]. The Black-Scholes model quickly attracted a great deal of attention from the communities of both academic researcher and engineer. In order to improve the efficiency of the model and also fit the practical market, a sequence of option pricing models were proposed and studied, for instance, jump-diffusion model ^[2,3], stochastic volatility model ^[4,5] and the fractional option pricing models based on Lévy process including of finite moment log stable (FMLS) ^[6], KoBol ^[7,8] and CGMY ^[9] model.

In order to solve the fractional option pricing models based on Lévy process, a number of numerical approaches were proposed and well studied in the past decades. Cont and Voltchkova ^[10] first presented a finite difference methods for solving the fractional European option pricing model driven by exponential Lévy process and studied the stability and convergence. Then, Cartea and del-Castillo-Negrete ^[11] rewrote FMLS, CGMY and KoBol option pricing models as a general fractional partial differential equation and studied the shifted Grünwald difference (SGD) formula for the numerical solution. Marom and Momoniat ^[12] further compared the numerical solutions of three fractional option pricing models based on Lévy process. Chen and Wang ^[13] developed a numerical scheme with second-order accuracy in both the spatial and time mesh size for pricing European and American option under a geometric Lévy process. Recently, Zhang et al. ^[14] constructed an implicit numerical scheme with second-order accuracy for FMLS model and used BiCGstab method to solve the discreted linear equations.

As we known, exotic options play important roles and are widely used in the practical finance market ^[15]. However, it is a great challenge to obtain the solution for traditional numerical methods since the non-smooth payoffs usually lead to serious degradation in the convergence of the numerical schemes and result in inaccurate and discontinuous solution near the strike or barrier. For instance, the well known second-order implicit schemes, Crank-Nicolson method, are prone to spurious oscillations unless the time step size is small enough. To overcome this difficulty, Wade et al. ^[16] studied the Padé schemes to smooth the Crank-Nicolson scheme to get fourth-order schemes for pricing barrier European option models, see also ^{[15,17,18,19]} for the Padé schemes for different exotic option models.

For the fractional exotic options under FMLS model, a class of fourth-order numerical schemes are presented and studied in this paper. We first discretize the fractional option pricing models with weighted and shifted Grünwald difference (WSGD) formula, which is of second-order accuracy in space direction. Then, by making use of the Padé schemes for the time direction, an L-stable and fourth-order accurate scheme is obtained. The convergence of the proposed numerical scheme is proved when the spatial discretization matrix is positive definite and has lower Hessenberg Toeplitz structure, without the assumption of the self-adjoint operator. The proposed method is adapted to be implemented on parallel processors by making use of partial fraction. Numerical experiments on digital option and barrier option are presented to verify the efficiency and accuracy of our numerical schemes.

The structure of this paper is organized as follows. The weighted and shifted Grünwald difference schemes for space discretization is presented in Section 2. Time stepping schemes on Padé approximation and the implementations are presented in Section 3. We then analyze and prove the convergence of the numerical scheme in Section 4. Numerical experiments are given to show the accuracy and efficiency of the proposed schemes in Section 5. Finally, conclusions are drawn in Section 6.

2.   The FMLS model and the discretization in space

Denote $S_t$ as the asset price at time $t$, the FMLS model ^[6] can be written as follow:

where $V(x, t)$ is the price of the option at the time $t$ before the expiry time $T$, $x = \ln S_t$, $_{- \infty}D _ { x } ^ { \alpha } (1 < \alpha < 2)$ is the left Riemann-Liouville derivative ^[20], $\nu = - \frac { 1 } { 2 } \sigma ^ { \alpha } \sec \left(\frac { \alpha \pi } { 2 } \right)$, $S_t$ is the price of underlying asset at time $t$, $r$ is the risk free interest rate and $\sigma$ is the volatility.

In order to solve the FMLS model numerically, we first transform (2.1) into a forward problem by using the transformation $t^* = T - t$, and then drop $*$ for simplicity of the notation. Then we truncate the interval of $x$ to a finite interval $\left[B_{ d }, B_{ u } \right]$, and consider the following FMLS model for pricing European option

where the left Riemann-Liouville derivative ^[20] is defined as

and the initial and boundary conditions are as follows:

Let $M$ and $N$ be the number of the uniform discrete points in the space and time direction respectively, $h = \left(B _ { u } - B _ { d } \right) / M$ and $\tau = T / N$ be the corresponding step length. Define $t _ { j } = j \tau (j = 0, 1, 2, \cdots, N), x _ { i } = B _ { d } + i h (i = 0, 1, 2, \cdots, M)$, then the discrete equation can be obtained. We discrete the first order derivative and $\alpha$ order left Riemann-Liouville fractional derivative by central difference scheme and weighted and shifted Grünwald difference (WSGD) scheme ^[21] respectively.

The second-order WSGD scheme was first proposed by Tian et al. ^[21], which is a more general and flexible approach and independent on the changed fractional order. It is further applied into the numerical solution of time fractional sub-diffusion equation ^[22,23], as well as fractional Black-Schole equation in the FMLS model ^[14]. Recently Liu et al. ^[24] developed a class of second-order $\theta$ schemes based on the WSGD formula for solving the nonlinear fractional cable equation.

Using WSGD scheme, the fractional derivative can be approximated as

where

Denote $V_i^j = V(x_i, t_j), \mathbf { V }_j = \left(V_1^j, V_2^j, \cdots, V_{M-1}^j \right)^T, i = 0, 1, 2, \cdots, M, j = 0, 1, 2, \cdots, N$, then the semidiscretization equation is given by

where

Denote $\zeta = \frac { \nu } { h ^ { \alpha } }$, $\xi = \frac { r - \nu } { 2 h }$, it leads to the following semi-equation

where $\mathbf { A } = r \mathbf { I }-\zeta \mathbf { B } -\xi \mathbf { C}$, $\mathbf { B } = \left[b _ { i j } \right] _ { (M - 1) \times (M - 1) }$ defined by

$\mathbf { C } = {\rm tridiag }\{ - 1, 0, 1 \}$ and ${f}_j = (0, 0, \cdots, 0, (\xi + \zeta \omega_0^{(\alpha)}) ({ B(t_{j+1})+B(t_{j})})$. Since both the matrices $\mathbf{B}$ and $\mathbf{C}$ are Toeplitz matrices, the matrix $\mathbf{A}$ is also a Toeplitz matrix ^[25].

3.   Time stepping schemes

In order to smooth the oscillations caused by the non-smooth payoff functions and improve the accuracy, we construct time stepping schemes with Padé approximation. The discretization (2.8) in space leads to the following system of initial value problem

in a Hilbert space $X$, where $v$ denotes the initial condition $v(x)$ in (2.3). We assume the resolvent set $\rho (\mathbf{A})$ (The points $\lambda$ for which $\lambda\mathbf{I}-\mathbf{A}$ has a bounded inverse in $X$ comprise the resolvent set $\rho (\mathbf{A})$ of $\mathbf{A}$) satisfies, for some $\gamma \in (0, \frac{\pi}{2})$ ^[26],

Also, assume there exists $C\geq 1$ such that

The exact solution of (3.1) satisfies the following recurrence formula

where $\tau = T/N$, $j = 0, 1, 2, \cdots, N-1$.

Consider now its discrete analogue of the form

where $\{s_i\}_{i = 1}^{m'}\subset [0, 1]$ are the the distinct numbers selected as integral points to approximate ${\mathbf { v}(t_{j+1})}$ in formula (3.4).

The time discretization scheme (3.5) is accurate of order $q$ in time which can be described as follow.

Lemma 3.1. ^[26] The time discretization scheme $(3.5)$ is accurate of order $q$ if

and, for $0\leq l\leq q$,

or, equivalently,

It is shown in ^[26] that for the case $m' = q$ ($m'$ is the number of quadrature points and $q$ is the accuracy of the scheme), the conditions of the Lemma 1 can be achieved by choosing the rational functions $R(z)$ satisfying (3.6), selecting distinct real numbers, by Gaussian Quadrature, $\{Q_i(z)\}_{i = 1}^q$, and finally solving the system

This system (3.9) is of Vandermonde type (whose determinant is not zero), which gives the rational functions $\{Q_i(z)\}_{i = 1}^q$ as linear combinations of the terms on the right hand side of (3.9).

For the case when the number of quadrature points $m'$ is less than the order of the scheme $q$, an alternative formula similar to (3.9) is given in ^[26]. The accuracy conditions are reformulated by defining

and requiring that

and a moment condition

on the quadrature points, with $p(s) = \prod\limits_{i = 1}^{m'}\left(s-s_{i}\right)$. The formula to obtain the rational functions $\{Q_i(z)\}_{i = 1}^q$ in ^[26] is

For the rest of this chapter, we will use Padé approximation as $R(z)$ above to constructe the high-order numerical scheme.

3.1.   Padé approximations

Let $P_{n, m}(z)$ and $Q_{n, m}(z)$ be two polynomials of degree $n$ and $m$ respectively, the $(n+m)$ th order rational Padé approximation of the exponential function $e^{-z}$ can be written as

where

and

The Padé approximation $R_{n, m}(z)$ to the exponential function $e^{-z}$ is of the order $(n+m)$.

Definition 3.1. The rational approximation $R_{n, m}(z)$ of $e^{-z}$ is said to be $A$-stable if $|R_{n, m}(z)| < 1$ whenever $\Re(z) < 0$ and L-stable if in addition $|R_{n, m}(z)| \rightarrow 0$ as $\Re(z) \rightarrow - \infty$.

It is known from ^[27] that $R_{n, m}(z) = e^{-z} +O\left(|z|^{m+n+1}\right)$ as $z \rightarrow 0$, and we consider the $L$-stable $(0, 2m)$-Padé approximations and $A$-stable $(m, m)$-Padé approximations ^[28] for the exponential function $e^{-z}$.

Here for practical purpose, we are particularly interested to the following $A$-stable and $L$-stable Padé approximation of $e^{-z}$ respectively:

and

Replace the matrix exponential $e^{-\tau \mathbf{A}}$ by $(n, m)$ Padé approximation $R_{n, m}(\tau \mathbf{A})$, the recurrence relation is approximated by

which is the fully discretization of (2.2). The $\{Q_{n, m}^{(i)}(z)\}_{i = 1}^2$ are rational functions, which have same denominator as those $R_{n, m}(z)$ and $\{s_i\}_{i = 1}^{2}$ are the Gaussian points.

Using the result of equation (3.11), we still have the same accuracy when we choose the corresponding Padé approximation

which is a linear system in $Q_{n, m}^{(i)}(z)$ and could be solved easily since the matrix of the coefficients on the left is of Vandermonde's type.

3.2.   A fourth order $L$-stable method

Consider the fourth order $L$-stable Padé approximation $R_{0, 4}(z)$ with $s_1 = \frac{3-\sqrt{3}}{6}$ and $s_2 = \frac{3+\sqrt{3}}{6}$, the system reduces to

Solving the Eq (3.14), it leads to the following fourth order schemes

where

3.3.   Partial fractional form of the schemes

Both the schemes discussed above require to take inverse of higher order matrix polynomial which can cause computational difficulty due to higher power of matrix $\mathbf{A}$.

For overcoming this difficulty, Khaliq et al. ^[29], Gallopoulos and Saad ^[30] and references therein developed these schemes in a partial fraction decomposition (with complex arithmetic) that allows efficient and accurate computations on serial or parallel machines.

The partial fraction form of the rational functions $R_{n, m}(z)$ and $\{Q_{n, m}^{(i)}(z) \}_{i = 1}^{2}$ requires us to consider two cases, $n < m$ and $n = m$ for subdiagonal and diagonal Padé schemes respectively.

If $n < m$, then we have

and for the case $n = m$, the partial fraction form for $R_{n, m}(z)$ and $Q_{n, m}^{(i)}(z)$ is given by Gallopoulos and Saad ^[30]

where $R_{n, m}(z)$ as well as $Q_{n, m}^{(i)}(z)$ have $q_1$ real and $2q_2$ complex pole $c_{{i}}$ with $q_1+2q_2 = m$, and $w_j = \frac{R_{n, m}(c_j)}{Q'_{n, m}(c_j)}$ and $w_{ij} = \frac{\mathcal{N}^{(i)}_{n, m}(c_j)}{\mathcal{D}^{'(i)}_{n, m}(c_j)}$.

The polynomial $\mathcal{N}^{(i)}_{n, m}(z)$ and $\mathcal{D}^{(i)}_{n, m}(z)$ are the numerator and denominator of the function $Q_{n, m}^{(i)}(z)$ respectively.

The poles and weights for $R_{n, m}(z)$ and $Q_{n, m}^{(i)}(z)$ are:

The algorithm becomes

where

which can be solved in parallel on two machines for speedup, or on a serial machine.

4.   Convergence analysis

In this section, we prove the convergence of the proposed scheme in the case that the spatial discretization matrix $\mathbf{A}$ is a lower Hessenberg Toeplitz matrix, without the assumption of a self-adjoint operator in ^[26].

We begin with the proof of the positive definiteness of matrix $\mathbf{A}$. The matrix $\mathbf{A}$ is positive definite if and only if its symmetric part $\mathbf{W} = (\mathbf{A}+\mathbf{A}^T)/2$ is positive definite ^[31], which means its eigenvalues are all positive.

Theorem 4.1. Assume the fractional parameter $\alpha$ satisfying $1 < \alpha < 2$, the matrix $\mathbf{A} = r \mathbf { I }-\zeta \mathbf { B } -\xi \mathbf { C}$ defined in $(2.8)$ as the following

is positive definite.

Proof. Consider now the matrix $\mathbf{W} = [{(-\zeta \mathbf { B } -\xi \mathbf { C})+(-\zeta \mathbf { B } -\xi \mathbf { C})^T}]/2$ defined by

Use Gerschgorin Disk Theorem and note that $w_{ii} = {\zeta (\alpha+2)(\alpha-1)}/{2} > 0$, we only need to prove it is row diagonally dominant and column diagonally dominant. It is clear that the $i$th and the $(M-i-1)$th rows are the same. Without loss of generality, we choose $1 \leq i \leq\lceil \frac{M-1}{2} \rceil$.

For $i = 1$, we have

For $i = 2, 3, \ldots, \lceil \frac{M-1}{2} \rceil$, using the properties in (2.5), we have

Therefore, the matrix $\mathbf{W}$ is row diagonal dominant. Because of its Toeplitz structure, it is column diagonally dominant as well, and thus the matrix $\mathbf{A}$ is positive definite.

Assume the function $R(\cdot)$ is the $L$-stable $(0, 2m)$-Padé approximation of order $q$ in (3.6), and $v$ is the initial condition in (3.1), then the convergence of the scheme is established in the following two theorems.

Theorem 4.2. Assume that $\mathbf{A}$ defined in $(3.1)$ and satisfying $(3.2)$ and $(3.3)$ is diagonalizable. For the $R(\cdot)$ and $q > 0$ in $(3.6)$, there exists a constant $C > 0$ such that for $n > 1$

Proof. Let $\mathbf{A}$ have eigenvalues $\{\lambda_i \}_{i = 1}^{M-1}$ and corresponding orthonormal eigenvectors $\{w_i \}_{i = 1}^{M-1}$. Suppose $v = \sum\limits_{j = 1}^{M-1}\alpha_jw_j$, then we have

and

It follows that

Using the identity $a^n-b^n = (a-b)\sum_{j = 0}^{n-1}a^jb^{n-j-1}$, it follows that

Without the confusion, we will reuse the constant $C$ from line to line. From Theorem 4.1 we know that $\mathbf{A}$ is positive definite, thus $\Re(\lambda_j) > 0, \; j = 1, 2\ldots, M-1$. For any integer $k\geq 0$ and $l > 0$, there exists a constant $C$ such that

We also find there exists $0 < c < 1$ such that $| R(z)|\leq e^{-cz}$ for (0, 4)-Padé scheme, thus we have the following bound

Therefore,

and the proof is complete.

Use the results of Theorem 4.2, and define the spaces $\dot{H}^{s} = \mathcal{D}\left(\mathbf{A}^{s / 2}\right)$ in ^[26] with the norm as following

we complete the proof of the convergence in the following theorem.

Theorem 4.3. Suppose that $f^{(l)}(t) \in \dot{H}^{2 q-2 l}$ for $l < q$ and $t\geq 0$, then there exists a constant $C$ such that

i.e., the time discretization scheme $(3.5)$ is accurate of order $q$.

Proof. The error $\mathcal{E}^{n} = \mathbf{v}^n-\mathbf{v}(t_n)$, for $n\geq 2$, can be written as:

where $E(t) = e^{-t \mathbf{A}}$,

The error term $\mathcal{E}_0^{n}$ can be approximated by the established result from Theorem 4.2 as follows:

After inserting the term $R^{n-j-1}(\tau \mathbf{A})\mathcal{I}_kf(t_j)$ in the error term $\mathcal{E}_q^{n}$ and rearranging its terms, it can be derived as

Following the approach given in ^[26], we have the following estimate for $\mathcal{E}_1^{n}$ and $\mathcal{E}_2^{n}$,

which is bounded by the right hand side of (4.5). Also

Since Eq (4.9) can be incorporated into the right hand side of (4.10), we obtain the following estimate for the main scheme:

Combining (4.7) and (4.11), it leads to

which completes the proof.

5.   Numerical experiments

In this section, numerical performance of different Padé schemes compared with the Crank-Nicolson scheme is given for the numerical solution of both fractional digital options and fractional barrier options.

Since the non-smooth payoffs of fractional exotic options usually result in inaccurate and discontinuous solution, or serious errors when estimating the hedging parameters, e.g., Delta, Vega and Gamma values, we compare the price of the options under different schemes, as well as the Delta ^[16,32,33], which is the rate of change of the option value with respect to the asset price and can be approximated in the following way:

The order of the numerical scheme is defined as

where $V_\ast(\cdot, 0)$ denotes the exact solution at $t = 0$ and $V_{\tau}(\cdot, 0)$ denotes the numerical solution with time step $\tau$ at $t = 0$. In our experiments, the exact solution $V_\ast(\cdot, 0)$ is approximated by the numerical solution using a dense mesh with $M = N = 8192$ and we set $\tau_2/\tau_1 = 2$ to obtain the convergence order.

5.1.   A fractional digital call option

Example 5.1. Consider a fractional digital call option pricing model as follows

where $\alpha = 1.5, \; r = 0.05, \; \sigma = 0.25, \; B_u = \ln{100}, \; B_d = \ln{0.1}, \; T = 1, \; K = 50$ and $\nu = -\frac{1}{2} \sigma^{\alpha} \sec \frac{\alpha \pi}{2}$ where the payoff function is

where we take the average of playoff at $x = \ln K$ from mathematical viewpoint to restore the discontinuity in the payoff ^[15].

In Figures 1–3, the surfaces of the price and the corresponding Delta value of the fractional European digital option are plotted for the Crank-Nicolson, (2, 2)-Padé and (0, 4)-Padé schemes respectively when $N = 16$ and $M = 4096$.

From Figures 1–3, it is observed that the second-order Crank–Nicolson method suffers from oscillations with non-smooth payoff function while the (0, 4)-Padé scheme can provide reliable and smooth option values and Delta value with little oscillation.

The reason for the oscillation phenomenon at the strike price is that the time discretization grids are coarse ^[33,34]. One possible remedy is reducing the discrete step size in time direction or local mesh refinement strategy ^[35,36], which would highly increase the computational time. Here, the (0, 4)-Padé scheme could obtain the best accuracy cheaply and smooth the oscillations with relative less discrete points.

In Tables 1–3, we list the error of the numerical solution in $L_2$ and $L_\infty$ norm, as well as the corresponding order for the Crank-Nicolson scheme, (2, 2)-Padé and (0, 4)-Padé respectively when $\tau$ is varying and $M = 8192$.

From Tables 1–3, it is seen that the numerical results of digital option with both (2, 2)-Padé scheme and (0, 4)-Padé scheme have fourth-order accuracy, which is much better than those of the Crank-Nicolson scheme.

In Figure 4, we plot the curves of option price and the corresponding Delta value versus the price of the asset when $t = 0$ for the Crank-Nicolson, (2, 2)-Padé and (0, 4)-Padé scheme respectively.

From Figure 4, it is further confirmed that the (0, 4)-Padé scheme can significantly reduce the oscillations of the solution near the strike price and smooth both the price of option and the Delta value, compared with the Crank-Nicolson and (2, 2)-Padé schemes. It is possibly because of (0, 4)-Padé scheme is a high order scheme, so that it can quickly converges to the exact solution with less time layer.

5.2.   A fractional barrier put option

Example 5.2. Consider a fractional barrier put option pricing model as follows

where $\alpha = 1.5, \; r = 0.05, \; \sigma = 0.25, \; B_u = \ln{100}, \; B_d = \ln{0.1}, \; T = 1, \; K = 50, \; E = 20$ and $\nu = -\frac{1}{2} \sigma^{\alpha} \sec \frac{\alpha \pi}{2}$. The payoff function is

In Figures 5–7, the price surfaces of the fractional barrier put option and the corresponding Delta value are plotted for the Crank-Nicolson, (2, 2)-Padé and (0, 4)-Padé schemes respectively when $N = 16$ and $M = 4096$.

From Figures 5–7, it is observed that the option price and Delta value of the Crank–Nicolson method still suffer from serious oscillations near the barrier while those of the (0, 4)-Padé scheme provide the best approximate results.

It is also seen that though the (2, 2)-Padé scheme makes use of the same interpolation points with the (0, 4)-Padé scheme, the price of option and Delta values of the (2, 2)-Padé scheme still oscillate near the barrier.

Moreover, in Figure 8, we plot the curves of option price and the corresponding Delta value versus the asset price when $t = 0$ for the Crank-Nicolson, (2, 2)-Padé and (0, 4)-Padé scheme respectively.

From Figure 8, it is further verified that the (0, 4)-Padé scheme is the best one, which can significantly reduce the oscillations of the option price and the Delta value near the barrier.

In Tables 4–6, the error of the numerical solution in $L_2$ and $L_\infty$ norm, as well as the corresponding order for the Crank-Nicolson, (2, 2)-Padé and (0, 4)-Padé schemes are listed respectively when $\tau$ is varying and $M = 8192$.

From Tables 4–6, it is observed that the price of the digital option with both (2, 2)-Padé scheme and (0, 4)-Padé scheme can achieve the fourth-order accuracy, while the results of the (0, 4)-Padé scheme are more accurate than those of the (2, 2)-Padé scheme.

5.3.   A fractional double barrier call option

Example 5.3. Consider a fractional double barrier call option pricing model as follows

where $\alpha = 1.5, \; r = 0.05, \; \sigma = 0.25, \; B_u = \ln{100}, \; B_d = \ln{0.1}, \; T = 1, \; K = 20, \; E_1 = 40, \; E_2 = 70$ and $\nu = -\frac{1}{2} \sigma^{\alpha} \sec \frac{\alpha \pi}{2}$. The payoff function is

In Figures 9–11, the price surfaces of the fractional double barrier call option and the corresponding Delta value are drawn for the Crank-Nicolson, (2, 2)-Padé and (0, 4)-Padé schemes respectively with $N = 16$ and $M = 2048$.

From Figures 9–11, it is observed that both the Crank–Nicolson and (2, 2)-Padé schemes suffer from spurious oscillations near the barrier while (0, 4)-Padé approximation provides reliable option values and smooth Delta value.

In Figure 12, we plot the curves of option price and the corresponding Delta value versus the asset price when $t = 0$ for Crank-Nicolson, (2, 2)-Padé and (0, 4)-Padé scheme respectively.

From Figure 12, it is further confirmed that the (0, 4)-Padé scheme can significantly reduce the oscillation of the price of the option near the barrier and smooth the corresponding Delta value.

In Tables 7–9, the error of the numerical solution in $L_2$ and $L_\infty$ norm, as well as the corresponding order for the Crank–Nicolson, (2, 2)-Padé and (0, 4)-Padé schemes are listed respectively when $\tau$ is varying and $M = 8192$.

From Tables 7–9, it is seen that the numerical price of the fractional double barrier call option for both (2, 2)-Padé scheme and (0, 4)-Padé scheme can achieve fourth-order accuracy while the numerical results of the Crank–Nicolson scheme only obtains the second-order accuracy. Among the three schemes, the (0, 4)-Padé scheme is the most accurate one, which requires less discrete points than the other two schemes to achieve the same accuracy.

6.   Conclusions

A class of fourth order Padé schemes for pricing fractional exotic options under FMLS model are proposed and studied, which make use of the 2nd-order weighted and shifted Grünwald difference scheme in space direction and the 4th-order Padé schemes in time direction. The convergence of the Padé schemes are proved in detailed under the FMLS model. Numerical experiments on fractional digital option and fractional barrier options are given to verify the 4th-order precision, and show that the (0, 4)-Padé scheme can significantly reduce the oscillations of the solution near the strike price or barrier, smooth the Delta value and save the computational time.

Acknowledgements

The authors are very grateful to the referees for their constructive comments and valuable suggestions, which greatly improved the original manuscript of the paper. This work is supported by the National Natural Science Foundation of China (No. 11971354).

Conflict of interest

The authors declare there is no conflicts of interest.