Coupling Technique of Haar Wavelet Transform and Variational Iteration Method for a Nonlinear Option Pricing Model

: Compared with the linear Black–Scholes model, nonlinear models are constructed through taking account of more practical factors, such as transaction cost, and so it is difﬁcult to ﬁnd an exact analytical solution. Combining the Haar wavelet integration method, which can transform the partial differential equation into the system of algebraic equations, the homotopy perturbation method, which can linearize the nonlinear problems, and the variational iteration method, which can solve the large system of algebraic equations efﬁciently, a novel numerical method for the nonlinear Black–Scholes model is proposed in this paper. Compared with the traditional methods, it has higher efﬁciency and calculation precision.


Introduction
In 1958, Modiglian and Miller used the no-arbitrage equilibrium analysis method in modern finance when studying how the financial activities of the company create value. Referring to the parameters of the analysis method, the original Black-Scholes model derives partial differential equations (PDEs) and simulates price changes over time [1]. In this model, there are too many unpractical assumptions, such as no-arbitrage, no transaction cost, no arbitrage opportunities, etc. Therefore, the application of the nonlinear Black-Scholes equation has gradually expanded to more fields, such as property value evaluation, option pricing, and stock prices, etc. Considering the above factors in the nonlinear Black-Scholes model, a more accurate fixed value can be calculated when more realistic parameters are put in. The nonlinear Black-Scholes equation is as follows [2]: where S(t) denotes the underlying asset, t ∈ (0, T), T denotes the expiry date, σ is the volatility (measures the standard deviation of the returns), r is the riskless interest rate, and σ denotes a non-constant volatility. It is almost impossible to find the exact analytical solution of the nonlinear Black-Scholes model. Elbeleze et al. applied the homotopy perturbation Sumudu transformation method (HPSTM) to the fractional Black-Scholes option pricing equation and successfully obtained approximate and numerical solutions. [3]. Rashidinia and Jamalzadeh proposed a modified B-spline collocation approach for pricing American-style Asian options [4]. Chen proposed a new operator splitting method for the American option problem; this method can solve the fractional Black-Scholes American option model effectively by decomposing the linear complementarity problem into a linear boundary value problem and an algebraic system [5]. After obtaining the optimal algebra, Mandal simplified the given fractional partial differential equations (FPDE) to fractional where ξ 1 = k/m, ξ 2 = (k + 0.5)/m, ξ 3 = (k + 1)/m. The integer m = 2j, (j = 0, 1, . . . , J) indicates the level of the wavelet, and J is the maximal level of resolution, while k = 0, 1, . . . , m − 1 is the translation parameter. The subscript of the Haar wavelet function h is defined as i = m + k + 1, and so i ∈ [2, 2 J+1 ]. It should be pointed out that the Haar scaling function corresponds to the subscript i = 1. Obviously, Haar function can be viewed as the Daubechies wavelet of order 1. For a long time, it has been thought that the Haar wavelets are valueless in the numerical calculations for nonlinear PDEs as they are not continuous. This shortcoming has been overcome by using regularization techniques with interpolation splines [9,10] to some extent, although this method makes the solution process more complicated. Another effective method is the Haar wavelet integration technique proposed by Chen and Hsiao [11,12]. In this method, the highest derivative function appearing in the PDEs is approximated by the Haar series, and the other derivatives are obtained through integrations. In order to solve the nonlinear Burgers equation with Dirichlet boundary conditions, Jiwari proposed an effective numerical scheme based on the Haar wavelet and quasi-linearization process [13]. Hariharan and his associates [14,15] used the Haar wavelet transform method to solve the equations similar to hyperbolic partial differential equations and one-dimensional convection-diffusion equations, and the parameter values can be in any region during the solving process. Compared with other wavelet numerical methods, the Haar wavelet transform method for solving equations is not only simple, but also lower in computational cost [16]. Therefore, the Haar integration method for PDEs has been developed and widely used in various branches of sciences and engineering [17][18][19][20][21]. As mentioned above, the homotopy perturbation method and variational iteration method are effective tools to solve the nonlinear problems, which have been developed and applied to solve various nonlinear problems by Ji-Huan He [22] and by others [23]. Unlike analytical perturbation methods, HPM does not depend on small parameters which are difficult to find, and so it has been widely used to solve various nonlinear problems [24].
In this paper, combining the Haar integration method, the variational iteration method, and the homotopy perturbation method, a coupling technique for solving the nonlinear Black-Scholes model is proposed. This method possesses both of the merits appearing in two methods. Compared with other methods, both the iteration time step and the precision are improved effectively.

Discretization format of Nonlinear Black-Scholes Model
In order to solve the nonlinear Black-Scholes Equation (1), it is necessary to perform a variable transformation as follows: Substituting Equation (2) into Equation (1), the following equation can be obtained: where D = 2r The initial condition is and the boundary conditions are According to Equation (2), the curve of the initial solution will have a sudden change when x = 0 at the same time, S = K. Obviously, when x = 0, the curve of the initial solution in the value space is smooth [25]. It can be seen that the adaptive numerical solution method can avoid the sudden change of the solution when x = 0, when using the Haar wavelet numerical method to solve the equation, the integration process is as follows: As mentioned above, ξ 1 = k/m, ξ 2 = (k + 0.5)/m, ξ 3 = (k + 1)/m. The integer m = 2j, (j = 0, 1, . . . , J) indicates the level of the wavelet, and J is the maximal level of resolution, and k = 0, 1, . . . , m − 1 is the translation parameter. Without loss of generality, the discrete point sequence of the variable τ ∈ [0, T] can be defined as The discrete points are usually defined as According to the Haar wavelet integration method [16], the two-order derivative term in the nonlinear Black-Scholes equation can be approximately represented as ∂ ∂τ where a s is taken as a constant approximately in the subinterval τ ∈ [τ s , τ s+1 ]. Integrating Equation (12) with respect to τ, it can be obtained that Integrating Equation (13) with respect to x from −∞ to x, it can be obtained that Similarly, According to the boundary condition (7), it is easy to understand that Obviously, it should be a constant c instead of the infinity in the boundary condition in order to get the numerical solution. As a matter of fact, it is big enough to meet the precision requirement when c > 100. That is, Substituting Equation (20) into Equation (15), we obtain Similarly, substituting Equation (20) into Equation (16), we obtain Mathematics 2021, 9, 1642

Solve the Large System of Algebraic Equations
There are many ways to solve the large system of algebraic equations, such as the Gaussian elimination method, finite-difference schemes, the Newton iteration method, and so on.

Finite-Difference Schemes
Although there are several ways to solve the problem, here is only the solution of the finite-difference method as an example, and an explicit form is given [2].
There are several numerical methods of solving the problem. This work's emphasis is on the finite-difference schemes, thus other methods will not be mentioned here.
To apply finite-difference schemes to the transformed problem subject to the conditions with the corresponding volatilities, we start by replacing x ∈ R and τ ∈ [0, T] by a bounded interval x ∈ [−R, R], R > 0. We discretize the new computational domain by a uniform grid ( We denote the approximate solution of the problem in x i at time τ n by U n i ≈ u(x i , τ n ) and treat the initial and boundary conditions in the following way: For a more appropriate treatment of the boundary conditions, so-called articifial boundary conditions [26] can be introduced to limit the unbounded spatial domain of the problem. We bear in mind that we say a scheme is of order (m, n) if its truncation error is of order o(k m + h n ).
To discretize the problem we introduce the following notation for the forward difference quotient with the spatial step size h: where we leave out the error term o(h) Similarly, the backward difference quotient with respect to the spatial variable is denoted as and the central difference quotient as For the second spatial derivative we introduce the standard difference quotient with the error term o(h 2 ). Note that central differences in the time variable are never used in practice because they always lead to bad numerical schemes that are inherently unstable.

Variational Iteration Method for System of Algebraic Equations
In order to improve efficiency, the variational iteration method is introduced to solve the wavelet coefficient.
Then, Equation (25) can be rewritten as the matrix format as follows Equation (35) can be rearranged and simplified as where The solution of the set of Equation (36) is the wavelet coefficient, and so the fast algorithm for Equation (14) is very important to solve the Black-Scholes model. There are many ways to solve the linear simultaneous equations, such as the partial pivoting, the total pivoting, the Jacobi iterative method, the Gauss-Seidel iterative method, and so on. The variational iterative method [27] proposed by him is a more accurate approximation algorithm compared with the classical iterative method. In this method, a simple identified method for the Lagrange multipliers is proposed by making the correction equations stationary. To solve the set of Equation (36) by means of the variational iteration method, 2M Lagrange multipliers should be introduced to correct the initial approximation solution  1 , α 1,2 , . . . , α 1,2M ) T and α 2 = (α 2,1 , α 2,2 , . . . , α 2,2M ) T are the Lagrange multipliers. They could be solved as follows As a matter of fact, this method is the same as the Newton iterative method, as the Lagrange multiplier α 2 = 0.

Numerical Experiments and Discussion
According to Equation (2), when x = 0 the curve of the initial solution will have a sudden change; at the same time, S = K. Obviously, when x = 0, the curve of the initial solution in the value space is smooth (Figure 1). Therefore, the multi-resolution method is suitable in solving the problem.
Mathematics 2021, 9, x FOR PEER REVIEW 9 of 16 solution in the value space is smooth ( Figure 1). Therefore, the multi-resolution method is suitable in solving the problem. In this section, we try to solve three nonlinear models to test the efficiency of the approach proposed in this paper.

Leland's Model
If the transaction cost is included in the option transaction, it means that the Black-Scholes arbitrage argument is invalid in the discrete place. Leland proposed that the option price is the strategy solution of the nonlinear Black-Scholes Equation (1) by adjusting the volatility in the Black-Scholes Equation [28]. In this section, we try to solve three nonlinear models to test the efficiency of the approach proposed in this paper.

Leland's Model
If the transaction cost is included in the option transaction, it means that the Black-Scholes arbitrage argument is invalid in the discrete place. Leland proposed that the option price is the strategy solution of the nonlinear Black-Scholes Equation (1) by adjusting the volatility in the Black-Scholes Equation [28].
where σ represents the original volatility and The Leland number Le is Among them, δt is the frequency of option transactions, the transaction cost is k | ∆ | S / 2, and k is the transaction cost (per unit dollar). ∆ is the number of assets purchased when ∆ > 0, and the price of the number of assets sold when ∆ < 0.S will increase linearly with the increase of the monetary value of the asset purchased or sold.
It is obvious that the solution of the Leland model (43) is different from the one of the Black-Scholes equation. It can be seen from Figure 2 that the option price curve volatility has been significantly deformed in the range of 50-150 option prices, and with the increase of transaction costs within this range, the volatility of option prices will also increase. This shows that this research method can solve the Leland model correctly.  Although the Haar wavelet is the simplest wavelet, it has many excellent properties, such as compact support, orthogonality, interpolation, and precise analytical expressions. All these can be employed to act as the perfect weight function in the numerical method for PDEs. The sparseness of the wavelet coefficients is one of the appearances which is showed in Figure 3. In solving the Leland model, there are 4097 discrete points that are considered. Most of the wavelet coefficients corresponding to these discrete points are near to zero, except 157 coefficients that are greater than 0.005. Obviously, this is helpful to improve the efficiency of solving the set of the Equations (36).  Although the Haar wavelet is the simplest wavelet, it has many excellent properties, such as compact support, orthogonality, interpolation, and precise analytical expressions. All these can be employed to act as the perfect weight function in the numerical method for PDEs. The sparseness of the wavelet coefficients is one of the appearances which is showed in Figure 3. In solving the Leland model, there are 4097 discrete points that are considered. Most of the wavelet coefficients corresponding to these discrete points are near to zero, except 157 coefficients that are greater than 0.005. Obviously, this is helpful to improve the efficiency of solving the set of the Equations (36). such as compact support, orthogonality, interpolation, and precise analytical expressions. All these can be employed to act as the perfect weight function in the numerical method for PDEs. The sparseness of the wavelet coefficients is one of the appearances which is showed in Figure 3. In solving the Leland model, there are 4097 discrete points that are considered. Most of the wavelet coefficients corresponding to these discrete points are near to zero, except 157 coefficients that are greater than 0.005. Obviously, this is helpful to improve the efficiency of solving the set of the Equations (36).

Barles-Soner Model
Employing an exponential utility function, Barles and Soner derived an complicated model and proved that as ε and k tend to zero, V is the unique solution of Equation (1), in which the modified volatility σ is defined as

Barles-Soner Model
Employing an exponential utility function, Barles and Soner derived an complicated model and proved that as ε and k tend to zero, V is the unique solution of Equation (1), in which the modified volatility σ is defined as In the following experiment, k = 0.05, riskless rate r = 0.1, σ = 0.2, strike price K = 100, and expiry date T = 1. The numerical solution of the model is showed in Figure 4. It is obvious that the difference between the Barles-Soner model and the linear model is not symmetric, but decreases closer to the expiry date. The maximal difference appeared around S ≈ 90 at one year to expiry. Compared with the Leland model, the Barles-Soner model provides a higher price; this is consistent with the known facts.  Figure  4. It is obvious that the difference between the Barles-Soner model and the linear model is not symmetric, but decreases closer to the expiry date. The maximal difference appeared around 9 0 S ≈ at one year to expiry. Compared with the Leland model, the Barles-Soner model provides a higher price; this is consistent with the known facts.

Risk Adjusted Pricing Methodology
Risk adjusted pricing methodology is a model proposed by M.Jandačka and D. Ševčovič, and is model for pricing derivative securities in the presence of both transaction costs as well as the risk from unprotected portfolio. Therefore, this model should provide a higher price than the Barles and Soner model, in theory. In this model, the modified

Risk Adjusted Pricing Methodology
Risk adjusted pricing methodology is a model proposed by M.Jandačka and D. Ševčovič, and is model for pricing derivative securities in the presence of both transaction costs as well as the risk from unprotected portfolio. Therefore, this model should provide a higher price than the Barles and Soner model, in theory. In this model, the modified volatility is defined as [29] where M ≥ 0 is the transaction cost measure and C ≥ 0 is the measure of risk premium. In the following experiment, riskless rate r = 0.1, σ = 0.2, strike price K = 100, and expiry date T = 1. The solution difference between the Black-Scholes model and the risk-adjusted pricing method is shown in Figure 5. Similar to the other two models, the numerical results illustrate the function of the model, which is consistent with option prices in practice and in theory.

Comparison of the algorithm precision
In this study, in order to illustrate the advantages of the Haar wavelet method, the linear Black-Scholes equation is used in this section to illustrate the analysis formula and is represented as C is the purchase price when the budget is bullish, K is the asset execution price, r represents the risk-free interest rate, T is the maturity, the interval period of the σ index change rate, and 1 ( ) N d conforms to a normal distribution. The comparison between the Haar wavelet method and the finite-difference method is shown in Figure 6.

Comparison of the Algorithm Precision
In this study, in order to illustrate the advantages of the Haar wavelet method, the linear Black-Scholes equation is used in this section to illustrate the analysis formula and is represented as C is the purchase price when the budget is bullish, K is the asset execution price, r represents the risk-free interest rate, T is the maturity, the interval period of the σ index change rate, and N(d 1 ) conforms to a normal distribution. The comparison between the Haar wavelet method and the finite-difference method is shown in Figure 6.  Figure 6 shows that the numerical solution of the Haar wavelet numerical method has not changed significantly as the time step ∆τ increases, along with the amount of the space discretization points compared with the finite-difference method. The maximal error appeared around the point S = 300 as ∆τ = σ 2 80 , M = 4096. With the decrease of the parameter ∆τ, the numerical precision obtained by the finite-difference method becomes better accordingly, and the maximal error appeared around the point S = 95. Unfortunately, if the ∆τ decreases to σ 2 400 , the finite-difference method would be invalid even the amount of the space discretization points increases to M = 8192. This relates to the condition number of the matrix from the finite-difference method [30,31]. It is no doubt that the choice of L can change the condition number of the system of algebraic equations discretized by the wavelet interpolation operator or the finite-difference method [32]. Therefore, the choice of L should take the condition number into account. In fact, if the condition number Cond(A) = 10k, then you may lose up to k digits of accuracy on top of what would be lost to the numerical method, due to loss of precision from arithmetic methods [33]. According to the general rule of thumb, the choice should follow the rule as follows:

Analysis of Experimental Results
Aiming at the problem that it is difficult to find exact solutions to large algebraic equations, we used the method proposed in this paper to calculate the solutions of the Leland model, the Barles-Soner model, and the risk-adjusted pricing method model, respectively, to test their effectiveness. The results show that when the volatility of the option price curve fluctuates within a certain range, as the transaction cost increases within this range, the volatility of the option price will also increase. This research method can correctly solve these three models. The disadvantage is that, in order to solve the Barles-Soner model, although the research method can be applied to a series of higher-priced options, due to the difference between the Barles-Soner model and the linear model, large errors will occur.

Robustness of the Algorithms
Robustness is an important indicator of the evaluation algorithm. The wavelet numerical method is not sensitive to the time step, and it can improve the computational efficiency, compared with the finite-difference method. The variational iterative method is an efficient approach to solve the set of equations, which is a generalization of the Newton iterative algorithm in itself. The combination of the Haar integral method and variational iteration method is very effective for solving various nonlinear PDE problems. In this research, the multi-scale Haar wavelet numerical method is used to solve the nonlinear Black-Scholes equation, which solves the space jump problem in the equation.

Conclusions
In the work of this paper, the coupling technology of the Haar wavelet transform and the variational iteration method, a new numerical solution method, is used to solve the problem of a nonlinear approximate solution of option pricing. Among them, the method based on Haar wavelet integration transforms partial differential equations into algebraic equations. This algorithm can transform the nonlinear approximate solution of option pricing into the approximate solution of algebraic equations. This conversion solution process makes the derivative function near the boundary of the defined domain continuous, which can eliminate the boundary effect, so that the accuracy of the numerical solution in the entire defined domain is improved.
On the other hand, in order to improve the efficiency of solving large-scale algebraic equations, another method is introduced into our work that the variational iteration method is used to calculate the results of wavelet coefficients. In order to be effective in the proposed algorithm, we analyzed the numerical application in detail. The accuracy of the simulation results of a variety of Black-Scholes model solutions is compared, including the comparison of the simulation results of the algorithm in this research with the existing Leland's model, Barles-Soner model, and other algorithm models.
The presented numerical results and comparisons show that the proposed multiscale Haar wavelet numerical method can effectively solve the space jump problem in the equation when solving, and then achieve the purpose of solving various nonlinear Black-Scholes equations, which proves the method proposed in this paper strengthens robustness, efficiency and calculation accuracy.