1. Introduction
In financial markets, the well-known Black–Scholes model has been widely accepted, but the assumption of the ideal conditions in this model, such as liquidity, less friction and completeness, is clearly unrealistic. Therefore, in recent years, the linear Black–Scholes model has been replaced by a nonlinear Black–Scholes equation with a nonlinear volatility function, where transaction costs, market liquidity or volatility uncertainty are not neglected. The market liquidity studied in this paper is an issue of very high concern in financial risk management. Therefore, illiquid markets and large trader effects have been modelled by various researchers. Among them, Backstein and Howison [Reference Bakstein and Howison6] developed a parametrized model for liquidity effects arising from the trading in assets. Bank and Baum [Reference Bank and Baum7] introduced a general continuous-time model for an illiquid financial market. In the study [Reference Bordag and Chmakova12], the authors used families of explicit solutions to test numerical schemes solving a nonlinear Black–Scholes equation. Cetin et al. [Reference Cetin, Jarrow and Protter14] extended the classical approach by formulating a new model that considered illiquidities. A general nonlinear model of illiquid markets with feedback effects was considered in the study by Federov and Dyshaev [Reference Federov and Dyshaev28]. Another remarkable study was presented by Frey [Reference Frey29]. In that study, standard derivative pricing theory, which is the assumption of agents acting as price takers on the market for the underlying asset, was investigated. They characterized the solution to the hedge problem in terms of a nonlinear partial differential equation, and presented results on the existence and uniqueness of this equation. Gulen and Sari [Reference Gulen and Sari33] defined a nonclassical numerical method to capture the behaviour of the nonlinear option pricing model in illiquid markets where the implementation of a dynamic hedging strategy affects the price of the underlying asset. Jarrow [Reference Jarrow36] analysed market manipulation trading strategies by large traders in a securities market. Platen and Schweizer [Reference Platen and Schweizer44] developed a diffusion model for stock prices explicitly incorporating the technical demand induced by hedging strategies starting from a microeconomic equilibrium approach. Another work by Schönbucher and Wilmott [Reference Schönbucher and Wilmott45] analysed the influence of dynamic trading strategies on the prices in financial markets. Sircar and Papanicolaou [Reference Sircar and Papanicolaou49] studied the nonlinear partial differential equation for the price of the derivative by perturbation methods, by numerical methods which are easy to use and can be implemented efficiently and by analytical methods.
How the hedging strategy affects the price of the underlying security was examined in the studies of Frey [Reference Frey30] and Frey and Patie [Reference Frey, Patie, Sandman and Schönbucher31]. Liu and Yong [Reference Liu and Yong41] also discussed the problem of the replication of a European contingent claim with maturity T and payoff
$f(S)$
for the stock price nonlinear Black–Scholes equation
$$ \begin{align} &\frac{\partial V}{\partial t}+\frac{\sigma^2S^2}{2\Big(1-\lambda(S,t)S\frac{\partial ^2 V}{\partial S^2}\Big)}\frac{\partial ^2 V}{\partial S^2}+rS\frac{\partial V}{\partial S}+rS\frac{\partial V}{\partial S}-rV=0, \quad (S,t) \in (0,\infty) \times (0,T], \end{align} $$
$$ \begin{align} &V(S,T)=f(S), \quad S \in \Omega:=(0,+\infty), \end{align} $$
where V, S, T,
$\sigma>0 $
and
$r \geq 0$
stand for the price of the option, the stock price, the maturity date, the asset volatility and the interest rate, respectively. Here,
$f(S)$
and
$\lambda (S,t)$
describe a continuous piecewise linear payoff function and the price impact factor, respectively. The existence and uniqueness of classical solution to the problem is studied in the work of Liu and Yong [Reference Liu and Yong41] where the conditions ensured are:
-
(i)
$f(e^x)$
is Lipschitz continuous; -
(ii)
$e^{-\beta \sqrt {1+x^2}}f(e^x)$
is bounded for some
$\beta \geq 0$
.
These conditions are valid when the payoff function of the European contingent claim is a continuous piecewise linear function [Reference Company, Jodar and Pintos17, Reference Company, Jodar and Pintos18]. Note that
$\lambda (S,t)$
is the price impact factor influence of the trader involved in the hedging strategies; some regularity conditions given in the work of Liu and Yong [Reference Liu and Yong41] are valid for this parameter [Reference Company, Jodar and Pintos17, Reference Company, Jodar and Pintos18]. The price impact function of the trader
$\lambda (S,t)$
demonstrates the hedging strategies: as a trader buys, the stock price goes up and vice versa. This is defined as follows:
$$ \begin{align*} \lambda(S,t)= \begin{cases} \dfrac{\gamma}{S}(1-e^{-\beta(T-t)}) & \underline{S} \leq S \leq \overline{S},\\ 0 & \text{otherwise,}\\ \end{cases}\end{align*} $$
where the constant price impact coefficient
$\gamma>0$
measures the price impact per traded share, and
$\underline {S}$
and
$\overline {S}$
respectively represent the lower and upper limit of the stock price within which there is a price impact. More details about the topic can be found in the work of Liu and Yong [Reference Liu and Yong41].
Using the change
$\tau =T-t$
,
$u(S,\tau )=V(S,t)$
, equations (1.1)–(1.2) transform to
$$ \begin{align*} &\frac{\partial u}{\partial \tau}-\frac{\sigma^2 S^2}{2\Big(1-\lambda(S,T-\tau)S\frac{\partial ^2 u}{\partial S^2}\Big)^2}\frac{\partial ^2 u}{\partial S^2}-rS\frac{\partial u}{\partial S}+ru=0, \quad (S,\tau) \in (0,+\infty) \in (0,T] ,\\ & u(S,0)=f(S), \quad 0<S<+\infty. \end{align*} $$
Studies on analytical solutions of the nonlinear Black–Scholes equation were executed for some special cases by various researchers [Reference Bordag, Sarychev, Shiryaev, Guerra and Grossinho10–Reference Bordag, Frey and Ehrhardt13]. Later, some authors continued to be interested in analytical solutions of the nonlinear Black–Scholes equation [Reference Duris, Tan, Lai and Ševčovič23, Reference Edeki, Ugbebor and Owoloko24, Reference Esekon, Onyango and Ongoti27, Reference Goxiola, Chavez and Santiago32, Reference Yan53]. However, since there is no exact solution for the case of a call or put terminal payoff, efficient numerical techniques have been needed to solve the nonlinear Black–Scholes equation.
In the following years, many researchers focused on the studies of nonlinear option pricing models in illiquid markets, which is an important issue in financial risk management. In this context, while some authors [Reference Ankudinova and Ehrhardt2, Reference Arenas, Gonzales-Parra and Caraballo3, Reference Company, Jodar and Pintos17, Reference Company, Jodar and Pintos18, Reference Ehrhardt and Valkov25, Reference Guo and Wang34, Reference Heider35, Reference Koleva and Vulkov39] were interested in capturing the behaviour of the model by a numerical scheme, some others [Reference Cuomo, Sica and Toraldo20, Reference Kim, Bae and Koo38, Reference Schönbucher and Wilmott45] took into account the implications of certain terms and approaches to possible financial outcomes. In addition, many authors [Reference Ankudinova and Ehrhardt2, Reference During, Fournier and Jüngel22, Reference Duris, Tan, Lai and Ševčovič23] successfully applied linearization techniques with various numerical schemes due to coping with analysing the nonlinear terms. To tackle the nonlinearity, Bordag and Chmakova [Reference Bordag and Chmakova12] reduced equation (1.1) to an ordinary differential equation in some special cases using Lie group theory. Also, Frey and Patie [Reference Frey, Patie, Sandman and Schönbucher31] solved a smooth version instead of the model in equations (1.1)–(1.2) due to strong nonlinearity. In the studies [Reference Company, Jodar and Pintos17–Reference Company, Jodar, Ponsoda and Ballester19], the nonlinear option pricing model was analysed by a finite difference scheme and its stability properties were investigated. Yet, the applicability of implicit numerical schemes to various nonlinear Black–Scholes models including transaction cost, market liquidity and volatility uncertainty was discussed by Heider [Reference Heider35]. At the same time, Guo and Wang [Reference Guo and Wang34] focused on the numerical solution of the nonlinear Black–Scholes equation, which models illiquid markets, based on locally one-dimensional (LOD) methods. In addition to all this, Koleva and Vulkov [Reference Koleva and Vulkov39] presented a wide class of nonlinear option models in the illiquid markets. Arenas et al. [Reference Arenas, Gonzales-Parra and Caraballo3] put a nonstandard finite difference scheme on the agenda in solving the nonlinear Black–Scholes equation in illiquid markets. Ehrhardt and Valkov [Reference Ehrhardt and Valkov25] investigated numerical solutions of two nonlinear Black–Scholes equations in illiquid markets. Although the nonlinear models of interest were analysed in the above studies with some success, since an exact solution does not exist for these problems, there is still a need to develop an effective approach to get more accurate solutions.
At this point, in this study, an effective alternative approach has been developed based on linearizing the option pricing model, followed by discretizing it with a numerical scheme. In this context, the Fréchet derivative is applied to the vanilla call option for an illiquid market and then the Newton iteration technique has been considered to produce the results. At the same time, the linearization technique helps to avoid considering the convergence issues of the Newton iteration applied to the nonlinear algebraic system containing many unknowns at each time step if an implicit method is used in time discretization [Reference Erdogan, Sari and Kocak26]. Thus, the technique allows us to use a large time step size when using an implicit time discretization scheme. The Fréchet derivative helps to enhance the convergence order of the proposed iterative scheme. After the original equation is linearized, the differential quadrature method (DQM), where approximations of the spatial derivatives are based on a polynomial of high degree in space and the second-order accurate trapezoidal rule in time, has been combined to obtain highly accurate solutions of the equation. The DQM, which is a discretization method using a considerably small number of grid points to solve various problems accurately, was introduced in the early 1970s by Bellman et al. [Reference Bellman, Kashef and Casti8]. The hybrid method has been seen to provide very accurate solutions with relatively little computational effort and very little storage requirement. To the best the authors’ knowledge, there is no study on the proposed linearization technique combined with DQM applied to the nonlinear European option problems.
2. Differential quadrature method
The DQM was proposed for the first time by Bellman et al. [Reference Bellman, Kashef and Casti8] as an efficient discretization technique to solve nonlinear partial differential equations. Later, the technique was successfully developed to investigate the solutions of many problems in various fields of science [Reference Bert and Malik9, Reference Chen, Lv and Bian16, Reference O’Mahoney43, Reference Shu46, Reference Shu and Richards48, Reference Xionghua and Zhihong52, Reference Ycel54]. In the study [Reference Shu46], numerical techniques based on differential quadrature were developed for the solution of partial differential and integral equations, and incompressible viscous flows were simulated using these techniques. Shu and Richards [Reference Shu and Richards48] applied generalized differential quadrature to solve the two-dimensional incompressible Navier–Stokes equations in the vorticity-stream-function formulation. Bert and Malik [Reference Bert and Malik9] presented a review of the differential quadrature method, which should be of general interest to the computational mechanics community. A combined method with the differential quadrature and Taylor series was applied to the two-dimensional inverse heat conduction problem by O’Mahoney [Reference O’Mahoney43]. Chen et al. [Reference Chen, Lv and Bian16] presented a method of state-space-based differential quadrature for free vibration of generally laminated beams. The polynomial-based differential quadrature (PDQ) and the Fourier expansion-based differential quadrature (FDQ) methods were applied to obtain eigenvalues of the Sturm–Liouville problem by Yucel [Reference Ycel54]. In addition to these, Xionghua and Zhihong [Reference Xionghua and Zhihong52] applied the differential quadrature method to solve American option problem and used it to overcome the difficulty in determining the optimal exercise boundary of American options.
In this technique, the partial derivative of a function with respect to a variable can be expressed by the sum of weighted function values at all grid points in that direction. The weighting coefficients do not change to any special problem and depend on choosing the grid selection that affects an important role in the accuracy of the solution. The point selection may or may not be evenly spaced. Equally spaced grid points can be considered to be easy and convenient to work with. However, to obtain a more accurate solution, the frequently used Chebyshev–Gauss–Lobatto points are preferred [Reference Ycel54]. For a domain specified by
$a \leq x \leq b$
and discretized by a set of unequally spaced points, the coordinate of any point i can be determined by
$$ \begin{align*} x_i=a+\frac{1}{2}\bigg(1-\cos\bigg(\frac{i-1}{N-1}\pi\bigg)\bigg)(b-a). \end{align*} $$
The values of the function
$u(x,t)$
at any time on the above grid points are given by
$u(x_i,t)$
,
$i=1,2,\ldots ,N$
. Here, N is the number of grid points. The differential quadrature discretizations of the first- and second-order spatial derivatives are respectively given by
$$ \begin{align*} \frac{\partial u}{\partial x}\bigg|_{x=x_i}&=\sum_{j=1}^{N}a_{ij}u(x_j,t), \quad i=1,2,\ldots,N,\\ \frac{\partial ^2 u}{\partial x^2}\bigg|_{x=x_i}&=\sum_{j=1}^{N}b_{ij}u(x_j,t),\quad i=1,2,\ldots,N, \end{align*} $$
where
$a_{ij}$
and
$b_{ij}$
are the weighting coefficients of the first- and second-order derivatives, respectively [Reference Shu47]. The weighting coefficients are determined by using many kinds of test functions such as the polynomials, the sine and cosine functions [Reference Striz, Wang and Bert51], the spline function [Reference Kashef and Bellman37] and various orthogonal polynomials [Reference Chang, Tsai and Lin15, Reference Mansell, Merryfield and Shizgal42]. The Lagrange interpolation function [Reference Shu and Richards48] is widely used, because it has no restriction on the choice of the points to avoid the ill-conditioning Vandermonde matrix in the calculation of the weighting coefficients,
$$ \begin{align*} a_{ij}&=\frac{1}{x_j-x_i}\prod_{k=1,k\ne i,j}^{N}\frac{x_i-x_k}{x_j-x_i}, \quad i\ne j, \\ a_{ii}&=-\sum_{j=1,j\ne i}^{N}a_{ij}. \end{align*} $$
The weighting coefficients of the second-order derivative can be considered as
$B=A^2$
, where
$A=(a_{ij})$
is the weighting coefficient matrix of the first derivative [Reference Xionghua and Zhihong52].
3. Linearization and discretization process
First, the Black–Scholes equations (1.1)–(1.2) have been linearized by using the Newton iteration approach and the Fréchet derivative approximation. For this, as
$\phi :U \rightarrow W$
at
$u \in U$
, the initial value problems in equations (1.1)–(1.2) can respectively be written as follows:
$$ \begin{align*} &\phi(u)=\frac{\partial u}{\partial \tau}-\frac{\sigma^2 S^2}{2\Big(1-\lambda(S,T-\tau)S\frac{\partial ^2 u}{\partial S^2}\Big)}\frac{\partial ^2 u}{\partial S^2}+rS\frac{\partial u}{\partial S}-ru=0, \quad (S,t) \in (0,\infty) \times (0,T], \\ &u(S,T)=f(S), \quad S \in \Omega:= (0,+\infty). \end{align*} $$
The solution of the operator equation,
$\phi (u)=0$
, can be approximated by the Newton method
$$ \begin{align} u_{k+1}=u_{k}+\theta_{k}, \end{align} $$
where k is the iteration and
$\theta _k$
is the refinement variable for correcting function
$u_k$
. The refinement variable
$\theta _k$
is obtained by solving the following equation:
$$ \begin{align} \phi'(u_k)\theta_k=-\phi(u_k). \end{align} $$
This approach is an important issue in functional analysis framework. The convergence criteria of the Newton iteration in a Banach space were studied by Ambrosetti and Giovanni [Reference Ambrosetti and Giovanni1]. The convergence conditions of the Newton iteration in infinite-dimensional Banach spaces were investigated in the work of Argvros [Reference Argvros4]. In addition, the convergence analysis of the Newton method combining a Fréchet derivative is studied by Korkut et al. [Reference Korkut, Kaymak and Tanoglu40]. As for this work, to the best knowledge of the authors, for the first time, the proposed linearization approach has been discussed to cope with the nonlinearities in option pricing problems.
Solving this problem to obtain
$\theta _k$
, equation (3.1) is used and iterated. The Kantorovich theorem guarantees that the Newton method converges for sufficiently suitable initial conditions [Reference Atkinson and Han5, Reference Deuflhard21].
With the use of the Fréchet derivative [Reference Erdogan, Sari and Kocak26],
$$ \begin{align} \phi'(u_k)(\theta_k)=\frac{\partial }{\partial \epsilon}\phi(u_k+\epsilon \theta_k )|_{\epsilon=0}, \end{align} $$
and the DQM for spatial discretization systems in equations (3.2)–(3.3) are then respectively expressed as follows:
$$ \begin{align} \frac{\partial \theta_k}{\partial \tau}&-\frac{1}{2}\frac{\sigma^2S^2}{(1-\lambda(S,T-\tau) SBV)^3}(1+\lambda(S,T-t)SBV)B\theta_k-rSA\theta_k+r\theta_k \nonumber \\[3pt]&\!\! =-\frac{\partial{u_k}}{\partial \tau}+\frac{1}{2}\frac{\sigma^2 S^2}{(1-\lambda(S,T-\tau)SBu_k)^2}Bu_k+rSAu_k-ru_k, \\[3pt]\theta_k(S,0)\,&\!\!=0,\nonumber \end{align} $$
where
$\phi :U \rightarrow W$
is an operator at
$u \in U$
and
$V,W$
are Banach spaces, and U is an open subset of V.
Next, the quadratic accurate trapezoidal rule for time discretization was implemented. The stability properties of the method can be found in the reference [Reference Smith50]. The time coordinate is discretized with uniformly spread M grids. The step length in the t-direction is denoted by
$dt={T}/{M}$
, in which T is the maturity time of the contract. After the linearization, and applying the DQM and time discretization, equation (3.4) can be expressed as
$$ \begin{align*} &\frac{\theta^{n+1}-\theta^{n}}{\Delta \tau}-\frac{1}{2}\frac{\sigma^2S^2}{(1-\lambda(S,T-\tau)SB({u^{n+1}+u^{n}})/{2})^3}\bigg(1+\lambda(S,T-\tau) S B\frac{u^{n+1}+u^{n}}{2}\bigg)\\[5pt]&\qquad \times B\frac{\theta^{n+1}+\theta^{n}}{2} -rSA\frac{\theta^{n+1}+\theta^{n}}{2}+r\frac{\theta^{n+1}+\theta^{n}}{2} \\[5pt]&\quad=-\frac{u^{n+1}+u^{n}}{\Delta \tau}+\frac{1}{2}\frac{\sigma^2S^2}{(1-\lambda(S,T-\tau) S B(u^{n+1}+u^{n})/{2})^2}B\frac{u^{n+1}+u^{n}}{2} \\[5pt]&\qquad +rSA\frac{u^{n+1}+u^{n+1}}{2}-r\frac{u^{n+1}+u^{n+1}}{2}. \end{align*} $$
The pseudo-code of the presented algorithm is given in Algorithm 1.
4. Results and analysis
In this section, in illiquid markets, the numerical solution of the Liu and Yong model [Reference Liu and Yong41] with the combined method discussed in the previous sections is presented and compared with the existing ones in the literature. All computations have been performed in double precision using MATLAB
$^{\circledR }$
2019.
Consider the European vanilla call option for an illiquid market with
$E=50, r=0.06, \sigma =0.4, T=1, \underline {S}=20, \overline {S}=80, \beta =100, \gamma =1$
taken from the literature [Reference Company, Jodar and Pintos18]. The initial boundary condition
$u(S,0)=\max \{S-E,0\}$
is the suitable value for
$u_0$
. After the hybridization of the Newton iteration with the Fréchet derivative is applied, the linearized model is discretized by the DQM in space and a quadratic trapezoid rule in time. The obtained linear algebraic equations are solved numerically with iteration in each time value for a tolerance of
$1\times 10^{-10}$
. In the Newton method, the derivatives are represented by the Fréchet derivatives instead of using the usual Jacobian matrices. The finite difference method (FDM) used by Company et al. [Reference Company, Jodar and Pintos18] has been recalculated here for the European vanilla call option prices as a reference solution.
Table 1 The results and relative errors for
$N=16$
.
Table 2 The results and relative errors for
$N=32$
.
Table 3 The results and relative errors for
$N=64$
.
The solutions obtained by the proposed method and the FDM, and relative errors
$\varepsilon $
defined by
$$ \begin{align*} \varepsilon_i=\bigg| \frac{u^{\text{FDM}}_i-u^{\text{prop}}_i}{u^{\text{FDM}}_i}\bigg|, \end{align*} $$
are presented in the related tables and figures. Here,
$u^{\text {FDM}}$
and
$u^{\text {prop}}$
indicate the solutions at the ith grid points obtained by the FDM and the proposed method, respectively, and are given in Tables 1, 2 and 3 for
$\delta t=0.01,0.001,0.0001$
and
$N=16,32,64$
, respectively. In addition, the spline interpolation is applied to the FDM solutions for obtaining the values on the related grid points. As seen from the results, the proposed method can generate highly accurate solutions even for relatively large spatial and time elements. Figures 1–3 illustrate the results in Tables 1–3, respectively. The figures show that the behaviours of the model with the proposed method and FDM are in very good agreement. Moreover, Figure 4 displays the response of the model for different spatial elements.
Figure 1 Hedge cost for the European vanilla call option for
$N=16$
,
$\Delta t=0.01,0.001,0.0001$
at maturity time
$T=1 $
.
Figure 2 Hedge cost for the European vanilla call option for
$N=32$
,
$\Delta t=0.01,0.001,0.0001$
at maturity time
$T=1 $
.
Figure 3 Hedge cost for the European vanilla call option for
$N=64$
,
$\Delta t=0.01,0.001,0.0001$
at maturity time
$T=1 $
.
Figure 4 Hedge cost for the European vanilla call option for
$M=19$
at maturity time
$T=1 $
.
Figure 5 Hedge costs with the parameter
$\gamma $
for
$N=64$
and
$\Delta t=0.0001$
.
Table 4 displays the solutions generated by the proposed method and the FDM for different spatial and time elements. As seen from the values, the proposed method allows us to use relatively large spatial and time step sizes compared with the FDM. However, even when a few iterations are used, the technique reaches high accurate solutions. In addition, for
$\Delta t=7.0671\times 10^{-4}$
, while the FDM has spurious oscillations, the proposed method works very well.
Figure 5 illustrates the variation of the option price with different values of the parameter
$\gamma $
simulating the illiquidity influence in the price of the option. As seen, as the value of the parameter
$\gamma $
grows, the price grows to
$S=E$
and for
$S>E$
, for different
$\gamma $
values, the prices are close to each other.
Due to the effect on the nonlinearity of the problem as seen in equation (1.1), the numerical behaviours of the Delta
${\partial V}/{\partial S}$
and Gamma
${\partial ^2 V}/{\partial S^2}$
of the option for different
$\gamma $
values are exhibited in Figures 6 and 7. From Figure 6, it can be concluded that the hedge ratio is increasing in
$\delta $
for
$S<E$
, and decreasing in
$\delta $
for
$S>E$
. Additionally, Figure 7 illustrates the effect of
$\gamma $
on the variation of the Gamma and it is seen that the Gamma flattens out as illiquidity increases moving its peak more to smaller values of S.
5. Conclusions and recommendation
Since an exact solution for the illiquid European option does not exist, the need for capturing the behaviour of the model leads researchers to obtain effective methods. In this paper, an effective hybrid numerical approach has been introduced to determine the behaviours represented by the nonlinear Black–Scholes equation resulting from the nonlinear property of volatility known as the illiquid European option model. For this purpose, first, the derivative in the Newton iteration formula has been evaluated by the Fréchet derivative in capturing the numerical behaviour of the model. In the investigation of the behaviour of the problem with nonlinear volatility, called the illiquid European option model, the proposed approach based on both the temporal variability and the spatial variability has been found to be very easy to apply and computationally inexpensive. When compared with the literature, the hybrid method has been found to be very accurate and efficient. We believe that this study can be adapted to deal with the real-life problems representing economic behaviour and financial problems with volatility.
Table 4 The results for different spatial and time elements.
Figure 6 Variation of the Delta.
Figure 7 Variation of the Gamma.
