AN IMEX-BASED APPROACH FOR THE PRICING OF EQUITY WARRANTS UNDER FRACTIONAL BROWNIAN MOTION MODELS

Abstract

In this paper, the pricing of equity warrants under a class of fractional Brownian motion models is investigated numerically. By establishing a new nonlinear partial differential equation (PDE) system governing the price in terms of the observable stock price, we solve the pricing system effectively by a robust implicit-explicit numerical method. This is fundamentally different from the documented methods, which first solve the price with respect to the firm value analytically, by assuming that the volatility of the firm is constant, and then compute the price with respect to the stock price and estimate the firm volatility numerically. It is shown that the proposed method is stable in the maximum-norm sense. Furthermore, a sharp theoretical error estimate for the current method is provided, which is also verified numerically. Numerical examples suggest that the current method is efficient and can produce results that are, overall, closer to real market prices than other existing approaches. A great advantage of the current method is that it can be extended easily to price equity warrants under other complicated models.

1 Introduction

Warrants are derivatives that give holders the right, but not the obligation, to buy or sell the underlying at a certain date for a prescribed price. Warrants can be classified into American and European styles with call or put features, depending on when and how to exercise the contracts. According to the types of issuers, warrants can also be divided into covered and equity warrants. Covered warrants are usually issued by financial institutions or dealers. Since this kind of warrants will not raise the number of the firm’s stocks after the expiration date, the pricing of them is almost the same as that of ordinary options. In contrast, equity warrants are usually issued by the listed firms. When an equity warrant is exercised, the proceeds will be part of the firm’s equity, which affects the value of all other claims of the firm, including its outstanding warrants. It is this dilution effect that has made the pricing of equity warrants totally different from that of the covered warrants. Owing to the large trading volumes of warrants being traded around the world, it is important to ensure that the warrants, especially the equity warrants, can be priced both accurately and efficiently.

The pricing of equity warrants can be dated back to the mid-1970s when the famous Black–Scholes (B–S) formula was established. Black and Scholes [5] showed in their paper that their formula can be modified for price equity warrants. Their formulation was, however, inaccurate, because the differences between options and equity warrants, especially the dilution effect, were totally neglected. Fortunately, in the 1980s, several studies on the pricing of equity warrants were carried out. By noticing that the distribution of stock returns will be changed after the equity warrant is exercised, Galai and Schneller proposed a modified B–S model to incorporate the dilution effect [11]. Their work was further extended by a number of authors who focused on various ways of correction for dilution (see [13, 15, 18, 21] and the references therein). In some of the work mentioned above [13, 15], the authors required the value of the firm and its volatility to be known in advance. This is impossible and, moreover, in the case when there are warrants outstanding, the value of the firm is, in fact, a function of the value of the warrant. To overcome this difficulty, Schulz and Trautmann proposed a method to price equity warrants by using stock price and its volatility [18]. Ukhov also developed an algorithm for the pricing of equity warrants by using observable variables [21].

Owing to the limitations of the B–S assumption on the underlying asset, a number of authors concentrated on applying different fractional Brownian motion (FBM) models to the finance field. They pointed out that these models are able to capture the long-range dependence of the underlying. For example, Almani et al. introduced a new FBM with a two-variable Hurst exponent, and found that the model is more precisely matched to the real values of the rate of the stock price [2]. Cheng and Xu considered the pricing of vulnerable options under a mixed FBM model with jumps[9]. By introducing an FBM to the constant elasticity of variance (CEV) model, Araneda and Bertschinger proposed a sub-fractional CEV model and considered the pricing of options underneath [3]. Han et al. raised a stochastic volatility model driven by both an FBM and a standard Brownian motion (BM), and obtained an analytical solution for the European option price [12]. For a variety of applications of FBMs in the field of derivative pricing, interested readers can refer to [17, 23] and the references therein.

Note that, for most of the algorithms documented for the pricing of equity warrants, the authors derived a closed-form analytical expression for the warrant prices with respect to the firm value under the assumption that the volatility of the firm is constant. Then, they expressed the volatility of the firm in terms of the warrant price and stock volatility. However, from a partial differential equation (PDE) point of view, the warrant price satisfies a nonlinear PDE because the volatility is also a function of the warrant price, and it cannot be expressed explicitly as a function of the firm value. On the other hand, even if these documented methods can produce a good approximated price under certain parameter settings, the applications of these methods are still limited. Only if the closed-form expression of the warrant price with respect to the firm value can be found, can the methods be applied. In most cases, however, such an expression is difficult to derive, especially under rather complicated models.

In this paper, we consider the pricing of equity warrants under a class of FBM models numerically. For convenience, we adopt the so-called generalized mixed fractional Brownian motion (GMFBM) model, because this model includes most of the FBM models as its special cases. This model was introduced by Th $\ddot {{\mathrm a}}$ le [20], and further applied to the option pricing field by Chen et al. [7, 8]. The contribution of the current work mainly includes both practical and theoretical aspects. Practically, unlike most authors in the literature, we solve the warrant price in terms of the observable stock price directly by an implicit-explicit (IMEX) finite difference method. By comparison with the existing methods, the current method can not only produce results that are, overall, closer to real market prices, but also very promising to be extended to price warrants under other complicated models, under which the closed-form expression of European option prices cannot be easily derived. Theoretically, we show that the coefficient matrix associated with the current approach is an M-matrix, which ensures the stability of the method in the maximum-norm sense. Most remarkably, a sharp error estimate of the current method is obtained, which suggests that our method is first-order convergent in both the time and spatial directions. Numerical results also agree with the theoretical statement.

The rest of the paper is organized as follows. In Section 2, we derive the PDE system governing the price of equity warrants under the GMFBM model in terms of observable variables. In Section 3, we introduce the IMEX numerical method in detail and also derive a sharp error estimate for the current method. In Section 4, numerical experiments are conducted to test the theoretical error estimate, and useful discussions are also provided. Concluding remarks are given in Section 5.

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 [20], 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 [4, 7, 16]. 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 [4, 10]. Nevertheless, Chen and He showed in their recent work [7] that the market modelled under the GMFBM could still be arbitrage free under certain parameter settings (see [7] for details).

Under the assumptions specified by Xu et al. [23] 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 [21]. 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 [21] 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. [1], 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 [22] 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 [14]. 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 [6, 8]. 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 [8]. 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 [23] 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 [23] (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 [23]. 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.

5 Conclusion

This paper considers the numerical pricing of equity warrants under the GMFBM model. A nonlinear PDE system governing the equity warrant price in terms of observable variables is derived first, and then solved effectively by an IMEX finite difference method. Compared with the documented methods for the pricing of equity warrants, the current method is based on a more rigorous mathematical deduction and can produce results that are overall closer to real market prices. It is also theoretically shown that the coefficient matrix associated with the current method is an M-matrix, which ensures its stability in the maximum-norm sense. Most remarkably, a sharp error estimate for the current method is also provided, which suggests that the proposed method is first-order convergent in both the time and spatial directions. Based on the current work, at least two future research directions can be expected. First, it is promising to be able to extend the current method to price warrants under other complicated models, under which closed-form expressions of European options cannot be easily derived. Second, application of the current method to the pricing of other similar financial derivatives is also plausible.