Abstract

Transforming the nonlinear Black-Scholes equation into the diffusion PDE by introducing the log transform of S and can provide the most stable platform within which option prices can be evaluated. The space jump that appeared in the transformation model is suitable to be solved by the sparse grid approach. An adaptive sparse approximation solution of the nonlinear second-order PDEs was constructed using Faber-Schauder wavelet function and the corresponding multiscale analysis theory. First, we construct the multiscale wavelet interpolation operator based on the definition of interpolation wavelet theory. The operator can be used to discretize the weak solution function of the nonlinear second-order PDEs. Second, using the couple technique of the variational iteration method (VIM) and the precision integration method, the sparse approximation solution of the nonlinear partial differential equations can be obtained. The method is tested on three classical nonlinear option pricing models such as Leland model, Barles-Soner model, and risk adjusted pricing methodology. The solutions are compared with the finite difference method. The present results indicate that the method is competitive.

1. Introduction

Black and Scholes [1] described a mathematical framework for calculating the fair price of a European option in which they used a no-arbitrage argument to derive a partial differential equation which governs the evolution of the option price with respect to the time to expiry, , and the price of the underlying asset, ; that is, where denotes the underlying asset, , denotes the expiry date, is the volatility (measures the standard deviation of the returns), and is the riskless interest rate. The solution of the famous Black-Scholes equation provides both the price for a European option and a hedging portfolio that replicates the option assuming that the market is frictionless; thus there are no transaction costs, the interest rates for borrowing and lending money are equal, all parties have immediate access to any information, and all securities and credits are available at any time and any size. That is, all variables are perfectly divisible and may take any real number. Moreover, individual trading will not influence the price. Another unpractical assumption is that there are no arbitrage opportunities, meaning that there are no opportunities of instantly making a risk-free profit. In practice, the volatility relates to the time , the stock price , or the derivatives of the option price itself. Taking all these factors into the option pricing model, the nonlinear Black-Scholes can be obtained as follows [2]: where denotes a nonconstant volatility. In this paper, we will be concerned with several transaction cost models from the most relevant class of nonlinear Black-Scholes equations for European options with a constant drift and a variable volatility .

In contrast with the traditional linear Black-Scholes equation, it is very difficult to find the analytical solution of the nonlinear Black-Scholes model. Elbeleze et al. proposed an asymptotic analytical method for fractional Black-Scholes model using Sumudu transform [3]. Song and Wang solved the fractional Black-Scholes with the finite difference method [4]. Yan proposed the wavelet precise integration method (WPIM) for the nonlinear model [5]. In order to eliminate the boundary effect in WPIM, an interval wavelet was constructed by Liu [6] based on the restricted variational principle.

Since is lognormal distributed, the direct discretization mesh of the Black-Scholes PDE will not have equally spaced -steps, so, Schaeffer et al. [7] transform the backward Black-Scholes PDE (2) into the diffusion PDE by introducing the log transform of and . This provides the most stable platform within which option prices can be evaluated. Implementing the finite difference algorithms to the diffusion equation is much simpler than the Black-Scholes (2). The disadvantage of the new model is the equally spaced jumps in which need more refined -step to meet the precision requirement. Obviously, the sparse grid approach for the PDEs is the powerful tool to solve the diffusion equation.

The sparse representation of functions via a linear combination of a small number of basic functions has recently received a lot of attention in several mathematical fields such as approximation theory as well as signal and image processing [7]. Many physical equations contain multiscale phenomena, such as the homogenization problems from material science and chemistry and multiscale systems in biology, computational electrodynamics, fluid dynamics, and atmospheric and oceanic sciences.

The main source of difficulty in multiscale computation is that accurate simulation of the system requires all phenomena to be fully resolved [8]. From the perspective of mathematics, multiscale methods began with representation of a function using a global basis, such as Taylor series or Fourier series. More sophisticated bases such as wavelet have appeared in computational physics in recent decades. Comparing with the common wavelet function, the interpolation wavelet function possesses the interpolation characteristic which is helpful to improve the efficiency of the algorithm. The familiar interpolation wavelet functions include Haar wavelet [9], Shannon wavelet [10], Shannon-Gabor wavelet [11], and Faber-Schauder wavelet [12]. Besides, the autocorrelation function of the Daubechies scaling function is also taken as an interpolation function to construct wavelet numerical method for PDEs. The wavelet functions which have all of the excellent numerical properties such as the compact support, smoothness, orthogonality, interpolation, and exact analytical expression are the perfect weight function in the numerical method for PDEs. Unfortunately, there is not any such wavelet function that has been constructed up to now. Shannon wavelet does not have compact support property. Hoffman et al. have suggested using the Shannon-Gabor wavelet instead of Shannon wavelet in numerical method. In some ways this improves the approximation to a Dirac delta function. However, the presence of the Gaussian window destroys the orthogonal properties possessed by the Shannon wavelet, effectively worsening the approximation to a Dirac delta function. The autocorrelation function of the Daubechies wavelet has no exact analytical expression, which destroys the approximation precision to some extent. Haar wavelet is discontinuous although it has compact support and orthogonality, which cannot be employed to solve the derivative of the approximated function.

Our study in this work will prove that the Faber-Schauder’s compact support and linear expression can be employed to construct the efficiency algorithm using the variational principle, although it has no smoothness property. In our method, we construct the weak solution function of the nonlinear second-order PDEs with variational principle firstly, which can be discretized into nonlinear system of ODEs with the multiscale wavelet interpolation operator which has only first-order derivative. Then, the variational iteration method [5, 13, 14] was employed to solve the nonlinear system of algebraic equations which is obtained by the discretized to the nonlinear PDEs by the wavelet collocation method. The multiscale wavelet interpolation operator [15] and the soft threshold method were constructed to get the sparse input data. A key difference between the proposed method and other methods is that the proposed method directly resolves all the significant scales in the solution. By construct, the other methods directly resolve only the coarse scales of the solution, and they separately reconstruct the fine-scale solution (as well as its effect on the coarse scales). Another difference is the choice scheme of the sparse basis [16]. As the solution of the PDEs are enforced at every time step by simply applying soft threshold [7] to the coefficients of the basis approximation in our method, that is, the approximate solution resides on a sparse subspace of a basis, the input data were reduced with the optimization method, and the traditional sparse grid multilevel methods reduce information with the error tolerance [17], which is also be called hard threshold.

2. Faber-Schauder Wavelet Function

2.1. Definition of Faber-Schauder Wavelet

The Faber-Schauder wavelet function can be defined as

The basis of compactly supported wavelets of is formed by the dilation and translation of the single function :

The wavelet basis induces a multiresolution analysis on , that is, the decomposition of the Hilbert space into a chain of closed subspaces:(1)such that(2)and(3).By defining as a complement of in ,(4).The space is represented as a direct sum:(5).On each fixed scale , the function form a basis of , and .

Based on the multiscale analysis theory, the Faber-Schauder can be used as the basis in , to represent any function defined in space , where . On each fixed scale , the discrete points about the variable can be defined as , .

Then, the discrete form of the Faber-Schauder basis function is obtained as follows

2.2. Properties of Faber-Schauder Wavelet

It is easy to prove that Faber-Schauder function satisfies the following relations:(i)(ii)(iii)(iv)(v) has orthogonality property; that is,

Proof. First-order derivative of the Faber-Schauder scaling function at point is Obviously, is the same as Haar wavelet function, which possesses orthogonality. So, the proof is completed.

3. Construction of the Multiscale Wavelet Interpolation Operator

In this section, we try to construct a multi-scale interpolator, which is independent of the wavelet functions. Let be the interpolation scaling function, which can be any wavelet function with interpolation property. By dilation and translation, the scaling function series can be obtained as which constitute the scaling function subspaces in as As has interpolation property, that is, , the interpolator can be defined as And the corresponding wavelet function space is defined as where

According to the multiscale wavelet approximation theory, the approximation of the function , can be represented as the coefficients and are defined as respectively. These expressions show that the wavelet coefficient represents the approximation error at . In other words, the interpolation wavelet coefficient represents the local regularity of the approximated function.

In the following, we will discuss the construction of the multiscale interpolator based on the wavelet multiscale analysis theory. The interpolation operator can be defined as where is the interpolation function. According to the wavelet transform theory, function can be expressed approximately as where , and the interpolation wavelet transform coefficient can be represented as where , , , and the is the restriction operator which is defined as Suppose Substitute (19) into (17); we obtain and if, then Substitute the restriction operator (18) and the wavelet transform coefficient (19) into (16); the approximate expression of the solution function can be obtained as According to the definition of the interpolation operator (15), it is easy to obtain the expression of the interpolation operator: And the corresponding -order derivate of the interpolation operator is Substituting (23) and (24) into PDEs, the nonlinear PDEs can be changed into a nonlinear ODEs, and then the corresponding analytical solution can be obtained with VIM [13].

4. Sparse Grids Approach for Parabolic PDEs Based on Faber-Schauder Wavelet

As mentioned above, Faber-Schauder scaling function has no smoothness property, which means that the second order derivative of the Faber-Schauder scaling function does not exist. Therefore, we cannot take the Faber-Schauder scaling function as the weight function directly to approximate the solution of PDEs with second order derivative. In order to solve this problem, the weak solution form (variational equation) of the PDEs can be employed. It is well known that the solution of the variational equation equals one of the corresponding PDEs. Consider the parabolic PDEs as follows: where is the definition domain in plane and denotes the differential operator. For convenience to construct the variational equation, the parameter should be discretized as , where and . Then, can be approximated as Substituting the above equation into (2), we obtain Obviously, it is the initial-boundary elliptic PDEs. Using the virtual displacement theory, the variation equation can be obtained as where and is the Sobolev space.

According to the interpolation wavelet transform theory, the variables and can be approximated as The first-order derivatives are respectively. Substituting (30)–(32) into (16), the sparse method for the parabolic PDEs based on the Faber-Schauder wavelet will be obtained.

Sparse operation via optimization: soft threshold [7]. At a given time step, the problem of projecting the updated solution onto a sparse subset is equivalent to fitting a solution with corresponding coefficients . This can be written as a constrained least squares fit as follows: This constrained optimization problem is related to the unconstrained problem as follows: where is the vector of coefficients. In general, this can be computed for a nonorthonormal basis, which is equivalent to a basis pursuit problem with norm as a sparse regularizer. The resulting minimizer is a proximal solution that lies on a sparse subset of the original coefficient domain. We can apply the soft threshold on the coefficients directly in order to induce sparsity in this way.

5. Numerical Experiments and Discussion

The Black-Scholes model requires a continuous portfolio adjustment in order to hedge the position without any risk. In the presence of transaction costs it is likely that this adjustment easily becomes expensive, since an infinite number of transactions are needed. Thus, the hedger needs to find the balance between the transaction costs that are required to rebalance the portfolio and the implied costs of hedging errors. The assumption of continuous portfolio adjustment is awkward in the presence of nonzero transactions cost since continuous revision implies infinite trading. Up to now, there are 3 classical modified options replication strategy which have been using widely in recent years. The difference among them is the volatility appeared in (2).

In this section, we are taking these 3 nonlinear Black-Scholes equations (2) as the examples to test the numerical algorithm described in the previous section. In order to solve the problem, it is necessary to perform a variable transformation as follows: Substitute (35) into (2); the following diffusion equation can be obtained: where Initial condition is And the boundary conditions are According to the transformation relation (35), it is easy to understand that the point is corresponding to the strike price . Obviously, the initial solution curve is smooth except around the position , where a sharp steep wave happened (Figure 1). So, an adaptive sparse grid approach is suitable for this problem.

Using the difference coefficient to approximate the partial differential operator , the Burgers equation becomes According to the virtual displacement theory, the variational form of the nonlinear Black-Scholes equation can be represented as: Substituting (30)–(32) into (41), we obtain where and . This is a system of algebraic equations, which can be solved easily.

5.1. Numerical Results
5.1.1. Leland’s Model

Assume that the trade just can do at discrete times: Leland deduces that the option price is the solution of the nonlinear Black-Scholes equation (2) with the modified volatility [18]: where represents the original volatility and Le the Leland number given by where denotes the transaction frequency, the transaction cost is , denotes the round trip transaction cost per unit dollar of the transaction, and the number of assets bought or sold at price is proportional to the monetary value of the assets bought or sold.

In the following experiment, , , riskless rate , , strike price , and expiry date .

The numerical solution is shown in Figure 2. It illustrates that the wavelet sparse grid approach proposed in this paper can capture the steep slope appearing in the solution; that is, there are more grid points around the boundary and the position around the strike price. Obviously, the adaptability not only improves the approximate precision but also decreases the calculation amount. The grid points of the finite difference method are collocated evenly, so, 4097 (= 212) grid points are needed to get the same precision with the wavelet sparse grid approach. In other words, comparing with the finite difference method, the retained coefficients are less than 4% in the sparse solution of the nonlinear Black-Scholes equation. Absolutely, this is helpful to improve the efficiency of the algorithm.

The numerical solution of the Leland model (2) obtained by the wavelet sparse grid approach is different from the solution of the Black-Scholes equation (1). The difference is shown in Figure 3. The main difference appears around the expiry price, which illustrates the impact of transaction cost on the option price. With the increase of the transaction size around the expiry price, the option price volatility brought by the transaction cost becomes larger and larger. This shows that the proposed method can solve the Leland model correctly.

5.1.2. Barles-Soner Model [2]

Using an exponential utility function, a more complicated model was constructed by Barles and Soner as follows: The difference between Black-Scholes model and Barles- Soner model obtained by the wavelet sparse grid method is shown in Figure 4. Comparing with the Leland model, the Barles-Soner model is more close to the linear Black-Scholes model and the numerical result is consistent with it. For comparison, the values of the parameters are the same as those in the Leland model.

5.1.3. Risk Adjusted Pricing Methodology [2]

The objective of this model is to find the optimal time-lag between the transactions. This can minimize the sum of the rate of the transaction costs and the rate of the risk from an unprotected portfolio. In this model, the modified volatility is defined as where is the transaction cost measure and is the risk premium measure.

The difference of the solution between the Black-Scholes model and the risk adjusted pricing methodology is shown in Figure 5. Similar to the Leland model, the numerical result illustrates the function of the model, which is consistence with option price in practice.

5.2. Comparison of the Algorithm Precision and Efficiency

In this section, we are taking the linear Black-Scholes equation to illustrate the advantage of the wavelet sparse grid approach proposed in this paper. Its exact analytical solution is represented as where is the call price, is the underlying asset price, is the strike price, is the riskless rate, is the maturity, is the volatility, and express the normal distribution. The comparison between the wavelet sparse grid approach and the finite difference method is shown in Figure 6.

It is obvious that the error arises around the strike price, which expresses the nonlinear B-S model, and the numerical approaches are valid. While the call option price is going far away from the strike price, the error of the finite difference method becomes larger instead of close to the boundary condition, which is due to the discontinuousness of the boundary condition. And the error of the wavelet sparse grid approach is becoming smaller and smaller; this due to the adaptability of the proposed method. There are more of the sparse grid points collocate around the boundary, which is helpful to improve the precision near the boundary.

6. Conclusions

The space jump that appeared in the nonlinear Black-Scholes equation is suitable to be solved by the sparse grid approach. Although Faber-Schauder wavelet has no second-order derivative, combining the variational principle, we have proposed a method to resolve fully the solutions of multiscale PDEs while only retaining important nodes based on the multilevel interpolation wavelet theory and the optimization principle. The reduced dynamics created by the sparse projection property capture the true phenomena exhibited by the solution. This sparse projection amounts to a shrinkage of the coefficients of the updated solution at every time step. The multiscale interpolation operator and the soft threshold employed in this paper not only work for the Faber-Schauder wavelet, but also work for any interpolation wavelets. Comparing with the finite difference method, the retained coefficients are less than 4% in the sparse solution of the nonlinear Black-Scholes equation.

Conflict of Interests

The author declares that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (no. 41171337) and the National Key Technology R&D Program of China (no. 2012BAD35B02).