Olivier Bokanowski

AND Athena Picarelli

Hasnaa Zidani

DYNAMIC PROGRAMMING AND ERROR ESTIMATES FOR STOCHASTIC CONTROL PROBLEMS WITH MAXIMUM COST

Keywords: AMS subject classifications: 49J20, 49L25, 65M15, 35K55 Hamilton-Jacobi equations, oblique Neuman boundary condition, error estimate, viscosity notion, reachable sets under state constraints, lookback options, maximum cost

HAL is a

Introduction

Consider the controlled process X u t,x (•) solution of the stochastic differential equation dX(s) = b(s, X(s), u(s))ds + σ(s, X(s), u(s))dB (s) s

∈ [t, T ] X(t) = x (1.1)
(where u ∈ U denotes the control process). We are interested in an optimal control problem where the cost function is given by :

J(t, x, u) = E ψ X u t,x (T ), max s∈[t,T ] g(X u t,x (s)) , (1.2)
and the value function is defined by: v(t, x) := inf u∈U J(t, x, u). (1.3) This problem arises from the study of some path-dependent options in finance (lookback options). Another application of control problem with maximum running cost comes from the characterization of reachable sets for state-constrained controlled stochastic systems (see section 2.3).

The main contributions to the study of this kind of problems can be found in [START_REF] Barles | Optimal control on the L ∞ norm of a diffusion process[END_REF] and [START_REF] Barron | The Bellman equation for control of the running max of a diffusion and applications to lookback options[END_REF] (see also [START_REF] Heinricher | Optimal control of the running max[END_REF] for the stationary case in the framework of classical solutions). In these works the dynamic programming techniques are applied on the L p -approximation of the L ∞ -cost functional (1.2), using the approximation:

a ∨ b ≃ (a p + b p ) 1 p (for p → ∞),
for any a, b ≥ 0, where a ∨ b := max(a, b). Then a characterization of the auxiliary value function v (see definition (2.4)) as solution to a suitable Hamilton-Jacobi equation is obtained as limit for p → ∞. A fundamental hypothesis in order to apply this approach is the positivity of the functions involved.

In our work, a direct derivation of a Dynamic Programming Principle for the function v gives an alternative and natural way for dealing with the maximum running cost problems under less restrictive assumptions. By this way, the optimal control problem (1.3) is connected to the solution of a HJB equation with oblique derivative boundary conditions. Here, the boundary conditions have to be understood in the viscosity sense (see [START_REF] Ishii | A boundary value problem of the Dirichlet type for Hamilton-Jacobi equations[END_REF]).

The second part of the paper is devoted to numerical aspects. First, a general scheme is proposed to deal with HJB equations with oblique derivative boundary conditions. Then, a convergence result is proved by using the general framework of Barles-Souganidis [START_REF] Barles | Convergence of approximation schemes for fully nonlinear second order equations[END_REF] based on the monotonicity, stability, consistency of the scheme (a precise definition of these notions, in the case of HJB equations with oblique derivative boundary condition, is given in section 4.1). Secondly, we focus on the fully discrete semi-Lagrangian scheme and we study some properties of the numerical approximation. Semi-Lagrangian schemes for second order HJB equations (in the general form of Markov chain approximations) have been first considered in the book [START_REF] Kushner | Probability Methods for Approximations in Stochastic Control and for Elliptic Equations[END_REF]. In the framework of viscosity solutions, semi-Lagrangian schemes for first order HJB equations have been studied in [START_REF] Capuzzo-Dolcetta | On a discrete approximation of the Hamilton-Jacobi equation of dynamic programming[END_REF]. Extensions to the second order case can be found in [START_REF] Menaldi | Some estimates for finite difference approximations[END_REF][START_REF] Camilli | An approximation scheme for the optimal control of diffusion processes[END_REF][START_REF] Crandall | Convergent difference schemes for nonlinear parabolic equations and mean curvature motion[END_REF][START_REF] Debrabant | Semi-Lagrangian schemes for linear and fully non-linear diffusion equations[END_REF]. In our context, the numerical scheme that will be studied couples the classical semi-Lagrangian scheme with an additional projection step on the boundary. For this particular scheme, we derive a rate of convergence, generalizing in this way the results already known for general HJB equations without boundary conditions. More precisely, an error bound of the order of h 1/4 (where h is the time step) is obtained for Lipschitz (unbounded) cost function. The activity related to error estimate for numerical approximations of nonlinear second order degenerate equations dated back few years ago with [START_REF] Krylov | On the rate of convergence of finite-difference approximations for Bellman's equations[END_REF][START_REF] Krylov | On the rate of convergence of finite-difference approximations for Bellman's equations with variable coefficients[END_REF] where some useful techniques of "shaking coefficients" were introduced, and then mollifiers were used to get smooth sub-solutions. These ideas were later used extensively by many authors [START_REF] Barles | On the convergence rate of approximation schemes for Hamilton-Jacobi-Bellman equations[END_REF][START_REF] Barles | Error bounds for monotone approximation schemes for Hamilton-Jacobi-Bellman equations[END_REF][START_REF] Barles | Error bounds for monotone approximation schemes for parabolic Hamilton-Jacobi-Bellman equations[END_REF][START_REF] Bonnans | Error bounds for a Stochastic Impulse Control Problem[END_REF].

The paper is organized as follows: Section 2 introduces the problem and some motivations. In particular it will be shown how it is possible to connect a state-constrained reachability problem to an optimal control problem associated to a cost functional of the form (1.2). Section 3 is devoted to the characterization of the value function by the appropriate HJB equation. A strong comparison principle is also stated. In Section 4 the numerical approximation is discussed and a general convergence result is provided. The semi-Lagrangian scheme is presented in Section 4.2 and the properties of this scheme are investigated. Section 5 focuses on the study of error estimates. Finally some numerical tests are presented in Section 6 to analyse the relevance of the proposed scheme.

2. Formulation of the problem 2.1. Setting and preliminary results. Let (Ω, F, P) be a probability space, {F t , t ≥ 0} a filtration on F and B(•) a Brownian motion in R p (with p ≥ 1). Given T > 0 and 0 ≤ t ≤ T , the following system of controlled stochastic differential equations (SDEs) in R d (d ≥ 1) is considered dX(s) = b(s, X(s), u(s))ds + σ(s, X(s), u(s))dB(s)

s

∈ [t, T ] X(t) = x, (2.1)
where u ∈ U set of progressively measurable processes with values in U ⊂ R m (m ≥ 1), U is assumed to be a compact set. The following classical assumption will be considered on the coefficients b and σ: We recall some classical results on (2.1) (see [START_REF] Yong | Stochastic controls: Hamiltonian systems and HJB equations[END_REF]):

Theorem 2.1. Under assumption (H1) there exists a unique F t -adapted process, denoted by X u t,x (•), strong solution of (2.1). Moreover there exists a constant C 0 > 0, such that for any u ∈ U,

0 ≤ t ≤ t ′ ≤ T and x, x ′ ∈ R d E sup θ∈[t,t ′ ] X u t,x (θ) -X u t,x ′ (θ) 2 ≤ C 2 0 |x -x ′ | 2 , (2.2a) E sup θ∈[t,t ′ ] X u t,x (θ) -X u t ′ ,x (θ) 2 ≤ C 2 0 (1 + |x| 2 ) |t -t ′ |, (2.2b) E sup θ∈[t,t ′ ] X u t,x (θ) -x 2 ≤ C 2 0 (1 + |x| 2 ) |t -t ′ |. (2.2c)
Consider two functions ψ, g such that:

(H2) ψ : R d × R → R and g : R d → R are Lipschitz functions. Their Lipschitz constants will be denoted by L ψ and L g respectively.

In this paper we deal with optimal control problems of the form

v(t, x) := inf u∈U E ψ X u t,x (T ), max s∈[t,T ] g(X u t,x (s)) . (2.3)
Under assumptions (H1) and (H2), v is a Lipschitz continuous function in x and a 1 2 -Hölder continuous function in t (by the same arguments as in Proposition 2.2 below). In order to characterize the function v as solution of an HJB equation the main tool is the well-known optimality principle. The particular non-Markovian structure of the cost functional in (2.3) prohibits the direct use of the standard techniques. To avoid this difficulty it is classical to reformulate the problem by adding a new variable y ∈ R that, roughly speaking, keeps the information of the running maximum. For this reason, we introduce an auxiliary value function ϑ defined on [0, T ] × R d × R by:

ϑ(t, x, y) := inf u∈U E ψ X u t,x (T ), max s∈[t,T ] g(X u t,x (s)) ∨ y . (2.4)
The following property holds:

ϑ(t, x, g(x)) = v(t, x), (2.5)
so from now on, we will only work with the value function ϑ since the value of v can be obtained by (2.5) .

Proposition 2.2. Under assumptions (H1) and (H2), ϑ is a Lipschitz continuous function in (x, y) uniformly with respect to t, and a 1 2 -Hölder continuous function in t. More precisely, there exists L ϑ > 0 (L ϑ depends only on C 0 , L ψ , L g ) such that:

|ϑ(t, x, y) -ϑ(t ′ , x ′ , y ′ )| ≤ L ϑ (|x -x ′ | + |y -y ′ | + (1 + |x|)|t -t ′ | 1 2 ), for all (x, y), (x ′ , y ′ ) ∈ R d+1 and for any t ≤ t ′ ∈ [0, T ].
Proof. Let be t ≤ t ′ ≤ T and x, x ′ ∈ R d . Notice that the following property holds for the maximum

|(a ∨ b) -(c ∨ d)| ≤ |a -c| ∨ |b -d| .
Then, the inequalities (2.2) yield to:

ϑ(t, x, y) -ϑ(t, x ′ , y) ≤ Ksup u∈U E sup s∈[t,T ] X u t,x (s) -X u t,x ′ (s) ≤ KC 0 |x -x ′ | for K = L ψ (L g + 1)
. In a similar way, we obtain

ϑ(t, x, y) -ϑ(t ′ , x, y) ≤ L ψ sup u∈U E X u t,x (T ) -X u t ′ ,x (T ) + L g sup s∈[t,t ′ ] X u t,x (s) -x + + L g sup s∈[t ′ ,T ] X u t ′ ,x (s) -X u t ′ ,X u t,x (t ′ ) (s) ≤ K ′ (1 + |x|) t -t ′ 1 2
for a positive constant K ′ > 0 that depends only on L ψ , L g and C 0 . The L ψ -Lipschitz behavior in the variable y is immediate. We conclude then the result with

L ϑ = K ∨ K ′ ∨ L ψ .
  -∂ t ϑ + H(t, x, D x ϑ, D 2 x ϑ) = 0 in [0, T ) × D -∂ y ϑ = 0 on[0, T ) × ∂D ϑ(T, x, y) = ψ(x, y) in D (2.6)
for a suitable choice of the domain D ⊂ R d+1 . Theorem 3.4 states a comparison principle for (2.6), adapting the result of [START_REF] Dupuis | On oblique derivative problems for fully nonlinear secondorder elliptic partial differential equations on nonsmooth domains[END_REF] to parabolic equations in unbounded domains. Thanks to this result, ϑ is finally characterized as the unique continuous viscosity solution of the HJB equation.

Starting from Section 4 the numerical aspects will be discussed. A numerical scheme for the approximation of ϑ will be proposed. Theorem 4.1 gives a general convergence result under hypotheses of stability, consistency and monotonicity.

Error estimates for the semi-discrete semi-Lagrangian scheme are investigated in the Theorems 5.4 and 5.5, while the rate of convergence for the fully discrete scheme is given in Theorem 5.6.

Motivations.

Reachability for state constrained control systems. A first motivation comes from the reachability analysis for state-constrained stochastic controlled systems. Let C ⊆ R d be a target set, and let K ⊆ R d be a set of state constraints. Let d C (resp. d K ) denotes the Euclidean distance functions from C (resp. K). The aim is to characterize and compute the backward reachable set R C,K t defined by

R C,K t := x ∈ R d , ∀ε > 0, ∃u ε ∈ U such that E[d C (X uε t,x (T ))] ≤ ε and E[ max s∈[t,T ] d K (X uε t,x (s))] ≤ ε .
Remark 2.3. Straightforward calculations yield to the following formulation:

R C,K t ≡ x ∈ R d , ∀ε > 0, ∃u ε ∈ U, E[d C (X uε t,x (T )) ∨ max s∈[t,T ] d K (X uε t,x (s))] ≤ ε .
(2.7) Indeed, this is easily deduced from the following inequalities:

E d C (X uε t,x (T )) ∨ E max s∈[t,T ] d K (X uε t,x (s)) ≤ E d C (X uε t,x (T )) ∨ max s∈[t,T ] d K (X uε t,x (s)) ≤ E d C (X uε t,x (T )) + E max s∈[t,T ] d K (X uε t,x (s)) .
Many works in the literature have been devoted to the reachability analysis for general stochastic controlled systems, we refer to [START_REF] Aubin | Stochastic viability and invariance[END_REF][START_REF] Aubin | The viability theorem for stochastic differential inclusions[END_REF]. Stochastic target problems arising in finance have been also analysed [START_REF] Soner | Stochastic target problems, dynamic programming and viscosity solutions[END_REF][START_REF] Bouchard | Stochastic targets with mixed diffusion processes and viscosity solutions[END_REF]. Here, we suggest to characterize the reachable sets by using the level set approach introduced by Osher and Sethian in [START_REF] Osher | Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF] in the deterministic case. At the basis of this approach there is the idea of considering the set R C,K t as a level set of a continuous function defined as the value function of a suitable optimal control problem. Recently, the level-set approach has been extended to deal with general deterministic optimal control problems in presence of state constraints, see for example [START_REF] Bokanowski | Reachability and minimal times for state constrained nonlinear problems without any controllability assumption[END_REF][START_REF] Altarovici | A general Hamilton-Jacobi framework for non-linear state-constrained control problems[END_REF]. The level set approach can be also extended to the stochastic setting by considering the following optimal control problem:

v(t, x) := inf u∈U E d C (X u t,x (T )) ∨ max s∈[t,T ] d K (X u t,x (s)) . (2.8)
Then the following equivalence can easily be proved Theorem 2.4. Under assumption (H1), for every t ≥ 0, we have:

x ∈ R C,K t ⇐⇒ v(t, x) = 0. (2.9)
Proof. The result follows immediately from Remark 2.3.

In view of Theorem 2.4, the reachable set R C,K t can be characterized by means of the value function of a control problem with supremum cost in the form of (2.3), with functions g and ψ defined by: g(x) := d K (x) and ψ(x, y) := d C (x) ∨ y.

Lookback options in finance. Another interest for computing expectation of supremum cost functionals is the lookback options in Finance. The value of such an option is typically of the form

v(t, x) = E e -T t r(s)ds ψ X T , max θ∈[t,T ] g(X θ )) ,
where X θ (the "asset") is a solution of a one-dimensional SDE (2.1), g(x) = x, r(•) is the interest rate, and ψ is the payoff function. Here the option value depends not only on the value of the asset at time T but also on all the values taken between times t and T . A detailed description of this model can be found in [START_REF] Wilmott | Paul Wilmott on quantitative finance[END_REF]. Typical payoff functions are ψ(x, y) = y -x (lookback floating strike put), ψ(x, y) = max(y -E, 0) (fixed strike lookback call), ψ(x, y) = max(min(y, E) -x, 0) (lookback limited-risk put), etc., see [START_REF] Wilmott | Paul Wilmott on quantitative finance[END_REF][START_REF] Barles | Convergence of numerical schemes for parabolic equations arising in finance theory[END_REF] (see also [START_REF] Barles | Optimal control on the L ∞ norm of a diffusion process[END_REF] for other examples and related american lookback options).

The HJB equation

In all the sequel, we will use the abbreviation a.s. for almost surely. Also, almost every will be abbreviated as a.e.. The abbreviation usc (resp. lsc) will be used for upper semi-continuous (resp. lower semi-continuous).

3.1. Dynamic Programming. The aim of this section is to characterize the value function ϑ defined by (2.4) as the unique viscosity solution of an HJB equation. Let us define the process

Y u t,x,y (θ) := max s∈[t,θ] g(X u t,x (s)) ∨ y, (3.1)
then the optimal control problem (2.4) can be re-written as

ϑ(t, x, y) := inf u∈U E ψ X u t,x (T ), Y u t,x,y (T ) . (3.2)
Stochastic optimal control problems with running maximum cost in the viscosity solutions framework have been studied in [START_REF] Barles | Optimal control on the L ∞ norm of a diffusion process[END_REF][START_REF] Barron | The Bellman equation for control of the running max of a diffusion and applications to lookback options[END_REF]. The arguments developed in that papers are based on the approximation technique of the L ∞ -norm and they only apply if ψ and g are positive functions. Here, we derive the HJB equation directly without using any approximation. For this, the first step is the Bellman's principle that we state in the following theorem.

Theorem 3.1 (Dynamic Programming Principle). Under hypothesis (H1) and (H2), for all (t, x) ∈ [0, T )×R d and all family of stopping times {θ u , u ∈ U} independent of F t with values on [t, T ]:

ϑ(t, x, y) = inf u∈U E ϑ(θ u , X u t,x (θ u ), Y u t,x,y (θ u )) . (3.3)
A proof of Theorem 3.1 can be obtained by adapting the same arguments developed by Bouchard and Touzi in [START_REF] Bouchard | Weak dynamic programming principle for viscosity solutions[END_REF], thanks to the fact that the couple of variables (X u t,x (•), Y u t,x,y (•)) satisfies the following fundamental property

X u t,x (s) Y u t,x,y (s) = X u θ,X u t,x (θ) (s) Y u θ,X u t,x (θ),Y u t,x,y (θ) (s) a.s.
for any stopping time θ with t ≤ θ ≤ s ≤ T . In our case the proof is even simpler than the one of [START_REF] Bouchard | Weak dynamic programming principle for viscosity solutions[END_REF] thanks to the uniform continuity of ϑ (Proposition 2.2).

3.2.

Hamilton-Jacobi-Bellman equation. Theorem 3.1 is the main tool for proving next result that characterizes ϑ as a solution, in viscosity sense, of a HJB equation with oblique derivative boundary conditions. Set

D := (x, y) ∈ R d+1 |y ≥ g(x) = Epigraph(g), (3.4)
where D is the interior of D.

Theorem 3.2. Under assumptions (H1) and (H2), ϑ is a continuous viscosity solution of the following HJB equation

-∂ t ϑ + H(t, x, D x ϑ, D 2 x ϑ) = 0 in [0, T ) × D (3.5a) -∂ y ϑ = 0 on [0, T ) × ∂D (3.5b) ϑ(T, x, y) = ψ(x, y) in D (3.5c) with H(t, x, p, Q) := sup u∈U -b(t, x, u)p - 1 2 T r[σσ T ](t, x, u)Q . ( 3.6)
(In all the paper D x , D 2 x denote respectively the first and second order derivatives with respect to x.)

Before starting the proof we recall the definition of viscosity solution for problem (3.5) (see [START_REF] Crandall | User's guide to viscosity solutions of second order partial differential equations[END_REF] and the references therein for a complete discussion on weak boundary conditions). Definition 3.3. A usc function ϑ (resp. lsc function ϑ) on [0, T ] × D is a viscosity sub-solution (resp. super-solution) of (3.5), if for each function ϕ ∈ C 1,2 ([0, T ] × D), at each maximum (resp. minimum) point (t, x, y) of ϑ -ϕ (resp. ϑ -ϕ) the following inequalities hold

       -∂ t ϕ + H(t, x, D x ϕ, D 2 x ϕ) ≤ 0 in [0, T ) × D min -∂ t ϕ + H(t, x, D x ϕ, D 2 x ϕ), -∂ y ϕ ≤ 0 on [0, T ) × ∂D min ϑ -ψ, -∂ t ϕ + H(t, x, D x ϕ, D 2 x ϕ), -∂ y ϕ ≤ 0 on {T } × ∂D min ϑ -ψ, -∂ t ϕ + H(t, x, D x ϕ, D 2 x ϕ) ≤ 0 on {T } × D. (3.7) (resp.        -∂ t ϕ + H(t, x, D x ϕ, D 2 x ϕ) ≥ 0 in [0, T ) × D max -∂ t ϕ + H(t, x, D x ϕ, D 2 x ϕ), -∂ y ϕ ≥ 0 on [0, T ) × ∂D max ϑ -ψ, -∂ t ϕ + H(t, x, D x ϕ, D 2 x ϕ), -∂ y ϕ ≥ 0 on {T } × ∂D max ϑ -ψ, -∂ t ϕ + H(t, x, D x ϕ, D 2 x ϕ) ≥ 0 on {T } × D.) (3
.8) Finally a continuous function ϑ is a viscosity solution of (3.5) if it is both a sub-and a super-solution.

Proof of Theorem 3.2. First, from the definition of ϑ, we obtain easily that ϑ(T, x, y) = ψ(x, y). Now, we check that ϑ is a viscosity sub-solution.

Let ϕ ∈ C 1,2 ([0, T ] × D) such that ϑ -ϕ attains a maximum at point ( t, x, ȳ) ∈ [0, T ) × D.
Without loss of generality we can always assume that ( t, x, ȳ) is a strict local maximum point (let us say in a ball of radius r > 0 centered in ( t, x, ȳ)) and ϑ( t, x, ȳ) = ϕ( t, x, ȳ). Thanks to Theorem 3.1, for any u ∈ U and for any sufficiently small stopping time θ = θ u , we have:

ϕ( t, x, ȳ) = ϑ( t, x, ȳ) ≤ E[ϑ(θ, X u t,x (θ), Y u t,x,ȳ (θ))] ≤ E[ϕ(θ, X u t,x (θ), Y u t,x,ȳ (θ))]. (3.9)
Two cases will be considered depending on if the point (x, ȳ) belongs or not to the boundary of D.

-Case 1: g(x) < ȳ. Consider a constant control u(s) ≡ u ∈ U . From the continuity of g and the a.s. continuity of the sample paths it follows that, for a.e. ω ∈ Ω, there exists s(ω

) > t such that g(X u t,x (s)) < ȳ if s ∈ [ t, s(ω))
. Let be h > 0, and let θ be a the following stopping time:

θ := inf s > t|(X u t,x (s), Y u t,x,ȳ (s)) / ∈ B r (x, ȳ) ∧( t + h)∧inf s > t|g(X u t,x (s)) ≥ ȳ , (3.10)
(here and through the paper B r (x 0 ) denotes the ball of radius r > 0 centered at x 0 ). One can easily observe that a.s.

Y u t,x,ȳ ( θ) = ȳ, then by (3.9) ϕ( t, x, ȳ) ≤ E[ϕ( θ, X u t,x ( θ), ȳ)] ∀u ∈ U.
By applying the Ito's formula, and thanks to the smoothness of ϕ, we get:

E θ t -∂ t ϕ -b(s, X u t,x (s), u)D x ϕ - 1 2 T r[σσ T ](s, X u t,x (s), u)D 2 x ϕds ≤ 0.
Since the stopping times

inf s > t|(X u t,x (s), Y u t,x,ȳ (s)) / ∈ B r (x, ȳ)
and inf s > t | g(X u t,x (s)) ≥ ȳ are a.s. strictly greather than t, it follows that a.s. θt > 0. Therefore,

E 1 θ -t θ t -∂ t ϕ -b(s, X u t,x (s), u)D x ϕ - 1 2 T r[σσ T ](s, X u t,x (s), u)D 2 x ϕds ≤ 0.
By the dominate convergence theorem and the mean value theorem, letting h going to 0, it follows

-∂ t ϕ -b( t, x, u)D x ϕ - 1 2 T r[σσ T ]( t, x, u)D 2 x ϕ ≤ 0, ∀u ∈ U,
and finally:

-∂ t ϕ + sup u∈U -b( t, x, u)D x ϕ - 1 2 T r[σσ T ]( t, x, u)D 2 x ϕ ≤ 0.
-Case 2: g(x) = ȳ. Assume that -∂ y ϕ( t, x, ȳ) > 0, otherwise the conclusion is straightforward. Consider a constant control u(s) ≡ u ∈ U . Thanks to the continuity of the sample paths and the smoothness of ϕ, for a.e. ω there is a time s(ω) > t and η > 0 such that:

ϕ(s, X u t,x (s), y) ≤ ϕ(s, X u t,x (s), ȳ) ∀s ∈ [ t, s], y ∈ [ȳ, ȳ + η).
Let θ be the stopping time given by:

θ := inf s > t|(X u t,x (s), Y u t,x,ȳ (s)) / ∈ B r (x, ȳ) ∧ inf s > t|g(X u t,x (s)) / ∈ [ȳ, ȳ + η) ∧ inf s > t|∂ y ϕ(s, X u t,x (s), ȳ) ≥ 0 ∧ ( t + h) .
By (3.9) one has ϕ( t, x, ȳ) ≤ E[ϕ( θ, X u t,x ( θ), ȳ)], which implies:

-∂ t ϕ + sup u∈U -b( t, x, u)Dϕ - 1 2 T r[σσ T ]( t, x, u)D 2 ϕ ≤ 0.
In conclusion,

min -∂ t ϕ + sup u∈U -b( t, x, u)Dϕ - 1 2 T r[σσ T ]( t, x, u)D 2 ϕ , -∂ y ϕ ≤ 0,
and ϑ is a viscosity sub-solution of equation (3.5).

It remains to prove that ϑ is a viscosity super-solution of (3.5). Let

ϕ ∈ C 1,2 ([0, T ]×D) be such that ϑ-ϕ attains a minimum at point ( t, x, ȳ) ∈ [0, T ) × D.
Without loss of generality we can always assume that ( t, x, ȳ) is a strict local minimum point and ϑ( t, x, ȳ) = ϕ( t, x, ȳ). Thanks to the DPP, for any sufficiently small stopping time θ = θ u one has

ϕ( t, x, ȳ) = ϑ( t, x, ȳ) = inf u∈U E[ϑ(θ, X u t,x (θ), Y u t,x,ȳ (θ))] ≥ inf u∈U E[ϕ(θ, X u t,x (θ), Y u t,x,ȳ (θ))].
This can be treated by using again the same arguments as for the sub-solution proof.

Comparison principle.

The section is concluded with a comparison principle for equation (3.5). There is a large literature dealing with Neumann-type or oblique derivative boundary conditions. We refer to [START_REF] Lions | Neumann type boundary conditions for Hamilton-Jacobi equations[END_REF][START_REF] Barles | Fully nonlinear Neumann type boundary conditions for first-order Hamilton-Jacobi equations[END_REF] for the first order case and to [START_REF] Ishii | Viscosity solutions of fully nonlinear Second-order elliptic partial differential equations[END_REF][START_REF] Ishii | Fully nonlinear oblique derivative problems for nonlinear Second-order elliptic PDEs[END_REF] for second order equations. For dealing with this kind of problems some regularity assumption on the domain D is often required. In our case the definition of D is strictly connected with the choice of the function g that is, without further hypothesis, only Lipschitz continuous. The result below has the advantage of taking into account also this non smooth case.

Theorem 3.4. Assume that (H1)-(H2) are satisfied. Let ϑ (resp. ϑ) a usc (resp. lsc) viscosity sub-solution (resp. super-solution) of (3.5) satisfying the following growth condition (for C > 0):

ϑ(t, x, y) ≤ C(1 + |x| + |y|), (resp. ϑ(t, x, y) ≥ -C(1 + |x| + |y|). Then ϑ ≤ ϑ on [0, T ] × D.
Such a comparison principle is proved in [START_REF] Dupuis | On oblique derivative problems for fully nonlinear secondorder elliptic partial differential equations on nonsmooth domains[END_REF] for elliptic equations in bounded domains. The arguments extend the ones used in [START_REF] Ishii | Viscosity solutions of fully nonlinear Second-order elliptic partial differential equations[END_REF][START_REF] Ishii | Fully nonlinear oblique derivative problems for nonlinear Second-order elliptic PDEs[END_REF] for the case when the domain D has a smooth boundary. For convenience of the reader we give the main steps of the proof in Appendix A, and show how the arguments can be adapted for the parabolic case in unbounded domains. Assertions of Theorems 3.2 and 3.4 lead to the following result: Corollary 3.5. Under assumptions (H1) and (H2), the value function ϑ is the unique continuous, Lipschitz in (x,y), viscosity solution of equation

(3.5) on [0, T ] × D.
The uniqueness result is stated on D, which means that the solution to the HJB equation coincides with the value function ϑ on [0, T ] × D. Then ϑ is extended in a unique way on [0, T ] × R d × R by: ϑ(t, x, y) = ϑ(t, x, y ∨ g(x)).

Numerical approximation

Now, we focus on general approximation schemes for (3.5). Let N ≥ 1 be an integer (number of time steps), and let h := T N , and t n := nh, for n = 0, . . . , N . Let ∆x = (∆x 1 , . . . , ∆x d ) ∈ (R * + ) d be a mesh step in R d , ∆y > 0 be a mesh step in R, and ρ := (h, ∆x, ∆y) be a set of mesh steps (in time and space).

For a given ρ, consider the corresponding space grid

G ρ := (x i , y j ) = (i∆x, j∆y), for (i, j) ∈ Z d × Z .
where

for i ∈ Z d , i∆x := (i 1 ∆x 1 , • • • , i N ∆x N ).
For any x ∈ R d , let j x ∈ Z be the upper integer part of g(x) ∆y , i.e,

j x := min j ∈ Z, j∆y ≥ g(x) .
Consider a projection operator (along the direction e y ) Π Gρ defined as follows:

Π Gρ (Φ)(t n , x i , y j ) := Φ(t n , x i , y j ), ∀j ≥ j x i , i ∈ Z d , Φ(t n , x i , y jx i ), ∀j < j x i , i ∈ Z d , for any Φ ∈ C([0, T ] × R d × R, R).
We aim at computing an approximation W n i,j of the value function ϑ(t n , x i , y j ), on the grid G ρ . By W ρ , we will denote the interpolation of the values W n i,j at (t n , x i , y j ). The values W n i,j are defined recusively as follows:

General scheme (GS) 1)Define W N i,j := ψ(x i , y j ∨ g(x i )), ∀i, j.

2)For n = N, . . . , 1, the value W n-1 i,j is obtained as solution of:

S ρ (t n-1 , x i , y j , W n-1 i,j , Π Gρ (W ρ )) = 0, ∀i, j with g(x i ) ≤ y j , (4.1)
W n-1 i,j

:= W n-1 i,jx i , ∀i, j with g(x i ) > y j , (4.2)
where

S ρ : [0, T ] × D × R × C([0, T ] × R d × R, R) → R is a scheme operator.
Typical operator S ρ can be obtained by using an explicit or implicit finite difference method on (3.5a) (see [START_REF] Bonnans | Consistency of generalized finite difference schemes for the stochastic HJB equation[END_REF][START_REF] Bonnans | Numerical schemes for the two dimensional second-order HJB equation[END_REF][START_REF] Kushner | Numerical methods for stochastic control problems in continuous time[END_REF]), or a semi-Lagrangian (SL) scheme ( [START_REF] Menaldi | Some estimates for finite difference approximations[END_REF][START_REF] Camilli | An approximation scheme for the optimal control of diffusion processes[END_REF][START_REF] Debrabant | Semi-Lagrangian schemes for linear and fully non-linear diffusion equations[END_REF]). In (4.1), the value of W n-1 i,j may depend on the whole function W ρ (i.e., all the values W k i,j for k = 0, • • • , N ). Of course, in case of an explicit time discretization, the step (4.1) could be re-written as:

S ρ (t n-1 , x i , y j , W n-1 i,j , Π Gρ (W ρ (t n , •, •))) = 0, ∀i, j with g(x i
) ≤ y j . However, the formulation (4.1) is more general and includes different kind of time-discretization like Euler implicit scheme or θ-methods.

The main idea of the numerical method described here above is to mix the use of a standard scheme for (3.5a), together with a "projection step" on ∂D in order to approximate the homogeneous oblique derivative boundary condition (3.5b). Let us point out that a similar method was introduced in [START_REF] Barles | Convergence of numerical schemes for parabolic equations arising in finance theory[END_REF] for the case g(x) ≡ |x|. However, the general case with possibly nonlinear function g requires some further attention on the boundary.

4.1.

A general convergence result. In this part we closely follow the arguments of Barles and Souganidis [START_REF] Barles | Convergence of approximation schemes for fully nonlinear second order equations[END_REF], using consistency, stability and monotonicity arguments. For this, we assume that the following hypotheses are considered: (H3) (stability) For any ρ, there exists a solution W ρ of the scheme which has linear growth (uniformly in ρ), that is, there exists a constant C > 0, independent of ρ, such that

|W ρ (t, x, y)| ≤ C(1 + |x| + |y|), on [0, T ] × D; (4.3)
(H4) (consistency) The scheme S ρ is consistent with respect to equation

(3.5a) in [0, T ] × D, that is, for all (t, x, y) ∈ [0, T ] × D and for every Φ ∈ C 1,2 ([0, T ] × D), lim ρ→0 [0,T ]×D∋(s,ξ,γ)→(t,x,y) ζ→0 S ρ (s, ξ, γ, Φ(s, ξ, γ)+ζ, Φ +ζ) = -∂ t Φ+H(t, x, D x Φ, D 2 x Φ).
(H5) (monotonicity) For any ρ, for r ∈ R, (t, x, y) ∈ [0, T ]×D, S ρ (t, x, y, r, Φ) depends only on the values of Φ in a neighborhood B η(ρ) (t, x, y) of

(t, x, y), with η(ρ) ≥ 0 such that η(ρ) ρ→0 -→ 0. Moreover, for any Φ 1 , Φ 2 functions on [0, T ] × R d × R → R, such that Φ 1 ≥ Φ 2 on B η(ρ) ((t, x, y)), it holds: S ρ (t, x, y, r, Φ 1 ) ≤ S ρ (t, x, y, r, Φ 2 ).
Let us point out that the monotonicity and the consistency are required here only for the operator scheme S ρ that corresponds to the discretization of the equation -∂ t ϑ + H(t, x, Dϑ, D 2 ϑ) = 0.

Notice also that the monotonicity assumption (H5) is slightly different from the classical one usualy required for general HJB equations. The reason for this comes from the fact that the operator scheme S ρ is defined on D and may use some values of the function W ρ outside the domain D that are obtained by oblique projection through the projection step. Our monotonicity assumption requires local dependency on the values of Φ in a neighborhood of the point under consideration, then a comparison of the scheme values is requested for any two functions Φ 1 and Φ 2 such that Φ 1 ≥ Φ 2 only on this neighborhood. The proof of Theorem 4.1 highlights what this requirement is needed for.

In the next section, we will check that assumptions (H3)-(H5) are well satisfied by the semi-Lagrangian scheme. Now, we state the convergence result for any general scheme satisfying (H3)-(H5). We start by proving that W is a viscosity sub-solution of equation (3.5) (the proof that W is a viscosity super-solution is analogous).

Let ϕ be in C 1,2 ([0, T ] × D) and let ( t, x, ȳ) be a local maximum point for W -ϕ on [0, T ] × D. Without loss of generality we can assume that ( t, x, ȳ)

is a strict local maximum in the restricted ball B r ( t, x, ȳ) ∩ ([0, T ] × D) for a certain r > 0, (W -ϕ)( t, x, ȳ) = max Br∩([0,T ]×D) (W -ϕ) = 0,
and ϕ ≥ 2 sup ρ W ρ ∞ outside the ball B r ( t, x, ȳ). We first assume that t ∈ (0, T ). Then we claim that:

-∂ t ϕ+H( t, x, D x ϕ, D 2 x ϕ) ≤ 0 if (x, ȳ) ∈ D, (4.4a) or min -∂ t ϕ+H( t, x, D x ϕ, D 2 x ϕ),-∂ y ϕ ≤ 0 if (x, ȳ) ∈ ∂D.(4.4b)
Following [START_REF] Crandall | User's guide to viscosity solutions of second order partial differential equations[END_REF], there exists a sequence of mesh steps ρ k and a sequence

(t k-1 , x k , y k ) in (0, T ) × D such that: ρ k →0 and (t k-1 , x k , y k )→( t, x, ȳ) as k →+∞, and (t k-1 , x k , y k ) is a global maximum point of W ρ k -ϕ, with (W ρ k -ϕ)(t k-1 , x k , y k ) = max [0,T ]×D (W ρ k -ϕ) = δ k k→∞ -→ 0 (4.5) and W ρ k (t k-1 , x k , y k ) k→∞ -→ W ( t, x, ȳ).
-Case 1: assume that ȳ > g(x). Thus, (x, ȳ) is in the open set D and for k large enough,

(x k , y k ) ∈ D. By continuity of g, y > g(x) in B η(ρ k ) (t k-1 , x k , y k ) for ρ k small enough. Therefore Π Gρ (W ρ k ) = W ρ k ≤ ϕ + δ k on B η(ρ k ) (t k-1 , x k , y k ). On the other hand, W ρ k (t k-1 , x k , y k ) = ϕ(t k-1 , x k , y k ) + δ k .
Hence, thanks to the monotonicity of the scheme, it follows

0 = S ρ k (t k-1 , x k , y k , W ρ k (t k-1 , x k , y k ), Π Gρ (W ρ k )) ≥ S ρ k (t k-1 , x k , y k , ϕ(t k-1 , x k , y k ) + δ k , ϕ + δ k ).
Using the consistency of the scheme, we obtain in the limit when ρ k → 0,

-∂ t ϕ + H( t, x, D x ϕ, D 2 x ϕ) ≤ 0. (4.6)
We conclude that (4.4a) is satisfied when (x, ȳ) ∈ D.

-Case 2: when ȳ = g(x), (x k , y k ) can be also on ∂D and the scheme may involve values W n kℓ on some points (x k , y ℓ ) that are not in D. Here, we need to consider two sub-cases.

-Sub-Case 2.1: if -∂ y ϕ( t, x, ȳ) ≤ 0, then (4.4b) holds.

-Sub-Case 2.2: if -∂ y ϕ( t, x, ȳ) > 0, then there exists a neighborhood V of ( t, x, ȳ) where ∂ y ϕ(t, x, y) is well defined, and -∂ y ϕ(t, x, y) > 0 for every (t, x, y) ∈ V. Therefore, y ≤ y ′ =⇒ ϕ(t, x, y) ≥ ϕ(t, x, y ′ ) ∀(t, x, y), (t, x, y ′ ) ∈ V.

(4.7)

For k large enough, B η(ρ k ) (t k-1 , x k , y k ) ⊂ V. Let (t, x, y) ∈ B η(ρ k ) (t k-1 , x k , y k ). If y ≥ g(x), then Π Gρ (W ρ k (t, x, y)) = W ρ k (t, x, y) ≤ ϕ(t, x, y) + δ k . (4.8)
Otherwise if y < g(x), Π Gρ (W ρ k (t, x, y)) = W ρ k (t, x, y jx ), and we have

Π Gρ (W ρ k (t, x, y)) = W ρ k (t, x, y jx ) ≤ ϕ(t,
x, y jx ) + δ k (by using (4.5)) (4.9) ≤ ϕ(t, x, y) + δ k (by using (4.7)) (4.10)

(for (4.10) we also used the continuity of g that ensures that for k big enough (t, x, y jx ) still lays in V). We conclude

Π Gρ (W ρ k ) ≤ ϕ + δ k on B η(ρ k ) (t k-1 , x k , y k ).
Thus by monotonicity (H5), we get

0 = S ρ k (t k-1 , x k , y k , W ρ k (t k-1 , x k , y k ), Π Gρ (W ρ k )) ≥ S ρ k (t k-1 , x k , y k , ϕ(t k-1 , x k , y k ) + δ k , ϕ + δ k ).
Using the consistency of the scheme, when ρ k → 0, we get

-∂ t ϕ + H( t, x, D x ϕ, D 2 x ϕ) ≤ 0. (4.11)
As conclusion for case 2, (4.4b) is satisfied when (x, ȳ) ∈ ∂D.

Similar arguments can be used in order to treat the case t = T (classical arguments enable to treat also the case of the boundary t = 0). We conclude that W is a viscosity sub-solution of the equation (3.5).

As already mentioned, we can also show that W is a super-solution. Then by the comparison principle on [0, T ]×D (Theorem 3.4 and using the growth condition (4.3)), the inequality W ≥ W holds on [0, T ] × D. Furthermore since the reverse inequality W ≤ W is always true, we deduce that W = W = ϑ on [0, T ] × D. Hence the convergence result is proved.

4.2.

A semi-Lagrangian scheme. We first introduce some new notations: for any u ∈ U , the drift b(•, •, u) (resp. volatility σ(•, •, u)) will be simply denoted by b u (resp. σ u ).

Consider the operator T :

C(R d+1 ) → C([0, T ] × R d+1
) defined by:

T (φ)(t, x, y) = min u∈U 1 2p 2p k=1 φ x x+hb u (t, x)+ hp(-1) k σ u ⌊ k+1
2 ⌋ (t, x), y (4.12) where (σ u k ) k=1,...,p are the column vectors of σ u , and ⌊q⌋ denotes the integer part of q. The notation

[•] ≡ [•]
x stands for a monotone, P 1 interpolation operator on the x-grid (x i ), satisfying for every Lipschitz continuous function φ (with Lipschitz constant L φ ):

       (i) [φ] x (x i ) = φ(x i ), ∀i, (ii) |[φ] x (x) -φ(x)| ≤ L φ |∆x|, (iii) |[φ] x (x) -φ(x)| ≤ C|∆x| 2 D 2 x φ ∞ if φ ∈ C 2 (R d ), (iv) for any functions φ 1 , φ 2 : R d → R, φ 1 ≤ φ 2 ⇒ [φ 1 ] x ≤ [φ 2 ] x . (4.13)
The operator T corresponds to a discretisation of the equation (3.5a) by a semi-Lagrangian scheme (see [START_REF] Menaldi | Some estimates for finite difference approximations[END_REF][START_REF] Camilli | An approximation scheme for the optimal control of diffusion processes[END_REF][START_REF] Debrabant | Semi-Lagrangian schemes for linear and fully non-linear diffusion equations[END_REF]). Now, define an approximation method for the system (3.5) as follows:

Fully discrete scheme (FDS) 1)Initialisation step: for all i, j, set W N i,j = ψ(x i , y j ∨ g(x i )). 2)For n = N, . . . , 1:

-Step 1 Compute W n-1 2 i,j = T (W n •,• )(t n , x i , y j ), for all i, j (here W n •,• denotes the values {W n i,j | i, j ∈ Z d × Z}); -Step 2 Compute W n-1 i,j = Π G ρ (W n-1 2
), for all i, j. The FDS scheme is a particular GS method, with the following formula for the operator S ρ :

S ρ (t -h, x, y, r, Φ) := 1 h r -Π G ρ (T (Φ)(t, x, y)) . (4.14)
Remark 4.2. We point out that in the FDS, for j such that y j < g(x i ), the projection step fixes the value of W n-1 i,j as W n-1 i,jx i . Of course, other alternatives could be considered, for instance, one could consider the scheme:

1) Initialization step: for n = N , for all i, j, set W N i,j = ψ(x i , y j ∨ g(x i )).

2) For n = N, . . . , 1, W n-1 i,j = T ( W n
•,• )(t n , x i , y j ), for all i, j, where

T (φ)(t, x, y) = min u∈U 1 2p 2p k=1 φ x,y x+hb u (t,x)+ hp(-1) k σ u ⌊ k+1 2 ⌋ (t,x), y ∨ g(x) ,
with [•] x,y standing for an interpolation in both variables x and y.

All the convergence results stated in the paper also hold for this new scheme.

In the next subsection, we will study some properties of the approximated solution W ρ . Before this, we define also the following semi-discrete scheme where we consider only a time-discretization:

Semi-discrete scheme (SDS) 1)For n = N , for every (x, y) ∈ R d+1 , set V N (x, y) = ψ(x, y ∨ g(x)).

2)For n = N, . . . , 1, define V n-1 as the function defined by:

V n-1 (x, y) = T 0 (V n )(t n , x, y ∨ g(x)), (4.15)
where, T 0 is defined from

C(R d+1 ) into C([0, T ] × R d+1
) by:

T 0 (φ)(t, x, y) := min u∈U 1 2p 2p k=1 φ x + hb u (t, x) + (-1) k hpσ u ⌊ k+1 2 ⌋ (t, x), y .
Unlike the fully discrete scheme, no interpolation step is required in the SDS. Straightforward calculations lead to the following consistency estimate, for any φ ∈ C 2,4 ((0, T ) × R d × R):

1 h (φ(t -h, x, y) -T 0 (φ)(t, x, y)) --∂ t φ + H(t, x, D x φ, D 2 x φ) (4.16) ≤ K 1 max u∈U ∂ 2 tt φ ∞ + |b u (t, x)| 2 D 2 x φ ∞ + |b u (t, x)||σ u (t, x)| 2 D 3 x φ ∞ +|σ u (t, x)| 4 D 4 x φ ∞ h,
for any (x, y) ∈ R d , with K 1 a positive constant independent of φ. Moreover, T and T 0 satisfy the following properties :

Lemma 4.3. For every φ 1 , φ 2 ∈ C([0, T ] × R d × R), we have: (i) (T 0 (φ 1 ) -T 0 (φ 2 )) + ≤ (φ 1 -φ 2 ) + ; (ii) T φ 1 -T φ 2 ∞ ≤ φ 1 -φ 2 ∞ ; (iii) T 0 φ 1 -T 0 φ 2 ∞ ≤ φ 1 -φ 2 ∞ ; (iv) for any Lipschitz continuous function φ on R d × R, T φ -T 0 φ ∞ ≤ L φ ∆x, where L φ is a Lipschitz constant of φ.
Notice that assertion (i) corresponds to a discrete comparison principle. The proof of the lemma is straightforward and will be omitted.

4.3.

Convergence of the SL scheme. The aim of this subsection is to prove a convergence result for the scheme FDS. Theorem 4.4. Let assumptions (H1)-(H2) be satisfied. Let W ρ be defined by the fully discrete scheme. Assume also that

|(∆x,∆y)| h ≤ C (where C is a positive constant) (4.17)
Then the FDS is convergent, i.e. W ρ converges towards ϑ, as ρ tends to 0, uniformly on each compact set of [0, T ] × D.

The claim of this theorem would be a consequence of the general result of theorem 4.1, provided the monotonicity, consistency and stability properties of the FDS are established. The monotonicity and the consistency are easy to check. If ψ is a bounded function, the stability of the scheme would be immediate because in this case W n ∞ ≤ ψ ∞ . However, in the general case the function ψ is not necessarily bounded. In order to derive the stability property, we first prove some results on the semi-discrete scheme that we state in the following next two propositions. Proposition 4.5. If assumptions (H1) and (H2) are satisfied, the solution V of the semi-discrete scheme is Lipschitz continuous in (x, y) and 1 2 -Hölder continuous in t: there exists a constant L V ≥ 0, for every 0 ≤ n, m ≤ N , and every (x 1 , y 1 ), (x 2 , y 2 ) ∈ R d+1 such that

|V n (x 2 , y 2 ) -V m (x 1 , y 1 )| ≤ L V |x 2 -x 1 | + |y 2 -y 1 | + (1 + |x 1 |)|t n -t m | 1/2 .
Proof. Let (ξ n , ξ n+1 , . . . , ξ N -1 ), be independent random variables, with values in {0, . . . , 2p}, and such that P(ξ n = k) = 1 2p , ∀k = 1, . . . , 2p. Let u = (u n , u n+1 , . . . , u N -1 ) denotes a sequence of controls, with u i ∈ U , and let Z n k = Z n,x,u k , n ≤ k ≤ N , be the random walk such that:

Z n n = x, Z n k+1 = Z n k + hb u k (t k+1 , Z n k ) + √ hσ u k ξ k (t k+1 , Z n k ), k ≥ n (4.18)
where the notation σu

k (t, x) := (-1) k √ pσ u ⌊ k+1
2 ⌋ (t, x) is used. Notice that Z n,x,u k depends only on the controls (u n , . . . , u k-1 ). Direct calculations lead to the following expressions:

V n (x, y) = min un∈U 1 2p 2p k=1 V n+1 x + hb un (t n+1 , x) + √ hσ un k (t n , x), y ∨ g(x) = min un∈U E V n+1 Z n,x,u n+1 , y ∨ g(x) ,
and, in the same way,

V n (x, y) = min un∈U E min u n+1 ∈U E V n+2 (Z n+1,Z n,x,u n+1 ,u n+2 , y ∨ g(x) ∨ g(Z n,x,u n+1 )) , = min un∈U E min u n+1 ∈U E V n+2 (Z n,x,u n+2 , y ∨ g(x) ∨ g(Z n,x,u n+1 )) ,
and so on, and finally (also using V N (x, y) ≡ ψ(x, y ∨ g(x))):

V n (x, y) =min un E min u n+1 E • • • min u N -1 E V N Z n,x,u N , y∨ max i=n,...,N -1 g(Z n,x,u i ) • • • =min un E min u n+1 E • • • min u N -1 E ψ Z n,x,u N , y ∨ max i=n,...,N g(Z n,x,u i ) • • • .
Then, we have the equivalent representation formula:

V n (x, y) = min u=(un,...,u N -1 ) E ψ(Z n,x,u N , y ∨ max i=n,...,N g(Z n,x,u i )) . (4.19) By using |(x 2 ∨ y 2 ) -(x 1 ∨ y 1 )| ≤ |x 1 -x 2 | ∨ |y 1 -y 2 | ≤ |x 1 -x 2 | + |y 1 -y 2 |,
we obtain:

|V n (x 2 , y 2 ) -V n (x 1 , y 1 )| ≤ ≤ L ψ max u E |Z n,x 2 ,u N -Z n,x 1 ,u N | + L g max n≤i≤N |Z n,x 2 ,u i -Z n,x 1 ,u i |+L g |y 2 -y 1 | ≤ L ψ (L g + 1) max u E max n≤i≤N |Z n,x 2 ,u i -Z n,x 1 ,u i | + |y 2 -y 1 | , (4.20)
where max u denotes the maximum over (u n , u n+1 , . . . , u N -1 ) ∈ U N -n . By (H1), there exists a constant C 1 ≥ 0, that depends only on (T, L σ , L b ), such that for every 0 ≤ m, n ≤ N and every x, x 1 , x 2 ∈ R d , the following classical estimates hold:

E max n≤i≤N Z n,x 2 ,u i -Z n,x 1 ,u i ≤ C 1 |x 2 -x 1 | (4.21a) E max 0≤i≤p Z n,x,u m+i -Z n,x,u m ≤ C 1 |t m+p -t m | 1/2 (1 + |x|). (4.21b)
Combining (4.20) and (4.21a), for n = m, the Lipschitz property follows:

|V n (x 1 , y 1 ) -V n (x 2 , y 2 )| ≤ L ψ (L g + 1)(C 1 |x 2 -x 1 | + |y 2 -y 1 |).
On the other hand, for 0 ≤ n ≤ m ≤ N by using again (4.19) we have:

|V m (x, y) -V n (x, y)| ≤ C max u E max n≤i≤N Z n,x,u i -max m≤i≤N Z m,x,u i . (4.22)
for some constant C. Therefore, since Z m,x,u m = x:

|V m (x, y) -V n (x, y)| ≤ C max u E max n≤i≤m Z n,x,u i -x (4.23) +C max u E max m+1≤i≤N Z n,x,u i -max m+1≤i≤N Z m,x,u i .(4.24)
The right term of (4.23) is bounded by

C(1 + |x|)|t m -t n | 1/2
. In order to bound (4.24), by using that Z n,x,u i = Z m,Z n,x,u m ,u i (∀i ≥ m) and (4.21a) the following estimate is obtained:

E max m+1≤i≤N Z m,Z n,x,u m ,u i -max m+1≤i≤N Z m,x,u i ≤ C 1 E Z n,x,u m -x ≤ C 2 1 (1 + |x|)|t m -t n | 1/2 .
Hence it holds, for some constant C,

|V n (x, y) -V m (x, y)| ≤ C(1 + |x|)|t m -t n | 1/2 .
Together with the Lipschitz property the desired result follows.

Remark 4.6. It is not clear whether the solution W ρ obtained by the fully discrete scheme satisfies a Lipschitz continuity property or not. The main difficulty is that a representation formula as (4.19) is not guaranteed for W ρ . Proposition 4.7. If assumptions (H1) and (H2) are satisfied, there exists a constant C > 0 independent of ρ such that, ∀n = 0, . . . , N :

W n -V n ∞ ≤ C |(∆x, ∆y)| h . ( 4

.25)

Proof. Consider the operator Π :

C(R d × R) → C(R d × R) defined by Πφ(x, y) = φ(x, y ∨ g(x)
). With this notation, we have

V n-1 = Π(T 0 (V n )). On the other hand, W n-1 ij = Π G ρ (T (W n )(t n , x i , y j ))
. Therefore, by using Lemma 4.3, we get:

W n-1 -V n-1 ∞ ≤ Π G ρ (T (W n ))-Π(T 0 (V n )) ∞ ≤ Π G ρ (T (W n ))-Π G ρ (T 0 (V n )) ∞ + Π G ρ (T 0 (V n )) -Π(T 0 (V n )) ∞ ≤ T (W n )-T 0 (V n ) ∞ + Π G ρ (T 0 (V n ))-Π(T 0 (V n )) ∞ ≤ T (W n )-T (V n ) ∞ + T (V n )-T 0 (V n ) ∞ + L V ∆y ≤ W n -V n ∞ + L V |∆x| + L V ∆y.

By recursion, it follows:

W n -V n ∞ ≤ W N -V N ∞ + nL V (|∆x| + ∆y). Since V N (x, y) = ψ(x, y ∨ g(x)) on R d+1 , W N is
an interpolation of V N on the grid G ρ and n ≤ T h , we obtain the desired result. Now we turn to the proof of Theorem 4.4.

Proof of Theorem 4.4. Since V n is a Lipschitz function in the (x, y) variables, uniformly in n, then V n has a linear growth with constants independent of n. By Lemma 4.7, and by using assumption (4.17),

W n -V n ∞ ≤ C |(∆x, ∆y)| h ≤ K.
Hence W n has also a linear growth, with constants independent of n.

Moreover, the consistency of the scheme comes from the estimate (4.16) and using that |∆x| 2 h → 0 by assumption (4.17). Also, the monotonicity property (as defined in (H5)) is deduced from the monotonicity of the interpolations and the monotonicity of the operator T 0 .

In consequence, the convergence result follows from Theorem 4.1.

Remark 4.8. The convergence result of Theorem 4.4 is still valid in the case where ψ is a bounded Lipschitz function and |∆x| 2 h → 0.

Error bounds for the semi-Lagrangian scheme

In the sequel, for any function f defined from some set Q ⊂ (0, ∞) × R d+1 into either R, R d or some space of matrices, we set

[f ] 1 := sup (t,x) =(s,y) |f (t, x) -f (s, y)| (|x -y| + |t -s| 1/2 )
,

and |f | 1 := f ∞ + [f ] 1 .
In this section we aim to give an error bound for the fully discrete scheme FDS, as well as the semi-discrete SL scheme (4.15). For this aim a supplementary hypothesis will be required (together with (H1) and (H2)): (H1') There is a constant M 1 independent of u ∈ U such that:

|σ u | 1 + |b u | 1 ≤ M 1 .
The error estimates theory is based on some technique of "shaking coefficients" and regularization introduced by Krylov in [START_REF] Krylov | On the rate of convergence of finite-difference approximations for Bellman's equations[END_REF][START_REF] Krylov | On the rate of convergence of finite-difference approximations for Bellman's equations with variable coefficients[END_REF] and studied later by many authors [START_REF] Barles | On the convergence rate of approximation schemes for Hamilton-Jacobi-Bellman equations[END_REF][START_REF] Barles | Error bounds for monotone approximation schemes for Hamilton-Jacobi-Bellman equations[END_REF][START_REF] Barles | Error bounds for monotone approximation schemes for parabolic Hamilton-Jacobi-Bellman equations[END_REF][START_REF] Bonnans | Error estimates for stochastic differential games: the adverse stopping case[END_REF][START_REF] Bonnans | Error bounds for a Stochastic Impulse Control Problem[END_REF][START_REF] Debrabant | Semi-Lagrangian schemes for linear and fully non-linear diffusion equations[END_REF]. The main idea consists of regularizing the exact solution ϑ in order to obtain a smooth sub-solution ϑ ε (for an approximation parameter ε > 0) to the equation (3.5), then the upper-bound error estimate can be obtained by using the consistency estimate (4.16). The regularisation procedure takes advantage from the fact that the exact solution ϑ is Hölder continuous, which enables to obtain an estimate of |ϑ -ϑ ε | ∞ of order O(ε). The lower bound is obtained by similar arguments, but in this case we need to build a smooth sub-solution to the discrete equation. For this, a regularization of the numerical solution will be considered. However, as noticed earlier (remark 4.6), it is not clear if the solution W ρ of the fully discrete scheme is Hölder continuous or not, and then it is not clear how to build directly a regularization W ρ ε (with regularisation parameter ε) with an error |W ρ -W ρ ε | of order O(ε). For this reason, we will start by deriving an error estimate between the solution ϑ and the solution V of the semi-discrete scheme SDS, using the fact that V has Hölder continuity properties. Then the estimates for the FDS are derived as a consequence of Lemma 4.7.

For the regularization procedure, consider a smooth function µ : R d+2 → R, with µ ≥ 0, supported in (0, 1) × B 1 (0), with R R d µ(s, x) dxds = 1, and define µ ε as the following sequence of mollifiers:

µ ε (t, x, y) := 1 ε d+3 µ( t ε 2 , x ε , y ε ) in R × R d+1 .
5.1. The rate of convergence for the semi-discrete scheme. For any ε > 0, let E be the set of progressively measurable processes (α, χ) valued in

[-ε 2 , 0] × B(0, ε) ⊂ R × R d that is, E := prog. meas. process (α, χ) valued in E ,
where

E := (a, e) ∈ R × R d , -ε 2 ≤ a ≤ 0, |e| ≤ ε .
On other hand, let M := 2 1 + L 2 g and introduce g ε : R d → R defined by

g ε (x) := g(x) -M ε.
(

Finally, let us denote by D ε the set defined as follows:

D ε := (x, y) ∈ R d+1 , y > g ε (x) .
Remark 5.1. The choice of g ε is such that the following property holds:

(x, y) ∈ D =⇒ (x -e 1 , y -e 2 ) ∈ D ε ∀ (e 1 , e 2 ) ∈ R d+1 , |(e 1 , e 2 )| ≤ ε.
5.1.1. Upper bound. Now, we start by introducing a perturbed control problem (with "shaking coefficients"). For any ε > 0, consider the following value function

ϑ ε (t, x, y) := inf u∈U , (α,χ)∈E E ψ X u,(α,χ) t,x (T ), max s∈[t,T ] g ε (X u,(α,χ) t,x (s)) ∨ y , (5.2)
where X u,(α,χ) t,x (•) is the solution of the perturbed system of SDEs

dX(s) = b(s + α(s),X(s) + χ(s),u(s))ds+σ(s + α(s),X(s)+χ(s),u(s))dB(s) X(t) = x.
(5.3)

Remark 5.2. The functions σ and b are only defined for times t ∈ [0, T ], but they can be extended to times [-2ε 2 , T + 2ε 2 ] in such a way that assumption (H1') still holds.

Proposition 5.3. Under assumptions (H1)-(H2) and (H1') (extending eventually b and σ as prescribed in Remark 5.2), the following holds:

(i) ϑ ε is a Lipschitz continuous function in x and y and a 1 2 -Hölder continuous function in t. More precisely, for t, ≤ t ′ ∈ [0, T ], x, x ′ ∈ R d and y, y ′ ∈ R, we have:

|ϑ ε (t, x, y) -ϑ ε (t ′ , x ′ , y ′ )| ≤ L ϑ (|x -x ′ | + |y -y ′ | + (1 + |x|)|t -t ′ | 1 2 ). (ii) |ϑ(t, x, y) -ϑ ε (t, x, y)| ≤ Cε on [0, T ]×D,
where the constant C only depends on T and the Lipschitz constants of b, σ, g and ψ.

Proof. The Lipschitz and Hölder continuity follows by the same arguments as in Proposition 2.2. Thanks to the Lipschitz continuity of ψ and g and the choice of g ε one has

|ϑ(t, x, y) -ϑ ε (t, x, y)| ≤ K sup u∈U , (α,χ)∈E E sup s∈[t,T ] X u t,x (s) -X u,(α,χ) t,x (s)
ϑ(t n , x, y) -V n (x, y) ≤ C h 1 4
, for all n ≥ 0 (with nh ≤ T ) and (x, y) ∈ D.

Proof. The proof is splitted into three steps:

-Step 1. Let (ã, (ẽ 1 , ẽ2 )) ∈ R × R d+1 be such that -ε 2 ≤ ã ≤ 0 and |(ẽ 1 , ẽ2 )| ≤ ε. Let ( t, x, ȳ) ∈ (0, T ) × D, and define for ν > 0:

I ν ( t, x, ȳ) := (t, x, y)|t ∈ [ t -ν 2 , t], (x, y) ∈ B ν (x, ȳ) . We claim that ϑ ε (• -ã, • -ẽ1 , • -ẽ2 ) is a viscosity sub-solution of -∂ t ϑ + H(t, x, D x ϑ, D 2 x ϑ) ≤ 0 on I ε ( t, x, ȳ). (5.4)
To prove this claim, we notice first that I ε ( t, x, ȳ) ⊂ D ε . Moreover, from Theorem 3.2, for any ε > 0, ϑ ε is a viscosity solution of the following equation

     -∂ t ϑ ε + sup (a,e)∈E H(t + a, x + e, D x ϑ ε , D 2 x ϑ ε ) = 0 (-2ε 2 , T ) × D ε , -∂ y ϑ ε = 0 (-2ε 2 , T ) × ∂D ε , ϑ ε (T, x, y) = ψ(x, (g(x) -M ε) ∨ y). (5.5) Let ϕ ∈ C 2,4 ([-2ε 2 , T ] × D ε ) be such that ϑ ε (• -ã, • -ẽ1 , • -ẽ2 ) -ϕ achieves a local maximum at ( t, x, ỹ) on I ε ( t, x, ȳ). Clearly ( t -ã, x -ẽ1 , ỹ -ẽ2 ) is also a local maximum of ϑ ε -ϕ(• + ã, • + ẽ1 , • + ẽ2 ) on I 2ε (and I 2ε ⊂ D ε ).
Since ϑ ε is a viscosity solution of equation (5.5), we obtain:

-∂ t ϕ( t, x, ỹ) + sup (a,e 1 )∈E H( t -ã + a, x -ẽ1 + e 1 , D x ϕ, D 2 x ϕ) ≤ 0.
Taking (a, e 1 ) = (ã, ẽ1 ), we get the result.

-Step 2. Define following mollification ϑ ε :

ϑ ε (t, x, y) := (ϑ ε * µ ε )(t, x, y) = |(e 1 ,e 2 )|≤ε -ε 2 ≤a≤0
ϑ ε (t-a, x-e 1 , y-e 2 )µ(a, e 1 , e 2 )dade.

We recall the following properties of the mollifiers:

|ϑ ε (t, x, y) -ϑ ε (t, x, y)| ≤ [ϑ ε ] 1 ε, (5.6)
with

[w ǫ ] 1 ≤ CL ϑ , C ≥ 0. Moreover, for any i ≥ 1 or j ≥ 1, D i x ϑ ε ∞ ≤ CL ϑ ε 1-i , D j t ϑ ε ∞ ≤ CL ϑ ε 1-2j (5.7)
(where D i x denotes the i-th derivative with respect to x, and D j t the j-th derivative with respect to t). Since ϑ ε is a limit of convex combinations of ϑ ε (• -ã, • -ẽ1 , • -ẽ2 ), then ϑ ε satisfies in the viscosity sense

-∂ t ϑ ε + H(t, x, D x ϑ ε , D 2 x ϑ ε ) ≤ 0 in (0, T ) × D. (5.8) Taking into account that ϑ ε is in C ∞ ([0, T ] × R d+1
), we conclude that (5.8) holds in classical sense on [0, T ] × D. From the consistency estimate (4.16) along with (5.8) and by assumption (H1'), we get:

ϑ ε (t n-1 , x, y) -T 0 (ϑ ε )(t n , x, y) ≤ C h 2 ε 3 .
Combining these bounds with Lemma 4.3, we get

ϑ ε (t n-1 , x, y) -V n-1 (x, y) ≤ T 0 (ϑ ε )(t n , x, y) -T 0 (V n )(x, y) + C h 2 ε 3 ≤ ϑ ε (t n , •) -V n + ∞ + C h 2 ε 3 .
Therefore, by recursion, it comes

(ϑ ε (t n , •) -V n ) + ∞ ≤ (ϑ ε (0, •) -V N )) + ∞ + CT h ε 3 .
(5.9) -Step 3. By Proposition 5.3 together with the inequalities (5.6) and (5.9), we obtain for n ≥ 0,

(ϑ(t n , •) -V n ) + ∞ ≤ Cε + CT h ε 3 , n ≥ 0.
The choice ε 4 = h leads to (for n ≥ 0):

(ϑ(t n , •) -V n ) + ∞ ≤ C h 1 4 ,
which concludes the proof.

5.1.2. Lower bound. For obtaining the lower bound we will apply exactly the same techniques as used for the upper bound, reversing the role of the equation and the scheme. The key point is that the solution V of the semidiscrete scheme SDS is Lipschitz continuous. Then it is possible to use the techniques of "shaking coefficients" and regularisation to build a smooth sub-solution

V ε satisying V n ε -V n ∞ ≤ Cε and V n-1 ε (x, y) -T 0 (V n ε )(t n , x, y ∨ g(x)) ≤ 0 in [0, T ] × R d × R.
Then by consistency estimate and comparison principle, we conclude the following result:

Theorem 5.5. Let assumptions (H1)-(H2) and (H1') be satisfied. There exists a constant C > 0 such that

ϑ(t n , x, y) -V n (x, y) ≥ -C h 1 4
for all n ≥ 0 (with N h ≤ T ) and (x, y) ∈ D.

5.2.

The fully discrete scheme. The section is concluded with the following theorem that provides error estimates for the fully discrete semi-Lagrangian scheme. The result is a simple consequence of Theorems 5.4, 5.5 and Lemma 4.7.

Theorem 5.6. Let assumptions (H1)-(H2) and (H1') be satisfied. If W is the solution of the fully-discrete scheme (4.14) , there exists a constant C > 0 independent from n = 0, . . . , N such that

W n -ϑ(t n , •) ∞ ≤ C h 1 4 + |(∆x, ∆y)| h .
Proof. The result follows by using

W n -ϑ(t n , •) ∞ ≤ W n -V n ∞ + V n -ϑ(t n ,
•) ∞ . and thanks to Theorems 5.4, 5.5 and to Lemma 4.7.

Numerical tests

In this section we present some numerical results for reachability problems in presence of state constraints (as in Section 3). The dynamics is given by the following controlled SDE in R 2 :

   dX(s) = u(s) 1 0 ds + u(s)σ(X(s))dB s , s ≥ t, X(t) = x (6.1)
where B is a one-dimensional Brownian motion (p = 1), U = [0, 1] ⊂ R and the function σ(x) ∈ R 2 will vary depending on the example. The target set is

C = x ≡ (x 1 , x 2 ) ∈ R 2 , 0 ≤ x 1 ≤ 0.4, |x 2 | ≤ 0.
5 and the set of state-constraints is

K = R 2 \ x ≡ (x 1 , x 2 ) ∈ R 2 , -0.4 < x 1 < -0.2, |x 2 | < 0.1 .
Given a final time T , the aim is to approximate the set R C,K t defined by (2.7), for t = 0. We recall that the following auxiliary value is introduced:

ϑ(t, x, y) = inf u∈U E d C (X u t,x (T )) ∨ max s∈[t,T ] d K (X u t,x (s)) ∨ y
and that the characterization of the backward reachable set is given by

x ∈ R C,K t ⇐⇒ ϑ(t, x, 0) = 0. (6.2)
In all the following tests, t = 0, T = 0.5 are fixed and the computational domain for (x 1 , x 2 , y) is

(x 1 , x 2 , y) ∈ Ω = [-1, 1] × [-1, 1] × [0, 1].
The numerical scheme implemented is the semi-Lagrangian scheme (FDS). We denote by N t the number of time-steps, N x 1 = N x 2 and N y are the number of mesh steps for the space variables (x 1 , x 2 , y), and

h := T N t , ∆x i = 2 N x i , ∆y = 1 N y .
Different simulations show that the results are not very sensitive to the step discretization of the variable y. Unless otherwise precised, we set N y = 10. For all tests, since u ∈ [0, 1] and because the dynamics depends linearly on the control u it is sufficient to take only two controls (N u = 2), e.g., u ∈ {0, 1} for the discretization of the control variable.

The different figures (see e.g. Figure 1) represent points in the (x 1 , x 2 ) plane. The obstacle is represented in black, and the target in dark grey. Then, an arbitrary non negative threshold ǫ := 10 -5 is choosen, and our approximation of the reachable set is given by

R C,K t := (x 1 , x 2 ) ∈ R 2 , ϑ(0, x 1 , x 2 , ǫ) ≤ ǫ (6.3)
(plotted in green).

Remark 6.1. This procedure is not very sensitive to the choice of ǫ in the range 10 -2 -10 -6 .

Remark 6.2. We have also tested the scheme of Remark 4.2 and obtain similar results.

Test 1. In this test, there is no diffusion, that is: Test 2. In this test a non-zero volatility is considered:

σ(x) := 0 0 .
σ(x) := 0 5
In that case the backward reachable set reduces to the target set, see Figure 1(b). In fact for any point outside the target, as soon as u = 0, even if the drift steers the system in the direction of the target, there is always a non-zero probability to go to far way in the x 2 orthogonal direction and therefore to not reach it. for three different meshes. It also shows how the scheme converges. Notice that the points which are immediately above or below the obstacle are not in the reachable set since in presence of diffusion the state-constraint will be violated with a non-zero probability. Also, in Table 1, various error norms are given to study the convergence of the scheme. For a given n ≥ 1 we have chosen

N x 1 = N x 2 = n, N t = n/4 and N y = n/4.
Errors have been computed by taking n = 160 for the reference value, and a convergence of order greater than one is observed in this simple example (better than the expected order of convergence). In Figure 3 we have plotted the results obtained with three different meshes, using (N x 1 , N x 2 , N t ) = (n, n, n/4) for n ∈ {100, 200, 400}. Altough the first computation plotted in Figure 3 (left, with n = 100) is not very accurate, the other computations (with n = 200 and n = 400) clearly show the good convergence of the scheme. where D, ν :

R n → R n , H : [0, T ] × R n × R n × S n (S n
is the space of n × n symmetric matrices), and w 0 : R n → R satisfy the following properties: In the sequel, it will be usefull to use the concept of semijets: Proof. The result can be obtained adapting the arguments in [START_REF] Dupuis | On oblique derivative problems for fully nonlinear secondorder elliptic partial differential equations on nonsmooth domains[END_REF], thanks to the local compactness of the set D.

(P1) D is a locally compact subset of R n such that there exists b > 0 such that 0≤t≤b B tb (x + tν(x)) ⊂ D c , ∀x ∈ ∂D, 0≤t≤b B tb (x -tν(x)) ⊂ D, ∀x ∈ ∂D; (P2) H ∈ C(R × D × R n × S n ) and there is a neighborhood U of ∂D in D and a function ω 1 ∈ C([0, ∞]) satisfying ω 1 (0) = 0 such that ∀t ∈ [0, T ), x ∈ U, p, q ∈ R n , X, Y ∈ S d |H(t, x, p, X) -H(t, x, q, Y )| ≤ ω 1 (|p -q| + X -Y ); (P3) There is a function ω 2 ∈ C([0, ∞]) satisfying ω 2 (0) = 0 such that H(t, y, p, -Y ) -H(t, x, p, X) ≤ ω 2 (α |x -y| 2 + |x -y| (|p| + 1)) ∀α ≥ 1, t ∈ [0, T ), x, y ∈ D, p ∈ R n , X, Y ∈ S n such that -α I 0 0 I ≤ X 0 0 Y ≤ α I -I -I I . (P4) ν is C 2 (R n ) ∩ W 2,∞ (R n ) with |ν(x)| = 1,
Proof of Theorem A.2. We will prove the theorem for u and v sub-and super-solution of (A.1) with boundary condition respectively replaced by ν(x), Du + α and ν(x), Dv -α for a certain α > 0. It means that on (0, T ) × ∂D, u and v satisfy min ( ν(x), p + α , a + H(t, x, p, X)) ≤ 0 ∀(a, p, X) ∈ P 1,2,+ u(t, x) (A.9)

max ( ν(y), q -α , b + H(t, y, q, Y )) ≥ 0 ∀(b, q, Y ) ∈ P 1,2,-v(t, y). (A.10)
We observe that it is always possible to consider a sub-solution u γ such that

lim t→T u γ (t, x) = -∞ ∂ t (u γ ) + H(t, x, Du γ , D 2 u γ ) < 0 defining, for instance, u γ (t, x) := u(t, x) -γ T -t .
The desired comparison result is then obtained as a limit for γ → 0. Given δ > 0, let us define ρ δ := sup

(t,x)∈[0,T )×D u(t, x) -v(t, x) -2δ(1 + |x| 2 ) .
The growth condition on u and v implies that there exists (s, z)

∈ [0, T ) × D such that ρ δ = u(s, z) -v(s, z) -2δ(1 + |z| 2 ). Set ρ := lim δ→0 ρ δ . Assume ρ ≤ 0. Since u(t, x) -v(t, x) ≤ 2δ(1 + |x| 2 ) + ρ δ , for every (t, x) ∈ [0, T ) × D, we have that u(t, x) ≤ v(t, x) on [0, T ) × D.
Now, assume that ρ > 0 and let us show a contradiction. Thanks to (P1) and (P4), if z ∈ ∂D, then there exists δ ∈ (0, 1 2

) such that B tδ (x+tν(z)) ⊂ D c for x ∈ B δ (z) ∩ ∂D, t ∈ (0, 2δ]. Set Γ := t>0 B tδ (tν(z)) • . It comes that y -x / ∈ Γ, if x ∈ ∂D ∩ B δ (z) • , y ∈ D ∩ B δ (z) • . (A.11)
In what follows we will restrict our attention to the events on the set

B δ (z) ∩ D, then we can assume ν ∈ W 2,∞ (R n , R n ) and |ν(x)| = 1.
Thanks to Lemma A.3 we can define

Φ(t, x, y) := u(t, x) -v(t, y) -α|x -z| 2 -ϑ ε (x, y) -δ(1 + |x| 2 ) -δ(1 + |y| 2 ).
and thanks to the growth conditions and the semicontinuity of u and v we can state that there exists ( t, x, ȳ) ∈ [0, T ) × D maximum point for Φ. Thanks to (A.2) and (A.3) the following inequalities hold ρ δ -ε ≤ Φ(s, z, z) ≤ Φ( t, x, ȳ) (A.12)

≤ u( t, x) -v( t, ȳ) -α|x -z| 2 -θ |x -ȳ| 2 ε -δ(1 + |x| 2 ) -δ(1 + |ȳ| 2 )
and by classical arguments (see [START_REF] Crandall | User's guide to viscosity solutions of second order partial differential equations[END_REF]), extracting a subsequence if necessary, we can conclude that |x -ȳ| 2 ε ε→0 -→ 0

x, ȳ ε→0 -→ z.

If t = 0: ρ δ -ε ≤ Φ(0, x, ȳ) ≤ w 0 (x)-w 0 (ȳ)-α|x-z| 2 -θ |x -ȳ| 2 ε -δ(1+|x| 2 )-δ(1+|ȳ| 2 ) so since the right-hand term in the inequality tends to zero for ε, δ → 0 and lim δ,ε→0

(ρ δ -ε) > 0, one can assume that t > 0 for ε and δ small enough.

Let us define ũ(t, x) := u(t, x) -α|x -z| 2 -δ(1 + |x| 2 ) ṽ(t, x) := v(t, x) + δ(1 + |x| 2 ).

Since ( t, x, ȳ) is a maximum point for Φ and property (A.8) holds, for x and y sufficiently close to x and ȳ, we have the following inequalities ũ(t, x) -ṽ(t, y)

≤ ũ( t, x) -ṽ( t, ȳ) + ϑ ε (x, y) -ϑ ε (x, ȳ)

≤ ũ( t, x) -ṽ( t, ȳ) + p • (x -x) + q • (y -ȳ)+ (a, D x ϑ ε (x, ȳ), X) ∈ P 1,2,+ ũ( t, x) (-b, -D y ϑ ε (x, ȳ), -Ỹ ) ∈ P 1,2,-ṽ( t, ȳ).

Recalling the definition of ũ and ṽ one has so that the desired contradiction is obtained by taking the limit for ε, α → 0.

In order to prove (A. [START_REF] Bonnans | Consistency of generalized finite difference schemes for the stochastic HJB equation[END_REF]) and (A.17) we proceed as in [START_REF] Dupuis | On oblique derivative problems for fully nonlinear secondorder elliptic partial differential equations on nonsmooth domains[END_REF]. In fact for ε and δ sufficiently small If (a, p, X) ∈ P 1,2,+ u α (t, x) and (b, q, Y ) ∈ P 1,2,-v α (t, x), by property (P2) one has a + H(t, x, p, X) ≤ a + C + H(t, x, p + αDh, X + αD 2 h) -C + ω 1 (αM ), (A. [START_REF] Bouchard | Weak dynamic programming principle for viscosity solutions[END_REF] b + H(t, x, q, Y ) ≥ b -C + H(t, x, q -αDh, Y -αD 2 h) + C -ω 1 (αM ) (A. Observing that P 1,2,+ u(t, x) = P 1,2,+ u α (t, x) + C, αDh, αD 2 h P 1,2,-v(t, x) = P 1,2,-v α (t, x) -C, αDh, αD 2 h one has a + C, p + αDh(x), X + αD 2 h(x) ∈ P 1,2,+ u(x) b -C, q -αDh(x), Y -αD 2 h(x) ∈ P 1,2,-v(x) then, by the definition of viscosity sub-and super-solution a + C + H(t, x, p + αDh, X + αD 2 h) ≤ 0 on D min ν(x), p + αDh(x) , a + C + H(t, x, p + αDh, X + αD 2 h) ≤ 0 on ∂D and b -C + H(t, x, q -αDh, Y -αD 2 h) ≥ 0 on D max ν(x), p -αDh(x) , b -C + H(t, x, q -αDh, Y -αD 2 h) ≥ 0 on ∂D.

For α small enough, taking C = w 1 (αM ), we can finally conclude by inequalities (A.18),(A.20) and (A. [START_REF] Camilli | An approximation scheme for the optimal control of diffusion processes[END_REF]),(A.21)) that a + H(t, x, p, X) ≤ 0 on D min ( ν(x), p + α, a + H(t, x, p, X)) ≤ 0 on ∂D and +b + H(t, x, q, Y ) ≥ 0 on D max ( ν(x), p -α, b + H(t, x, q, Y )) ≥ 0 on ∂D. In other words for α small enough u α and v α are respectively sub-and supersolution of (A.1) with boundary conditions ν(x), Du α +α and ν(x), Dv α -α and since u α α→0 → u and v α α→0 → v, we can prove the comparison theorem for them instead of u and v.

  H1) b : [0, T ] × R d × U → R d and σ : [0, T ] × R d × U → R d×p are continuous functions. Moreover, there exist L b , L σ > 0 such that for every x, y ∈ R d , t ∈ [0, T ], u ∈ U : |b(t, x, u) -b(t, y, u)| ≤ L b |x -y|, |σ(t, x, u) -σ(t, y, u)| ≤ L σ |x -y|.

Theorem 4 . 1 .

 41 Assume (H1)-(H2). Let S ρ be a scheme satisfying (H3)-(H5). Then, when ρ tends to 0, W ρ converges to the unique viscosity solution of (3.5) uniformly on each compact subset of [0, T ] × D. Proof. Let us define W (t, x, y) := lim sup [0,T ]×D∋(s,ξ,γ)→(t,x,y) ρ→0 W ρ (s, ξ, γ) and W (t, x, y) := lim inf [0,T ]×D∋(s,ξ,γ)→(t,x,y) ρ→0 W ρ (s, ξ, γ).

Theorem 5 . 4 .

 54 +M ε and the result is obtained by classical estimates thanks to the hypothesis on b and σ. Let assumptions (H1)-(H2) and (H1') be satisfied. There exists a constant C ≥ 0 such that

Figure 1 (

 1 Figure1(a) shows the approximation obtained using (N x 1 , N x 2 , N t ) = (800, 800, 200) (and the represented level set is 10 -5 ).The boundary of R C,K t (resp. R C,K t ) is also represented by a black doted line (resp. red line). The result perfectly matches the expected solution.

Figure 1 .

 1 Figure 1. (a): Test 1, no diffusion, and (b): Test 2, with diffusion.

Test 3 .

 3 In this test the volatility is now given by σ(x) = 5 d Θ (x) 0 1 where d Θ denotes the distance function to the set Θ := (x 1 , x 2 ), |x 2 | ≥ 0.3 . Hence for any point (x 1 , x 2 ), if |x 2 | ≥ 0.3 the volatility vanishes. According to the drift term, the target is finally reached, see Figure 2. This Figure shows the approximation of R C,K t

Figure 2 .

 2 Figure 2.(Test 3) Vertical diffusion, using (N x 1 , N x 2 , N t ) = (n, n, n/4), for n ∈ {200, 400, 800}.

  n L ∞ -error order L 1 -error order L 2 -error order 10 016633 1.48 0.00345 2.87 0.01979 1.84 40 0.06746 1.30 0.00111 1.63 0.00668 1.56 80 0.02500 1.43 0.00024 2.20 0.00194 1.78 Test 4 (oblique diffusion) In this test the coefficient σ is now given by σ(x) = 5 d Θ (x

Figure 3 .

 3 Figure 3.(Test 4) Oblique diffusion. using (N x 1 , N x 2 , N t ) = (n, n, n/4) for n ∈ {100, 200, 400}.

and w 0

 0 is a Lipschitiz continuous function. It is clear that under assumption (H1)-(H2), equation (3.5) is a particular case of (A.1), with n = d + 1, D ≡ Epigraph(g), ν ≡ -e y , H(t, •) ≡ H(T -t, •). We notice also that under (H1)-(H2), the set defined by D = Epigraph(g) satisfies the property (P1), with b := 1/ 1 + L 2 g . The properties (P2)-(P4) are also satisfied under (H1)-(H2).

Definition A. 1 (

 1 Parabolic semijets). Let O be a locally compact subset of R n . We define the parabolic semijets of a function w : [0, T ] × O → R in Lemma A.4. If (P1) and (P4) hold, then there exists h ∈ C 2 (D) such that: h ≥ 0 on D and h = 0 in D \ U (where U is a neighborhood of ∂D as in Property (P2)) and ν(x), Dh(x) ≥ 1 x ∈ ∂D.

  |x -x| 3 + |y -ȳ| 3 ) ≤ ũ( t, x) -ṽ( t, ȳ) + p • (x -x) + q • (y -|x -x| 3 + |y -ȳ| 3 ) ≤ ũ( t, x) -ṽ( t, ȳ) + p • (x -x) + q • (yx) -(y -ȳ)| 2 + λ|x -x| 2 + λ|y -ȳ| 2 , with p = D x ϑ ε (x, ȳ), q = D y ϑ ε (x, ȳ) and λ = |x-ȳ| 2 ε + ε.As a consequence of the Crandall-Ishii Lemma it follows that there exist X, Ỹ ∈ S m and a, b ∈ R such that -

  (a, D x ϑ ε (x, ȳ) + 2α(x -z) + 2δ x, X + 2αI + 2δI) ∈ P 1,2,+ u( t, x) and (-b, -D y ϑ ε (x, ȳ) -2δ ȳ, -Ỹ -2δI) ∈ P 1,2,-v( t, ȳ), so just setting X := X + 2δI and Y := Y + 2δI we get(a, D x ϑ ε (x, ȳ) + 2α(x -z) + 2δ x, X + 2αI) ∈ P 1,2,+ u( t, x) and (-b, -D y ϑ ε (x, ȳ) -2δ ȳ, -Y ) ∈ P 1,2,-v( t, ȳ), 15) Let us assume that for ε small enough D x ϑ ε (x, ȳ) + 2α(x -z) + 2δ x, ν(x) + α > 0 (A.16) -D y ϑ ε (x, ȳ) -2δ ȳ, ν(ȳ) -α < 0 (A.17)then by the definition of viscosity sub-and super-solution of (A.1) we geta + H( t, x, D x ϑ ε (x, ȳ) + 2α(x -z) + 2δ x, X + 2αI) < -c< 0 and -b + H( t, ȳ, -D y ϑ ε (x, ȳ) -2δ ȳ, -Y ) ≥ 0 and from (P2), (P3) (A.14) and (A.15), it followsc < -(a + b) + H( t, ȳ, -D y ϑ ε (x, ȳ) -2δ ȳ, -Y ) -H( t, x, D x ϑ ε (x, ȳ) + 2α(x -z) + 2δ x, X + 2αI) ≤ H( t, ȳ, -D y ϑ ε , -Y + CλI + 2δI) -H( t, x, -D y ϑ ε , X -CλI -2δI) + ω 1 (|D x ϑ ε + D y ϑ ε | + 2α + Cλ + 2δ(2 + |x| + |ȳ|) + 2α|x -z|) ≤ ω 1 (|D x ϑ ε + D y ϑ ε | + 2α + Cλ + 2δ(2 + |x| + |ȳ|) + 2α|x -z|)+ + ω 2 (C |x -ȳ| 2 ε + |x -ȳ|(2δ|ȳ| + 1))

D 2 ⇒

 2 x ϑ ε (x, ȳ), ν(x) ≥ -α D x ϑ ε (x, ȳ), ν(x) + 2α(x -z), ν(x) + 2δ x, ν(x) + α > 0 D y ϑ ε (x, ȳ), ν(ȳ) ≥ -α 2 ⇒ -D y ϑ ε (x, ȳ), ν(ȳ) -2δ ȳ, ν(ȳ) -α < 0,then thanks to (A.4) and (A.5), (A.[START_REF] Bonnans | Consistency of generalized finite difference schemes for the stochastic HJB equation[END_REF]) and (A.17) are finally obtained.In order to conclude the proof we need to argue how it is possible to modify our boundary conditions.Let us assume that u and v are respectively viscosity sub-and super-solution of the original problem. Thanks to Lemma A.4 we can define u α (t, x) := u(t, x) -αh(x) -Ct and v α (t, x) := v(t, x) + αh(x) + Ct.

  [START_REF] Camilli | An approximation scheme for the optimal control of diffusion processes[END_REF] and ν(x), p = ν(x), p + αDh(x) -α ν(x), Dh(x) (A.20)≤ ν(x), p + αDh(x) -α ν(x), p ≥ ν(x), p -αDh(x) + α (A.21)where M := max(|Dh| + D 2 h ) (from the construction of h follows the boundness of Dh and D 2 h).

Table 1 .

 1 Test 3, convergence table.

May 1, This work was partially supported by the EU under the 7th Framework Programme Marie Curie Initial Training Network "FP7-PEOPLE-2010-ITN", SADCO project, GA number 264735-SADCO.

1

(t, x) ∈ (0, T ) × O by: P 1,2,+ w(t, x) := (a, p, X) ∈ R × R n × S n | as (0, T ) × O ∋ (s, y) → (t, x) w(s, y) ≤ w(t, x) + a(s -t) + p, y -x + 1 2 X(y -x), (y -x) +o(|s

Then we define the closures of the previous sets:

Theorem A.2. Assume (P1)-(P4) hold. Let u (resp. v) be a usc subsolution (resp. lsc super-solution) of (A.1) satisfying the following growth conditions

Before starting the proof we state some important preliminary results.

Lemma A.3. If (P1) and (P4) hold, then there exists a family {w ε } ε>0 of C 2 functions on R n × R n and positive constants θ, C such that

and

for ε > 0 and x, y ∈ R n .

Proof. The result is presented in [START_REF] Dupuis | On oblique derivative problems for fully nonlinear secondorder elliptic partial differential equations on nonsmooth domains[END_REF].


Axel Bacher
email: bacher@lipn.univ-paris13.fr

Directed and multi-directed animals on the King's lattice

We define the directed King's lattice to be the square lattice with diagonal (next nearest neighbor) bonds and with the preferred directions { , , , , }. We enumerate directed animals on this lattice using a bijection with Viennot's heaps of pieces. We also define and enumerate a superclass of directed animals, the elements of which are called multi-directed animals. This follows Bousquet-Mélou and Rechnitzer's work on the directed triangular and square lattices. Our final results show that directed and multi-directed animals asymptotically behave similarly to the ones on the triangular and square lattices.

Introduction

An animal on a lattice is a finite and connected set of vertices. The enumeration of animals (up to a translation) is a longstanding problem in statistical physics and combinatorics. The problem, however, is extremely difficult, and little progress has been made [START_REF] Klarner | A procedure for improving the upper bound for the number of n-ominoes[END_REF][START_REF] Guttmann | On the number of lattice animals embeddable in the square lattice[END_REF]. A more realistic goal, therefore, is to enumerate natural subclasses of animals.

The class of directed animals is one of the most classical of these subclasses. Directed animals have been enumerated in a variety of lattices; let us cite, nonexhaustively, the square and triangular lattices [START_REF] Nadal | Directed lattice animals in 2 dimensions: numerical and exact results[END_REF][START_REF] Hakim | Exact results for 2D directed animals on a strip of finite width[END_REF][START_REF] Dhar | Equivalence of the two-dimensional directed-site animal problem to Baxter's hard-square lattice-gas model[END_REF][START_REF] Gouyou-Beauchamps | Equivalence of the two-dimensional directed animal problem to a one-dimensional path problem[END_REF][START_REF] Bétréma | Modèles avec particules dures, animaux dirigés et séries en variables partiellement commutatives[END_REF], Bousquet-Mélou and Conway's lattices L n [START_REF] Bousquet-Mélou | Enumeration of directed animals on an infinite family of lattices[END_REF][START_REF] Corteel | Bijections for directed animals on infinite families of lattices[END_REF], and the "strange" or n-decorated lattices [START_REF] Conway | Directed animals on two-dimensional lattices[END_REF][START_REF] Bousquet-Mélou | New enumerative results on two-dimensional directed animals[END_REF] (Figure 1). Unsolved lattices include, notably, the honeycomb lattice [START_REF] Guttmann | Hexagonal lattice directed site animals[END_REF].

The class of multi-directed animals is a superclass of directed animals, first introduced by Klarner [START_REF] Klarner | Cell growth problems[END_REF] on the square and triangular lattices. Bousquet-Mélou and Rechnitzer [START_REF] Bousquet-Mélou | Lattice animals and heaps of dimers[END_REF] clarified Klarner's definition and introduced a variant class on the square lattice. Moreover, they gave closed expressions for the generating function of multi-directed animals and showed that it is not D-finite.

The goal of this paper is to enumerate directed and multi-directed animals on a new lattice. We call King's lattice the square lattice with added diagonal bonds, or next nearest neighbor bonds. We also consider the preferred orientations { , , , , } (Figure 1, right). Directed animals on the King's lattice are a superclass of directed animals on Bousquet-Mélou and Conway's lattice L 3 , which has arcs { , , } [START_REF] Bousquet-Mélou | Enumeration of directed animals on an infinite family of lattices[END_REF].

Several techniques have been used to enumerate directed animals on the various lattices. Among them are direct bijections with other combinatorial objects [START_REF] Gouyou-Beauchamps | Equivalence of the two-dimensional directed animal problem to a one-dimensional path problem[END_REF], comparison with gas models [START_REF] Dhar | Equivalence of the two-dimensional directed-site animal problem to Baxter's hard-square lattice-gas model[END_REF][START_REF] Bousquet-Mélou | New enumerative results on two-dimensional directed animals[END_REF][START_REF] Borgne | Directed animals and gas models revisited[END_REF][START_REF] Albenque | A note on the enumeration of directed animals via gas considerations[END_REF] and the use of Viennot's theory of heaps of pieces [START_REF] Viennot | Heaps of pieces. I. Basic definitions and combinatorial lemmas[END_REF][START_REF] Bétréma | Modèles avec particules dures, animaux dirigés et séries en variables partiellement commutatives[END_REF][START_REF] Corteel | Bijections for directed animals on infinite families of lattices[END_REF][START_REF] Bousquet-Mélou | Lattice animals and heaps of dimers[END_REF][START_REF] Viennot | Multi-directed animals, connected heaps of dimers and Lorentzian triangulations[END_REF]. In this paper, we use the last method; we show that directed animals on the King's lattice are in bijection with heaps of segments, already defined in [START_REF] Bousquet-Mélou | Empilements de segments et q-énumération de polyominos convexes dirigés[END_REF]. 2 Animals on the King's lattice and heaps of segments 2.1 Definitions Definition 1. We call segment a closed real interval of the form [i, j], where i and j are integers such that j > i. Two segments are called concurrent if they intersect, even by a point. A heap of segments is a finite sequence of segments, considered up to commutation of non-concurrent segments.

The heaps of segments described here are the same as in [START_REF] Bousquet-Mélou | Empilements de segments et q-énumération de polyominos convexes dirigés[END_REF], except that the segment reduced to a point is not allowed. More information on heaps of pieces in general can be found in [START_REF] Viennot | Heaps of pieces. I. Basic definitions and combinatorial lemmas[END_REF]. Graphically, a heap is built by dropping segments in succession; a segment either falls on the ground or on another segment concurrent to it. Examples are shown in Figures 2 and3.

Directed animals and pyramids of segments

Let A be an animal; we say that a site t of A is connected to another site s if there exists a directed path (i.e. respecting the preferred directions of the lattice) from s to t visiting only sites of A. We say that the animal A is directed of source s if every site t of A is connected to s. The source s is not unique; it may be any of the bottommost sites of A (see Figure 2, left). By convention, we call source of A the leftmost bottom site.

In Figure 2 is illustrated a bijection between directed animals and pyramids of segments (or heaps with only one segment lying on the ground). This bijection works identically to the classical bijection between directed animals on the square lattice and strict pyramids of dimers [START_REF] Viennot | Heaps of pieces. I. Basic definitions and combinatorial lemmas[END_REF][START_REF] Bétréma | Modèles avec particules dures, animaux dirigés et séries en variables partiellement commutatives[END_REF]. Fig. 2. Left: a directed animal on the King's lattice (represented, for clarity, as a polyomino on the dual lattice) with its source circled. Right: the pyramid of segments obtained by replacing each maximal sequence of ℓ consecutive sites by a segment of length ℓ.

Multi-directed animals and connected heaps of segments

Let A be an animal. For any abscissa i, we denote by b(i) the ordinate of the bottommost site of A at abscissa i (or b(i) = +∞ if there is no site of A at abscissa i). We call source of A a site that realizes a local minimum of b and keystone of A a site that realizes a local maximum. In case several consecutive sites realize a minimum or maximum, the source or keystone is the leftmost one (Figure 3, left). This is a purely arbitrary choice that does not alter the definition. Definition 2. Let A be an animal. The animal A is said multi-directed if it satisfies the two conditions:

for every site t of A, there exists a source s such that t is connected to s; for every keystone t of A, there exist two sources s ℓ and s r , to the left and to the right of t respectively, such that t is connected to both s ℓ and s r . Moreover, the directed paths connecting t to s ℓ and s r do not go through a keystone at the same height as t.

As a directed animal has only one source and no keystone, every directed animal is multi-directed. Multi-directed animals are in bijection with connected heaps of segments (or heaps without an empty column). A multi-directed animal and its corresponding heap are depicted in Figure 3. Definition 2 can be adapted in the directed square and triangular lattices; the animals thus defined are in bijection with connected heaps of dimers. Bousquet-Mélou and Rechnitzer also defined multi-directed animals in bijection with connected heaps of dimers in [START_REF] Bousquet-Mélou | Lattice animals and heaps of dimers[END_REF], in a slightly different way. Our definition of multidirected animals has the advantages of being more intrinsic and of having a vertical symmetry.

Enumeration of directed animals

In this section, we use the bijection with pyramids of segments to enumerate directed animals on the King's lattice. We call half-animal a directed animal with no site on the left side of its source. The associated pyramids are called half-pyramids. We adapt Bétréma and Penaud's methods [START_REF] Bétréma | Modèles avec particules dures, animaux dirigés et séries en variables partiellement commutatives[END_REF] to decompose the pyramids of segments, which yields the following result. Theorem 3. The generating functions S(t) and D(t) of half-animals and animals are:

S(t) = 1 -3t - √ 1 -6t + t 2 4t ; D(t) = 1 4 1 + t √ 1 -6t + t 2 -1 .
The decomposition of the half-pyramids is sketched in Figure 4. Interestingly, the generating function S(t) is already known in combinatorics. Its coefficients are the little Schröder numbers, A001003 in the OEIS [START_REF] Sloane | The on-line encyclopedia of integer sequences[END_REF]. The coefficients of D(t) also appear as A047781. This is remindful of the triangular lattice, where the half-animals are enumerated by the Catalan numbers [START_REF] Bétréma | Modèles avec particules dures, animaux dirigés et séries en variables partiellement commutatives[END_REF].

P = P 1 P 2 + P ′ 1 P 2
Fig. 4. Sketch of the two cases in the decomposition of half-pyramids. The generating function of the possible heaps P 1 and P 2 is 1 + S(t), while the generating function of the possible heaps P ′ 1 is S(t). This shows the identity S = t(1 + S) 2 + tS(1 + S), from which we derive the value of S(t).

Enumeration of multi-directed animals

In this section, we enumerate multi-directed animals, or, equivalently, connected heaps of segments. To do this, we adapt the Nordic decomposition, invented by Viennot to enumerate connected heaps of dimers [START_REF] Viennot | Multi-directed animals, connected heaps of dimers and Lorentzian triangulations[END_REF].

Asymptotic results

Finally, we derive asymptotic results from Theorems 3 and 4.

Theorem 5. Let D n and M n be the number of directed and multi-directed animals of area n, respectively. As n tends to infinity, we have:

D n ∼ κ 3 + √ 8 n n -1/2 ; M n ∼ λ µ n ,
with µ = 6.475.... The average width of directed animals grows like √ n, while the average width of multi-directed animals grows like n. Finally, the series M(t) is not D-finite.

The results on directed animals are a straightforward application of singularity analysis [START_REF] Flajolet | Analytic combinatorics[END_REF]Theorem VI.4]. The results on multi-directed animals are more involved. Similar results already exist on the square and triangular lattices, including the non-D-finiteness of the series M(t) [START_REF] Bousquet-Mélou | Lattice animals and heaps of dimers[END_REF].

Fig. 1 .

 1 Fig.1. Directed animals on a selection of lattices. From left to right: the square lattice, the triangular lattice, the lattice L 3 , and the King's lattice.

Fig. 3 .

 3 Fig.3. Left: a multi-directed animal with four sources (circled) and three keystones (boxed). The directed paths connecting one keystone, denoted by t, to the sources s ℓ and s r are shown. Right: the corresponding connected heap of segments, with has four minimal pieces (one for each source of the animal).

Theorem 4 .S( 1 +

 41 Let M = M(t) be the generating function of multi-directed animals. Let S = S(t), D = D(t) be the power series defined in Theorem 3, R = S +t(1+S) and Q = 2(1t)St. The generating function M is given by: S) k QR k 1 -QR k .


Axel Bacher

Gilles Schaeffer

Multivariate Lagrange inversion formula and the cycle lemma ⋆

, and the extension of the cycle lemma to more than 2 variables was elusive. The cycle lemma has found a lot of applications in combinatorics, so we expect our multivariate extension to be quite fruitful: as a first application we mention economical linear time exact random sampling for multispecies trees.

Introduction

For any power series g(x) with g(0) = 0, there exists a unique power series f (t) solution of the equation f = tg( f ). The Lagrange inversion formula says that the nth coefficient of f (t) is 1 n [x n-1 ]g(x) n . This formula is now known as a fundamental tool to derive tree enumeration results. The two simpler and most classical examples are:

if g(x) = exp(x) then f (t) is the exponential generating function of Cayley trees, so that the number of rooted Cayley trees with n nodes is equal to 1 n

x n-1

(n-1)! exp(nx) = n n-2 ; if g(x) = (1 + x) 2 then f (t) is the ordinary generating function of binary trees, so that the number of rooted binary trees with n nodes is equal to

1 n [x n-1 ](1 + x) 2n = 1 n 2n n-1 = 1 n+1 2n
n . These examples are in a sense generic: the bijection between doubly rooted Cayley trees and endofunctions underlies Labelle's proof of the Lagrange inversion formula [START_REF] Labelle | Une nouvelle démonstration combinatoire des formules d'inversion de Lagrange[END_REF], while the cyclic lemma used in [START_REF] Dvoretzky | A problem of arrangements[END_REF] to count ballot numbers underlies Raney's proof [START_REF] Raney | Functional composition patterns and power series reversion[END_REF].

Our interest is in the multivariate extension of the Lagrange inversion formula. Fix an integer k ≥ 1. Let bold letters denote k-dimensional vectors; write

x n = x n 1 1 • • • x n k k and x n-1 = x n 1 -1 1 • • • x n k -1 k
. Let h(x) and g 1 (x), . . . , g k (x) be power series in x such that for i = 1, . . . , k, g i (0) = 0. Again there is a unique familly of power series f(t) solution of the system of equations May 27, 2013 and the multivariate Lagrange inversion formula admits the two equivalent following formulations (among several others):

f i = t i g i ( f i ) for i = 1, . . . , k ⋆ Version of
[t n ]h(f(t)) = [x n-1 ]h(x)g(x) n det δ i, j - x i g j (x) ∂ g j (x) ∂ x i (1) = 1 n 1 • • • n k [x n-1 ] ∑ T ∂ T (h, g n 1 1 , . . . , g n k k ) (2)
where the sum is over oriented 0-rooted Cayley trees (non-plane trees with arcs going toward 0) with vertices {0, . . . , k} and the derivative ∂ G with respect to a directed graph G with vertex set V = {0, . . . , k} and edge set E is defined as

∂ G ( f 0 (x), . . . , f k (x)) = k ∏ j=0 ∏ (i, j)∈E ∂ ∂ x i f j (x)
Several variants of (1) are given in [START_REF] Gessel | A combinatorial proof of the multivariable Lagrange inversion formula[END_REF] but (2) appeared more recently, implicitely in [START_REF] Goulden | Multivariable Lagrange inversion, Gessel-Viennot cancellation, and the matrix tree theorem[END_REF] and explicitely in [START_REF] Bender | A multivariate Lagrange inversion formula for asymptotic calculations[END_REF]. As far as we know all combinatorial proofs of the multivariate Lagrange inversion extend Joyal and Labelle's approach [START_REF]Une théorie combinatoire des séries formelles[END_REF][START_REF] Labelle | Une nouvelle démonstration combinatoire des formules d'inversion de Lagrange[END_REF]: [START_REF] Gessel | A combinatorial proof of the multivariable Lagrange inversion formula[END_REF] proves another variant of (1), [START_REF] Goulden | Multivariable Lagrange inversion, Gessel-Viennot cancellation, and the matrix tree theorem[END_REF][START_REF] Bousquet | Two bijective proofs for the arborescent form of the Good-Lagrange formula and some applications to colored rooted trees and cacti[END_REF] prove [START_REF] Bernardi | Counting trees using symmetries[END_REF]. A completely different approach was recently proposed in [START_REF] Bernardi | Counting trees using symmetries[END_REF]. Chottin [START_REF] Chottin | Une démonstration combinatoire de la formule de Lagrange à deux variables[END_REF][START_REF] Chottin | Énumération d'arbres et formules d'inversion de séries formelles[END_REF] instead proposed a remarkable extension of Raney's strategy to prove the two variable Lagrange inversion formula. Yet he failed to move to three variables and the problem of proving the multivariate Lagrange inversion formula with the cycle lemma was considered as difficult. Apart from the theoretical interest of such a proof, an extension of the cycle lemma is desirable in view of its numerous applications, to tree and map enumeration [START_REF] Lothaire | Combinatorics on words[END_REF]Chapter 11], [START_REF] Lothaire | Applied Combinatorics on words[END_REF]Chapter 9], probability [START_REF] Pitman | Enumerations of trees and forests related to branching processes and random walks[END_REF], and random sampling [START_REF] Dershowitz | The cycle lemma and some applications[END_REF]. We present such an extension in this paper.

Generalized Cycle Lemma

Following [START_REF] Chottin | Énumération d'arbres et formules d'inversion de séries formelles[END_REF] and the modern accounts of Raney's proof, our combinatorial construction is in terms of encodings of rooted plane trees by sequences of nodes. To deal with the multivariate case, we introduce colored trees.

A colored tree is a plane, rooted tree in which all edges have a color in the set {1, . . . , k}. A colored bush is a colored tree that can have "pending" edges with no node attached. We call such edges free edges; an edge is occupied if a node is attached to it.

Let ℓ be an integer with 0 ≤ ℓ ≤ k. An ℓ-bush is a colored bush such that occupied edges have color more than ℓ and free edges have color at most ℓ. In particular, a k-bush only has a single node, while a 0-bush is a colored tree. Therefore, ℓ-bushes can be seen as intermediate objects between colored nodes and colored trees. Definition 1. Let 0 ≤ ℓ ≤ k. We denote by S ℓ the set of tuples of the form S = (S 0 , . . . , S ℓ , e 1 , . . . , e k ), where S 0 , . . . , S ℓ are sequences of ℓ-bushes and e i is an edge of color i in S (refered to as a marked edge), satisfying the following conditions.

1. The sequence S 0 has only one element; for i = 1, . . . , ℓ, the number of elements of the sequence S i is equal to the number of edges of color i in S. 2. Let T be the graph with vertices 0, . . . , ℓ and an arc i → j if the edge e j is in the sequence S i . The graph T is a 0-rooted Cayley tree.

According to the previous remark, the objects of S ℓ may also be seen as intermediates between two objects. If k = ℓ, then S 0 , . . . , S k are simply sequences of nodes. If ℓ = 0, the unique element of the sequence S 0 is a colored tree.

Theorem 2 (generalized Cycle Lemma).

There is a bijection between the sets S ℓ and S ℓ-1 that works by attaching the elements of the sequence S ℓ to the edges of color ℓ.

The actual description of the bijection is given in Section 3; an example is given in Figure 1. Observe that for ℓ = k = 1, the statement is equivalent to the standard Cycle Lemma: there is a bijection between p-uples of rooted plane trees with one pointed node (represented here as a unique tree with an extra node of degree p at the root) and pairs formed of a node of degree p having a marked free edge and a sequence of nodes such that the number of nodes in the sequence equals the total number of free edges. Iterating the bijection yields a bijection between S 0 and S k . Given a Cayley tree T with vertices {0, . . . , k}, the generating function of the associated subset of

S k is hg n 1 1 • • • g n k
k with the proper derivatives to mark edges (taking a derivative of g n j j with respect to x i amounts to marking an edge of color i in the sequence S j ). It can thus be seen that the formula ( 2) is a corollary of the generalized Cycle Lemma.

The bijection

We now describe without proof our bijection between the sets S ℓ and S ℓ-1 , for 1 ≤ ℓ ≤ k, illustrated in Figure 1. This bijection uses Prüfer codes of Cayley trees [START_REF] Prüfer | Neuer beweis eines satzes über permutationen[END_REF]. The Prüfer code of a 0-rooted Cayley tree T with vertices 0, . . . , s is a sequence p 1 , . . . , p s with p s = 0; there exists a permutation σ such that the parent vertex of j in T is p σ ( j) for j = 1, . . . , s. Moreover, every sequence corresponds to a unique tree.

Let S = (S 0 , . . . , S ℓ , e 1 , . . . , e k ) be an element of S ℓ ; let T be the associated Cayley tree. In the following, we call ℓ-edge an edge of color ℓ. Since we are dealing with plane trees, there is a natural order on the set of ℓ-edges of S given by depth-first traversals of the ℓ-bushes.

Let r be the number of ℓ-edges in the sequences S 0 , . . . , S ℓ-1 . Since T is a Cayley tree, it does not contain the arc ℓ → ℓ, which means that the marked edge e ℓ is not in the sequence S ℓ . Thus, we have r ≥ 1. Let u be the root of the first element of S ℓ . By Definition 1, the sequence S ℓ has exactly r more elements than it has ℓ-edges. We combine the elements of S ℓ using the Cycle Lemma; we denote by b 1 , . . . , b r the resulting ℓ -1-bushes, so that the node u is in b 1 .

First case: the vertex ℓ is a leaf of T . In this case, we attach the bushes b 1 , . . . , b r to the remaining free ℓ-edges in order, with a cyclic permutation chosen so that the bush b 1 is attached to the marked edge e ℓ .

Second case: the vertex ℓ is not a leaf of T . In the Cayley tree T , all arcs going toward ℓ correspond to marked edges in the sequence S ℓ . Let m 1 , . . . , m s be the ℓ -1-bushes that we just formed containing at least one marked edge.

We break up the Cayley tree T in the following manner. For i = 1, . . . , s, let T i be the forest composed of the colors of the marked edges in m i and their descendants in T . Let T 0 be the tree composed of the remaining vertices of T with the vertex ℓ deleted; in other words, the vertices of T 0 are exactly the nondescendants of ℓ. We also assume that the order of the m i 's was chosen so that the T i 's are ordered according to their lowest label. Now, attach the bushes b 1 , . . . , b r to the free ℓ-edges in order, with a cyclic permutation chosen so that b s is attached to the marked edge e ℓ . For j = 1, . . . , s, let p j = i if m j is attached to an edge in the sequence with an index in T i . As the edge e ℓ is in the sequence corresponding to the parent of ℓ in T , which is in T 0 , we have p s = 0.

We can therefore regard the sequence p 1 , . . . , p s as the Prüfer code of a 0rooted Cayley tree T with labels {0, . . . , s}. Let σ be the permutation associated to this tree and swap the bushes m 1 , . . . , m s according to the permutation σ . Let S ′ 0 , . . . , S ′ ℓ-1 be the sequences resulting from this procedure. As no free ℓ-edges remain, every element of the sequences S ′ 0 , . . . , S ′ ℓ-1 is an ℓ -1-bush. Finally, let e ′ ℓ be the parent edge of the node u defined at the beginning. Let T ′ be the graph associated with the marked edges e 1 , . . . , e ℓ-1 . In the first case above, this is the tree T with the leaf ℓ deleted; in the second, it is a compound of the tree T 0 and the forests T 1 , . . . , T s , arranged according to the tree T ; it is therefore a Cayley tree. This shows that the tuple (S ′ 0 , . . . , S ′ ℓ-1 , e 1 , . . . , e ′ ℓ , . . . , e k ) is an element of S ℓ-1 .

We conclude with some comments. We use Cayley trees and Prüfer codes in a manner that may seem needlessly complicated; however, in the first stage, we see how important the condition that the graph T contains no edge ℓ → ℓ is, which is implied by the fact that T is a Cayley tree. The construction aims at ensuring that the graph T ′ describing the marked edges e 1 , . . . , e ℓ-1 remains a Cayley tree.

Fig. 1 .

 1 Fig.1. An example with k = 3, a tree of S 0 and the corresponding elements of S 1 , S 2 and S 3 (from left to right), with their associated Cayley trees (bottom). Marked edges are represented as double lines.


Olivier Mareschal

Sebastien Loiseau

Aurelien Fougerat

Laurie Valbin

Gaelle Bazin Lissorgues

Sébastien Saez

Christophe Dolabdjian

Rachid Bouregba

Gilles Poullain

Piezoelectric Aluminum Nitride Resonator for Oscillator

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.