Nonparametric model calibration for derivatives

Consistently ﬁtting vanilla option surfaces is an important issue in derivative modelling. In this paper, we consider three diﬀerent models: local and stochastic volatility, local correlation, hybrid local volatility with stochastic rates, and address their exact, nonparametric calibration. This calibration process requires solving a nonlinear partial integro-diﬀerential equation. A modiﬁed alternating direction implicit algorithm is used, and its theoretical and numerical analysis is performed.


Introduction
One of the most important challenge for real-life applications of a model to derivatives trading is the issue of calibration. Similar to common situations in many areas of physics and engineering, once a model has been suggested, its parameters have to be estimated using external data. In the case of derivative modelling, those data are the liquid (tradable) options, generally called vanilla products. It is well-known since the pioneering work of [1] and its celebrated extension by [2] that the knowledge of market data such as the prices of vanilla options across all strikes and maturities is equivalent to the knowledge of the risk-neutral marginals of the underlying stock distribution. Here, we are interested in applications of this result to three different cases. The first one is quite classic, and is the starting point of our work: local and stochastic volatility models. Such models are very useful in practice, since they offer both the flexibility and realistic dynamics of stochastic volatility, and the exact calibration properties of local volatility. The problem of calibrating local and stochastic volatility models has been dealt with for a while now, for instance by [3]. However, practitioners seem to agree that the stability of its resolution becomes uncertain when the volatility's volatility is too large. We shall see that indeed some kind of instability appears, and offer explanations to the phenomenon. The second case we focus on is the correlation between assets. Empirical measures give a certain set of results. However, when modelling a basket on multiple underlyings, a problem occurs. If one uses local volatility models for each underlying and correlates their brownian motions using the empirical correlation, the basket obtained will not reproduce the vanillas quoted on the market. This raises significant issues when hedging products on multiple underlyings. One of the solution for this problem is the known 'local correlation' approach: the correlation matrix for the n underlyings is deformed using a parameter, function of the time and the basket level. Here, we use that approach to obtain a calibration equation for the basket, relatively similar to the one appearing in the local and stochastic volatility model, and then numerically solve said equation in a two-underlyings framework. The last topic we shall be interested in are interest rates, we study a hybrid model: local volatility with stochastic rates. Using a partial differential equation approach similar to the local correlation and the local and stochastic volatility, we write a calibration equation for the vanillas of this hybrid model, solve it and verify the accuracy of the fit. The general form of our calibration equations is nonlinear partial and integro-differential. For their resolution, we chose to adapt the alternating direction implicit scheme (very efficient to solve classic linear second order parabolic equations, [4]). Being in a nonlinear non-local framework, many questions arise. Is it relevant to use ADI algorithms to solve the equations stemming from our calibration problems? How should we deal with the nonlocal term? Is the finite difference scheme we chose consistent, and what is the order of the truncation error? Can we detect an instability in certain cases? Is it possible to explain it? The aim of this work is to address and at least partially answer these questions. The paper is organized as follows. In Section 2, we quickly present the case of the local and stochastic volatility model. Section 3 is devoted to the local correlation, its calibration equation and the fit we obtain in the case of a basket on two underlyings. In Section 4, we do the same thing in the stochastic rates frame. Finally, in Section 5 and 6, we adress the questions raised previously concerning the ADI algorithm used for the resolution, Section 7 is a brief conclusion.
2 Local and Stochastic Volatility models 2

.1 Partial Integro-Differential Equation for the calibration of LSV models
The diffusion model is assumed to be the following (S t , t ≥ 0) is the stock price process and (y t , t ≥ 0) the stochastic component of the volatility. The function b simply transforms that factor into a proper volatility. a is the local volatility part of the model, exactly as in Dupire's formula, its value shall be specified depending on the aimed vanillas. ξ is the volatility of the volatility factor (commonly called 'vovol') and µ is a drift term. W 1 and W 2 are one-dimensional standard brownian motions with correlation ρ.
Let us now consider a surface of vanilla prices C(T, K) and the corresponding Local Volatility σ D . Under suitable regularity and ellipticity assumptions, the following proposition can be proved Proposition 1. The diffusion model defined by (1-2) has a density p(t, S, y) with respect to Lebesgue's measure. Moreover, if the model fits the surface of vanillas C(T, K) then necessarily a 2 (t, S) = σ 2 D (t, S) R p(t, S, y)dy R b 2 (y)p(t, S, y)dy Proof. The exact assumptions and the existence proof can be found in [5]. Here, the main concern is the calibration result. Let us assume that the model fits exactly the surface C. Letting (S 0 , y 0 ) denote the initial state of the system, the joint density p(t, S, y) of the couple (S t , y t ) verifies Kolmogorov forward equation Applying Fubini, let q(t, S) = R p(t, S, y)dy be the first marginal density of our couple. It is possible to integrate the previous equation and obtain In the case of a local volatility model (b = 1 and a = σ D ), the density q D of the Spot process solves the equation The vanillas of the LSV model being perfectly fitted, we have q = q D . Identifying the terms in the two last formulas gives Using this proposition, and reintroducing the value of a in Kolmogorov forward equation, the joint density p(t, S, y) is then solution of the nonlinear partial integro-differential equation There is thus equivalence between the existence of a model (1-2) that calibrates the vanillas C and the existence of a solution p to the pide (4).
Remark 2.1. The quotient R b 2 pdy R pdy is nothing but the conditional expectation of the volatility squared, knowing the spot process. This result is not original in itself (by applying the theorem from [6] for instance), the partial differential equation method however is unusual, and will be used on the other models as well.
The theoretical study of equation (4)-(5) can be found in [7]. Existence of solutions is proved under strong assumptions (especially on b, which must be sufficiently close to a constant). The general resolution remains an open problem.

Numerical results
It seems to be well-known among practitioners, that instabilities occur in their calibration when the volatility's volatility (in the notations, function ξ) is too large. This seems to confirm the theoretical limitations met trying to prove the global existence of a solution: when the function b oscillates too much (a change of scale in the factor y t clearly shows the equivalence between a b that moves a lot and a large ξ), the resolution of the equation is not guaranteed anymore. To assess these statements, we considered our problem from a practical viewpoint. In this section, the calibration that stems from solving the partial differential equation (4) is studied, for two stochastic volatility models: a lognormal one and a 'Cox-Ingersoll-Ross' process. The details of the algorithm used for the resolution and a study of the instabilities will be treated later, in sections 5 and 6.

Lognormal volatility
Starting with a simple mean reverting model for the volatility factor, the function b is chosen as an exponential with a 2 (t, S) = σ 2 D (t, S) R p(t, S, y)dy R exp(2y)p(t, S, y)dy Equation (4) is solved using the functions we just chose and the local volatility σ D associated to the Eu-roSTOXX 50 implied volatility surface of 2009/04/02. Once function p, density of the couple (S t , y t ), is found, we compute the vanilla prices for different strikes and maturities using this density. To have a point of comparison, the same prices are also computed with the local volatility σ D , both of them are then compared to the targeted prices (column TP). Let us then plot the error between the original vanillas and the ones obtained with the model. The calibration is quite efficient, the errors are equivalent to the ones of the local volatility model.

Cox-Ingersoll-Ross process
We also focus on the calibration of a model inspired from the interest rates framework: the volatility is assumed to follow a CIR process.
Detailled properties of this process are described by [8]. In particular, as long as 2κα > γ 2 , y t is strictly positive a.e. Once again, equation (4) is solved with this stochastic volatility. 5

Application to the 'Local Correlation' model
In this chapter, we are interested in the calibration of a market with n stocks and a basket on those stocks. The purpose is to define a diffusion model for those underlyings that is able to reproduce their implied volatility surface as well as the one of the basket. The notations are the following, let (S i t ) 1≤i≤n denote the n stocks involved in our problem. The basket's value is given by where the set (w i ) 1≤i≤n stands for the weights of the different underlyings. They are assumed to be constant in the rest of our work. Let us also fix n + 1 surfaces of vanillas (C i (T, K)) 1≤i≤n and C B (T, K).

Inconsistencies between stock and basket options
The naive approach to solve this problem is simply to consider n local volatility models The functions σ i are easily determined to fit the surfaces (C i (T, K)) 1≤i≤n with this diffusion. The correlation matrix ρ = (ρ ij ) 1≤i,j≤n associated to the standard brownian motions W i t of each underlying can be estimated with historical data. The model is now entirely defined. By equation (10) of B t , the vanilla prices for the basket are completely determined and are equal to However, there is no particular reason for the surface computed in this framework to be equal to C B (T, K). In fact, the skew of the basket is more pronounced on the market than in a model with constant correlations between the underlyings, [9].

'Local Correlation' model
In the manner of B. Dupire who decided to let the volatility depend on the level of the spot process, a degree of freedom is added to our model by distorting the matrix of correlation with a function of B t . This method appears for instance in [5,10] where the basket level induces some feedback on the values of the different underlyings. In our context, the new correlation matrix is taken as a linear combination of ρ and the constant matrix with only 1 as coefficients. This matrixρ is equal tõ We shall see while writing the calibration equation that λ has to be chosen as a function of the time and of B t . The matrixρ ij can be seen as an analogous of Dupire's local volatility, a 'Local Correlation' so to speak. Assuming that ρ is a proper correlation matrix (definite, positive...), and that the coefficients of the diffusion have a sufficient regularity, our diffusion model then possesses a density in the more general case of a matrix ρ function of the couple (t, B). It is also possible to write a condition for the vanillas of the model to be fitted Proposition 2. The diffusion model defined by (11) with a correlation function of the couple (t, B), has a density p(t, S 1 , .., S n ) with respect to Lebesgue's measure. Furthermore, this model calibrates the surface Proof. The existence of the transition density p(t, S 1 , .., S n ) stems from the assumptions made on the regularity of the coefficients, and on ρ, for more details see [5]. We can now write the calibration problem for the vanillas of the basket B t . The density just defined satisfies Kolmogorov forward equation To ease the problem, it is useful to change the coordinates (S 1 , ..., S n ) into (B, S 2 , ..., S n ) with S 1 defined by (14). After computations, the equation becomes .., S n ) are defined above. Integrating the equation against the variables (S 2 , ..., S n ), and writing q = 1 w1 pdS 2 ..dS n , the density of the marginal law of B satisfies Comparing this equation to Dupire's equation for the local volatility σ B ∂q ∂t − 1 2 if our model reproduces the vanillas σ D , then the following equality must be verified Reciprocally, the condition we just wrote is clearly sufficient for the options to be calibrated Remark 3.1. Let us note that this condition is written as an equality between two functions of the time and of B. The other variables are no longer represented.
Assuming that condition (13) is not verified, the model defined by (11) does not fit the vanillas of the basket B t . It has to be enriched to solve the calibration problem. Our choice is to distort the correlation matrix. The new matrixρ is described by (12). Hence, let Θ denote the matrix Θ ij = 1 for all 1 ≤ i, j ≤ n. We also notice that, the trace of ρ being equal to n, its smallest eigenvalue K ρ is smaller than 1.
to be stricly positive. The diagonal coefficients ofρ are still 1. As for the other terms, thanks to the first term in relation (15), they still belong to the interval ] − 1, 1[. We introduce the new correlation matrix in condition (13), this gives is a function of B and t. It is now possible to use this value to write a pide on the densityp. Any solution of the following equation is a density that calibrates the vanillas of the basket where L 1 is linear and verifies is the nonlinear part of the equation Remark 3.2. The operator L 1 + L 2 stems from a change of coordinates on a uniformly elliptic operator. It is also elliptic, uniformly on any domain where theβ i are bounded away from 0 by a strictly positive constant.
Furthermore, the initial condition is where S i 0 is the market value at instant 0 of the i-th stock. Applying this initial condition to (16), the initial value of λ is equal to For a theoretical study of the calibration equation, we refer the reader to [5].

Resolution of the equation for the calibration of a basket
This subsection focuses on the results of the calibration for a two-underlyings basket. Let us consider two assets, both of them are assumed to generate the following implied volatility surface Distorting this theoretical surface by a factor of 0.9, and thus making the prices of the basket inconsistent with the prices of the underlyings, the calibration algorithm is applied. Solving the partial integro-differential equation (17)  Here are some results for other tests. We calibrate a model with the following parameters ρ 12 = 0, w 1 = 0.3 and w 2 = 0.7, the targeted surface is the theoretical one distorted with two factors: first 0.95 and second 1.05. The results are rather satisfactory, especially at the money. To keep the computations to a reasonnable duration, a sparse initial surface was used. This explains why the calibration is not better far from the money, the fitting method is nontheless valid. Now follows an outlook of the values taken by the new correlationρ at different maturities when the theoretical surface is distorted by factors 0.95 and 1.05. The parameters are: ρ 12 = 0.5, w 1 = 0.7 and w 2 = 0.3. As expected, the Local Correlation and the distorsion factor evolve in the same direction. The underlyings must be more correlated when the implied volatility of the basket is higher, and reciprocally. Furthermore, it appears that in the case of the 0.95 distorsion, the correlation has to violently decrease for high values of B: the two underlyings must be anti-correlated when they are both large. As for the influence of the maturity, let us first state that in the computation of λ, when the denominator is smaller than 10 −6 , we chose not to change the correlation, to avoid numerical errors. It appears that, as long as B is in a zone where λ was actually computed, the framework chosen to test the calibration actually generates a local correlation constant in time.

Application to stochastic interest rates
This section is dedicated to hybrid local volatility models with stochastic rates. The interest rate is assumed to be stochastic and to follow a diffusion equation. The volatility depends on the level of the spot process exactly as in a local volatility model. The idea is to compute its exact value for the vanillas in this model to be calibrated.

Calibration of the hybrid local volatility model
The risk-neutral diffusion of the model is written as The two brownian motions are correlated with a constant correlation denoted by ρ. Classic regularity and ellipticity assumptions are made on the coefficients of the diffusion (described in [5]) to get the Proposition 3. The diffusion model defined above has a transition density with respect to Lebesgue's measure. The value of σ that fits its vanillas is given by where p(t, S, y) is the density of the couple (S t , y t ), and r(t) a deterministic curve of rates used in the computation of Dupire's local volatility σ D .
Proof. The existence of the density p(t, S, y) stems from the assumptions on the coefficients. Let us prove formula (20). The function p solves the forward parabolic equation with the initial condition p(0, S, y) = δ S0,y0 . As previously, the equation is integrated with respect to y, writing q(t, S) = R p(t, S, y)dy r(t, y)p(t, S, y)dy) + R r(t, y)p(t, S, y)dy = 0 This equation needs to be matched with Both of them can be written as where the second line stems from a simple integration by parts. Reintroducing this into the previous equations, we get In order to calibrate the vanillas, all that remains to be done is match the marginal density q with q D , giving the necessary condition

qds q
Replacing q by R p(t, S, y)dy completes the proof. The calibration equation for the vanillas of our hybrid model is thus with σ given by formula (20). Using similar technics to the other cases, an existence result can be obtained under certain assumptions, but this is not the scope of this paper. It is however noteworthy that one of the necessary hypothesis is the small variation of function r with respect to the deterministic curve r.

Numerical calibration
In this section, the theoretical results above are applied to calibrate a given diffusion model. Assuming the instantaneous rate to obey a Vasicek model (or in other words an Ornstein-Uhlenbeck process), the diffusion equations become The ADI algorithm described in section 5 is applied to equation (21) with the coefficients associated to this diffusion. The initial condition is p 0 (S, r) = δ S0,r0 . As in the two previous sections, this partial integro-differential equation is solved with a variable change for the spot process x = ln(S). The grid chosen is [−10σ √ t, 10σ To assess the quality of the calibration, call and put options on the spot process are computed by integration on the grid and compared to the targeted prices (target columns

Algorithm for the resolution of the calibration equation
In the previous sections, we wrote the calibration equations and gave graphs for their efficiency. In this one, we describe the algorithm used to solve them: a classic alternating direction implicit scheme. The nonlinear term is handled using a forward induction at first, and then a predictor-corrector method. The strong feature of an ADI scheme is its convergence rate in time and space: O(δx 2 ) + O(δt 2 ). The nonlinearity of the equation challenges this assertion. It is however possible to prove that it remains true in this case too. 13

Alternating Direction Implicit scheme
The calibration equation is a second-order parabolic equation. One of the most efficient method to solve such equations is a finite-difference approximation with alternating direction methods. For more informations on the subject, the reader can look into [4], numerous articles have also been published, in particular [11,12,13,14]. Let us now consider the following equation where f, g, α, β and γ are functions of t, x and y, ρ a constant and I(p) the quotient of integrals x, y)dy p(t, x, y)dy (24) Remark 5.1. We restricted ourselves to two-dimensional equations since they cover all the concrete examples studied previously. But the computations that follow are true in the general case. The cost in time however becomes an issue in higher dimensions.
The domain for the numerical resolution is ]0, T [×]x * , x * [×]y * , y * [. The first step is to take care of the initial condition. Instead of the Dirac, the initial condition p 0 is chosen as a gaussian distribution with very small variance. It obviously approximates our initial condition. It also verifies (on any bounded domain) the properties of regularity and strict positivity required in the theoretical study of the equation (though in the present section, we are only interested in its numerical resolution). Let p(t, x, y) = p 0 (x, y) if x = x * , x * or y = y * , y * be the boundary conditions.
The algorithm is based upon a predictor-corrector approach. Let ∆x, ∆y and ∆t be increments of the variables x, y and t, where ∆x = x * −x * I , ∆y = y * −y * J and ∆t = T N with I, J and N integers. The sets of points in the x,y,t-plane is given by x n = x * + i ∆x x * −x * , y j = y * + j ∆y y * −y * and t n = n ∆t T , for 0 ≤ i ≤ I, 0 ≤ j ≤ J and 0 ≤ n ≤ N . We construct four sequences p n , p * n , q n and q * n of space-dependent functions with n between 0 and N.
A classic alternating direction implicit scheme functions as follows: let us define the initial functions and then by induction ) + δ 2 xy (ρf n g n I(p n )p n ) + δ 2 y ( ) + δ 2 xy (ρf n g n I(p n )p n ) + δ 2 y ( where f n (i, j) designates f (n∆t, x i , y j ) (the same thing being true for the other coefficients of the equation).
δ is a difference operator for the space derivatives. For instance, with 1 ≤ i ≤ I − 1 and 1 ≤ j ≤ J − 1 4∆x∆y The equations (25) and (26) form two tridiagonal systems that can be solved very efficiently. A recursion formula can be computed on the functions q n defined above And eventually ) + δ 2 xy (ρf n g n I(p n )p n ) + δ 2 y ( In the litterature concerning alternating direction implicit schemes, the functions q n are often called the predicted value of the solution p, a second 'corrector' step usually follows. In our case, the equation being nonlinear, we start by chosing p n+1 = q n+1 and study the finite-difference approximation that results. Proof. Let p be a classic solution of equation (22) on a bounded domain, with p |t=0 = p 0 . We assume that p is strictly positive and sufficiently differentiable for all the quantities in the sequel to be properly defined. All the derivatives that appear are bounded as a consequence of the regularity assumptions on p. Writing p n i,j = p n (x i , y j ) = p(t n , x i , y j ), a simple Taylor expansion with remainder gives As for the space derivatives, clearly where the different constants θ are between 0 and ∆t, ∆x or ∆y depending on the context, they may change from one formula to another. Let E denote the truncation error for scheme (27), we have Applying the previous expansions and using the fact that p verifies equation (22) e 3 enables us to compensate for both the nondiagonal term and the nonlocal term I(p n ) that cannot be computed implicitely e 41 and e 42 are the terms corresponding to the first order space-derivatives At last, e 5 is the correction term stemming from the alternating direction implicit scheme Thanks to the regularity of p and of the coefficients of (22), the upper bounds e 1 ≤ K∆t 2 , e 21 +e 41 ≤ K∆x 2 , e 23 + e 42 ≤ K∆y 2 and e 22 ≤ K(∆y 2 + ∆x 3 ∆y ) are obtained. It is also easily proven that e 5 ≤ K∆t 2 . The term that prevents us from getting an error in O(∆t 2 ) is e 3 . All we have is e 3 ≤ K∆t. This concludes the proof.
Remark 5.2. In a case with no I(p) term, the equation is a classic linear and parabolic one. In that case, when the off-diagonal 1 term is absent (ρ = 0 for instance), the previous scheme has an error in O(∆t 2 ). To obtain such an error in the general case, a second 'corrector' step is generally used: the predicted value q n+1 is introduced as an approximation of p n+1 in the cross-derivatives. Here, we use it in the nonlocal term too.
The correction step is the following Here too, one can compute p * Eventually, a Crank-Nicholson like formula appears Let us study the consistency of this new scheme Proof. To prove the consistency, let q n+1 be defined as ) + δ 2 xy (ρf n g n I(p n )p n ) + δ 2 y ( The computations are almost identical to the previous proposition. This time the error is equal to − δ 2 xy (ρ (f g) n+1 I(q n+1 )q n+1 + (f g) n I(p n )p n 2 ) − δ 2 y ( g 2 n+1 p n+1 + g 2 n p n 4 ) + δ x ( α n+1 p n+1 + α n p n 2 ) + δ y ( β n+1 p n+1 + β n p n 2 ) + γ n+1 p n+1 + γ n p n 2 Using the same decomposition, the errors e * 1 , e * 23 , e * 41 and e * 42 do not change. e * 21 is slightly different but still in O(∆x 2 ). As for e * 22 and e * 5 , we now have which also verify e * 22 ≤ K(∆y 2 + ∆x 3 ∆y ) and e * 5 ≤ K∆t 2 . The real difference can be seen in The important feature of the predictor is that the difference q n+1 − p n+1 is O(∆t 2 ). This gives e * 3 ≤ K(∆t 2 ) and concludes the proof.

Time Convergence Rate of the modified ADI algorithm
In this brief section, we compare the convergence of the algorithm with the theoretical rates computed in the previous part. To do so, the calibrated value of 1-year at-the-money vanillas is computed for different number N of time steps. We then plot the error between this price and the targeted value against N. The next graph is obtained with the one-step predictor algorithm The error is clearly in O(∆t) as was proved in Proposition 4. We conduct the same experiment with this time both the predictor and the corrector steps. This time too, numerical experiments seem to agree with theory. The error appears to be in O(∆t 2 ). The predictor/corrector scheme serves its purpose.

Instabilities of the solutions: numerical explosion for 'oscillating' volatilities
In the previous sections, we described the algorithm used to solve the different calibration equations concerned by our work, and also mentioned the theoretical research performed on the subject. Though a partial existence result for equation (4) was found, it was still impossible to prove it in the general case: a strongly variable function b. This last section is devoted to the local and stochastic volatility model, from that point of view. The numerical resolution of the calibration is performed for strongly variable functions b.
At first, let us go back to the equation for the calibration of a local and stochastic volatility model. For the sake of simplicity, the interest rate is assumed to be zero. The equation is the following The existence of a solution to this equation (on a bounded domain, with regularized boundary conditions) was previously stated under certain assumptions on b, the essential one being: there exists a constant b * such that |b − b(y 0 )| ≤ b * for a given y 0 . Using the algorithm described in the previous section, we are now interested in the behavior of the numerical solution of (28) when the function b violates the assumption. The model chosen for this study is the mean-reverting volatility already expounded b(y) = exp(y) α(t, y) = γ β(t, y) = κ (δ − y) with γ, κ and δ three strictly positive constants. The values chosen in section 2 for the different parameters of the model and of the algorithm led us to a satisfactory calibration. Using them once again, we plot on the left (Figure 1 (a)) the density p(T, x, y) for T = 1 year and x = ln(S/S 0 ) close to 0. As expected, the solution p is perfectly smooth. Now, bouncing on the idea of strongly variable functions b, another density is plotted on the right (Figure 1 (b)) with this time a function b equal to b(y) = exp(10y). The solution is not smooth anymore. On the contrary, some kind of instability seems to occur.
Remark 6.1. To check that no other numerical effects are involved in the instability, the adjoint equation of (4) was also studied. From a theoretical viewpoint, it admits a solution without any restrictive assumption on b, see [15,5]. We make the same test and plot its numerical solution for b(y) = exp(Cy) with C = 10 and C = 15 On both these graphics (Figure 2), no sign of instability can be seen. With higher values of C, for instance C = 20, the function p computed numerically takes meaningless values (10 80 ), but at no time does it start to oscillate.

Conclusion
Using methods inspired from local and stochastic volatility models, we were able to write calibration equations for two other cases: the so-called local correlation model and a hybrid local volatility and stochastic rates model. Their numerical resolution, based upon an alternating direction implicit scheme, produces a satisfactory fit under certain assumptions (confirmed by the theoretical difficulties met when studying them). When those hypothesis are not verified, an instability occurs. Explaining it brought us to consider Hadamard stability and a certain class of integro-differential equations. Unfortunately, the criterion we found can not be applied in the case of the local and stochastic volatility model, it remains an open problem. As far as the ADI algorithm is concerned, we managed to adapt it to deal with the nonlinearity of our equation. Its consistency was also proved, with a convergence rate in time in O(∆t 2 ). A result confirmed numerically.