{'fraction_non_alphanumeric': 0.03813182369591589, 'fraction_numerical': 0.013141932875050546, 'mean_word_length': 4.63091985428051, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Comparison of the analytical approximation formula and Newton's method for solving a class of nonlinear Black-Scholes parabolic equations

Market illiquidity, feedback effects, presence of transaction costs, risk from unprotected portfolio and other nonlinear effects in PDE based option pricing models can be described by solutions to the generalized Black-Scholes parabolic equation with a diffusion term nonlinearly depending on the option price itself. Different linearization techniques such as Newton's method and analytic asymptotic approximation formula are adopted and compared for a wide class of nonlinear Black-Scholes equations including, in particular, the market illiquidity model and the risk-adjusted pricing model. Accuracy and time complexity of both numerical methods are compared. Furthermore, market quotes data was used to calibrate model parameters.


Introduction
According to the classical theory due to Black, Scholes and Merton an option in a stylized and idealized financial market can be priced by a solution V = V (S, t) to the linear Black-Scholes parabolic equation: where r > 0 is the interest rate of a zero-coupon bond, q ≥ 0 is the dividend yield rate andσ > 0 is a constant historical volatility of the underlying asset price process {S t , t ≥ 0}, which is assumed to follow a stochastic differential equation of the geometric Brownian motion with a drift r − q (cf. [22,18,25]). The linear Black-Scholes equation with a constant volatilityσ has been derived under several restrictive assumptions, for example, zero transaction costs, perfectly replicated portfolio, frictionless, market completeness, etc. In this paper, the main goal is to investigate and compare two numerical approximation methods for solving a class of nonlinear generalizations of the linear Black-Scholes equation (1) in which the volatility is assumed to be a function of the underlying asset price S and Gamma of the option (the Greek Gamma is the second derivative ∂ 2 S V ), i.e. σ = σ(∂ 2 S V, S) .
The motivation for solving the nonlinear Black-Scholes equation (1) with the volatility function σ of the form (3) arises from more realistic option pricing models in which one can take into account nontrivial transaction costs, market feedbacks, risk from unprotected portfolio and other effects. In the last decades, some of the restrictive assumption of the classical Black-Scholes theory [4] have been relaxed in order to model, for instance, presence of constant transaction costs (see e.g. Leland [23], Hoggard et al. [16]), non-constant transaction costs (see e.g. Amster et al. [1],Ševčovič andŽitňanská [26]), uncertain volatility model (cf. Avellaneda and Paras [2]), feedback and illiquid market effects due to large traders choosing given stock-trading strategies (cf. Frey [12], Frey and Patie [13], Frey and Stremme [14], Schönbucher and Wilmott [24]), imperfect replication and investor's preferences (cf. Barles and Soner [3]), risk from an unprotected portfolio (cf. Kratka [21], Jandačka andŠevčovič [19]). Efficient techniques and fast computational methods for pricing derivative securities is a practical task in financial quotes markets. Therefore, realistic PDE based option models including, in particular, nonlinear generalizations of the Black-Scholes equation have to be solved in a fast and efficient way. However, in most important cases there is no explicit formula except for some special cases with non-standard pay-off diagrams (cf. Bordag [6]). This is the reason why numerical methods for solving nonlinear Black-Scholes equation have to be developed and analyzed. In this paper, attention is focused on a class of nonlinear Black-Scholes equations. In particular, the nonlinear volatility model developed by Frey et al. [12,11,13,14] and the risk-adjusted pricing methodology model proposed and investigated by Kratka [21] and Jandačka andŠevčovič [19,25] are the main concern of this work. In a series of papers [12,11,13,14] Frey et al. considered a model in which the price of an underlying asset is affected by specific hedging strategies due to a large trader. Supposing that a large trader uses a given stock-holding strategy α t and the underlying stock price process satisfies the following SDE: where µ is a drift parameter, σ > 0 is the volatility of the process and 0 ≤ ρ <ρ is the so-called market liquidity parameter. It is worth noting that the quantity 1/(ρS t ) measures the size of the change in the stock-holding position of the large trader. Notice that if α t ≡ 0 or ρ = 0, the stock price S t follows the geometric Brownian motion. In [12] Frey (see also [13,14]) showed that the option price is then a solution to a nonlinear volatility Black-Scholes equation of the form: for 0 ≤ S < ∞ and 0 ≤ t < T where T is the maturity time. The nonlinear volatility function σ is given by whereσ is a constant historical volatility. A solution V = V (S, t) is subject to the terminal pay-off condition describing call or put option with expiration price E > 0, i.e.
Another nonlinear model was proposed by Kratka [21]. It was further generalized and analyzed by Jandačka andŠevčovič in [19,25]. The model is constructed following the classical Leland approach for modeling transaction costs (cf. [23]) in which the time between consecutive portfolio rearrangements is subject to optimization with respect to the risk arising from an unprotected portfolio. In this risk-adjusted pricing methodology (RAPM) model the nonlinear volatility function has the form: Construction of explicit solutions to equation (5) with the nonlinear volatility function as the one defined in (8) were recently provided by Bordag and Frey [5] (see also [6]). Several invariant solutions were constructed by means of the invariant Lie group theory. These invariant solutions depend on various parameters restricting the class of solutions. In particular, not every pay-off diagram can be considered. In general, there is no exact pricing formula for the case of a call or put terminal payoffs. Hence efficient numerical techniques for solving such nonlinear Black-Scholes equations are required.
A numerical method proposed and investigated by Jandačka andŠevčovič [19] is based on the transformation H = S∂ 2 S V, x = ln(S/E), τ = T − t, which transforms equation (5) with σ = σ(S∂ 2 S V ) into a porous media type of quasilinear parabolic equation: where β(H) = 1 2 σ 2 (H)H is an increasing function. For instance, in the case of the volatility function given by (6) one obtains β(H) =σ [25] for details). In the recent paper [26],Ševčovič andŽitňanská investigated the nonlinear equation (9) in the context of modeling variable transaction costs. The existence of classical Hölder smooth solutions was proved and useful bounds for the solution were derived.
The transformation technique developed in [19] allows for construction of a semiimplicit finite volume based numerical scheme for solving (9). There are other approaches dealing mainly with the nonlinear equation (5) for the option price rather than for its transformation H = S∂ 2 S V . Another method using quasilinearization technique for solving the fully nonlinear parabolic equation (5) was proposed and analyzed by Koleva and Vulkov [20]. A consistent monotone explicit finite difference numerical scheme was analyzed by Company et al. in the context of the Frey and Patie model (5) with a nonlinear volatility function given by (6). In [10] Ehrhardt and Valkov derived an unconditionally stable explicit numerical scheme for solving the same problem and provided necessary numerical analysis of the scheme.
In this paper, two numerical approximation methods based on the asymptotic perturbation analysis and the Newton linearization technique are developed. These methods are used to solve a wide class of nonlinear Black-Scholes equations. The first method is the asymptotic perturbation method which is based on asymptotic expansion of the solution into power series in a small model parameter. The first order expansion then corresponds to an explicit analytic approximation formula requiring only one-dimensional numerical integration which can be computed in a fast and efficient way. The second method is based on Newton's iterative method for solving the corresponding nonlinear problem in each temporal discretization level. It is applicable to a rather general nonlinear case not restricted by any specific types of equations and boundary and terminal conditions. In [15] Heider used Newton's iterative method for solving equation (5) with four types of nonlinear volatilities and different finite different schemes. Note that different variants of Newton's linearization and their implementation are also discussed and compared in this paper.
The paper is organized as follows. In section 2 an explicit analytic approximation formula for solving a general class of nonlinear volatility models is derived. In section 3 an algorithm utilizing Newton's method for solving equation (5) is described and analyzed. Several comparisons of both methods are discussed in section 4. Examples of solution to the Frey and Patie model and RAPM model are presented. Finally, section 5 contains an example of model calibration to real market quotes data.

Analytic approximation formula based asymptotic perturbation analysis
In this section, an analytic approximation formula for pricing European call or put options with a nonlinear volatility is derived. Typically this paper considers a wide class of nonlinear volatility functions taking the following form: The powers γ, δ, the parameter ε as well as the function A(T −t) depend on the chosen nonlinear volatility model. For example, in the case of the Frey and Patie model with the nonlinear volatility function given by (6), σ(∂ 2 S V, S) =σ/(1 − ρS∂ 2 S V ) ≈ σ(1 + ρS∂ 2 S V ) and the parameters are ε = ρ, γ = 1, δ = 2, A(T − t) =σ 2 , and the small model parameter ε can be identified with 0 < ρ 1.
For the RAPM model with the nonlinear volatility function given by (8) the parameters can be identified as follows: and the small parameter ε is identified with 0 < µ 1. Equation (5) can now be rewritten as whereV is the prescribed pay-off diagram. The problem is to seek the option price in the form of an asymptotic expansion in terms of a small parameter (cf. [17]). More precisely, where the leading term V 0 ≡ V BS is simply a solution to the linear Black-Scholes model. The aim here is to derive an asymptotic approximation formula obtained from the first two terms in the asymptotic expansion, i.e.
In order to obtain an explicit formula for the second term V 1 in the expansion, equation (10) is first approximated as follows: where L 0 is a linear and L 1 is a nonlinear differential operator in V , Hence the first order approximation of the equation L(V, ε) = 0 reads as follows: satisfying the initial condition: Equation (14) with the initial condition (15) can be separated into a system of equations in powers of ε, i.e.
O(ε 0 ) : The solution V 0 can be obtained by solving the linear Black-Scholes equation. The second equation for V 1 is a non-homogeneous PDE with zero initial condition.
Therefore, equation (16) can be solved once the value of V 0 (S, t) is evaluated to In order to solve equation (16) one adopts the usual transformation (see e.g. where Equation (16) is thus transformed to as follows: The termd 1 corresponds to d 1 after transformation (17). It is given bỹ Finally, equation (19) can be simplified to as the equation below (20) satisfying the growth condition |u(x, τ )| ≤ M e b|x| 2 for all x ∈ R, τ ∈ [0, T ] where M, b are some constants. Then u(x, τ ) is given by the formula: is a constant depending on the model parameters and the functions Q(τ, ξ) and Λ(τ, ξ) are defined as follows: The proof of this theorem is a straightforward application of the variation of constants formula and can be found in the Appendix. As a consequence of the previous theorem an explicit expression for the first order approximation of the option price can be obtained.
where K is a constant given by and The analytic approximation of the option price V (S, t) can then be evaluated by using equation (12).
3 Implicit finite difference scheme using Newton's method A standard way of solving equation (5) numerically is to use implicit temporal discretization in combination with a finite difference method for approximating the derivatives. Note the volatility term appearing in (6) and (8) is nonlinear and at each time level an iterative technique is to be applied. The method of frozen coefficient technique is commonly applied to handle the nonlinearity though sometimes it converges slowly without proper initial guess. To obtain a better convergence rate, Newton's method has to be employed in combination with a temporal implicit discretization scheme. Newton's method is a linearization technique with many variants and each takes different implementation. In this section two approaches are discussed. The first approach (denoted by NM1) addresses the root-finding problem of the nonlinear system derived from an implicit scheme in which calculation of the Jacobian matrix is used to update the approximate solution. The second approach (NM2) linearizes the original equation in which a correction term is to be solved and used to update the approximate solution.

Newton's Method (NM1)
Using standard finite difference notations and the transformation τ = T − t an implicit finite difference scheme which replaces equation (5) reads as follows: The volatility function σ as given by (6) may be discretised as Here S i = (i − 1)∆S, i = 1, · · · , M , and n = 1, · · · , N − 1, where M and N are the numbers of grid points for spatial and temporal discretization respectively. The above equation can be simplified as follows: where H(V n+1 ) is an M × M tridiagonal matrix whose elements nonlinearly depend on V n+1 . Introducing the mapping turns the original problem to the construction of a solution V n+1 of the equation G(V n+1 ) = 0 at each time level. Newton's method is applied to solve the root-finding problem which requires the Jacobian matrix of the function G to be computed. An initial guess chosen as the solution V from the previous time level usually reduces the number of Newton's iterations.
turns the original problem to the construction of a solution V n+1 of the equation G(V n+1 ) = 0 at each time level. Newton's method is applied to solve the rootfinding problem which requires the Jacobian matrix of the function G to be computed. An initial guess chosen as the solution V from the previous time level usually reduces the number of Newton's iterations.

Algorithm 1: Newton's Method (NM1)
Input: initial guess V n+1 , tol, initial condition In order to compute the Jacobian matrix efficiently, a decomposition of the matrix H(V n+1 ) may be adopted as below Note that H 1 and H 2 are constant tridiagonal matrices. By using this decomposition, the Jacobian matrix of G becomes k ) 2 is treated as a row vector. This decomposition simplifies the computing of the Jacobian matrix in terms of the nonlinear volatility. Each gradient ∇(σ n+1 k ) 2 can be obtained by either deriving the analytic expression for σ or by using a finite difference approximation of the spatial derivatives.
In order to compute the Jacobian matrix efficiently, a decomposition of the matrix H(V n+1 ) may be adopted as below Note that H 1 and H 2 are constant tridiagonal matrices. By using this decomposition, the Jacobian matrix of G becomes k ) 2 is treated as a row vector. This decomposition simplifies the computing of the Jacobian matrix in terms of the nonlinear volatility. Each gradient ∇(σ n+1 k ) 2 can be obtained by either deriving the analytic expression for σ or by using a finite difference approximation of the spatial derivatives.

Waveform-Newton's Method (NM2)
The second approach of applying Newton's linearization is to consider a smooth function F representing the nonlinear Black-Scholes equation, i.e.
Here V τ , V S , V SS abbreviate the partial derivatives ∂ τ , ∂ S V, ∂ 2 S V respectively. The linearization of the function F at (V * τ , V * S , V * SS , V * ) in direction (v τ , v S , v SS , v) reads as follows: where D 2 represents all higher order terms and the partial derivatives are evaluated at (V * τ , V * S , V * SS , V * ).
Equation (26) transforms equation (5) into a linear partial differential equation of the correction term v with zero boundary and initial conditions. This equation can be solved easily because all coefficients of equation (26) are determined. Similar to the first approach (NM1), these coefficients can be evaluated either by the analytic expression for σ or by a finite difference approximation. Eventually, the problem becomes Again, an initial guess can be set to the solution from previous time level in the algorithm. The main difference between algorithm NM1 and NM2 is the linearization error O(D 2 ). Figure 1 illustrates this error which can be reduced by refining the mesh using more grid points. Both approaches can approximate to the same value with ∆S and ∆t small enough and can be easily applied to different nonlinear volatilities models as well as different types of options.

Numerical experiments
In this section a comparison is made of two different numerical approximation methods for computing prices of European call options based on Newton's methods (NM1, NM2) and the analytic asymptotic approximation formula developed in section 2. In the asymptotic approximation formula, the Frey and Patie model (6) and the RAPM model (8) are characterized by the following parameters: (ε, δ, γ) = (ρ, 2, 1) and (ε, δ, γ) = (µ, 4/3, 1), respectively. For the finite difference Newton's Common model parameters were chosen as:σ = 0.4, E = 100, r = 0.03, q = 0, S min = 0, S max = 300, T = 1/12 and a transformation τ = T − t was used. The tolerance for Newton's iterations was set as tol = 10 −8 . The initial guess in Newton's methods at the first time level was chosen as the constant value of 1. In the subsequent temporal levels the initial guess was taken from the approximate solution at the previous time level. The fast and robust Thomas algorithm for tridiagonal solver was used in Newton's method. Calculation of integrals for the asymptotic formula was done by using the built-in Matlab function integral.

Comparison of numerical methods with explicit invariant solution
In order to ensure all the numerical solvers mentioned in section 3 are accurate, the explicit invariant solutions for the Frey and Patie model derived by Bordag in [5, (86), (87)] with parameters c = −0.05, d 1 = 0, d 2 = 30 were computed and taken as reference solutions for evaluating experimental order of convergence. The boundary conditions and initial conditions were generated from these invariant solutions.
The table containing the experimental order of convergence (or convergence ratio) is constructed from the convergence rate of the error defined as follows: Here the error Err is defined as τ ) is the solution from numerical solver, andV (S, τ ) is from

Comparison of accuracy of Newton's method and asymptotic analytic formula
In Figure 2 errors between different methods were plotted in order to analyze the changes of the numerical approximation with respect to different model parameter ρ in the Frey and Patie model and µ in the RAPM model. The error V (S, τ ) −Ṽ (S, τ ) / Ṽ (S, τ ) for S ∈ [0.5E, 1.5E] was computed with the l ∞ maximum norm where V was calculated from Newton's method andṼ was evaluated by the asymptotic formula.
The difference between Newton's method and the asymptotic formula can be reduced by taking smaller values of model parameters as shown in the Figure 2. When ρ and µ become larger, the difference increases. Notice that in the asymptotic formula, higher order terms such us O(ρ 2 ) and O(µ 2 ) are ignored. These terms can not be neglected when considering larger values of the model parameters.

Time complexity comparison of Newton's method and analytic asymptotic formula
The comparison of time complexity is based on the implementation under the same Matlab computing environment in order to ensure fair comparison. Since the CPU time depends on the software implementation, the comparison is chosen to be based The model parameters were chosen as ρ = 0.005 and µ = 0.005. For all Newton based methods the spatial variable S was stored in a vectorized form in order to speed up computation. The ratio of grid sizes was taken as ∆S/∆τ = 3600 and (∆τ ) n+1 /(∆τ ) n = 0.5. Tables 3, 4, 5 and 6 show the computation times and the values of EOTC. NM1,2(Fo) corresponds to computing the analytic form of the Jacobian matrix and the coefficients. NM1,2(Nu) corresponds to using a finite difference approximation of the Jacobian matrix and the coefficients. Abbreviation ,Asym' stands for results computed by means of the analytic approximation formula derived in section 2, and ,Frozen' stands for results obtained by the frozen coefficient method instead of the Newton one.
The results from evaluating the computational complexity and the experimental order of time complexity shows that the analytic approximation formula has the advantage when considering smaller time steps ∆τ . Hence it can be successfully adopted for model calibration using high frequency data. When all the numerical methods converge, Newton's method seems to have worse performance when compared to frozen coefficients methods as can be seen from Figure 3 which shows the number of iterates for the example with grid points M = N = 200 and ρ = 0.01. Clearly, for the first few time levels, the method of frozen coefficients requires higher   Figure 4, it is important to notice that the option price increases for asset prices close to E when these parameters are increasing. In fact ρ and µ can be calibrated using market data to observe whether the market of underlying asset has high or low liquidity.
In the calibration experiments the parameter ρ for the Frey and Patie model was calibrated by using the call option time series from Apple Inc. (AAPL) in NASDAQ quotes market. Bisection method was used in the search algorithm as described in Algorithm ??. The parameters in the calibration process were fixed as r = 0.01, E = 106, q = 0, andσ = σ impl was computed as the implied volatility from the market quotes prices. As for the solution method (,Solver') both analytic approximation formula and Newton based methods were used. Table 7 shows similar calibration results for both methods when the parameter ρ is not large using these market data. This means that the analytic approximation formula has a benefit of

Conclusion
In this paper two different linearization numerical methods for solving the nonlinear Black-Scholes equation are proposed and analyzed. Numerical results are compared in their accuracy and time complexity for the Frey and Patie illiquid market model and the risk-adjusted pricing methodology model. It turns out that the analytic approximation formula is more suitable for computation when the model parameters are sufficiently small. In particular it can be applied in calibrating parameters using market data efficiently as it is a time consuming process for a full temporalspatial finite difference approximation scheme based on Newton's method. On the other hand, the analytic approximation formula becomes restrictive as the error increases when the parameters become larger. The Newton's method is easy to implement and suits various types of nonlinear Black-Scholes equations. There are different approaches to implement Newton's method and two of them are discussed in this paper. Although time complexity is a general problem, it can be improved by combining other techniques or by using the so-called Newton-like methods to approximate the Jacobian matrix in order to reduce the number of iterates. Both techniques in fact can be extended to solve other types of nonlinear option pricing models, and the resulting numerical solutions may also be considered as a benchmark solution when exact solutions do not exist.
is given by the variation of constant formula and is given by The solution of equation (20) can be written Let us consider the change of variables in (22) and introduce the function: because β = −σ 2 2 α 2 − r (see (18)). In order to simplify further notation and let EXP denote the power of the exponential function, i.e.
Consider the function Λ(τ, ξ) defined in (22). Then the inner integral can be calculated as follows: Now let us consider the case δ = 1. Since Substituting the terms P, Q(τ, ξ), Λ(τ, ξ) yields the form of the solution u(x, τ ) as stated in Theorem 2.1.