Xiaolu Tan
email: xiaolu.tan@polytechnique.edu

A splitting method for fully nonlinear degenerate parabolic PDEs

Keywords: Numerical scheme, nonlinear degenerate PDE, splitting method, viscosity solution AMS MSC 2010: Primary 65C05, secondary 49L25

 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 simulationregression 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 [START_REF] Bonnans | A fast algorithm for the two dimensional HJB equation of stochastic control[END_REF] 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 [START_REF] Debrabant | Semi-Lagrangian schemes for linear and fully nonlinear diffusion equations[END_REF]. It can be easily constructed to be monotone, but they need next to use a finite 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 [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF] 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. [START_REF] Cheridito | Second-order backward stochastic differential equations and fully nonlinear parabolic PDEs[END_REF] and Soner et al. [START_REF] Soner | Wellposedness of second order backward SDEs, Probability Theory and Related Fields[END_REF]. 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 [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF] 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. [START_REF] Bardou | When are swing options bang-bang?[END_REF], [START_REF] Bonnans | Energy contracts management by stochastic programming techniques[END_REF] and [START_REF] Chen | A Semi-Lagrangian Approach for Natural Gas Storage Valuation and Optimal Operation[END_REF]), 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. [START_REF] Dai | Guaranteed minimum withdrawal benefit in variable annuities[END_REF] 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 [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF] 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 [START_REF] Debrabant | Semi-Lagrangian schemes for linear and fully nonlinear diffusion equations[END_REF], or Chen and Forsyth [START_REF] Chen | A Semi-Lagrangian Approach for Natural Gas Storage Valuation and Optimal Operation[END_REF]. 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 [START_REF] Gobet | A regression-based Monte Carlo method to solve backward stochastic differential equations[END_REF] and also in Bouchard and Warin [START_REF] Bouchard | Monte-Carlo valorisation of American options: facts and new algorithms to improve existing methods[END_REF], 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 [START_REF] Chen | A Semi-Lagrangian Approach for Natural Gas Storage Valuation and Optimal Operation[END_REF] for example) in every dimension. In particular, it permits to treat problems in high dimensions (up to 5 dimensions in [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF] and up to 6 dimensions in [START_REF] Bouchard | Monte-Carlo valorisation of American options: facts and new algorithms to improve existing methods[END_REF]). In our context, we shall provide a four dimensional numerical example.

The rest of the paper is organized as follows. In Section 2, we introduce a degenerate PDE and a splitting scheme which combines the probabilistic scheme in [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF] and semi-Lagrangian scheme. Then we provide a local uniform convergence result as well as a rate of convergence, where the main idea is to adapt the viscosity solution technique proposed in Barles and Souganidis [START_REF] Barles | Convergence of approximation schemes for fully nonlinear second order equations[END_REF] and Barles and Jakobsen [START_REF] Barles | Error bounds for monotone approximation schemes for parabolic Hamilton-Jacobi-Bellman equations[END_REF]. 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.

Notation:

Let |η| := η 1 + • • • + η d for η ∈ N d . Given T ∈ R + and d, d ∈ N, we denote Q T := [0, T ) × R d × R d , Q T := [0, T ] × R d × R d and C 0,1 (Q T ) := ϕ : Q T → R such that |ϕ| 1 < ∞ ,
where

|ϕ| 0 := sup Q T |ϕ(t, x, y)| and |ϕ| 1 := |ϕ| 0 + sup Q T ×Q T |ϕ(t, x, y) -ϕ(t , x , y )| |x -x | + |y -y | + |t -t | 1 2
.

In this paper, the constant C is used in many inequalities, its value may vary from line to line.

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

Let T ∈ R + , µ : [0, T ] × R d → R d and σ : [0, T ] × R d → S d
be continuous, denote a(t, x) := σ(t, x)σ(t, x) T , we define a linear operator L X on the smooth functions ϕ :

Q T → R by L X ϕ(t, x, y) := ∂ t ϕ(t, x, y) + µ(t, x) • D x ϕ(t, x, y) + 1 2 a(t, x) • D 2 xx ϕ(t, x, y).
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:

dX t = µ(t, X t ) dt + σ(t, X t ) dW t , (2.1)
where W = (W t ) 0≤t≤T is a d-dimensional standard Brownian motion.

Given a nonlinear function

F : (t, x, y, r, p, Γ) ∈ R + × R d × R d × R × R d × S d → F (t, x, y, r, p, Γ) ∈ R,
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 Γ.

Next, we give the degenerate part which involves with the partial gradient with respect to y. Given functions

l α,β , c α,β , f α,β i , g α,β j α∈A, β∈B, 1≤i≤d, 1≤j≤d
defined on Q T with index space A and B, we denote f α,β := (f α,β i ) 1≤i≤d and g α,β := (g α,β j ) 1≤j≤d , and define the Lagrangian L α,β by L α,β ϕ(t, x, y) := l α,β (t, x, y) + c α,β (t, x, y)ϕ(t, x, y) Finally, let us introduce the degenerate fully nonlinear parabolic PDE which will be considered throughout the paper:

+ f α,β (t, x, y) • D x ϕ(t, x, y) + g α,β (t, x, y) • D y ϕ(t,
-L X v -F (•, v, D x v, D 2 xx v) -H(•, v, D x v, D y v) (t, x, y) = 0, on Q T , (2.2)
with terminal condition v(T, x, y) = Φ(x, y).

(2.

3)

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:

-L X v(t, x, y) -F (•, v, D x v, D 2 xx v)(t, x, y) = 0 (2.4)
and

-∂ t v(t, x, y) -H(•, v, D x v, D y v)(t, x, y) = 0, (2.5)
then to solve them separately. Equation (2.4) is nonlinear and non-degenerate for every fixed y, then it can be treated by the probabilistic scheme proposed in Fahim et al. [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF].

Equation (2.5) is a first order Hamilton-Jacobi-Bellman-Isaacs (HJBI) equation, we shall solve it by semi-Lagrangian scheme. Then, combining the two schemes sequentially, we get the splitting scheme.

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 [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF], we define Xt,x h by the Euler scheme of the diffusion process X in (2.1):

Xt,x

h := x + µ(t, x) h + σ(t, x) • (W t+h -W t ), ∀(t, x) ∈ [0, T ] × R d .
Let v h denote the numerical solution, then the probabilistic scheme of [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF] for equation (2.4) is given by

v h (t n , x, y) = T h [v h ](t n , x, y) := E v h (t n+1 , Xtn,x h , y) + hF (t n , x, y, ED h v h (t n , x, y)), (2.6)
where

ED h ϕ(t n , x, y) := E ϕ(t n+1 , Xtn,x h , y)H tn,x,h i (∆W n+1 ) : i = 0, 1, 2 ,
with ∆W n+1 := W tn+1 -W tn and the Hermite polynomials are defined by H t,x,h

0 (w) := 1, H t,x,h 1 (w) := σ T (t, x) -1 w h and H t,x,h 2 (w) := σ T (t, x) -1 ww T -hI d h 2 σ(t, x) -1 .
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

E ϕ t n+1 , Xtn,x h , y H tn,x,h i (∆W n+1 ) = E D i x i ϕ t n+1 , Xtn,x h , 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. [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF].

EJP 18 (2013), paper 15.

The second PDE (2.5) is a first order HJBI equation, its semi-Lagrangian scheme is given by v h (t n , x, y) = S h [v h ](t n , x, y) := inf α∈A sup β∈B hl α,β (t n , x, y) + hc α,β (t n , x, y)v h (t n+1 , x, y) + v h t n+1 , x + hf α,β (t n , x, y), y + hg α,β (t n , x, y) .

(2.7) Remark 2.2. The semi-Lagrangian scheme S h is deduced intuitively from the discrete version of equation (2.5):

v h (t n+1 , x, y) -v h (t n , x, y) h + inf α∈A sup β∈B l α,β (t n , x, y) + c α,β (t n , x, y)v h (t n+1 , x, y) + v h (t n+1 , x + hf α,β (t n , x, y), y + hg α,β (t n , x, y)) -v h (t n+1 , x, y) h = 0.
Finally, we are ready to introduce the splitting scheme S h

v h (t n+ 1 2 , x, y) := T h [v h ](t n , x, y) with T h defined in (2.6), (2.9)
and

v h (t n , x, y) = S h • T h [v](t n , x, y) := inf α∈A sup β∈B h l α,β (t n , x, y) + h c α,β (t n , x, y) v h (t n+ 1 2 , x, y) + 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. [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF]. The first one is the local uniform convergence in the context of Barles and Souganidis [START_REF] Barles | Convergence of approximation schemes for fully nonlinear second order equations[END_REF], and the second is a rate of convergence.

We first recall that an upper semicontinuous (resp., lower semicontinuous) function v (resp. v) on Q T is called a viscosity subsolution (resp., supersolution) of (2.2) if, for any (t, x, y) ∈ Q T and any smooth function ϕ satisfying

0 = (v -ϕ)(t, x, y) = max Q T (v -ϕ) resp., 0 = (v -ϕ)(t, x, y) = min Q T (v -ϕ) , we have -L X ϕ -F (t, x, y, ϕ, D x ϕ, D
v on Q T satisfying v(T, •) ≤ v(T, •), we have v ≤ v.
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

a -1 • F Γ ≤ 1 on R × R d × R d × R × R d × S d . (2.11) (iv) F p ∈ Image(F Γ ) and F T p F -1 Γ F p ∞ < +∞.
|l α,β | 0 + |c α,β | 0 + |f α,β i | 0 + |g α,β j | 0 < ∞. Assumption M : F r -1 4 F T p F -1 Γ F p ≥ 0 and c α,β ≥ 0 for every α ∈ A, β ∈ B.
Remark 2.5. Assumption M is imposed to guarantee the monotonicity of the splitting scheme S h • T h . However, it is not crucial as soon as Assumptions F and H hold true.

In fact, as discussed in Remark 3.13 of [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF] 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 [START_REF] Barles | Error bounds for monotone approximation schemes for parabolic Hamilton-Jacobi-Bellman equations[END_REF]. 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 [START_REF] Barles | Error bounds for monotone approximation schemes for parabolic Hamilton-Jacobi-Bellman equations[END_REF]. Assumption HJB : Assumptions F and M hold and F is a concave Hamiltonian, i.e.

µ • p + 1 2 a • Γ + F (t, x, y, r, p, Γ) = inf γ∈C L γ (t, x, y, r, p, Γ), with L γ (t, x, y, r, p, Γ) := l γ (t, x, y) + c γ (t, x, y)r + f γ (t, x, y) • p + 1 2 a γ (t, x, y) • Γ.
A splitting method for nonlinear degenerate PDEs And B = {β} is a singleton, hence H is also a concave Hamiltonian, so that it can be written as

H(t, x, y, r, p, q) = inf α∈A l α (t, x, y) + c α (t, x, y)r + f α (t, x, y) • p + g α (t, x, y) • q
Moreover, the functions l, c, f , g and σ satisfy that sup α∈A,γ∈C

|l α + l γ | 1 + |c α + c γ | 1 + |f α + f γ | 1 + |g α | 1 + |σ γ | 1 < ∞
Assumption HJB+ : Assumption HJB holds true, and for any δ > 0, there exists a finite set {α i , γ i } I δ i=1 such that for any (α, γ)

∈ A × C : inf 1≤i≤I δ |l α -l αi | 0 + |c α -c αi | 0 + |f α -f αi | 0 + |σ α -σ αi | 0 ≤ δ, and inf 1≤i≤I δ |l γ -l γi | 0 + |c γ -c γi | 0 + |f γ -f γi | 0 + |g γ -g γi | 0 ≤ δ.
Theorem 2.7. Suppose that the terminal condition function Φ is bounded and Lipschitzcontinuous. Then there is a constant C such that (i) under Assumption HJB, we have v -v h ≤ Ch 1 4 , (ii) under Assumption HJB+, we have -Ch

1 10 ≤ v -v h ≤ Ch 1 4
, where 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. [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF]. 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 [START_REF] Barles | Convergence of approximation schemes for fully nonlinear second order equations[END_REF]: the monotonicity, the consistency of the scheme and the stability of the numerical solutions. Moreover, as discussed in Remark 3.2 of [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF], we need also to show that lim inf

(t ,x ,y ,h)→(T,x,y,0) v h (t , x , y ) ≥ Φ(x, y) and lim sup (t ,x ,y ,h)→(T,x,y,0) v h (t , x , y ) ≤ Φ(x, y).
(2.12) Remark 2.9. By the definition of the numerical scheme S h •T h in (2.10), the numerical solution v h is only defined on the time grid (t n ) 0≤n≤n product R d × R d . However, we can use linear interpolation method to extend v h on the whole space Q T .

Proposition 2.10. Let Assumptions F, H and M hold true, then for two functions ϕ and ψ defined on Q T with exponential growth, we have

ϕ ≤ ψ =⇒ S h • T h [ϕ] (t, x, y) ≤ S h • T h [ψ] (t, x, y).
Proof. By Lemma 3.12 and Remark 3.13 of [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF],

ϕ ≤ ψ implies that T h [ϕ](t, x, y) ≤ T h [ψ](t, x, y). Then since c α,β ≥ 0 according to Assumption M, it follows immediately by (2.10) that S h • T h [ϕ](t, x, y) ≤ S h • T h [ψ](t, x, y).
We first define a consistency error function, then prove that our splitting scheme

S h • T h is consistent.
EJP 18 (2013), paper 15.

Definition 2.11. Given a smooth function ϕ defined on Q T , the consistency error function of scheme S h • T h is given by

Λ ϕ h (•) := ϕ(•) -S h • T h [ϕ](•) h + L X ϕ(•) + F (•, ϕ, D x ϕ, D 2 xx ϕ) + H(•, ϕ, D x ϕ, D y ϕ).
(2.13)

The scheme S h • T h is said consistent if Λ ϕ+c h (t , x , y ) → 0 as (c, h, t , x , y ) → (0, 0, t, x, y), (2.14)
for every (t, x, y) ∈ Q T and every smooth function ϕ with bounded derivatives.

Proposition 2.12. Let Assumptions F, H and M hold true, then the scheme S h • T h is consistent. In addition, if µ and σ are uniformly bounded, then the consistency error function Λ ϕ h is uniformly bounded by h E(ϕ), where

E(ϕ) := C 1 + |∂ tt ϕ| 0 + 2 i=0 |∂ t D i z i ϕ| 0 + 4 i=0 |D i z i ϕ| 0 with z := (x, y) ∈ R d+d ,
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

Λ ϕ h (•) 0 ≤ h E(ϕ).
(2.15)

First, let us denote

L Xt,x ϕ(t , x , y) := ∂ t ϕ(t , x , y) + µ(t, x) • D x ϕ(t , x , y) + 1 2 a(t, x) • D 2 xx ϕ(t , x , y),
then by Itô's formula,

E h (t, x, y, ϕ) := T h [ϕ](t, x, y) -ϕ(t, x, y) = h L X ϕ(•) + F (•, ϕ, D x ϕ, D 2 xx ϕ) (t, x, y) + h 2 1 h 2 E t+h t u t L Xt,x L Xt,x ϕ(s, Xt,x s , y) ds du (2.16) + h 2 1 h F (•, ED h ϕ)(t, x, y) -F (•, ϕ, Dϕ, D 2 xx ϕ)(t, x, y) . Denote E 1 (t, x, y, ϕ) := L X ϕ(t, x, y) + F (•, ϕ, D x ϕ, D 2 xx ϕ)(t, x, y)
and by E 2 (t, x, y, ϕ) the last two terms of the above equality (2.16) divided by h 2 , then E h (t, x, y, ϕ) can rewritten as

E h (t, x, y, ϕ) = h E 1 (t, x, y, ϕ) + h 2 E 2 (t, x, y, ϕ).
Clearly, by the boundedness of µ and σ, together with Assumption F, there is a constant

C independent of h such that E 2 (•, ϕ) 0 ≤ C 1 + |∂ tt ϕ| 0 + 2 i=0 |∂ t D i x i ϕ| 0 + 4 i=0 |D i x i ϕ| 0 ,
A splitting method for nonlinear degenerate PDEs and moreover, E 1 is Lipschitz in z := (x, y) with coefficient

L E1 ≤ C 1 + |∂ t D z ϕ| 0 + |D z ϕ| 0 + |D 2 zz ϕ| 0 + |D 3 zzz ϕ| 0 .
By simplifying c α,β (t, x, y), l α,β (t, x, y), f α,β (t, x, y), g α,β (t, x, y) into (c α,β , l α,β , f α,β , g α,β ),

we deduce that where E 3 (t, x, y, ϕ) is defined by the last equality of (2.17), and it satisfies

1 h S h [(ϕ + E h (•, ϕ))](t, x, y) -ϕ(t, x, y) -E h (t, x, y, ϕ) = 1 h inf α∈A sup β∈B hl α,β + hc α,β ϕ(t, x, y) + ϕ(t, x + f α,β h, y + g α,β h) -ϕ(t, x, y) + hc α,β E h (t, x, y, ϕ) + E h (t, x + f α,β h, y + g α,β h) -E h (t, x, y, ϕ) = inf α∈A sup β∈B l α,β + c α,β ϕ(t, x, y) + (f α,β • D x ϕ + g α,β • D y ϕ)(t, x, y) + 1 h ϕ(t, x + f α,β h, y + g α,β h) -ϕ(t, x, y) -(f α,β D x ϕ + g α,β D y ϕ)(t, x, y) + c α,β E h (t, x, y) + 1 h E h (t, x + f α,β h, y + g α,β h, ϕ) -E h (t, x, y, ϕ) =: H(•, ϕ, D x ϕ, D y ϕ)(t, x, y) + hE 3 (t, x, y, ϕ), (2.17
|E 3 (t, x, y, ϕ)| ≤ C |D 2 zz ϕ| 0 + 1 h E h (t, x, y, ϕ) + 2|E 2 (t, x, y, ϕ)| + L E1 ≤ E(ϕ).
Combining the estimations (2.16) and (2.17), and by (2.13) as well as the equality

ϕ(t, x, y) -S h • T h [ϕ](t, x, y) h = ϕ(t, x, y) -T h [ϕ](t, x, y) h + ϕ(t, x, y) + E h (t, x, y, ϕ) -S h [ϕ + E h (•, ϕ)](t,
Φ be L ∞ -bounded, then (v h ) h is L ∞ -bounded, uniformly in h for h small enough.
Proof. Suppose that |v h (t n+1 , •)| 0 ≤ C n+1 , then from Lemma 3.14 of [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF], there exists a constant C independent of h such that

v h t n+ 1 2 , • 0 ≤ C n+1 (1 + hC) + hC. It follows from (2.10) that when h < C -1 , |v h (t n , •)| 0 ≤ (1 + hC)(C n+1 (1 + hC) + hC) + hC ≤ (1 + 3hC)C n+1 + 3hC. Therefore, |v h (t n , •)| 0 ≤ C e C T
for some constant C (independent of h) from the discrete Gronwall inequality. 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 lim (t ,x ,y ,h)→(T,x,y,0) v h (t , x , y ) = Φ(x, y).

Proposition 2.14. Let Assumptions F, H and M hold true, and Φ be Lipschitz and uniformly bounded. Then (v h ) h is Lipschitz in (x, y), uniformly in h.

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 [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF].

Suppose that v h (t n+1 , •) is Lipschitz with coefficient L n+1 , then by the proof of Lemma 3.16 of [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF], the function v h (t

n+ 1 2 , •) = T h [v h ](t n , •) is Lipschitz in x with coeffi- cient L n+1 ((1 + Ch) 1/2 + Ch) + Ch; moreover, v h (t n+ 1 2 , •) is Lipschitz in y with coefficient L n+1 (1 + Ch) by Lemma 3.14 of [13]. It follows that v h (t n+ 1 2 , •) is Lipschitz in (x, y) with coefficient L n+ 1 2 ≤ L n+1 ((1 + Ch) 1/2 + Ch) + Ch. Next, we can easily verify by (2.10) that v h (t n , •) is Lipschitz in (x, y) with coefficient L n ≤ L n+ 1 2 (1 + Ch) + Ch.
Therefore, the proof is concluded by the discrete Gronwall inequality.

We can also prove that v h is 1/2-Hölder in t as was done in Lemma 3.17 of [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF] 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).

Proposition 2.15. Let Assumptions F, H and M hold true, and Φ be Lipschitz and uniformly bounded. Then |v h (t n , x, y)

-Φ(x, y)| ≤ C √ T -t n .
Proof. We first introduce vh as the numerical solution of (2.4) computed by scheme T h , i.e. vh (T,

•) := Φ(•) and vh (t n , •) := T h [v h ](t n , •).
Clearly, by Lemmas 3.14 and 3.17 of [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF], (v h ) h>0 is uniformly bounded and satisfies

|v h (t n , •) -Φ(•)| ≤ C(T -t n ) 1/2 , uniformly in h. (2.18)
We claim that

|v h (t n , x, y) -v h (t n , x, y)| ≤ C(T -t n ). (2.19)
Then by (2.18), we conclude the proof. Thus it is enough to prove the claim (2.19).

We first recall that by Assumption F and (2.6), for a constant c ∈ R, we have

T h [v h + c](t, x, y) ≤ T h [v h ](t, x, y) + c + hF r |c|. Suppose that for L large enough, |v h (t n+1 , x, y) -v h (t n+1 , x, y)| ≤ L(T -t n+1 ).
It follows by the monotonicity of T h and the uniform boundedness of v h and vh that

|v h (t n , x, y) -v h (t n+ 1 2 , x, y)| ≤ L(T -t n+1 ) + Ch.
And hence by (2.10),

|v h (t n , x, y) -v h (t n , x, y)| ≤ L(T -t n+1 ) + 2Ch ≤ L(T -t n ), which confirms (2.19).
We remark finally that with Propositions 2.10, 2.12, 2.13, 2.14 and 2.15 together with Theorem 2.1 of Barles and Souganidis [START_REF] Barles | Convergence of approximation schemes for fully nonlinear second order equations[END_REF], Theorem 2.6 holds true.

Proof for rate of convergence

As in [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF], 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 [START_REF] Barles | Error bounds for monotone approximation schemes for parabolic Hamilton-Jacobi-Bellman equations[END_REF]. We first derive some technical Lemmas similar to that in [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF].

S h • T h [ϕ + δ](t, x, y) ≤ S h • T h [ψ](t, x, y) + δ(t) -h(b -λc), ∀t ≤ T -h and x ∈ R d .
EJP 18 (2013), paper 15.

A splitting method for nonlinear degenerate PDEs Proof. First, from the proof of Lemma 3.21 in [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF], we have

T h [ϕ + δ](t, x, y) ≤ T h [ϕ](t, x, y) + (1 + hλ 1 ) δ(t + h).
It follows by the definition of the splitting scheme S h • T h in (2.10) that

S h • T h [ϕ + δ](t, x, y) ≤ S h • T h [ϕ](t, x, y) + (1 + hλ 1 )(1 + hλ 2 ) δ(t + h).
By the monotonicity of the splitting scheme S h • T h , we get

S h • T h [ϕ + δ](t, x, y) ≤ S h • T h [ψ](t, x, y) + δ(t) + ζ(t), where ζ(t) := (1 + hλ)δ(t + h) -δ(t).
Finally, using exactly the same arguments as in the proof of Lemma 3.5 of [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF], it follows that

ζ(t) ≤ -h(b -λc),
which concludes the proof. Proposition 2.17. Let Assumptions F, H and M hold true, h ≤ 1 and ϕ, ψ be two bounded functions defined on Q T satisfying

1 h ϕ -S h • T h [ϕ] ≤ g 1 and 1 h (ψ -S h • T h [ψ]) ≥ g 2 , on Q T
for some bounded functions g 1 and g 2 . Then for every n = 0,

• • • , N, (ϕ -ψ)(t n , x, y) ≤ e λ(T -tn) |(ϕ -ψ) + (T, •)| 0 + (T -h)e λ(T -tn) |(g 1 -g 2 ) + | 0 , with some constant λ ≥ |F r | ∞ + sup α,β |c α,β | 0 + |F r | ∞ sup α,β |c α,β | 0 .
Proof. With Lemma 2.16, the proof is exactly the same as in Proposition 3.20 of [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF].

Note that we replace β by λ in our proposition. Now, we are ready to give the Proof of Theorem 2.7 (i). First, under Assumption HJB, we can rewrite the original PDE (2.2) as a standard HJB

-∂ t v - inf α∈A,γ∈C (l α + l γ ) + (c α + c γ )v + (f α + f γ ) • D x v + g α • D y v + 1 2 (σ γ σ γT ) • D 2 xx v = 0.
With Assumption HJB and the Lipschitz terminal condition, it satisfies a comparison result and admits a unique viscosity solution in C 0,1 (Q T ) (see e.g. Proposition 2.1 of [START_REF] Barles | Error bounds for monotone approximation schemes for parabolic Hamilton-Jacobi-Bellman equations[END_REF]). Then by the shaking coefficients method, we can construct a bounded subsolution

v ε ∈ C 0,1 (Q T ) such that v -ε ≤ v ε ≤ v. Let ρ ∈ C ∞ c (Q T ) be a positive function supported in (t, x, y) : t ∈ [0, 1], |x| ≤ 1, |y| ≤ 1
with unit mass, and define

w ε (t, x, y) := v ε * ρ ε , where ρ ε (t, x, y) := 1 ε d+d +2 ρ t ε 2 , x ε , y ε .
Then w ε is a smooth subsolution of (2.2) and satisfies |w ε -v| ≤ 2ε. Moreover, since v ε ∈ C 0,1 (Q T ) is uniformly Lipschitz in (x, y) and 1/2-Hölder in t, it follows that

w ε ∈ C ∞ , and ∂ η0 t D η1+η2 x η 1 y η 2 w ε ≤ Cε 1-2η0-|η1|-|η2| , ∀(η 0 , η 1 , η 2 ) ∈ N 1+d+d \ {0}. (2.20)
A splitting method for nonlinear degenerate PDEs Now, let us consider the consistency error function Λ w ε h (t, x, y) defined in (2.13). By Proposition 2.12 and (2.20), it follows that there exists a constant C independent of ε and h for 0 ≤ h ≤ 1 such that

|Λ w ε h | 0 ≤ R(h, ε) := Chε -3 .
(2.21) Moreover, since w ε is a subsolution of equation (2.2), it follows by the definition of Λ w ε h in (2.13) that

w ε ≤ S h • T h [w ε ] + Ch 2 ε -3 .
Finally, by Proposition 2.17, we get

w ε -v h ≤ C (ε + hε -3 ), and v -v h = v -w ε + w ε -v h ≤ C(ε + hε -3 )
and it follows by a minimization technique on ε that

v -v h ≤ C inf ε>0 ε + hε -3 ≤ C h 1 4 .
(2.22)

Proof of Theorem 2.7 (ii) : Under Assumption HJB+, we can apply the switching system method of Barles and Jakobsen [START_REF] Barles | Error bounds for monotone approximation schemes for parabolic Hamilton-Jacobi-Bellman equations[END_REF] which constructs a smooth supersolution closed to viscosity solution to PDE (2.2) and provides the lower bound:

v -v h ≥ -inf ε>0 Cε 1 3 + R(h, ε) = -C h 1 10 , (2.23)
where R(h, ε) is defined in (2.21).

Basis projection and simulation-regression method

To get an implementable scheme, we need to specify how to compute the expecta-

tions E ϕ(t n+1 , Xtn,x h , y)H tn,x,h i (∆W n+1 ) i=0,1,2
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 [START_REF] Bouchard | Discrete-time approximation and Monte-Carlo simulation of backward stochastic differential equations[END_REF], function basis regression [START_REF] Gobet | A regression-based Monte Carlo method to solve backward stochastic differential equations[END_REF] and cubature method [START_REF] Crisan | Solving Backward Stochastic Differential Equations using the Cubature Method[END_REF], 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 approximation X of X

Xtn+1 := Xtn + µ(t n , Xtn )h + σ(t n , Xtn )∆W n+1 , (3.1)
where ∆W n+1 := W tn+1 -W tn . Then with simulations of process X as well as W , one can estimate the conditional expectations

E ϕ(t n+1 , Xtn+1 , y) H tn, Xtn ,h i (∆W n+1 ) Xtn i=0,1,2 .
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 [START_REF] Longstaff | Valuing American options by simulation: a simple leastsquare approach[END_REF] in American option pricing context or Gobet, Lemor and Warin [START_REF] Gobet | A regression-based Monte Carlo method to solve backward stochastic differential equations[END_REF] in BSDE context to estimate the conditional expectations

E ϕ(t n+1 , Xtn+1 , Y )H tn, Xtn ,h i (∆W n+1 ) Xtn , Y i=0,1,2 , (3.2)
with which we shall make the splitting scheme S h • T h implementable.

Remark 3.1. (i)

The distribution of Y may be chosen arbitrarily according to the concrete context. (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.

Basis projection scheme and simulation-regression method

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 λi := arg min λ E ϕ(t n+1 , Xtn+1 , Y )H tn, Xtn ,h i (∆W n+1 ) - K k=1 λ k e k ( Xtn , Y ) 2 , (3.3)
then the projected approximation of (3.2) is denoted by

Ẽ ϕ(t n+1 , Xtn+1 , Y )H tn, Xtn ,h i (∆W n+1 ) Xtn , Y := K k=1 λi k e k ( Xtn , Y ). (3.4)

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 [START_REF] Bouchard | Monte-Carlo valorisation of American options: facts and new algorithms to improve existing methods[END_REF] 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 • Th . Concretely, it is defined as follows: + Th [ṽ h ] t n , x + f α,β (t n , x, y)h, y + g α,β (t n , x, y)h . (3.5)

Th [ṽ h ](t n , x, y) := Ẽ ṽh (t n+1 , Xtn,x h , y) + hF (•, ẼDṽ h (•))(t n , x,

Simulation-regression method

Next, we propose to use a simulation-regression method to approximate λ. We still suppose that d = d = 1 for simplicity. Let ( Xm tn ) 0≤n≤N , (∆W m n ) 0<n≤N , Y m 1≤m≤M be M independent simulations of X, ∆W and Y , where X 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:

λi,M = arg min λ M m=1 ϕ(t n+1 , Xm tn+1 , Y m )H tn, Xm tn ,h i (∆W m n+1 ) - K k=1 λ k e k ( Xm tn , Y m ) 2 . (3.6)
A raw regression estimation of the conditional expectations (3.2) from these M samples is given by

ĒM ϕ(t n+1 , Xtn+1 , Y )H tn, Xtn ,h i (∆W n+1 ) Xtn , Y := K k=1 λi,M k e k ( Xtn , Y ), i = 0, 1 , 2.
:= Γ i ( Xtn , Y ) ∨ ĒM ϕ(t n+1 , Xtn+1 , Y )H tn, Xtn ,h i (∆W n+1 ) Xtn , Y ∧ Γ i ( Xtn , Y ).
Remark 3.3. As observed in Bouchard and Touzi [START_REF] Bouchard | Discrete-time approximation and Monte-Carlo simulation of backward stochastic differential equations[END_REF], the truncation method is an important technique to obtain a L p -convergence. By Lemma (2.15), we can choose Γ 0 (x, y) = Γ 0 (x, y) and Γ 0 (x, y) = -Γ 0 (x, y) with a function Γ 0 satisfying Γ 0 (x, y) ≤ Φ(x, y) + C T -t n for some constant C.

(3.9) Remark 3.4. In Gobet et al. [START_REF] Gobet | A regression-based Monte Carlo method to solve backward stochastic differential equations[END_REF], the authors propose the following minimization problem in place of (3.6):

min λ 0 ,λ 1 M m=1 ϕ(t n+1 , Xm tn+1 , Y m ) - K k=1 λ 0 k e k ( Xm tn , Y m ) - K k=1 λ 1 k e k ( Xm tn , Y m )∆W m n+1 2 ,
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)

vh (t n , x, y) = S h • TM h [v h ](t n , x, y) := inf α∈A sup β∈B hl α,β (t n , x, y) + hc α,β (t n , x, y) TM h [v h ](t n , x, y) (3.10) + TM h [v h ] t n , x + f α,β (t n , 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 • Th and S h • TM h , we can no longer use the same arguments as in Fahim et al. [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF], since there is no uniform convergence property in L p for the Monte-Carlo error ( ÊM -E)(R) as in the Assumption E of [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF]. 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 [START_REF] Barles | Convergence of approximation schemes for fully nonlinear second order equations[END_REF],

we can not hope to obtain a probabilistic L 2 (Ω)-convergence as in Gobet et al. [START_REF] Gobet | A regression-based Monte Carlo method to solve backward stochastic differential equations[END_REF]. 3), 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, Ẽ ϕ(t n+1 , Xtn+1 , Y )H tn, Xtn ,h i (∆W n+1 ) Xtn , Y (3.11)
= K k=1 E ϕ(t n+1 , Xtn+1 , Y )H tn, Xtn ,h i (∆W n+1 ) ( Xtn , Y ) ∈ 1 B k 1 B k ( Xtn , Y ).
Let us use (e k ) 1≤k≤K = (1 B k ) 1≤k≤K as projection basis in (3.3) and (3.6), where (B k ) 1≤k≤K is a partition of D. Given a bounded function ϕ on D, a process X and a random variable Y , we shall consider the random variables of the form

R i (ϕ) := ϕ(t n+1 , Xtn+1 , Y ) H tn, Xtn ,h i (∆W n+1 ), i = 0, 1, 2,
Γ i (x, y) 2 + Γ i (x, y) 2 ≤ C Γ(x, y) 2 h -i , for some function Γ(x, y). Then for every (x, y) ∈ B k , E ( ÊM -Ẽ) 2 [R i (ϕ) | Xtn = x, Y = y] ≤ C 1 M h -i |ϕ| 2 0 + Γ 2 (x, y) P(( Xtn , Y ) ∈ B k ) . (3.13)
The proof is almost the same as that of Theorem 5.1 of Bouchard and Touzi [START_REF] Bouchard | Discrete-time approximation and Monte-Carlo simulation of backward stochastic differential equations[END_REF], 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(( Xtn , Y ) ∈ B k ) ≈ Cδ d+d , where C depends on the density of ( Xtn , Y ). As the total volume of D is fixed and finite, let

Ĉ(δ) := sup N,n,k,x,y C 1 M h -i b 2 + Γ 2 (x, y) P(( Xtn , Y ) ∈ B k ) , (3.14)
it follows that Ĉ(δ) ≈ Cδ -(d+d ) . Now, let us give a local uniform convergence as well as a rate of convergence for the simulation-regression scheme S h • TM h .

Theorem 3.9. Let Assumptions F, H and M hold true, F be uniformly bounded, Φ be bounded and Lipschitz continuous, and the PDE (2.2) satisfy a comparison result for bounded viscosity solutions. In addition, given a time step h, there is a D-partition hypercubes (B h k ) 1≤k≤K h with edge δ h such that δ h h -1 → 0 as h → 0. Let the truncation function Γ 0 satisfies (3.9), and we use hypercubes (1 B h k ) 1≤k≤K h as projection basis func- tions and with sample number (3.14). Then there exists a function v, such that vh → v locally uniformly, a.s. where vh is the numerical solution of scheme S h • TM h defined in (3.10) with terminal condition Φ. Moreover, v is the unique bounded viscosity solution of (2.2) and (2.3). Theorem 3.10. Let Assumption HJB+ hold, Φ be bounded and Lipschitz continuous, and assume that we use hypercubes (1 B h k ) 1≤k≤K h as projection basis functions whose longest edge satisfies δ h ≤ Ch 11 10 , and we choose simulation number M = M h such that

M = M h such that Ĉ(δ h ) h -2 M -1 h → 0, where Ĉ(δ h ) is defined in
lim sup h→0 h -1 20 -2 Ĉ(δ) M -1 h < ∞.
Then there is a constant C > 0, s.t.

v -vh L 2 (Ω) ≤ Ch 1 10
, where vh is the numerical solution of scheme S h • TM h in (3.10) with terminal condition Φ and v is the unique bounded viscosity solution of (2.2) and (2.3).

Some analysis on the basis projection scheme S h • Th

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

Proposition 3.11. Let (B k ) 1≤k≤K be a partition of domain D, and the three projections (i = 0, 1, 2) of (3.3) use the same hypercubes (1 B k ) 1≤k≤K as projection function basis.

Then under Assumptions F, H and M,

• i) The basis projection scheme S h • Th is monotone.

• ii) If the terminal condition Φ is uniformly bounded, then the numerical solution ṽh of scheme S h • Th in (3.5) is uniformly bounded for h small enough.

Proof. In view of Remark 3.7, we replace the conditional expectations in S h • T h by the new conditional expectations (3.11), and then get the projection scheme S h • Th . Therefore, all the arguments still hold in the proof of Lemma 3.2 and 3.3 of [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF] for Th , so do Propositions 2.10 and 2.13. Therefore, Proposition 3.11 holds true.

Similar to the consistency error function Λ ϕ h for scheme S h • T h defined in (2.13), we define the consistency error function Λϕ h for scheme S h • Th by with E(ϕ) defined in Proposition 2.12.

Λϕ h (•) := ϕ(•) -S h • Th [ϕ](•) h + L X ϕ(•) + F (•, ϕ, D x ϕ, D 2 xx ϕ) + H(•, ϕ, D x ϕ, D y ϕ). ( 3
Proof. In view of Remark 3.7, the error caused by conditional expectation on hypercube is bounded by Cδ|D i+1 z i+1 ϕ| 0 for D i z i ϕ. Thus we get immediately the new consistency error Ẽ(ϕ) with Proposition 2.12. Proposition 3.13. Suppose that the three projections in (3.3) use the same D-partition hypercubes as projection function basis, then Lemma 2.16 and Proposition 2.17 hold true if we replace the scheme S h • T h by S h • Th .

Proof.

With the Proposition 3.11 and under Assumptions F, H and M, we see that all the arguments are still true in the proofs of Lemma 2.16 and Proposition 2.17 for scheme S h • Th , in view of Remark 3.7. So we get the same results for the basis projection scheme S h • Th .

The proof for convergence results of scheme S h • TM h

To prove Theorem 3.9, we shall mimic the proof of Theorem 4.1 in [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF], which uses the arguments of [START_REF] Barles | Convergence of approximation schemes for fully nonlinear second order equations[END_REF] 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 that v * (T, x, y) = v * (T, x, y) = Φ(x, y). Then it is enough to prove that v * and v * 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. A splitting method for nonlinear degenerate PDEs by uniform boundedness of vh and manipulation on ϕ, there is a sequence (t k , x k , y k , h k ) → (t 0 , x 0 , y 0 , 0) such that vh k (t k , x k , y k ) → v * (t 0 , x 0 , y 0 ) and

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 ), EJP 18
C k := min(v h k -ϕ) = (v h k -ϕ)(t k , x k , y k ) → 0.
From the monotonicity of scheme S h • Th , it follows that

S h • Th [ vh k ] ≥ S h • Th [ ϕ + C k ],
and hence

0 = vh k (t k , x k , y k ) -S h • TM h [v h k ](t k , x k , y k ) = vh k (t k , x k , y k ) -S h • Th [v h k ](t k , x k , y k ) + h k R k ≤ ϕ(t k , x k , y k ) + C k -S h • Th [ϕ + C k ](t k , x k , y k ) + h k R k , where R k := h -1 k (S h k • TM h k -S h k • Th k )[v h k ](t k , x k , y k ). We claim that
R k → 0 P-a.s. along some subsequence.

(3.16)

Then, from the consistence of scheme S h • Th in Proposition 3.12,

-L X ϕ -F (•, ϕ, D x ϕ, D 2 xx ϕ) -H(•, ϕ, D x ϕ, D y ϕ) (t 0 , x 0 , y 0 ) ≥ 0,
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 • TM h and S h • Th , and the boundedness of c α,β ,

E|R k | 2 ≤ (1 + Ch k ) 2 1 h 2 k E Th k [v h k ] -TM h k [v h k ] 2 (t k , x k , y k ) ≤ C (1 + Ch k ) 2 1 h 2 k ( Ẽ -ÊM ) 2 R 0 (v h k ) + h k R 1 (v h k ) + h k R 2 (v h k ) ,
where R i (v h k ) is defined in (3.12). And therefore by Lemma 3.8

E|R k | 2 ≤ C (1 + Ch k ) 2 1 h 2 k Ĉ(δ) + h 2 k 1 h k Ĉ(δ) + 1 h 2 k Ĉ(δ) 1 M h k ≤ C h -2 k Ĉ(δ) M -1 h k ,
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.

Proof of Theorem 3.10. With Proposition 3.13, we can proceed as in the proof of Theorem 2.7. Then there is a subsolution w h of (2.2) such that

v ≤ w h + C and w h -ṽh ≤ C( h -3 + h -1 δ + δ -2 ). Moreover, since h -1 vh -S h • Th [v h ] ≥ -R h [v h ], where R h [ϕ] := 1 h (S h • Th -S h • TM h )[ϕ] ,
it follows from Proposition 3.13 that ṽh -vh ≤ C|R h [v h ]|, and then

v -vh = v -ṽh + ṽh -vh ≤ C + h -3 + h -1 δ + δ -2 + R h [v h ] .
Similarly, we have the other side of the error boundary and get v -vh 2 ≤ C (

1 3 + h -3 + h -1 δ + δ -2 ) 2 + R h [v h ] 2 .
(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

E v -vh 2 ≤ C h 1 20 + 1 M h 1 h 2 Ĉ(δ) ≤ C h 1 20 .

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 dS t = rS t dt + σ t S t dW t , where r is the constant interest rate, σ is the volatility process which is bounded between the lower volatility σ and the upper volatility σ. Denote A t := t 0 S u du, an Asian option is an option with payoff g(S T , A T ) at maturity T , whose pricing equation is

∂ t v + rsD s v + 1 2 max σ≤σ≤σ σ 2 s 2 D 2 ss v + sD a v -rv (t, s, a) = 0, (4.1)
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 :

∂ t v + rsD s v + 1 2 σ 2 0 s s D 2 ss v + 1 2 max σ≤σ≤σ (σ 2 -σ 2 0 )s 2 D 2 ss v + sD a v -rv = 0. (4.2)
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%.

∂ 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.

Figure 2: The price of Asian option with payoff (A -K 1 ) + -(A -K 2 ) + in UVM with HW IR and in BS model with HW IR. In case that ∆t = 0.005, it takes 309.4 seconds for the splitting method using 5 × 10 5 simulations.

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 [START_REF] Porchet | Problems of Valuation and Organization in Energy Markets[END_REF].

A hydropower plant manages a dam, which is filled by rain precipitations with nonnegative rate A t , which follows equation dA t = µ a A t dt + σ a A t dW 1 t .

Denote by B t the volume of water in the dam, then dB t = (A t -q t ) dt, 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 A splitting method for nonlinear degenerate PDEs 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 (a -q)D b u + qsD x u = 0.

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 ( Xtn , Y )) 1≤k≤K as defined in (3.3)

  )

Lemma 2 . 16 . 3 +

 2163 Let Assumptions F, H and M hold true and h ≤ 1, define λ 1 := |F r | ∞ , λ 2 := sup α,β |c α,β | 0 , λ := λ 1 +λ 2 +λ 1 λ 2 . Then, for every (a, b, c) ∈ R , and every bounded function ϕ ≤ ψ defined on Q T , with function δ(t) := e λ(T -t) (a + b(T -t)) + c, we have

(3. 7 )

 7 Then with a priori upper bounds Γ i ( Xtn , Y ) and lower bounds Γ i ( Xtn , Y ), we define the regression estimation of (3.2):ÊM ϕ(t n+1 , Xtn+1 , Y )H tn, Xtn ,h i (∆W n+1 ) Xtn , Y(3.8)

Definition 3 . 5 .Remark 3 . 6 .Remark 3 . 7 .

 353637 However, we can get a convergence result if we choose the local hypercubes as function basis. Let us restrict the numerical resolution on [0, T ] × D instead of Q T , where D ⊂ R d+d is a bounded domain. Clearly, we need to assume that the boundary conditions on the domain D c := R d+d \ D are available for scheme S h • TM h . Given a domain D ⊆ R d+d , a class of hypercube sets (B k ) 1≤k≤K is called a partition of D whenever ∪ K k=1 B k = D and B i ∩ B j = ∅. The simplest examples of partition of D is the uniform partition. With uniform interval [x k , x k ) and [y k , y k ), B k are of the form [x k , x k ) × [y k , y k ). Recently, Bouchard and Warin [6] proposed a partition based on the simulations. They first sort all the simulations and then divide the space in a non-uniform way such that they have the same number of simulation particles in every hypercube B k . If we use hypercubes (1 B k ) 1≤k≤K as basis function in the projections (3.

( 3

 3 .12) and then give an estimation for the regression error ( ÊM -Ẽ) [R i (ϕ) | Xtn = x, Y = y]. EJP 18 (2013), paper 15.

Lemma 3 . 8 .

 38 Suppose that the a priori estimations used in (3.8) satisfy

  (2013), paper 15.

Figure 1 :

 1 Figure 1: The comparison of some numerical methods for pricing Asian option with payoff (A -K 1 ) + -(A -K 2 ) + in UVM, with parameters S 0 = 100, K 1 = 90, K 2 = 110, T = 1, r = 0.05, σ = 0.18, σ = 0.22 and σ 0 = 0.2. When ∆t = 0.005, a single computation takes 3.74 seconds for finite difference method, and 131.1 seconds for our splitting method using 5 × 10 5 simulations.

  , and λi,M k be simulated regression estimations of λi k with M simulations of X, Y as defined in(3.6). Then for (x, y)∈ B k , Ẽ R i (ϕ) | Xtn = x, Y = y = λi k and ÊM R i (ϕ) | Xtn = x, Y = y = Γ i (x, y) ∨ λi,M k ∧ Γ i (x, y). i (ϕ)e k ( Xtn , Y )] E[e 2 k ( Xtn , Y )]and λi,Mk = E M [R i (ϕ)e k ( Xtn , Y )] E M [e 2 k ( Xtn , Y )],where E M is the empirical expectation defined as follows: given M simulations (U m ) 1≤m≤M of random variable U , E M [U ] := 1 M M m=1 U m .

Figure 3 :Proof of Lemma 3 . 8 .

 338 Figure 3: Solution of optimal management for a hydropower plant, with σ x = 1 and σ x = 1.5. Using 2 × 10 6 simulations, the splitting scheme takes 639.2 seconds for a single calculation when ∆t = 0.0333.

  ee]) 2 Var(Re) + λ2 Var(ee) + 2Γ 2 Var(ee) .

( 5 . 1 )

 51 When e = 1 B k , we have E[e 2 ( Xtn , Y )] = E[e( Xtn , Y )] = P(( Xtn , Y ) ∈ B k ), E[eR i ] ≤ C |ϕ| 0 h i/2 E[e] and λ ≤ C |ϕ| 0 h i/2, and then it follows by (5.1) that (3.13) holds true.

Remark 2.4. Assumption F is

  almost the same as the Assumption F in[START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF], here we just add a variable y in the nonlinear operator F .

	Assumption H : The coefficients in Hamiltonian H are all uniformly bounded, i.e.
	sup
	(α,β)∈A×B, 1≤i≤d, 1≤j≤d

  and denote the new numerical splitting scheme by S h • TM

h , which is

TM h [v h ](t n , x, y) := ÊM vh (t n+1 , Xtn,x h , y) + h F (•, ÊM Dv h (•))(t n , x,

y), and ÊM D h ϕ(t n , x, y) = ÊM ϕ(t n+1 , Xtn,x h , y)H tn,x,h i (∆W n+1 ) : i = 0, 1, 2 , so that S h • TM h is defined by

  .15) Let δ denote the longest edge of hypercubes (B k ) 1≤k≤K , then the projection error for a Lipschitz continuous function is proportional to δ. Moreover, if we use hypercubes (1 B k ) 1≤k≤K as projection function basis, then under Assumptions F, H and M, the consistency error function Λϕ h is uniformly bounded by h Ẽ(ϕ), where Ẽ(ϕ) := E(ϕ) + Ch -1 δ |D z ϕ| 0 + h|D 2 zz ϕ| 0 + h|D 3 zzz ϕ| 0 , for z := (x, y),

	Proposition 3.12.

  in a stochastic context. Given vh the numerical solution of scheme S h • TM

	Proof of Theorem 3.9. h , we denote
	v * (t, x, y) :=	lim inf	vh (t , x , y ), v * (t, x, y) :=	lim sup	vh (t , x , y ).
		(t ,x ,y ,h)→(t,x,y,0)		(t ,x ,y ,h)→(t,x,y,0)

4.1.2 Asian option pricing in UVM with Hull-White interest rate: a three- dimensional 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 dynamicdr t = b(θ t -r t )dt + σ r dB t ,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

  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 uAs in the examples in Section 5.2 of[START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF], we truncate the optimization on θ and rewrite the equation in form of (2.2).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

		+ 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
								sa u +	1 2	σ 2 x D 2 xx u
	+	max -n≤θ≤n	θ 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 -	1 2	σ 2 x D 2 xx u
	+ max q

q (a -q)D b u + qsD x u = 0.

ejp.ejpecp.org

Acknowledgments. The author is grateful to J. Frédéric Bonnans, Damien Lamberton and Nizar Touzi for fruitful discussions and two anonymous referees for helpful comments and suggestions. This research is supported by the Chair Financial Risks of the Risk Foundation sponsored by Société Générale, the Chair Derivatives of the Future sponsored by the Fédération Bancaire Française, and the Chair Finance and Sustainable Development sponsored by EDF and CA-CIB.