A splitting method for fully nonlinear degenerate parabolic PDEs

Motivated by applications in Asian option pricing, optimal commodity trading etc., we propose a splitting scheme for fully nonlinear degenerate parabolic PDEs. The splitting scheme generalizes the probabilistic scheme of Fahim, Touzi and Warin [13] to the degenerate case. General convergence as well as rate of convergence are obtained under reasonable conditions. In particular, it can be used for a class of Hamilton-Jacobi-Bellman equations, which characterize the value functions of stochastic control problems or stochastic differential games. We also provide a simulation- regression method to make the splitting scheme implementable. Finally, we give some numerical tests in an Asian option pricing problem and an optimal hydropower management problem.


Introduction
Numerical methods for parabolic partial differential equations (PDEs) are largely developed in the literature, on finite difference scheme, finites elements scheme, semi-Lagrangian scheme, Monte-Carlo method, etc. For nonlinear PDEs, and especially in high dimensional cases, the numerical resolution becomes a big challenge.
A typical kind of nonlinear parabolic PDEs is the Hamilton-Jacobi-Bellman (HJB) equation, which characterizes the solution of the optimal control problems. In this context, for finite difference method, one can only use the explicit scheme, since the implicit scheme needs to invert too many matrices. In the one dimensional case, the explicit finite difference scheme can be easily constructed and the monotonicity is guaranteed by the CFL condition. In high dimensional cases, Bonnans and Zidani [4] propose a numerical algorithm to construct a monotone scheme. Another numerical method for general HJB equations is the semi-Lagrangian scheme proposed in Debrabant and Jakobsen [12]. It can be easily constructed to be monotone, but they need next to use a finite * CMAP, École Polytechnique, Paris. E-mail: xiaolu.tan@polytechnique.edu difference grid as well as an interpolation method to make it implementable. It hence can be viewed as a finite difference scheme.
Generally speaking, finite difference and semi-Lagrangian schemes are easily implemented and perform quite well in low dimensional cases; and in high dimensional cases, the Monte-Carlo method is preferred. Recently, Fahim, Touzi and Warin [13] proposed a probabilistic method for nonlinear parabolic PDEs, which is closely related to the second order backward stochastic differential equation (2BSDE) developed in Cheridito et al. [9] and Soner et al. [18]. With simulations of a diffusion process, they propose the estimations of the value function and its derivatives by conditional expectations, by which they can approximate the nonlinear part of the PDE and then get a convergent scheme. However, their scheme can only be applied in the non-degenerate cases.
We want to generalize the probabilistic scheme of Fahim, Touzi and Warin [13] to the degenerate case, motivated by its applications in finance. For example, in Asian option pricing problems, we must consider the cumulative average stock prices A t ; for lookback options, we consider also the historical maximum and/or minimum stock prices M t , m t . They are all degenerate variables without a diffusion generator, and hence the pricing equation turns to be a degenerate parabolic equation. In some optimal commodity trading models(see e.g. [1], [7] and [8]), the storage amount of commodities is an important state variable, and the optimization problem induces a PDE which degenerates on storage amount variable. In life insurance, Dai et al. [11] proposed a financial pricing model for a Variable Annuities product Guaranteed Minimum Withdrawal Benefit (GMWB). In their model, the price of GMWB depends on two variables: the reference account and the guaranteed account, where the latter degenerates and the pricing equation is a degenerate parabolic PDE.
For these degenerate PDEs, the degenerate part is separable. Therefore, a natural solution is the splitting scheme. Our idea is to use the probabilistic scheme to treat the non-degenerate part, and use the semi-Lagrangian scheme to solve the degenerate part, and by combining the two methods, we get a splitting scheme. In particular, it generalizes the probabilistic scheme of Fahim, Touzi and Warin [13] to the degenerate case.
Another contribution of the paper is to propose a simulation-regression technique to make the semi-Lagrangian scheme implementable, in place of the classical finite difference method together with interpolation technique as used in Debrabant and Jakobsen [12], or Chen and Forsyth [8]. In the simulation-regression method, we can use global polynomials, or local hypercubes or local polynomials etc. as regression function basis. The global polynomial method means to approximate a function with some polynomials on the whole space, while the local basis method means to discretize first the space into local rectangles, and then to approximate the corresponding function with some polynomials on every local rectangle. As illustrated in Gobet, Lemor and Warin [14] and also in Bouchard and Warin [6], the local hypercubes and local polynomials basis method are very efficient in concrete cases. Moreover, they show that in practice, it is enough to choose a small number (about five or six) of discretization points in every dimension for the local basis method, while for finite difference method, one needs many more discretization points (more than 50 points in [8] for example) in every dimension. In particular, it permits to treat problems in high dimensions (up to 5 dimensions in [13] and up to 6 dimensions in [6]). In our context, we shall provide a four dimensional numerical example.
proposed in Barles and Souganidis [3] and Barles and Jakobsen [2]. In Section 3, we propose a simulation-regression technique to approximate the conditional expectations used in the splitting scheme, making the scheme implementable. We shall also discuss the choices of function basis used in the regression and then provide some convergence results for this implementable scheme. Finally, Section 4 provides some experimental examples.

The degenerate PDE and splitting scheme
In this section, we first introduce a nonlinear parabolic PDE which has a separable degenerate part. We next propose a splitting scheme, and for which we provide a local uniform convergence result of the splitting scheme when the PDE satisfies a comparison result for bounded viscosity solutions, as well as a rate of convergence when the nonlinear part of the PDE is a concave Hamiltonian.

A degenerate nonlinear PDE
x) T , we define a linear operator L X on the smooth functions ϕ : We say that L X is a linear operator associated to the diffusion process X = (X t ) 0≤t≤T defined by the stochastic differential equation: Given a nonlinear function we then get a nonlinear operator F (t, x, y, ϕ, D x ϕ, D 2 xx ϕ) on ϕ. We denote by F p and F Γ the derivative of function F w.r.t. p and Γ.
Finally, let us introduce the degenerate fully nonlinear parabolic PDE which will be considered throughout the paper: The PDE (2.2) is composed by three separable parts: the linear part L X , the nonlinear part F , and the first order degenerate part H.

A splitting scheme
As observed above, the three parts in PDE (2.2) are separable, we can then propose a splitting numerical scheme to solve it. The idea is to split (2.2) into the following two equations: Let us first give a time discrete grid (t n ) n=0,··· ,N with t n := nh, where h := T /N for N ∈ N. As in [13], we defineX t,x h by the Euler scheme of the diffusion process X in Let v h denote the numerical solution, then the probabilistic scheme of [13] for equation Remark 2.1. The scheme T h is well defined as soon as Det(σ(t, x)) = 0 for each (t, x) ∈ [0, T ) × R d . When ϕ is smooth, by integration by parts, one can verify that , y , i = 0, 1, 2.
For more details on this fact and of the probabilistic scheme T h of (2.6), we refer to Fahim et al. [13].
Finally, we are ready to introduce the splitting scheme S h • T h for the original PDE + v h t n+ 1 2 , x + f α,β (t n , x, y)h, y + g α,β (t n , x, y)h . (2.10) Clearly, when Det(σ(t, x)) = 0 for every (t, x) ∈ [0, T ) × R d , the scheme S h • T h is well defined and it gives a unique numerical solution v h .

The convergence results
We shall provide two convergence results for the splitting scheme S h • T h in (2.10), similar to Fahim et al. [13]. The first one is the local uniform convergence in the context of Barles and Souganidis [3], and the second is a rate of convergence.
Let us now give some assumptions on the equation (2.2), and then provide a first convergence result.
Assumption F : (i) The diffusion coefficients µ and σ are Lipschitz in x and continuous in t, σσ T (t, x) > 0 for all (t, x) ∈ [0, T ] × R d and T 0 σσ T (t, 0) + µ(t, 0) dt < ∞. (ii) The nonlinear operator F is uniformly Lipschitz in (x, y, r, p, Γ), continuous in t and sup (t,x,y)∈Q T |F (t, x, y, 0, 0, 0)| < ∞. (iii) F is elliptic and satisfies  It is clear that Assumptions F and H hold true for a class of HJB equations as well as a class of HJBI equations which characterize the value functions of the stochastic differential game problems. We next provide a rate of convergence in case that F and H are both concave Hamiltonians, i.e. when the nonlinear equation (2.2) is a HJB equation. We shall use the arguments developed by Barles and Jakobsen [2]. The following stronger assumptions implies that the nonlinear PDE (2.2) satisfies a comparison result for bounded functions, and has a unique bounded viscosity solution given a bounded and Lipschitz continuous function Φ, see e.g. Proposition 2.1 of [2]. Assumption HJB : Assumptions F and M hold and F is a concave Hamiltonian, i.e.
And B = {β} is a singleton, hence H is also a concave Hamiltonian, so that it can be written as Moreover, the functions l, c, f , g and σ satisfy that sup α∈A,γ∈C Assumption HJB+ : Assumption HJB holds true, and for any δ > 0, there exists a finite v is the unique bounded viscosity solution of (2.2) introduced in Theorem 2.6.
Remark 2.8. The above convergence rate is the same as that obtained in Fahim et al. [13]. It may not be the best rate in general. However, to the best of our knowledge, it is the optimal rate that we can prove in this stochastic control problem context so far.

Proof of local uniform convergence
To prove the local uniform convergence in Theorem 2.6, we shall verify the criteria proposed in Theorem 2.1 of Barles and Souganidis [3]: the monotonicity, the consistency of the scheme and the stability of the numerical solutions. Moreover, as discussed in Remark 3.2 of [13], we need also to show that lim inf Proof. By Lemma 3.12 and Remark 3.13 of [13] x, y). Then since c α,β ≥ 0 according to Assumption M, it follows immediately by We first define a consistency error function, then prove that our splitting scheme Definition 2.11. Given a smooth function ϕ defined on Q T , the consistency error function of scheme S h • T h is given by for every (t, x, y) ∈ Q T and every smooth function ϕ with bounded derivatives.
for a constant C independent of ϕ and h.

Proof.
For every (t, x, y) ∈ Q T , the value Λ ϕ h (t, x, y) is independent of the value of (µ(t,x), σ(t,x)) when (t,x) = (t, x). Hence we can always change the value of µ and σ outside the neighborhood of (t, x) without influence on the definition of consistency in (2.14). Therefore, without loss of generality, we can just suppose that µ and σ are uniformly bounded and show that for every smooth function ϕ with bounded derivatives of any order, the consistency error function Λ ϕ h defined in (2.13) satisfies x, y) and by E 2 (t, x, y, ϕ) the last two terms of the above equality (2.16) Clearly, by the boundedness of µ and σ, together with Assumption F, there is a constant C independent of h such that x, y, ϕ) is defined by the last equality of (2.17), and it satisfies Combining the estimations (2.16) and (2.17), and by (2.13) as well as the equality Proof. Suppose that |v h (t n+1 , ·)| 0 ≤ C n+1 , then from Lemma 3.14 of [13], there exists We have shown in the above the monotonicity, consistency and stability of scheme S h • T h , the rest is to confirm (2.12). In fact, we will provide a little stronger property of (v h ) h>0 which implies that Proposition 2.14. Let Assumptions F, H and M hold true, and Φ be Lipschitz and uniformly bounded.
Proof. To prove the that v h is Lipschitz in (x, y), we shall use the discrete Gronwall inequality as in the proof of Lemma 3.16 of [13].
We can also prove that v h is 1/2−Hölder in t as was done in Lemma 3.17 of [13] for their numerical solution. However, to avoid the heavy calculation in their proof, we shall give a weaker result which is enough to guarantee the condition (2.12).
Proof. We first introducev h as the numerical solution of (2.4) computed by scheme ·). Clearly, by Lemmas 3.14 and 3.17 of [13], (v h ) h>0 is uniformly bounded and satisfies It follows by the monotonicity of T h and the uniform boundedness of v h andv h that |v h (t n , x, y) − v h (t n+ 1 2 , x, y)| ≤ L(T − t n+1 ) + Ch.

Proof for rate of convergence
As in [13], our arguments to prove the rate of convergence in Theorem 2.7 are based on Krylov's shaking coefficient method, and our analysis stays in the context of Barles and Jakobsen [2]. We first derive some technical Lemmas similar to that in [13].
Proof. First, from the proof of Lemma 3.21 in [13], we have It follows by the definition of the splitting scheme S h • T h in (2.10) that By the monotonicity of the splitting scheme S h • T h , we get Finally, using exactly the same arguments as in the proof of Lemma 3.5 of [13], it follows that which concludes the proof.
for some bounded functions g 1 and g 2 . Then for every n = 0, · · · , N, Proof. With Lemma 2.16, the proof is exactly the same as in Proposition 3.20 of [13].
Note that we replace β by λ in our proposition.

Basis projection and simulation-regression method
To get an implementable scheme, we need to specify how to compute the expecta- in the splitting scheme S h • T h . When analytic closed formulas are not available in the concrete examples, we usually use Monte-Carlo simulation-regression method to estimate them. Some estimations were discussed in recent works, e.g. Malliavin estimations [5], function basis regression [14] and cubature method [10], etc.
All of these methods need the simulations of X. Given a discrete time grid (t n ) 0≤n≤N , where t n := n h and h := T /N , we define a Euler approximationX of X X tn+1 :=X tn + µ(t n ,X tn )h + σ(t n ,X tn )∆W n+1 , However, these methods are usually discussed in a non-degenerate context, in other words, they can be used for a given fixed y, which is not appropriate for the implementation of our splitting scheme S h • T h .
One solution is to discretize the space of Y into a discrete grid (y i ) i∈I , and then for each fixed y i , we simulate the diffusion process X and get estimations of the conditional expectations for all x with every fixed y i , then use the interpolation method to get the estimation of theses expectations for all x and y. This is a combination of finite difference method and Monte-Carlo method, which may lose the advantages of Monte-Carlo method in high dimensional cases. Therefore, we propose to simulate the diffusion process X with Euler scheme and to simulate Y with a continuous probability distribution (e.g. normal distribution, uniform distribution, etc.) independent of X. And then we use a regression method like in Longstaff and Schwartz [15] in American option pricing context or Gobet, Lemor and Warin [14] in BSDE context to estimate the conditional expectations with which we shall make the splitting scheme S h • T h implementable. (ii) In practice, if we choose local hypercubes or local polynomials as functions basis for the regression method, we still need to discretize the space. However, as discussed in the introduction, this discretization can be coarse in practice, which permits to keep the advantage of the simulation-regression method in high-dimensional cases (see also the numerical examples in Section 4).
In the following, we first give a basis projection scheme as well as a similationregression method to estimate the regression coefficient. Then we discuss the convergence of Monte-Carlo errors in our context.

The basis projection scheme
To compute the conditional expectations (3.2), we first project them on a functional space spanned by the basis functions (e k (x, y)) 1≤k≤K , where K ∈ N ∪ {+∞}. We recall that H t,x,h 2 is a matrix of dimension d×d, H t,x,h 1 is a vector of dimension d and H t,x,h 0 = 1. In order to simplify the presentation, we shall suppose that d = d = 1. All of the results can be easily extended to the case d > 1 and/or d > 1. Let then the projected approximation of (3.2) is denoted bỹ

Remark 3.2.
There are several choices for function basis (e k (x, y)) 1≤k≤K , for example global polynomials, local hypercubes or local polynomials, we refer to Bouchard and Warin [6] for some interesting discussions.
We replace the conditional expectations (3.2) in scheme S h • T h by their projected approximations (3.4), and denote the new splitting scheme by S h •T h . Concretely, it is defined as follows:

Simulation-regression method
Next, we propose to use a simulation-regression method to approximateλ. We still suppose that d = d = 1 for simplicity. Let (X m tn ) 0≤n≤N , (∆W m n ) 0<n≤N , Y m 1≤m≤M be M independent simulations ofX, ∆W and Y , whereX is defined in (3.1), the regression method with function basis (e k (x, y)) 1≤k≤K is to get the solution of the least square problem: A raw regression estimation of the conditional expectations (3.2) from these M samples is given bȳ Then with a priori upper bounds Γ i (X tn , Y ) and lower bounds Γ i (X tn , Y ), we define the regression estimation of (3.2):
(3.9) Remark 3.4. In Gobet et al. [14], the authors propose the following minimization problem in place of (3.6): which gives also a good estimation forλ i by the fact that ∆W n+1 is independent of the σ−field generated by Y, W 0 , ∆W 1 , · · · , ∆W n .
We replace the conditional expectations (3.2) in scheme S h • T h by their regression estimations (3.8) x, y)h, y + g α,β (t n , x, y)h .

The convergence results of simulation-regression scheme
To get a convergence result of schemes S h •T h and S h •T M h , we can no longer use the same arguments as in Fahim et al. [13], since there is no uniform convergence property in L p for the Monte-Carlo error (Ê M − E)(R) as in the Assumption E of [13]. To see this, let us consider the extreme case where the equation is totally degenerate (i.e. d = 0 and d > 0), and then we need to approximate an arbitrary bounded function in a functional space with finite number of basis functions, which does not give a uniform convergence. Also, since we are in the viscosity solution analysis context of Barles and Souganidis [3], we can not hope to obtain a probabilistic L 2 (Ω)−convergence as in Gobet et al. [14]. , where (B k ) 1≤k≤K is a partition of D ⊆ R d+d , then the projection approximation is equivalent to taking another conditional expectation on the σ-field generated by (X tn , Y ) ∈ B k 1≤k≤K , in other words, Then for every (x, y) ∈ B k , .
The proof is almost the same as that of Theorem 5.1 of Bouchard and Touzi [5], we report it in Appendix for completeness. Let ϕ be bounded by constant b, δ denote the longest edge of the hypercubes (B k ) 1≤k≤K , then the volume of B k is of order δ d+d , and P((X tn , Y ) ∈ B k ) ≈ Cδ d+d , where C depends on the density of (X tn , Y ). As the total volume of D is fixed and finite, let

Some analysis on the basis projection scheme S h •T h
In preparation of the proof for Theorems 3.9 and 3.10, we give some analysis on the scheme S h •T h . In general, we shall show that if we use the local hypercubes as projection function basis, then S h •T h is still monotone, consistent and stable.

The proof for convergence results of scheme S h •T M h
To prove Theorem 3.9, we shall mimic the proof of Theorem 4.1 in [13], which uses the arguments of [3] in a stochastic context.

Proof of Theorem 3.9. Givenv h the numerical solution of scheme
First, it is clear by the truncation function (3.9) as well as the boundedness of F that |v(t n , x, y) − Φ(x, y)| ≤ C(T − t n ) for some constant C, which implies thatv * (T, x, y) = v * (T, x, y) = Φ(x, y). Then it is enough to prove thatv * andv * are respectively viscosity supersolution and subsolution of (2.2) to conclude the proof with the comparison assumption. Here we shall only prove the supersolution property, since the subsolution property holds true with the same kind of argument.
Given (t 0 , x 0 , y 0 ) ∈ Q T and a test function ϕ ∈ C ∞ c (Q T ) such that 0 = min(v * − ϕ) = (v * − ϕ)(t 0 , x 0 , y 0 ), by uniform boundedness ofv h and manipulation on ϕ, there is a sequence (t k , x k , y k , h k ) → (t 0 , x 0 , y 0 , 0) such thatv h k (t k , x k , y k ) →v * (t 0 , x 0 , y 0 ) and From the monotonicity of scheme S h •T h , it follows that . We claim that R k → 0 P-a.s. along some subsequence. (3.16) Then, from the consistence of scheme S h •T h in Proposition 3.12, which is the required supersolution property. Therefore, it is enough to justify the claim (3.16) to conclude the proof. Indeed, by the definition of splitting scheme S h •T M h and S h •T h , and the boundedness of c α,β , (3.12). And therefore by Lemma 3.8 which turns to 0 by assumptions of the theorem. Further, the L 2 −convergence implies that R k → 0 in probability and hence it admits a subsequence which converges to 0 almost surely. We then proved the claim (3.16) and hence conclude the proof of the theorem.
Similarly, we have the other side of the error boundary and get (3.17) Finally, it is enough to take expectations on both sides of (3.17) and maximize the right side on for h = h 3 10 , which implies that

Numerical examples
We provide here some numerical examples, one is from Asian option pricing problem and the other is from an optimal management problem for a hydropower plant. In every numerical example, we use local polynomial basis for the simulation regression method. The space is discretized to get the local basis, where we use 5 discretization points for every dimension, i.e. the space is divided into 6 parts along every dimension. The polynomials are of second order degree, i.e. they are of the form a 0 + a 1 x + a 2 x 2 in the one-dimensional case. Further, we give also the computation time of each numerical example using a computer with 2.4GHz CPU and 4G memory.

Asian option pricing
Our first example is to price Asian option in an uncertain volatility model (UVM), whose pricing equation is a degenerate and nonlinear PDE. Then we also consider the problem in UVM with Hull-White interest rate.

Asian option pricing in UVM: a two-dimensional case
We consider an uncertain volatility model with risky asset S t given by where r is the constant interest rate, σ is the volatility process which is bounded between the lower volatility σ and the upper volatility σ. Denote an Asian option is an option with payoff g(S T , A T ) at maturity T , whose pricing equation with terminal condition v(T, s, a) = g(s, a).
To implement the splitting scheme, we rewrite (4.1) in form of the equation (2.2) with some constant σ 0 : Further we consider a call spread type payoff g(S, A) = (A − K 1 ) + − (A − K 2 ) + . With 50 independent computations for every time discretization using the splitting scheme, we get the mean value as well as its standard deviation. Moreover, as comparison, we implemented the Crank-Nicolson finite difference scheme of equation (4.1) with parameters ∆S = 1 and ∆A = 0.25. The results of our splitting scheme and Crank-Nicolson scheme for different ∆t are given in Figure 1. We notice that the standard deviation of the splitting method price from 50 independent computations is less than 1% of the mean value and the relative difference between the two schemes are less than 0.3%.

Asian option pricing in UVM with Hull-White interest rate: a threedimensional case
We can also consider the uncertain volatility model with a stochastic interest rate, e.g. Hull-White interest rate (HW-IR). In HW-IR model, the interest rate has dynamic where θ t is determined by the current interst rate curve, b is the drawback force coefficient and B = (B t ) t≥0 is another Brownian motion with correlation ρ to Brownian motion W which generates the dynamics of risky asset S. Then the value function v(t, s, a, r) := E e − T t rsds g(S T , A T ) S t = s, A t = a, r t = r solves the pricing equation ∂ t v + rsD s v + b(θ t − r)D r v + 1 2 (σ r ) 2 D 2 rr v − rv + max σ≤σ≤σ ρσσ r sD 2 rs v + 1 2 σ 2 s 2 D 2 ss v + sD a v (t, s, a, r) = 0, with terminal condition v(T, s, a, r) = g(s, a). Let S 0 = 100, K 1 = 90, K 2 = 110, T = 1, σ = 0.15, σ = 0.25, r 0 = 0.02, b = 0.01, σ r = 0.01, ρ = 0.2 and interest rate curve is f t = 0.02, ∀t > 0. As in (4.2), we rewrite the pricing equation in form of (2.2) with constant σ 0 . For g(S T , A T ) = (A T − K 1 ) + − (A T − K 2 ) + , we implement our splitting method with different constants σ 0 , and take the mean value of 50 independent computations. The results are given in figure 2. We notice that the solution seems to be close to 11.51, and when the time discretization ∆t is large, the numerical solution underestimates the value. Another phenomena is that when σ 0 is larger (e.g. σ 0 = 0.25), the performance of the numerical solutions seems more stable.

Optimal management of hydropower plant: A four-dimensional case
Let us consider an optimal management problem for a hydropower plant, which generalizes a little the model in Chapiter 2 of the thesis of Arnaud Porchet [16].
A hydropower plant manages a dam, which is filled by rain precipitations with nonnegative rate A t , which follows equation Denote by B t the volume of water in the dam, then where q t represents the water flow sent at time t to generate electricity. It makes a profit T 0 q t S t dt in period [0, T ], where S t represents the market electricity price, which follows dynamics dS t = µ s S t dt + σ s S t dW 2 t .
At the same time, the power station invests in electricity market with money θ t , then the total revenue of the power station X t follows equation dX t = θ t S t dS t + q t S t dt = θ t µ s dt + θ t σ s dW 2 t + q t S t dt.
The power station optimizes its expected utility EU (X T ) on the strategy (q t ) 0≤t≤T and (θ t ) 0≤t≤T . Formally, we get a Bellman equation u t + µ s sD s u + 1 2 σ 2 s s 2 D 2 ss u + µ a aD a u + 1 2 σ 2 a a 2 D 2 aa u + ρσ s σ a saD 2 sa u + max θ θ s µD x u + 1 2 θ 2 σ 2 s D 2 xx u + θρaσ a σ s D 2 ax u + θσ 2 s sD 2 sx u + max q (a − q)D b u + qsD x u = 0. u t + µ s sD s u + 1 2 σ 2 s s 2 D 2 ss u + µ a aD a u + 1 2 σ 2 a a 2 D 2 aa u + ρσ s σ a saD 2 sa u + Let µ a = 0, σ a = 0.2, µ s = 0, σ s = 0.2, ρ = 0, n = 5 and the utility function is given by U (x) := −e −ρx with ρ = 0.2. Using the different choices of σ x , we report the numerical result in Figure 3. We notice that the solution seems converge to the value −0.66, and when σ x is chosen larger (e.g. σ x = 1.2), the numerical solution is more stable w.r.t. the time step ∆t.

Appendix
We give here the proof for Lemma 3.8. Let (λ i k ) 1≤k≤K be the projection coefficients of R i (ϕ) on basis (e k (X tn , Y )) 1≤k≤K as defined in