2 Equity warrants under the GMFBM model

In this section, the formulation of the pricing of equity warrants is considered. In particular, a nonlinear pricing system in terms of the observable variables governing the price of equity warrants is established. To include most of the FBM models as special cases, the GMFBM model, which is a linear combination of a countable number of BMs and FBMs [Reference Thäle20], is adopted in the current work. For the completeness of the paper, an introduction to this model is provided in Appendix A.

Consider a European equity warrant with a call feature, issued by a firm for its own stocks. Suppose the firm has only two forms of financing, namely, N shares of stocks and M outstanding equity warrants. Each warrant gives its holder the right but not the obligation to receive k shares of stocks with payment X at the expiry time T. Let $w(V,t)$ be the current value of each equity warrant on a firm with value V.

Under the GMFBM model, we assume that $V_t$ , with t being the current time, satisfies the following stochastic partial differential equation (SDE),

(2.1) $$ \begin{align} dV_t=\mu V_tdt + \sigma_{V} V_t\diamond dZ_t, \end{align} $$

where $\diamond $ is the so-called Wick product, $\mu $ is the expected rate of return of the firm value and $\sigma _V$ is the standard deviation of the firm value. There are two main reasons why the Wick product instead of the classical It $\widehat {o}$ theory is used in (2.1). First, the mean of the stochastic integral defined based on the Wick product is zero, which is useful for both theoretical development and practical applications. Second, the It $\widehat {o}$ theory fails to apply under the GMFBM because $V_t$ is neither a Markov nor a semimartingale, unless under some particular parameter settings, as pointed out in many previous studies [Reference Björk and Hult4, Reference Chen and He7, Reference Lin16]. We further remark that when the Wick product is adopted, a new concept “wickbitrage” is defined, which is, however, not identical to the concept of “arbitrage” in the intuitive sense [Reference Björk and Hult4, Reference Cheridito10]. Nevertheless, Chen and He showed in their recent work [Reference Chen and He7] that the market modelled under the GMFBM could still be arbitrage free under certain parameter settings (see [Reference Chen and He7] for details).

Under the assumptions specified by Xu et al. [Reference Xu, Xiao, Zhang and Zhang23] for an equity warrant, except that the firm value now follows a GMFBM rather than a mixed fractional Brownian motion (MFBM), it can be deduced that w satisfies the PDE

(2.2) $$ \begin{align} \frac{{\partial} w}{{\partial} t}+\sigma^2_VV^2\sum^{I}_{i=1}H_i\alpha_it^{2H_i-1}\frac{{\partial}^2 w}{{\partial} V^2}+rV\frac{{\partial} w}{{\partial} V}-rw=0. \end{align} $$

Considering the dilution effect, the terminal value of the price should be

$$ \begin{align*} w(V,T)=\frac{1}{N+Mk}\max(kV-NX,0),\end{align*} $$

as stated by Ukhov [Reference Ukhov21]. On the other hand, the boundary conditions along the V-direction are similar to those of European call options. When the firm value becomes extremely large, the warrant will definitely be exercised at the expiry, and its current value should be

$$ \begin{align*} \lim_{V{\rightarrow}\infty}w(V,\tau)=\frac{kV-NXe^{-r(T-t)}}{N+Mk}.\end{align*} $$

On the contrary, when $V{\rightarrow } 0$ , which implies that the firm has no value, we have $ \lim _{V{\rightarrow } 0}w(V,t)=0$ . Furthermore, by considering the dilution effect, Ukhov pointed out in [Reference Ukhov21] that the unobservable variables, namely, V and $\sigma _V$ , can be estimated by observable variables through

$$ \begin{align*}V=Mw+NS \quad \text{and}\quad \sigma_V=\frac{S\sigma_S}{V({\partial} S/{\partial} V)}, \end{align*} $$

where S is the stock price of the firm and $\sigma _S$ is its volatility. From the expressions of V and $\sigma _V$ , it is clear that $\sigma _V$ is a function of w rather than a constant, and thus (2.2) is a nonlinear PDE.

Note that, in most of the algorithms documented for the pricing of equity warrants, the authors first derived closed-form analytical expression of w with respect to V by assuming that $\sigma _V$ is a constant. Then, by solving a nonlinear system, they computed $w(V,\tau )$ using the observable S and estimated $\sigma _V$ . However, from a PDE point of view, the governing equation (2.2) is nonlinear and cannot be solved analytically first. This inspires us to solve the warrant price in terms of the observable stock price directly, which is also one of the innovations of the current work.

Now, let $\tilde {w}(S,\tau )=w(V,t)$ , where $\tau =T-t$ . According to the chain rule and the relationship among V, w and $S,$

$$ \begin{align*} \frac{{\partial} w}{{\partial} t}=-\frac{{\partial} \tilde{w}}{{\partial} \tau},\quad \frac{{\partial} w}{{\partial} V}=\frac{{\partial} \tilde{w}}{{\partial} S}\frac{1}{N+M({\partial} \tilde{w}/{\partial} S)},\quad \frac{{\partial}^2 w}{{\partial} V^2}=\frac{{\partial}^2}{\tilde{w}}{{\partial} S^2}\Big(\frac{{\partial} S}{{\partial} V}\Big)^2\frac{1}{1+(M/N)({\partial} \tilde{w}/{\partial} S)}.\end{align*} $$

Substituting the above partial derivatives into (2.2) as well as the boundary conditions mentioned above, and dropping all the tildes for simplicity, we have the following nonlinear PDE system:

(2.3) $$ \begin{align} \begin{cases} Lw=0\\ w(S,0)=\max(kS-X,0)\\ w(0,\tau)=0\\ \lim_{S{\rightarrow}\infty}w(S,\tau)=kS-Xe^{-r(T-t)}, \end{cases} \end{align} $$

where $Lw$ is defined as

$$ \begin{align*} Lw =\frac{{\partial} w}{{\partial} \tau}-A(w,S,\tau)\frac{{\partial}^2 w}{{\partial} S^2}-B(w,S)\frac{{\partial} w}{{\partial} S}+rw\end{align*} $$

with

$$ \begin{align*} A(w,S,\tau)=\sigma^2_SS^2 \sum^{I}_{i}H_i\alpha_i(T-\tau)^{2H_i-1}\frac{1}{1+(M/N)({\partial} w/{\partial} S)}\end{align*} $$

and

$$ \begin{align*} B(w,S)=r(Mw+NS)\frac{1}{N+M({\partial} w/{\partial} S)}.\end{align*} $$

According to Agliardi et al. [Reference Agliardi, Popivanov and Slavona1], there exists a unique solution to (2.3). Moreover, the solution $w(S,\tau )\in L^2(H^1(\Omega ), (0,T)]\cap C^0 [L^2(\Omega ),[0,T]]$ , where $\Omega =[0,+\infty ).$ Financially, due to the relationship between the firm's stocks and equity warrants, it is reasonable to infer that the price of the equity warrant with a call feature will increase when the underlying stock price becomes larger. Therefore, we have $A(w,S,\tau )>0$ and $B(w,S)>0$ . It should also be remarked that, after the transformation of variables, the dilution effect has now become the nonlinear terms contained in the governing PDE. It is these nonlinear terms that have made the pricing of equity warrants very difficult, even numerically. In the next section, a hybrid finite difference method will be constructed to solve (2.3) effectively.

3 Numerical scheme

As mentioned in the last section, the PDE system governing the price of an equity warrant is nonlinear. It is usually impossible to derive closed-form analytical solutions to nonlinear PDE systems, and numerical approaches are preferred. Therefore, proposing an efficient numerical method to solve for (2.3) is the aim of this section.

3.1 The IMEX method

To begin with, we truncate the domain $[0,\infty )$ into a finite one as $[0,S_{\max })$ . According to Wilmott et al.’s [Reference Wilmott, Dewynne and Howison22] estimate that the upper bound of the stock price is typically three or four times the strike price, it is reasonable to set $ S_{\max }=4X/k$ . Similar to the pricing of option derivatives, such a truncation of the domain will only lead to a negligible error in the computed warrant prices [Reference Kangro and Nicolaidesn14]. On the other hand, to ensure the stability of the discrete scheme, we use uniform meshes with P elements and Q elements on the spatial interval $[0, S_{\max }]$ and the time interval $[0,T]$ , respectively. Clearly, the mesh sizes $\Delta S$ and $\Delta \tau $ are, respectively,

$$ \begin{align*}\Delta S=\frac{S_{\max}}{P} \quad \text{and}\quad \Delta \tau=\frac{T}{Q}.\end{align*} $$

As mentioned previously, the current PDE system is nonlinear, and difficult to be solved numerically. To tackle its nonlinearity, we use the implicit Euler scheme for the time direction and approximate the coefficients of the PDE explicitly. In addition, the upwind scheme is adopted to approximate the first-order spatial derivative to avoid the numerical instability resulting from very small t. Now, denoting the value on an arbitrage grid point $(p,q)$ as $W^q_p=w(p\Delta S,q\Delta t)$ with $p=0,1\ldots P$ and $q=0,1\ldots Q$ , the finite difference system written on that particular grid point can be summarized as

(3.1) $$ \begin{align} \begin{cases} L_{P,Q}W^{q}_p=0, & 1\leq p\leq P-1, \quad 1\leq q\leq Q-1 \\[.1cm] W^0_{p}=\max(kp\Delta S-X,0), &1\leq p\leq P-1\\[.1cm] W^q_{0}=0, &1\leq q\leq Q\\[.1cm] W^q_{P}=k S_{\max}-Xe^{-r(Q-q)\Delta\tau}, &1\leq q\leq Q, \\ \end{cases} \end{align} $$

where

(3.2) $$ \begin{align} L_{P,Q}W^{q}_p&=\frac{W^{q+1}_p-W^{q}_p}{\Delta\tau}-A^q_p\frac{W^{q+1}_{p+1}-2W^{q+1}_{p}+W^{q+1}_{p-1}}{(\Delta S)^2}\nonumber\\&\quad-B^q_p\frac{W^{q+1}_{p+1}-W^{q+1}_{p}}{\Delta S}+rW^{q+1}_p,\end{align} $$

with

$$ \begin{align*} A^q_p=\frac{N\sigma^2_sS^2_p\sum^I_{i=1}H_i\alpha_i(T-\Delta\tau q)^{2H_i-1}}{N+M(W^q_{p+1}-W^q_{p-1})/2\Delta S}>0\end{align*} $$

and

$$ \begin{align*} B^q_p=\frac{r(MW^q_p+Np\Delta S)}{N+M(W^q_{p+1}-W^q_{p-1})/2\Delta S}>0.\end{align*} $$

We remark that, by approximating the coefficients explicitly, the discrete system (3.1) becomes linear when solving for the warrant price for the $(q+1)$ th time step. Such a linear discrete system can be solved directly step by step until the Qth time step is reached and the desired warrant price can be obtained. In the following subsection, an error estimation for the current method is provided.

3.2 Error estimation

An error analysis is usually an indispensable part of designing any robust numerical approach. In this subsection, we have managed to provide an error estimation for the current IMEX method used to solve for the price of equity warrants under the GMFBM model. Our proof is based on the discrete maximum principle, which is a useful tool for estimating the computational error [Reference Cen and Chen6, Reference Chen, Yan, Lian and Zhang8]. For simplicity, the following notations is adopted.

$$ \begin{align*} \begin{array}{ll} \Omega_P=\{p\in \mathbb{Z}\,\mid 0\leq p\leq P\}, &{\partial}\Omega_P=\{0,P\},\\[3pt] \Omega_Q=\{q\in \mathbb{Z}\,\mid 0\leq q\leq Q\}, &{\partial}\Omega_Q=\{0\},\\[3pt] \Omega =\Omega_P\times \Omega_Q, &{\partial}\Omega={\partial}\Omega_P\times {\partial}\Omega_Q. \\[.2cm] \end{array}\end{align*} $$

First, we will show, through the following lemma, that our differential operator $L_{P,Q}$ defined by (3.2) satisfies a discrete maximum principle.

Lemma 3.1. (Discrete maximum principle) Suppose that $\{U^{q}_{p}\}$ are a set of discrete values satisfying $U_{p}^{q}\geq 0$ for $(p,q)\in {\partial }\Omega $ , and $L_{P,Q}U^{q}_{p}\geq 0$ for $(p,q)\in \Omega \setminus {\partial } \Omega $ . Then $U^{q}_{p}\geq 0$ holds true for $(p,q)\in \Omega $ .

Proof. Let $R_{(P-1)\times (P-1)}$ be the matrix associated with the discrete operator $L_{P,Q}$ . For $q\in \Omega _Q\setminus {\partial }\Omega _Q$ ,

$$ \begin{align*} \begin{array}{ll} R_{p,p-1}=-\dfrac{\Delta\tau}{(\Delta S)^2}A^q_p, &2\leq p\leq P-1,\\ R_{p,p}=1+r\Delta\tau+\dfrac{2\Delta\tau}{(\Delta S)^2}A^q_p+\dfrac{\Delta\tau}{\Delta S}B^q_p, &1\leq p\leq P-1,\\ R_{p,p+1}=-\dfrac{\Delta\tau}{(\Delta S)^2}A^q_p-\dfrac{\Delta\tau}{\Delta S}B^q_p, &1\leq p\leq P-2. \end{array} \end{align*} $$

A direct calculation shows that

$$ \begin{align*}R_{p,p-1}+R_{p,p}+R_{p,p+1}=1+r\Delta\tau> 0.\end{align*} $$

Therefore, $R_{(P-1)\times (P-1)}$ is diagonally dominant and has nonpositive off diagonal entries. Hence, this matrix is an irreducible M-matrix. The result of our lemma can then be obtained.

Next, by using the Taylor expansion theory, the following truncation error estimate can be obtained.

Lemma 3.2. Let $w(S,\tau )$ be a smooth and monotonically increasing function defined on $0\leq S\leq S_{{\mathrm \max }}$ and $0\leq \tau \leq T$ . Then the estimate for the truncation error is

$$ \begin{align*} |L_{P,Q}w(S_p,\tau_{q})-Lw(S_p, \tau_q)| \leq \widetilde{\widetilde{C}}(\Delta\tau+\Delta S), \end{align*} $$

where $ \widetilde {\widetilde {C}}$ is a positive constant independent of the mesh.

Proof. First, according to the definitions of A, B, $A^q_p$ and $B^q_p$ , it is not difficult to show that these are nonnegative and bounded by a constant C that is independent of the mesh. Therefore,

$$ \begin{align*} |A(w(S_p,\tau_{q}), S_p,\tau_q)-A^q_p|= C_1|A(w(S_p,\tau_{q}),S_p,\tau_{q})A^q_p |(\Delta S)^2\leq C_1 C^2(\Delta S)^2\end{align*} $$

and

$$ \begin{align*}|B(w(S_p,\tau_{q}), S_p)-B^q_p|= C_1|B(w(S_p,\tau_{q}),S_p)B^q_p|(\Delta S)^2\leq C_1 C^2(\Delta S)^2,\end{align*} $$

where

$$ \begin{align*} C_1=\frac{1}{6}M\bigg|\frac{{\partial}^3 w}{{\partial} S^3}\bigg||_{(S_{\xi},\tau_q)}| \quad \text{with } S_{\xi}\in [S_{(p-1)},S_{(p+1)}]. \end{align*} $$

According to the Taylor expansion theory,

$$ \begin{align*} &|L_{P,Q}w(S_p,\tau_{q})- Lw(S_p, \tau_q) | \\&\quad\leq |A(w(S_p,\tau_{q}),S_p,\tau_q)-A^q_p|\cdot\bigg|\frac{{\partial}^2 w}{{\partial} S^2}\bigg|_{(S_p,\tau_q)}+\frac{1}{12}(\Delta S)^2\frac{{\partial}^4 w}{{\partial} S^4}\bigg|_{(S_1,\tau_q)}\bigg| \\&\qquad +|B(w(S_p,\tau_{q}),S_p)-B^q_p|\cdot\bigg|\frac{{\partial} w}{{\partial} S}\bigg|_{(S_p,\tau_q)}+\frac{1}{2}\Delta S\frac{{\partial}^2 w}{{\partial} S^2}\bigg|_{(S_2,\tau_q)}\bigg| \\&\qquad +\frac{1}{2}\Delta\tau\bigg|\frac{{\partial}^2 w}{{\partial}\tau^2}\bigg|_{(S_p,\tau_1)}\bigg|+|A(w(S_p,\tau_{q}),S_p,\tau_q)|\cdot\bigg|\frac{1}{12}(\Delta S)^2\frac{{\partial}^4 w}{{\partial} S^4}\bigg|_{(S_1,\tau_q)}\bigg| \\&\qquad +|B(w(S_p,\tau_{q}),S_p)|\cdot\bigg|\frac{1}{2}\Delta S\frac{{\partial}^2 w}{{\partial} S^2}\bigg|_{(S_2,\tau_q)}\bigg| \\&\quad\leq \bigg(4C_1\widetilde{C}C^2+\frac{1}{12}\widetilde{C}C\bigg)(\Delta S)^2+\frac{1}{2}\widetilde{C}\Delta\tau+\frac{1}{2}C\widetilde{C}\Delta S\leq\widetilde{\widetilde{C}}(\Delta S+\Delta\tau), \end{align*} $$

where

$$ \begin{align*}S_1\in[S_{(p-1)},S_{(p+1)}], \quad S_2\in[S_{p},S_{(p+1)}], \quad \tau_1\in[\tau_{q},\tau_{(q+1)}], \end{align*} $$

$$ \begin{align*} \widetilde{C}=\max_{(p,q)\in\Omega}\bigg\{\bigg|\frac{{\partial}^2 w}{{\partial} S^2}\bigg|,\bigg|\frac{{\partial}^4 w}{{\partial} S^4}\bigg|,\bigg|\frac{{\partial}^2 w}{{\partial} \tau^2}\bigg| \bigg\} \quad \text{and} \quad \widetilde{\widetilde{C}}=\max \bigg\{\frac{1}{2}\widetilde{C},4C_1C^2\widetilde{C}+\frac{1}{12}\widetilde{C}C\bigg\}.\\[-41pt] \end{align*} $$

Now, let $w(S,\tau )$ be the solution of (2.3) and let $W^{q}_{p}$ be the solution of (3.1). We have the following error estimate for the current IMEX method.

Theorem 3.3. For the current numerical method, the error estimate

$$ \begin{align*}|w(S_{p},\tau_{q})-W_{p}^{q}|\leq \overline{C}(\Delta \tau+\Delta S), \end{align*} $$

holds true, where $(p,q)\in \Omega $ , and $\overline {C}$ is a constant independent of $\Delta \tau $ and $\Delta S$ .

Proof. We define a function $F^q_p$ on $\Omega $ as

$$ \begin{align*} F_{p}^{q}=\frac{\widetilde{\widetilde{C}}}{r}[2\Delta\tau+2\Delta S]>0. \end{align*} $$

A direct calculation shows that $L_{P,Q}F^q_p=2\widetilde {\widetilde {C}}[\Delta \tau +\Delta S]$ . Therefore, for $(p,q)\in \Omega $ ,

$$ \begin{align*} L_{P,Q}[w(S_{p},\tau_{q})-W_{p}^{q}-F^q_p]&= L_{P,Q}[w(S_{p},\tau_{q})]- L_{P,Q}[F^q_p] \\ &= L_{P,Q}[w(S_{p},\tau_{q})]-L[w(S_{p},\tau_{q})]- L_{P,Q}[F^q_p] \\ &\leq |L_{P,Q}[w(S_{p},\tau_{q})]-L[w(S_{p},\tau_{q})]|- L_{P,Q}[F^q_p] \\ &=-\widetilde{\widetilde{C}}[\Delta\tau+\Delta S]<0. \end{align*} $$

On the other hand, for $(p,q)\in \Omega $ ,

$$ \begin{align*} L_{P,Q}[w(S_{p},\tau_{q})-W_{p}^{q}+F^q_p]&= L_{P,Q}[w(S_{p},\tau_{q})]+L_{P,Q}[F^q_p] \\ &= L_{P,Q}[w(S_{p},\tau_{q})]-L[w(S_{p},\tau_{q})]+L_{P,Q}[F^q_p] \\ &\leq -|L_{P,Q}[w(S_{p},\tau_{q})]-L[w(S_{p},\tau_{q})]|+L_{P,Q}[F^q_p] \\ &=\widetilde{\widetilde{C}}[\Delta\tau+\Delta S]\geq 0. \end{align*} $$

Now, we consider the nodes $(p,q)\in {\partial }\Omega $ . Since the boundary conditions of the current problem is of Dirichlet type, it is not difficult to show that, for $(p,q)\in {\partial }\Omega $ ,

$$ \begin{align*} w(S_{p},\tau_{q})-W_{p}^{q}+F^q_p>0 \quad \text{and}\quad w(S_{p},\tau_{q})-W_{p}^{q}-F^q_p<0.\end{align*} $$

Therefore, according to the discrete maximum principle established in Lemma 3.1,

$$ \begin{align*}\max_{(p,q)\in \Omega}|w(S_{p},\tau_{q})-W_{p}^{q}|\leq \max_{(p,q)\in \Omega}|F^q_p|\leq \overline{C}(\Delta\tau+\Delta S),\end{align*} $$

where $\overline {C}$ is a constant independent of the mesh. This completes the proof of the current theorem.

4 Numerical examples and discussions

In this section, numerical results will be provided together with some useful discussions. According to the issues to be addressed, three numerical experiments regarding the convergence, accuracy and efficiency of the current method are conducted in this section.

Numerical experiment 1: In this numerical experiment, we investigate the error estimate and convergence rate of the current method. The parameters used in this example are $I=1$ , $\alpha _1=1$ , $ H_1=1/2$ , $N=100$ , $M=50$ , $k=1$ , $\sigma _S=0.75$ , $r=4.48\%$ and $T-t=3$ (years).

The method we adopt to obtain the convergence rate is quite standard [Reference Chen, Yan, Lian and Zhang8]. To obtain the convergence rate in the $\tau $ direction, we fix the grid size in the S direction to be fairly small, that is, $ \Delta S=S_{\max }/2000$ , and vary the number of time intervals from $40$ to $640$ . The difference $e^{P,Q}$ reported in Table 1 is measured by the discrete maximum norm associated with P uniform elements in the S direction and Q uniform elements in the $\tau $ direction, and is defined as

$$ \begin{align*} e^{P,Q}=\max_{(p,q)\in\Omega}|W^{q}_p-w_{\text{exact}}(S_p,\tau_q)|, \end{align*} $$

where $w_{\text {exact}}$ refers to the approximated solution calculated with very fine grid sizes, namely, $ \Delta S=S_{\max }/2000$ and $ \Delta \tau =T/4000$ . Since there is no analytical solution available for equity warrants under the GMFBM model, $w_{\text {exact}}$ is adopted as the benchmark solution in all numerical experiments in this section. The rate appearing in this table is then calculated from

$$ \begin{align*} R^{P,Q}=\log_2\frac{e^{P,Q}}{e^{P,2Q}}.\end{align*} $$

From Table 1, it is clear that, for fixed sufficiently large P value, the rate is very close to $1$ , indicating that our method is indeed first-order convergent in the time direction. Similarly, when we fix the time step size to be $ \Delta \tau =T/4000$ , and gradually increase the grid numbers in the S direction, we find that the rate approaches $1$ , as shown in Table 2. Therefore, a first-order convergence is also achieved in the spatial direction, which agrees well with the theoretical result shown in Theorem 3.3.

Table 1 Convergence rate in the $\tau $ direction.

Table 2 Convergence rate in the S direction.

Numerical experiment 2: In this example, we compare our results with those listed in the literature. This example consists of five warrants considered in [Reference Xu, Xiao, Zhang and Zhang23] from August 2005 to May 2008, with basic information shown in Table 3. The comparison of our results (denoted by the superscript 2), the corresponding values computed with the method proposed in [Reference Xu, Xiao, Zhang and Zhang23] (denoted by the superscript 1) and their actual market prices (denoted by the subscript “Act”) are displayed in Table 4. In addition, in this table, the subscript “B–S” means that the values are determined under the B–S model, whereas the subscript “FBM” means that the values are determined under the FBM model. Note that both the B–S model and the FBM model are special cases of the GMFBM model. From this table, one can clearly observe that our method produces results that are overall closer to real warrant prices than those determined by the algorithm proposed in [Reference Xu, Xiao, Zhang and Zhang23]. Moreover, when our method is applied to the FBM model and the B–S model, the FBM prices are overall closer to the real market prices than the corresponding B–S prices. This not only demonstrates the reliability of the current method, but also suggests that the FBM is a better choice than the BM to describe the evolution of the firm value.

Table 3 Basic information for five warrants. Note that other parameters are $k=1$ , $r_1=2.25\%$ and $r_2=4.14\%$ , where $r_1$ and $r_2$ are risk-free interest rates for one and two years, respectively.

Table 4 Comparison of the different approaches with actual market prices.

Numerical experiment 3: In this example, we examine the efficiency of the current method. Three FBM models are adopted, with the following parameters.

$$ \begin{align*} I=1, &\quad \alpha_1=1, \quad H_1=1/2 \quad \text{for the B--S model}, \\ I=1,& \quad \alpha_1=1, \quad H_1=0.628 \quad \text{for the FBM model, } \\ I=2, & \quad \alpha_1=1, \quad \alpha_2=0.3, H_1=1/2, H_2=0.88 \quad \text{for the MFBM model}. \end{align*} $$

Other parameters are $N=100$ , $M=50$ , $k=1$ , $\sigma _S=0.75$ , $r=4.48\%$ and $T-t=3$ (years). Furthermore, the “error” in this experiment is defined as

$$ \begin{align*} \text{error}=\frac{\parallel\! w(S,t)-w_{\text{exact}}(S,t)\!\parallel_2}{\parallel\! w_{\text{exact}}(S,t)\!\parallel_2}. \end{align*} $$

As shown clearly in Figure 1, the accuracy measured by “error” varies inversely with the efficiency measured by the “CPU time”. When the same level of accuracy is required, the computational time will increase if more FBMs are involved. Most impressively, one can observe that, for the current method, a high computational efficiency can be achieved while a satisfactory accuracy can still be maintained. For example, this method can produce a result within one second with error less than 0.02%. This level of accuracy and efficiency would certainly meet the practical needs of market practitioners.

Figure 1 Accuracy versus efficiency.