Applying the Barycentric Jacobi Spectral Method to Price Options with Transaction Costs in a Fractional Black-Scholes Framework

The aim of this paper is to show how options with transaction costs under fractional, mixed Brownian-fractional, and subdiffusive fractional Black-Scholes models can be efficiently computed by using the barycentric Jacobi spectral method. The reliability of the barycentric Jacobi spectral method for space (asset) direction discretization is demonstrated by solving partial differential equations (PDEs) arising from pricing European options with transaction costs under these models. The discretization of these PDEs in time relies on the implicit Runge-Kutta Radau IIA method. We conducted various numerical experiments and compared our numerical results with existing analytical solutions. It was found that the proposed method is efficient, highly accurate and reliable, and is an alternative to some existing numerical methods for pricing financial options.


Introduction
Modeling of financial derivatives has been of great interest in the past three decades.Numerous mathematical models have been developed from the classical Black-Scholes [1] framework to help investors in their decision making process.These models are based on an arbitrage argument, i.e., by continuously adjusting a portfolio consisting of a stock and a risk-free bond, an investor can exactly replicate the return to any option on the stock.
In the presence of transaction costs in capital markets, the presence of an arbitrage argument [2] is no longer valid [3], since perfect hedging is impossible.Due to the infinite variation of the geometric Brownian motion, the continuous replication policy incurs an infinite amount of transaction costs over any trading interval no matter how small it might be.Therefore, in recent years, one observes many generalizations of the Black-Scholes model to deal with the problem of option pricing and hedging with transaction costs.This leads to the Black-Scholes model but with an adjusted volatility.Leland [3] was the first to examine option replication in a discrete time setting, and proposed a modified replicating strategy, which depended upon the level of transaction costs, the revision interval, the option to be replicated and the environment.Subsequently, several authors proposed new models, but all in discrete time [4,5].However, these models based on the diffusion process known as geometric Brownian motion (GBM) have severe shortcomings, for example, long-range correlations, heavy-tailed and skewed marginal distributions, lack of scale invariance and periods of constant values, to enumerate these only.Recently, Wang [6] obtained a European call option pricing formula using a mean self-financing delta hedging argument in a discrete time setting and then showed how scaling and long-range dependence impacted dramatically on the pricing of options using the fractional and multifractional Black-Scholes model with transaction costs.
In continuous time, a lot of efforts have been done in order to alleviate the problem of an infinite amount of transaction costs over any trading interval when the asset process follows a geometric Brownian motion.Magdziarz [7] et al. introduced a subdiffusive geometric Brownian process which captures the subdiffusive characteristics of financial markets.They showed that their model is arbitrage-free.The same idea was used later by Wang et al. [8], Hui and Yun-Xiu [9] to obtain a Black-Scholes equation with transaction costs in subdiffusive fractional Brownian motion regime.However, closed-form solutions of these PDEs in finance are generally rare.Therefore, one has to consider numerical methods to obtain solutions.
In this paper, we are concerned with the application of the barycentric Jacobi interpolation to value options with transaction costs under fractional [6], mixed Brownian-fractional [10] and subdiffusive fractional [8] Black-Scholes models.Barycentric spectral methods were introduced by Baltensperger et al. [11] to solve boundary value problems.Recently, these methods have emerged in the field of finance as a promising tool to solve option pricing problems.Pindza and Patidar [12] proposed an accurate method, namely the barycentric Chebyshev spectral method, to price options in illiquid markets.Ngounda et al. [13] used the barycentric Chebyshev domain decomposition method to provide fast and accurate results for pricing European options with jumps, which was later extended and applied to Heston's volatility model (see [14]).Most of the work on barycentric spectral methods has been based on the use of uniform or Chebyshev grids.Recently, Wang et al. [15] computed explicit barycentric weights for Jacobi polynomial interpolation in the roots or extrema of classical orthogonal polynomials in terms of the nodes and weights of the corresponding Gaussian quadrature rule.Hence, we investigate the utility of this new barycentric interpolation in the field of finance.The semi-discretization of the PDE in time is realized by using a 7-stage 13th-order fully implicit Runge-Kutta Radau IIA method with adaptive time stepping [16].
This paper is structured as follows.Section 2 reviews the option pricing formulation with transaction costs under fractional, mixed Brownian-fractional and subdiffusive fractional processes.In Section 3, we introduce the Jacobi barycentric spectral method, semi-discretize the PDE in the asset space and propose a conformal map in order to improve the accuracy of our method.In Section 4, we perform numerous experiments in order to advocate the utility of the barycentric spectral method.Finally, Section 5 gives a summary and scope for further research.

Pricing with Transaction Costs under Fractional, Mixed Brownian-Fractional and Subdiffusive Fractional Processes
, , , , be a complete probability space carrying a fractional Brownian motion , where  is the set of all possible outcomes of the experiment known as the sample space, is the set of all events, is a real world probality, t is a natural filtration, t a risky underlying asset price process.Assume that the price of the underlying stock at time satisfies a fractional Black-Scholes model where  , and 0 are constants.Assume that the portfolio is revised every small fixed time step 0,1  ; transaction costs are proportional to the value of the transaction in the underlying.Let denote the round trip transaction cost per unit dollar of transaction.Then under all the assumptions given by Wang [6], the Black-Scholes equation with transaction costs assuming factional Brownian motion can be written as where the modified volatility is given by sign .
Here is the fractional Leland number [6].

 
Le H If we assume that the price of the underlying stock at time satisfies a mixed Brownian-fractional Brownian model then under all the assumptions given by Wang [10], the Black-Scholes equation with transaction costs in mixed Brownian-fractional Brownian motion can be written as where the modified volatility is given by and .
Wang et al. [8] obtained the modified volatility corresponding to the continuous Black-Scholes equation with transaction cost in subdiffusive regime as where the modified volatility is given by Here     represents the gamma function evaluated at  and

 
sign  is the second derivative of the option value with respect to the asset.
The difference between American and European call and put options is made by the initial and boundary conditions.In this work we focus exclusively on the European call option.Such options have the following initial and boundary conditions , where . We want to solve the three PDEs (2), ( 5) and ( 8) subjected to modified volatilities (3), ( 6) and ( 9), respectively.Before moving to the applications of the barycentric Jacobi spectral method to solve these problems, it is worthwhile to consider some preliminaries of this method.

Barycentric Jacobi Spectral Collocation Method
In this section, we turn our attention to the problem of barycentric Jacobi interpolation and present a fast algorithm for the efficient computation of the interpolation formula using Jacobi-Gauss-Lobatto (JGL) points.This interpolation is realized by a class of the Lagrange form of the interpolating polynomial, as follows.Let j x , be a set of distinct nodes.Then the polynomial of degree that interpolates the function at these points is given by [17] where the Lagrange polynomial corresponding to the node j  j x has the property with .Generally, the Lagrange form of the interpolating polynomial ( 11) is not advocated for numerical computations.In particular, it requires , 0,1, , additions and multiplications for each evaluation of N and every time a node j   p x x is modified or added, all Lagrange fundamental polynomials have to be recalculated.However, with slight modifications, the Lagrange formula is indeed of great practical use.This has been noted by several authors, including Henrici [18] and Werner [19].Berrut and Trefethen [20] modified the Lagrange polynomial through barycentric interpolation and proposed an improved Lagrange formula.Following [20], we define   , the numerator of in (11) divided by In addition, if we define the barycentric weight by i.e.,   Consequently, the Lagrange formula (11) becomes In particular, if   1 u x  , we obtain Dividing ( 16) by (17), we get the barycentric formula for N p as This is the most used form of Lagrange interpolation in practice and admits operations.In order to obtain good approximations via interpolation, the choice of interpolation nodes and barycentric weights is particularly important.For certain particular sets of points, such as equidistant points as well as Chebyshev points, the barycentric weights j  can be computed analytically [20].Recently, Wang et al. [15] computed explicit barycentric weights for polynomial interpolation in the roots or extrema of classical orthogonal polynomials in terms of the nodes and weights of the corresponding Gaussian quadrature rule.
In this paper, we are interested in the barycentric Jacobi interpolation.The Jacobi-Gauss-Lobatto quadrature rule is defined by where the Jacobi-Gauss-Lobatto points   0 is the Jacobi polynomial of degree .The following result gives an analytical formula of barycentric weights for the Jacobi-Gauss-Lobatto points.
and let be the corresponding j w weights of the interpolatory quadrature rule with weight function     . Then the simplified bary-centric weights are for Jacobi-Gauss-Lobatto points, the simplified barycentric weights j  are given by The proof of Theorem 3.1 can be found in [15].The computation of entries of the first and second order differentiation matrices and where and , is given as in [11] by where ., 0,1, , i j N  

Conformal Mappings for High Resolution of Non-Smooth Initial Conditions
It is well-known that the solution of Equation ( 8) is very sensitive to localization errors when is in the vicinity of , because the second derivative of the payoff does not exist at this point.Therefore, to increase accuracy it would be reasonable to use an adaptive mesh with high concentration of the mesh points around  , while a rarefied mesh could be used far away from this area.If we assume that , where both and are finite, then without further loss of generality we may assume that the interval of integration is In this paper, we use the conformal map g given in [21] by , where   and   determine the magnitude of the region of rapid change and the location, respectively.Here, represents the Jacobi-Gauss-Lobatto points. y In Figure 1, we plot the new grid obtained from the original Jacobi grid by using transformations (25) and ( 24 , where is the strike price.In addition, we show the transformation function around for , (purple) and (red).We observe that when  grid distribution one.When   increases we accommodate more grid points around the strike price.
A significant advantage of the barycentric Jacobi spectral method is that it eliminates tedious computations of transformed derivatives using the chain rule as is usually the case in other spectral collocation methods.

Application to the Black-Scholes PDE
We discretize the Black-Scholes PDEs (2), ( 5) and ( 8) in the asset (space) direction by means of barycentric Jacobi spectral method.Let

 
x g y  be the transformed Jacobi-Gauss-Lobatto points, the first step is to transform that better suits the option at hand using , the PDE (8) together with its initial and boundary conditions yields , 0 , 11) into (28) yields the following system of nonlinear ODEs .
In order to write (29) in matrix form, we introduce the following matrix and vector notations

OPEN ACCESS JMF
where 1 identity matrix.Consequently, (29) can be expressed as an initial value problem of the form where We use a 7-stage 13th-order fully implicit Runge-Kutta Radau IIA method with step size control to integrate the system of ODEs (31).The method is B-stable and stiffly accurate.Details of the method can be found in [16].

Numerical Results and Discussions
We apply the barycentric Jacobi spectral method to value Black-Scholes equations with transaction costs under fractional (FBS), mixed Brownian-fractional (MBS) and subdiffusive fractional (SBS) processes.To show the efficiency of the present method we report the root mean square norm error   2 L of the solution computed with grid points, given by and the maximal norm error   where   j u S and are the benchmark and computed values of the solution at .

  exact j u S u j
The analytic solutions of Black-Scholes equations with transaction costs under FBS, MBS and SBS regime are possible if we assume that the Greek S  is always positive for the above mentioned models as in the case of standard Black-Scholes model, then for European call options under SBS [6] regime is known, and expressed as where T t    and is the cumulative probability distribution function for a standardized normal variable If we assume the same sign behaviour for  in the MBS model, then the analytic solution is given by where with the modified volatility given by Here In the case of FBS model, the closed-form solution is given as with In order to test the accuracy of the method, we present a comparison against the above exact solutions.In Figure 2 (top left) we plot together numerical and exact solutions for comparison purposes.We select numerical values of the parameters to be , 100   .We choose the tolerance of the 7-stage 13th-order fully implicit Runge- Kutta Radau IIA method [16] to be , so that the error is dominated by the spatial error.Clearly, it is observed that all numerical solutions are in good agreement with exact ones.
To see how good our numerical approach approximates exact solutions, we plot the absolute error, i.e., absolute distance between the exact solution and the numerical approximation for all three models.This shows very good accuracy for our method.

Effect of Changes in N and S max
To determine the convergence of the discretization scheme, we solve the problem by keeping some parameters fixed and varying others.We fist investigate the effect of the truncation domain on the errors by varying between 120 and 500 and keeping other parameters fixed as 3 show the -norm and 2 -norm errors, respectively.We notice that the errors remain bounded and do not vary significantly in term of the truncation domain.This is an important result, since we can accurately solve the option pricing problem on a small truncated domain, which will result in better efficiency.

L  L
In the next experiment, we investigate the convergence of the barycentric Jacobi spectral method.We vary the number of collocation points between 20 and 300, with parameters   ; and we plot the results in Figure 4 (top left and right).
For all three models, our method converges rapidly to the exact solutions.This is usually known as spectral or exponential convergence.The efficiency of our approach is measured in Figure 4 (bottom left and right) by the CPU time.The results on the efficiency of our method are very satisfactory for all three methods.An accuracy of can be obtained in less than a second.

Convergence of the Method
We explore the effect of changes in time, as well as grid-stretching on the accuracy of the model, keeping

Table 1.
We notice that when the grid stretching parameter   increases, the accuracy of our method is improved.By accommodating more grid points around the strike price , we can overcome the poor convergence of naive application of numerical methods when pricing options.We also observe that our method is highly accurate (even for long maturity options).

K
We also explore the effect of changes in time, as well as Jacobi parameters  and  on the accuracy of Once again the barycentric Jacobi spectral method achieves very good accuracy for all three models (even for long term maturity).The Jacobi parameters chosen here do not impact the accuracy of our methodology significantly.It would be useful to investigate how to choose these parameters in an optimal way, however, this is beyond the scope of this paper.

Conclusion
The work of Leland [3] on the discrete option pricing model and replication with transaction cost has paved the way to develop the continuous version.We exploited the continuous version by Wang et al. [8]   spectral-based solutions to the model.In practice, option pricing problems are solved numerically since analytical solutions rarely exist.We acknowledge that many studies have been conducted in the domain of finance and both numerical and analytical solutions have been investigated.However, the barycentric Jacobi spectral method has just been proven to have a better accuracy and has not been studied in the field of PDEs in finance.Furthermore, this method obtains solutions with greater accuracy than the usual well-known and studied numerical schemes.The barycentric Jacobi spectral method has full differentiation matrices, which require more computational effort than sparse differentiation matrix methods.In order to improve the efficiency of the rational Jacobi spectral method, we are currently investigating the domain decomposition algorithm which yields block diagonal matrices and we expect to have a greater accuracy, and at least less computational time and memory.
) together.The new grid in y


 decreases the grid approaches a Jacobi-Gauss-Lobatto

Figure 1 .
Figure 1.Asset direction grids S obtained from the mapped Jacobi grids x.

8 10 Figure 2 .Figure 3 .
Figure 2. Solutions of the European call options with transaction costs under the three models (top left), absolute error (top right), the Delta (bottom left) and Gamma (bottom right) of the three models with , 100 K  0.2   , , , , , , 0.03 r  1 T  200 N 

Figure 4 .
Figure 4. Convergence of L  and 2 L -norms and efficiency of the barycentric Jacobi spectral method with parame- ters , 100 K  0.2  

Table 1 . Norm infinity and norm relative of errors for the European call options with transaction costs under subdif- fusive regime (SBS), fractional Brownian (FBS) and multifractional Brownian (MBS) motions.
to construct