AN INVESTIGATION ON NONLINEAR OPTION PRICING BEHAVIOURS THROUGH A NEW FRÉCHET DERIVATIVE-BASED QUADRATURE APPROACH

Abstract

Complicated option pricing models attract much attention in financial industries, as they produce relatively better accurate values by taking into account more realistic assumptions such as market liquidity, uncertain volatility and so forth. We propose a new hybrid method to accurately explore the behaviour of the nonlinear pricing model in illiquid markets, which is important in financial risk management. Our method is based on the Newton iteration technique and the Fréchet derivative to linearize the model. The linearized equation is then discretized by a differential quadrature method in space and a quadratic trapezoid rule in time. It is observed through computations that the accurate solutions for the model emerge using very few grid points and time elements, compared with the finite difference method in the literature. Furthermore, this method also helps to avoid consideration of the convergence issues of the Newton approach applied to the nonlinear algebraic system containing many unknowns at each time step if an implicit method is used in time discretization. It is important to note that the Fréchet derivative supports to enhance the convergence order of the proposed iterative scheme.

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 [6] developed a parametrized model for liquidity effects arising from the trading in assets. Bank and Baum [7] introduced a general continuous-time model for an illiquid financial market. In the study [12], the authors used families of explicit solutions to test numerical schemes solving a nonlinear Black–Scholes equation. Cetin et al. [14] 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 [28]. Another remarkable study was presented by Frey [29]. 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 [33] 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 [36] analysed market manipulation trading strategies by large traders in a securities market. Platen and Schweizer [44] 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 [45] analysed the influence of dynamic trading strategies on the prices in financial markets. Sircar and Papanicolaou [49] 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 [30] and Frey and Patie [31]. Liu and Yong [41] 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

(1.1) $$ \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} $$

(1.2) $$ \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 [41] where the conditions ensured are:

1. (i) $f(e^x)$ is Lipschitz continuous;

2. (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 [17, 18]. 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 [41] are valid for this parameter [17, 18]. 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 [41].

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 [10–13]. Later, some authors continued to be interested in analytical solutions of the nonlinear Black–Scholes equation [23, 24, 27, 32, 53]. 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 [2, 3, 17, 18, 25, 34, 35, 39] were interested in capturing the behaviour of the model by a numerical scheme, some others [20, 38, 45] took into account the implications of certain terms and approaches to possible financial outcomes. In addition, many authors [2, 22, 23] successfully applied linearization techniques with various numerical schemes due to coping with analysing the nonlinear terms. To tackle the nonlinearity, Bordag and Chmakova [12] reduced equation (1.1) to an ordinary differential equation in some special cases using Lie group theory. Also, Frey and Patie [31] solved a smooth version instead of the model in equations (1.1)–(1.2) due to strong nonlinearity. In the studies [17–19], 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 [35]. At the same time, Guo and Wang [34] 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 [39] presented a wide class of nonlinear option models in the illiquid markets. Arenas et al. [3] put a nonstandard finite difference scheme on the agenda in solving the nonlinear Black–Scholes equation in illiquid markets. Ehrhardt and Valkov [25] 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 [26]. 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. [8]. 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. [8] 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 [9, 16, 43, 46, 48, 52, 54]. In the study [46], 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 [48] applied generalized differential quadrature to solve the two-dimensional incompressible Navier–Stokes equations in the vorticity-stream-function formulation. Bert and Malik [9] 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 [43]. Chen et al. [16] 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 [54]. In addition to these, Xionghua and Zhihong [52] 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 [54]. 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 [47]. The weighting coefficients are determined by using many kinds of test functions such as the polynomials, the sine and cosine functions [51], the spline function [37] and various orthogonal polynomials [15, 42]. The Lagrange interpolation function [48] 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 [52].

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

(3.1) $$ \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:

(3.2) $$ \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 [1]. The convergence conditions of the Newton iteration in infinite-dimensional Banach spaces were investigated in the work of Argvros [4]. In addition, the convergence analysis of the Newton method combining a Fréchet derivative is studied by Korkut et al. [40]. 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 [5, 21].

With the use of the Fréchet derivative [26],

(3.3) $$ \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:

(3.4) $$ \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 [50]. 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 [41] 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 [18]. 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. [18] 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.