Madalina Deaconu
email: madalina.deaconu@inria.fr

Antoine Lejay
email: antoine.lejay@iecn.u-nancy.fr

SIMULATION OF DIFFUSIONS BY MEANS OF IMPORTANCE SAMPLING PARADIGM

Keywords: 1991 Mathematics Subject Classification. Primary 60C05, secondary 65C, 65M, 68U20 Stochastic Differential Equations, Monte Carlo methods Random walk on squares, Random walk on rectangles Variance reduction, Simulation of rare events, Dirichlet/Neumann problems

The aim of this paper is to introduce a new Monte Carlo method based on importance sampling techniques for the simulation of stochastic differential equations. The main idea is here to combine random walk on squares or rectangles methods with importance sampling techniques.

The first interest of this approach is that the weights can be easily computed from the density of the one-dimensional Brownian motion. Compared to the Euler scheme this method allows to obtain a more accurate approximation of diffusions when one has to consider complex boundary conditions. The method provides also an interesting alternative to perform variance reduction techniques and to simulate rare events.

Introduction

Monte Carlo methods are sometimes the unique alternative in order to solve numerically partially differential equations (PDE) involving an operator of the form

L = 1 2 d i,j=1
a i,j (.)

∂ 2 ∂x i ∂x j + d i=1 b i (.) ∂ ∂x i .
The operator L is the infinitesimal generator associated to the solution of the stochastic differential equation (SDE)

(1)

X t = X 0 + t 0 σ(X s ) dB s + t 0 b(X s ) ds with σσ * = a.
It is well known that, for T > 0 fixed, the solution on the cylinder [0, T ] × D, of the parabolic PDE

      
∂u(t, x) ∂t + Lu(t, x) = 0, u(T, x) = g(x) for x ∈ D, u(t, x) = φ(t, x) for (t, x) ∈ [0, T ] × ∂D can be written as

u(t, x) = E t,x [g(X T ); T ≤ τ ] + E t,x [φ(τ, X τ ); τ < T ],
where τ stands for the first exit time of X from the domain D. E t,x means that the process X is starting from x at time t. Thus, an approximation of u(t, x) can be obtained by averaging g(X T )1 T ≤τ and φ(τ, X τ )1 τ <T over a large number of realizations of paths of X. Elliptic PDE may be considered as well.

A large spectra of methods has been already proposed in order to simulate X, see for example the books of P. Kloeden and E. Platen [START_REF] Kloeden | Numerical solution of stochastic differential equations[END_REF] and of G.N. Milstein and M.V. Tretyakov [START_REF] Milstein | Stochastic numerics for mathematical physics[END_REF]. Most of these methods are extensions of the Euler scheme, which provides a very efficient way to simulate [START_REF] Arouna | Adaptative Monte Carlo method, a variance reduction technique[END_REF] in the whole space. This method becomes harder to set up in a bounded domain, either with an absorbing or a reflecting boundary condition. Nevertheless some refinements have been proposed (see for example [START_REF] Bossy | A symmetrized Euler scheme for an efficient approximation of reflected diffusions[END_REF][START_REF] Gobet | Weak approximation of killed diffusion using Euler schemes[END_REF][START_REF] Gobet | Efficient schemes for the weak approximation of reflected diffusions[END_REF][START_REF] Jansons | Exponential timestepping with boundary test for stochastic differential equations[END_REF][START_REF] Pettersson | Approximations for stochastic differential equations with reflecting convex boundaries[END_REF][START_REF]Euler's approximations of solutions of SDEs with reflecting boundary[END_REF]). To improve the quality of the simulation or to speed it up, variance reduction techniques can be considered; see for example [START_REF] Arouna | Adaptative Monte Carlo method, a variance reduction technique[END_REF][START_REF] Arouna | Algorithmes stochastiques et méthodes de Monte Carlo[END_REF][START_REF] Bardou | Contrôle dynamique des erreurs de simulation et d'estimation de processus de diffusion[END_REF][START_REF] Kebaier | Réduction de Variance et discrétisation d'équations différentielles stochastiques. Théorèmes limites presque sûres pour les martingales quasi-continues à gauche[END_REF][START_REF] Heath | A variance reduction technique based on integral representations[END_REF][START_REF] Kohatsu-Higa | Variance reduction methods for simulation of densities on Wiener space[END_REF][START_REF] Newton | Variance reduction for simulated diffusions[END_REF][START_REF] Zou | Robust variance reduction for random walk methods[END_REF]. This list is not intended to be exhaustive.

In the simplest situation, for a = Id and b = 0, the underlying diffusion process is the Brownian motion. E. Muller proposed in 1956 a very simple scheme to solve a Dirichlet boundary value problem, this method is called the random walk on spheres method [START_REF] Muller | Some continuous Monte Carlo methods for the Dirichlet problem[END_REF]. The idea is to simulate successively, for the Brownian motion, the first exit position from the largest sphere included in the domain and centered in the starting point. This exit position becomes the new starting point and the procedure is iterated until the exit point is close enough to the boundary. Nevertheless, simulating the exit time from a sphere is numerically costly. In [START_REF] Mil'shtejn | An algorithm for random walks over small ellipsoids for solving the general Dirichlet problem[END_REF], G.N. Milstein and N.F. Rybkina proposed to use this scheme for solving [START_REF] Arouna | Adaptative Monte Carlo method, a variance reduction technique[END_REF] by freezing locally the value of the coefficients. In a first approach, spheres (that become ellipsoids) were used. Later on [START_REF] Mil'shteȋn | Solution of the first boundary value problem for equations of parabolic type by means of the integration of stochastic differential equations[END_REF] (see also the book [START_REF] Milstein | Stochastic numerics for mathematical physics[END_REF]), G.N. Milstein and M.V. Tretyakov used time-space parallelepipeds with a cubic space basis. For this last approach, it is easier to keep track of the time but the involved random variables are costly to simulate. In order to overcome these difficulties, one may think to use tabulated values. This is memory consuming as the random variables to simulate depend on one or two parameters. The method of random walk on squares was also independently developed in the PhD thesis of O. Faure [START_REF] Faure | Simulation du mouvement brownien et des diffusions[END_REF]. For the Brownian motion, this method is still a good alternative to the random walk on spheres (see [START_REF] Campillo | A Monte Carlo method without grid for a fractured porous domain model[END_REF] for an application in geophysics).

In [START_REF] Deaconu | A random walk on rectangles algorithm[END_REF], we have proposed a scheme for simulating the exact exit time and position from a rectangle for the Brownian motion starting from any point inside this rectangle. Compared to the random walk on spheres method, this method has the following advantages :

• It can be used whatever the dimension and, as for the random walk on squares, a constant drift term may be added. • The rectangles can be chosen prior to any simulation, and not dynamically.

There is no need to consider smaller and smaller spheres or squares when the particle is near the boundary. • The method can be also adapted and used for the simulation of diffusion processes killed on some part of the boundary. The method we propose here is based on the idea to simulate the first exit time and position from a parallelepiped by using an importance sampling technique (see for example [START_REF] Fishman | Monte Carlo. Concepts, algorithms, and applications[END_REF][START_REF] Glasserman | Monte Carlo methods in financial engineering[END_REF]). The exit time and position from a parallelepiped for a Brownian motion with locally frozen coefficients is chosen arbitrarily, and a weight is computed at each simulation. By repeating this procedure, we get the density on the boundary or at a given time of the particles, by weighting the simulated paths. As we will see, the weights are rather easily deduced from the density of the one-dimensional Brownian motion killed when it exits from [-1, 1]. All involved expressions are numerically easy to implement. This new algorithm is slower than the Euler scheme for smooth coefficients, but it is faster than the random walk on squares [START_REF] Campillo | A Monte Carlo method without grid for a fractured porous domain model[END_REF][START_REF] Milstein | Stochastic numerics for mathematical physics[END_REF] and the random walk on rectangles [START_REF] Deaconu | A random walk on rectangles algorithm[END_REF]. It can be used to simulate the Brownian motion as well as solutions of stochastic differential equations for specific complex situations as: (a) complex geometries (the boundary conditions are correctly taken into account); (b) fast estimation of the exit time of a domain for the Brownian motion (only few rectangles are needed); (c) variance reduction; (d) simulation of rare events. This algorithm could be relevant for many domains: finance, physics, biology, geophysics, etc. It may also be used locally (for example, it can be mixed with the Euler scheme and used when the particle is close to the boundary) or combined with other algorithms, such as population Monte Carlo methods (see Section 4.5).

We conclude this article with numerical simulations illustrating various examples. It has to be noted that choosing "good" distributions for the exit time and position from a rectangle is not an easy task in order to reduce the variance. We then plan to study in the future how to construct algorithms that minimize the variance, as in [START_REF] Arouna | Adaptative Monte Carlo method, a variance reduction technique[END_REF][START_REF] Bardou | Contrôle dynamique des erreurs de simulation et d'estimation de processus de diffusion[END_REF]. We have to consider for this a high dimension optimization problem. Outline. In Section 2, we present the importance sampling technique applied to the exit time and position for a (drifted) Brownian motion from a rectangle. In Section 3, we recall briefly some results about the density of the one-dimensional Brownian motion with different boundary conditions. The explicit expressions are given in Appendix A. In Section 4, we present our algorithm and compute its weak error. Four test cases are presented in Section 5, we compare also our algorithm with other methods in this last section.

2. Algorithm for the exit time and position from a right time-space parallelepiped by using an importance sampling method

The aim of this part is to give a clear presentation of our method. In order to avoid ambiguous notations we consider in this section the situation of a twodimensional space domain. The results can be easily generalized to higher space dimension.

We are looking for an accurate approximation of the exit time and position from a right time-space parallelepiped which is a geometric figure in the 3-dimensional space.

For

L 1 , L 2 > 0 given let R be the rectangle [-L 1 , L 1 ] × [-L 2 , L 2 ]
. The rectangle R is the space basis of the right time-space parallelepiped R T = [0, T ] × R for a fixed T > 0. We can also consider R ∞ = R + × R, and set in this case T = +∞.

For T < +∞, the right time-space parallelepiped R T has six sides which are denoted by

S 0,1 = {T } × R, S 0,-1 = {0} × R, S 1,η = [0, T ] × [-L 1 , L 1 ] × {ηL 2 } for η ∈ {-1, 1}, S 2,η = [0, T ] × {ηL 1 } × [-L 2 , L 2 ] for η ∈ {-1, 1}.
In other words, each side of R T is labelled by a couple (i, η) ∈ {0, 1, 2} × {-1, 1}. For i ∈ {1, 2} the side S i,η is perpendicular to the unit vector in the i-th direction. For i = 0, the side S 0,-1 corresponds to the rectangular initial basis while the side S 0,1 corresponds to the top of the time-space parallelepiped R T for T < +∞ (See Figure 1).

From now on, we shall identify each side with the corresponding (i, η)-indices.

We consider a time-homogeneous diffusion process (X t ) t≥0 living in R. On each side of R, the process X may be reflected or absorbed. Moreover, if T < +∞, the process is stopped at time T . We can thus identify the sides of R with the sides S i,η of R T for i ∈ {1, 2} and η ∈ {-1, 1}. We denote by R the subset of {1, 2} × {-1, 1} that contains the indices of the sides on which a Neumann boundary condition holds (possibly, R = ∅). On this set the diffusion is reflected. Let us denote by D the subset of {1, 2} × {-1, 1} that contains the indices of the sides on which a Dirichlet boundary condition holds. On this set the diffusion is killed. Finally let us set A = D if T = +∞ and A = {(0, 1)} ∪ D if T < +∞. With these notations the time-space process t → (t, X t ) is absorbed when hitting one of the sides S i,η with (i, η) ∈ A.

x 1 x 2 t Side S 0,1 Side S 1,1 Side S 1,-1 x 1 x 2 t Side S 0,-1 Side S 2,1 Side S 2,-1 Figure 1. Convention for the sides of R T = [0, T ] × R. Let B = (B 1 , B 2
) be a two-dimensional Brownian motion and µ = (µ 1 , µ 2 ) a vector of R 2 . For i ∈ {1, 2}, we set

γ i,η = 1 if (i, η) ∈ R (reflection), 0 if (i, η) ∈ A (absorption).
We consider the two-dimensional diffusion process (X, P x ) x∈R whose coordinates are, for

x = (x 1 , x 2 ) ∈ R, (2)
X i t = x i + B i t + µ i t + γ i,1 L i t (X i ) -γ i,-1 -L i t (X i ), P x -a.s.,
where ±L i t (X i ) stands for the symmetric local time of X i at ±L i , respectively. We define

τ 0 = T , τ i = inf{t > 0 | |X i t | > L i } for i ∈ {1, 2} and τ = min i∈{0,1,2} τ i .
In addition, we set J = argmin i∈{0,1,2} τ i . With this notation, unless J = 0, the J-th component of X is the first to exit from the domain. For J ∈ {1, 2}, let us define ε = X J τ J /L J ∈ {-1, 1}. For J = 0 we set ε = 1, in this case X has not reached the sides of D before time T .

The couple (J, ε) labels the side in A of the parallelepiped R T = [0, T ] × R that the diffusion X hits first. Note that with our convention, the sides on which the process is reflected cannot be reached, so that τ i = +∞ if X i is reflected both at -L i and L i .

We are interested in computing E x [f (τ, X τ )] by a Monte Carlo method for a bounded, measurable function f , where τ is defined as above.

Instead of simulating (τ, X τ ), we will simulate some random variables according to the following procedure. The aim is to simulate (J, ε, τ, X τ ) by using an importance sampling technique. In order to do this we choose a probability P x which is absolutely continuous with respect to P x and we draw a realization of (J, ε, τ, X τ ). Let us set

α i,η = P x [(J, ε) = (i, η)]
for (i, η) ∈ A. For (i, η) ∈ A let k i,η denote the density under P x of (τ, X τ ) given {(τ, X τ ) ∈ S i,η }.

In order to simplify notations let us consider an underlying probability space (Ω, F, P x ) rich enough. Let Z be a random variable on this space, with distribution P x . Let A be a measurable event on this space. We suppose that, conditionally on A, Z has a density p(•|A) with respect to the Lebesgue measure. Let us introduce the following convention

P x [Z = z; A] = p(z|A)P x [A].
That is, for B a measurable event of (Ω, F, P x ),

P x [{Z ∈ B} ∩ A] = B p(z|A)P x [A] dz = B P x [Z = z; A] dz.
Consider now the following notations : Let (i, η) ∈ A. For i ∈ {1, 2} set j = 3 -i. Then for any θ > 0 and z ∈ S i,η , we define

(3) M i,η (θ, z) = P x [τ i = θ; X i τ i = ηL i ]P x [X j θ = z j ; τ j > θ] α i,η k i,η (θ, z) ,
where k i,η is the {X τ ∈ S i,η }-conditional density under P x of (τ, X τ ). If T < +∞, we define

(4) M 0,1 (T, z) = 1 α 0,1 k 0,1 (T, z) j∈{1,2} P x [X j T = z j ; τ j > T ],
where k i,η is the {X τ ∈ S i,η }-conditional density under P x of (τ, X τ ). We call M i,η weight.

Proposition 1. The weights M i,η defined in (3) and (4) satisfy

E x [f (τ, X τ )] = E x [M J,ε (τ, X τ )f (τ, X τ )],
for any measurable and bounded function f on ∂R T .

Before proving this proposition let us introduce the algorithm. The algorithm is described as follows : Algorithm 1. Let x be fixed in R.

(1) Draw a realization (J, ε) of (J, ε) ∈ A under P x .

(2) Draw a realization of the exit time and exit position (τ , X τ ) according to the density k J,ε on S J,ε .

(3) Compute the value of M J,ε (τ , X τ ) by

E x [M J,ε (τ, X τ )f (τ, X τ )] = E x [f (τ, X τ )].
We call M J,ε (τ , X τ ) weight.

If {(J (i) , ε (i) , τ (i) , X (i)
τ , w (i) )} i=1,...,N are N independent realizations of the random variables (J, ε, τ, X τ , M J,ε (τ, X τ )) constructed as above, by the law of large numbers we have

E x [f (τ, X τ )] = lim N →∞ 1 N N i=1 w (i) f (τ (i) , X (i) τ (i) ).
The main feature of our approach is that the weights M J,ε (τ, X τ ) can be easily evaluated.

Remark 1. In order to evaluate M i,η with (3) and ( 4), there is no need to know P x [(J, ε) = (i, η)]. It is important to notice that M i,η depends only on the onedimensional distributions of the drifted Brownian motion.

Proof of the Proposition 1. We want to prove that :

E x [f (τ, X τ )] = E x [M J,ε (τ, X τ )f (τ, X τ )],
for any measurable and bounded function f on ∂R T .

We remark first that if

p i,η = P x [(J, ε) = (i, η)] for (i, η) in A, then (5) E x [f (τ, X τ )] = (i,η)∈A p i,η α i,η E x [M i,η (τ, X τ )f (τ, X τ ) | (J, ε) = (i, η)].
Furthermore, for (i, η) ∈ D, if i = 2 set j = 1 and z = (z 1 , ηL 2 ) else, if i = 1 set j = 2 and z = (ηL 1 , z 2 ).

E x [f (τ, X τ ) | (J, ε) = (i, η)] = [0,T ]×[-L j ,L j ] f (θ, z)P x [(τ i , X j τ i ) = (θ, z j ) | (J, ε) = (i, η)] dθ dz j where P x [(τ i , X j τ i ) = (θ, z j ) | (J, ε) = (i, η)] is the {(J, ε) = (i, η)}-conditional den- sity of (τ i , X j τ i ) with respect to dt dz j . Hence E x [f (τ, X τ ) | (J, ε) = (i, η)] = E x f (τ, X τ )M i,η (τ, X τ ) (J, ε) = (i, η) ,
where

M i,η (θ, z) = P x [(τ i , X j τ i ) = (θ, z j ) | (J, ε) = (i, η)] k i,η (τ, X τ ) .
Let us note that M i,η (θ, z) = M i,η (θ, z)p i,η /α i,η . With [START_REF] Bossy | A symmetrized Euler scheme for an efficient approximation of reflected diffusions[END_REF], we can deduce that

E x [f (τ, X τ )] = E x [f (τ, X τ )M J,ε (τ, X τ )].
Indeed, it suffices to remark that for (i, η) ∈ D,

M i,η (θ, z) = 1 α i,η k i,η (θ, z) P x [(τ i , X j τ i ) = (θ, z j ); (J, ε) = (i, η)] = 1 α i,η k i,η (θ, z) P x [(τ i , X j τ i ) = (θ, z j ); X i τ i = ηL i , τ j > θ].
The independence of the coordinates of X leads to the desired equality. If T < +∞, similar computations imply that for z ∈ R,

M 0,1 (T, z) = 1 α 0,1 k 0,1 (T, z) P x X T = z; min i∈{1,2} τ i > T
and the conclusion also holds.

Let us evaluate these probabilities. For i ∈ {1, 2}, let p i (t, x 1 , x 2 ) be the solution of ( 6)

         ∂p i (t, x 1 , x 2 ) ∂t = 1 2 ∂ 2 p i ∂x 2 2 (t, x 1 , x 2 ) + µ i ∂p i ∂x 2 (t, x 1 , x 2 ) for (t, x 1 , x 2 ) ∈ R + × (-L i , L i ) 2 , p i (t, x 1 , x 2 ) --→ t 0 δ x 1 (x 2 ) for (x 1 , x 2 ) ∈ (-L i , L i ) 2 ,
with the following boundary conditions (b.c.)

p i (t, x 1 , -L i ) = 0 (Dirichlet b.c.) if (i, -1) ∈ A, ∂p i ∂x 2 (t, x 1 , -L i ) = 0 (Neumann b.c.) if (i, -1) ∈ R, p i (t, x 1 , L i ) = 0 (Dirichlet b.c.) if (i, 1) ∈ A, ∂p i ∂x 2 (t, x 1 , L i ) = 0 (Neumann b.c.) if (i, 1) ∈ R.
Thus, p i denotes the density of the drifted Brownian motion X i with possibly some reflection at the endpoints of (-L i , L i ), and killed when it exits from this interval by an endpoint where no reflection holds. For f a bounded measurable function from [-L i , L i ] to R, we have

E x 1 [f (X i t ); t < τ i ] = L i -L i p i (t, x 1 , x 2 )f (x 2 ) dx 2 for x 1 ∈ [-L i , L i ], where P x 1 is the distribution of X i with X i 0 = x 1 ∈ [-L i , L i ].
Let us note that the distribution of the marginal X i of X under P (x 1 ,x 2 ) depends only on x i .

We introduce the scale function Φ i,+ of X i defined by,

for x 2 ∈ [-L i , L i ], Φ i,+ (x 2 ) =      e 2µ i L i -e -2µ i x 2 e 2µ i L i -e -2µ i L i if µ i = 0, x 2 + L i 2L i if µ i = 0. The function Φ i,+ (x 2 ) has been normalized such that Φ i,+ (L i ) = 1. Let us note that Φ i,+ (x i ) = P x i [X i τ i = L i ] if Dirichlet boundary conditions hold at both endpoints -L i and L i . We also set Φ i,-(x 2 ) = 1 -Φ i,+ (x 2 ).
If Dirichlet boundary conditions hold both at -L i and L i , then we set for t > 0 and (

x 1 , x 2 ) ∈ [-L i , L i ] 2 , p i,± (t, x 1 , x 2 ) = p i (t, x 1 , x 2 ) Φ i,± (x 2 ) Φ i,± (x 1
) .

Via a Doob transform, for a bounded and measurable function f ,

E x 1 [f (X i t ); t < τ i | X i τ i = ±L i ] = L i -L i p i,± (t, x 1 , x 2 )f (x 2 ) dx 2 .
Let us set for

x 1 ∈ (-L i , L i ), q i (t, x 1 ) = - L i -L i ∂p i ∂t (t, x 1 , x 2 )f (x 2 ) dx 2 (7)
and

q i,± (t, x 1 ) = - L i -L i ∂p i,± ∂t (t, x 1 , x 2 )f (x 2 ) dx 2 . ( 8
)
We can easily deduce that

P x 1 [τ i ≤ t] = t 0 q i (s, x 1 ) ds and P x 1 [τ i ≤ t | X i τ i = ±L i ] = t 0 q i,± (s, x 1 ) ds.
In other words, q i (t, x 1 ) (respectively q i,± (t, x 1 )) is the density of the first exit time from

[-L i , L i ] for X i (respectively the first exit time from [-L i , L i ] for X i given {X i τ i = ±L i }).
Thanks to these expressions, M 0,1 (T, z) and M i,η (θ, z) are easily computed, since

P x i [X i θ = z i ; τ i > T ] = p i (θ, x i , z i ), P x i [τ i = θ; X i θ = ±L i ] = q i,± (θ, x i )Φ i,± (x i ) if (i, -1) ∈ A and (i, 1) ∈ A, P x i [τ i = θ; X i θ = L i ] = q i (θ, x i ) if (i, -1) ∈ R and (i, 1) ∈ A, P x i [τ i = θ; X i θ = -L i ] = q i (θ, x i ) if (i, -1) ∈ A and (i, 1) ∈ R.

Analytical expressions for the densities

In order to compute p i (t, x 1 , x 2 ) together with q i (t, x 1 ) and q i,± (t, x 1 ) by ( 7) and ( 8), one has to solve the equation [START_REF] Cérou | Genetic genealogical model in rare event analysis[END_REF]. By using a scaling principle, we may assume that L i = 1, as

p i (t, x 1 , x 2 ) = 1 L i p t L 2 i , x 1 L i , x 2 L i ; L i µ ,
where p(t, x 1 , x 2 ; δ) is solution to [START_REF] Cérou | Genetic genealogical model in rare event analysis[END_REF] with L i = 1 and a convective term µ i equal to δ.

There are basically two ways to obtain p(t, x 1 , x 2 ; δ). The first one is based on the spectral expansion of 1 2 + δ∇, since this operator may be reduced to a selfadjoint one with respect to the scalar product induced by the measure exp(-2δx 1 ). The second one is the method of images when δ = 0.

If δ = 0, the case of a Dirichlet boundary condition at both endpoints may be treated by using a simple transform that reduces the problem to δ = 0.

For the case of Neumann boundary condition at both endpoints, one can invert term by term the Laplace transform of a series for the Green function.

In the case of a mixed boundary condition, the previous method gives rise to series that cannot be used in practice, so only the spectral expansion should be used. In addition, the first eigenvalues have to be computed numerically.

As the formula are standard in most of the cases, we give the relevant expressions in Appendix A.

General domain

As stated before, we aim to solve by a Monte Carlo method a parabolic or an elliptic PDE. The idea is to represent the domain as the union of time-space parallelepipeds and to simulate the successive exit times and positions from these parallelepipeds. Attention has to be paid while doing this decomposition in order to control the error at each simulation step.

4.1.

From parallelepipeds to right parallelepipeds. Consider herein the notations of Section 2. Let us study first the parabolic PDE with constant coefficients λ, c and µ = (µ i ) i=1,...,d on the rectangle R T :

(9)                  ∂v(t, x) ∂t + 1 2 d i=1 ∂ 2 v(t, x) ∂x 2 i + d i=1 µ i ∂v(t, x) ∂x i + cv(t, x) = λ on R T , ∂v(t, x) ∂x i = 0 for x ∈ S i,η if (i, η) ∈ R, v(t, x) = φ(t, x) for x ∈ S i,η if (i, η) ∈ A, v(T, x) = g(x) if T < +∞.
We assume that a classical solution to this problem exists, which is for example the case if φ and g are continuous and bounded. Let X be the diffusion process whose components are given by [START_REF] Arouna | Algorithmes stochastiques et méthodes de Monte Carlo[END_REF]. Then it follows from the Itô formula applied to X that, for t ∈ [0, T ],

v(t, x) = E x [e c(τ -t) φ(τ -t, X τ -t )); τ < T -t] + E x [e c(T -t) g(X T -t ); τ = T -t] + E x λ τ -t 0 e c(τ -t-s) ds ,
where τ is as above the first exit time from R T .

Let us remark that is σ is an invertible

d × d-matrix then the function u(t, x) = v(t, σ -1 x) is solution to (10)                              ∂u(t, x) ∂t + 1 2 d i,j=1 [σσ * ] i,j ∂ 2 u(t, x) ∂x i ∂x j + d i=1 [µσ * ] i ∂u(t, x) ∂x i + cu(t, x) = λ on [0, T ] × σR, σ j,i ∂u(t, x) ∂x j = 0 for x ∈ σS i,η if (i, η) ∈ R, u(t, x) = φ(t, σ -1 x) for x ∈ σS i,η if (i, η) ∈ A, u(T, x) = g(σ -1 x) if T < +∞.
If n i is the unit vector orthogonal to the side σS i,η , then n i = (σ * ) -1 e i , where e i is the unit vector in the i-th direction. It follows that σσ * n i = σe i and thus

for x ∈ σS i,±1 , [σσ * ]n i • ∇u(t, x) = σ j,i ∂u(t, x) ∂x j ,
which means that a Neumann boundary condition in the co-normal direction holds in [START_REF] Del Moral | Genealogical particle analysis of rare events[END_REF]

on σS i,η if (i, η) ∈ R.
We can thus solve [START_REF] Del Moral | Genealogical particle analysis of rare events[END_REF] by reducing the problem to (9) and use a Monte Carlo method in order to compute the values of u(t, x).

: Q → R d , c : Q → R and f : Q → R. For all (t, x) ∈ Q, we denote by σ(t, x) a d × d-symmetric matrix such that σ(t, x)σ * (t, x) = a(t, x).
We set

L = 1 2 d i,j=1 a i,j (t, x) ∂ 2 ∂x i ∂x j + d i=1 b i (t, x) ∂ ∂x i .
Let us introduce the hypotheses needed to ensure the convergence of our algorithm. To set up a Monte Carlo numerical scheme, one needs three inter-connected ingredients :

(i) The existence and the uniqueness of a solution u to the following PDE ( 11)

           ∂u(t, x) ∂t + Lu(t, x) + c(t, x)u(t, x) + f (t, x) = 0 on [0, T ] × D, u(T, x) = g(x), x ∈ D, u(t, x) = φ(t, x) on Γ d ⊂ [0, T ) × ∂D, ∂ n u(t, x) = 0 on Γ n ⊂ [0, T ) × ∂D,
where ∂ n denotes the co-normal derivative along the lateral surface. Γ d (respectively Γ n ) are subsets of [0, T ) × ∂D on which a Dirichlet (respectively Neumann) boundary condition holds. (ii) The existence of a solution to the diffusion process associated to L. Note that since the simulation involves distributions and not stochastic integrals, we do not need strong existence for the associated SDE. (iii) The solution u can be expressed in terms of the probabilistic representation ( 12)

u(t, x) = E t,x exp τ t c(s, X s ) ds φ(τ, X τ )1 τ <T +E t,x exp T t c(s, X s ) ds g(X T )1 τ >T +E t,x τ ∧T t exp s t c(r, X r ) dr f (s, X s ) ds ,
where τ is the first exit time from [0, +∞) × D by a point of Γ d . Notation 1. We denote by P the set of time-space parallelepipeds P such that there exist 0 ≤ s < t ≤ T , L 1 , . . . , L d and x ∈ R d such that

P = [s, t] × (x + σ([-L 1 , L 1 ] × • • • × [-L d , L d ])), where σ is a d × d-matrix. Possibly t = +∞ (if T = +∞).
The assumptions that have to be done are the following: (H1) There exists a subset

P D of P such that Q = ∪ P ∈P D P . Besides, if P = [s, t] × U ∈ P for a parallelepiped U , then for all r ∈ [s, t), [r, t] × U ∈ P.
In other words, one can truncate the parallelepipeds in time. (H2) There exist Γ n , Γ d contained in ∂Q = [0, T ]×∂D and some subsets P n , P d of (H3) The differential operator L is the generator of a continuous diffusion process X that is reflected at Γ n and killed when hitting Γ d ∪ {T } × D. The probabilistic representation of the solution given by ( 12) holds (see for example [START_REF] Lions | Stochastic differential equations with reflecting boundary conditions[END_REF] for existence results of such reflected process, and [START_REF] Stroock | Multidimensional diffusion processes[END_REF] if there are no reflections). (H4) There exists an unique solution u of class C 1,2 on [0, T ) × D to [START_REF] Faure | Simulation du mouvement brownien et des diffusions[END_REF] which is continuous on [0, T ] × D. (H5) For a right parallelepiped R and a matrix σ let P = [s, t] × (x + σR) ∈ P D .

I such that Γ n ⊂ ∪ P ∈Pn ∂P , Γ d ⊂ ∪ P ∈P d ∂P . The closure of Γ n ∪ Γ d is equal to [0, T ]×∂D
We associate to P a vector b ∈ R d , two constants c, f and we construct the differential operator

L = 1 2 d k,l=1 a k, ∂ 2 ∂x k ∂x + d k=1 b k ∂ ∂ x k with a = σ σ * .
Fix δ > 0. We assume that the solution u to (11) satisfies for any y in the interior of x + σR,

E s,y τ s e b c(r-s) ∂u ∂t + Lu + cu -f (r, X r ) dr ≤ δ,
where X is the diffusion process generated by L and τ is its first exit time from P .

Remark 2. If T = +∞ and, the coefficients are time-homogeneous and

Γ d = [0, ∞) × γ d , Γ n = [0, ∞) × γ n , then v(x) = u(0, x) is solution to the elliptic PDE (13)      Lv(x) + c(x)v(x) = f (x) on D, v(x) = φ(x) on γ d ⊂ ∂D, ∂ n v(x) = 0 on γ n ⊂ ∂D.
Thus, by solving the parabolic PDE [START_REF] Faure | Simulation du mouvement brownien et des diffusions[END_REF], we may also solve the elliptic PDE [START_REF] Glynn | Importance sampling for stochastic simulations[END_REF].

We will thus focus only on [START_REF] Faure | Simulation du mouvement brownien et des diffusions[END_REF].

Remark 3. The result of the existence of a stochastic process reflected on some part of the boundary of [0, T ) × D is deduced from the existence of a stochastic process reflected on the lateral boundary [0, T )×D, which is killed when it hits Γ n .

4.3. The algorithm and its weak error. In order to simplify the notations, if T < +∞, we denote the final condition g of ( 11) by φ(T, x). Given (t, x) ∈ Q, the solution u(t, x) of ( 11) is computed by the Feynman-Kac formula. For this, we have to simulate the diffusion process X up to its first exit time τ from Q. We suppose here that the particle cannot exit by a part of boundary where a Neumann boundary condition holds. Let u be the solution of [START_REF] Faure | Simulation du mouvement brownien et des diffusions[END_REF]. Let us introduce the following notations

for s ≥ t,        Y s = 1 + s t c(r, X r )Y r dr = exp s t c(r, X r ) dr , Z s = s t f (r, X r )Y r dr.
Then u(t, x) is then given by ( 14)

u(t, x) = E t,x [φ(τ, X τ )Y τ + Z τ ].
We construct now the algorithm that approximates ( 14) by a Monte Carlo method.

Algorithm 2. Assume that we start initially at the point (t, x) ∈ Q and fix a number N of particles.

(1) For i = 1, . . . , N do (A) Set (θ 0 , Ξ 0 , Y 0 , Z 0 , W 0 ) = (t, x, 1, 0, 1) and k = 0. (B) Repeat (a) Choose an element P (k) ∈ P D of the form

P (k) = [θ k , s] × U , U ⊂ R d such that (θ k , Ξ k )
Y k+1 = Y k exp(c (k) (θ k+1 -θ k )) Z k+1 = Z k + f (k) θ k+1 θ k exp(c (k) s) ds. (d) If Ξ k+1 ∈ Γ d or θ k+1 = T , then exit from the loop. (e) Increase k. (C) Set (θ (i) , Ξ (i) , Y (i) , Z (i) , W (i) ) = (θ k+1 , Ξ k+1 , Y k+1 , Z k+1 , W k ). (2) Return (15) u(t, x) = 1 N N i=1 W (i) φ(θ (i) , Ξ (i) )Y (i) + W (i) Z (i) .
We denote from now on by P x the distribution of the Markov chain

Λ k = (θ k , Ξ k ), k ≥ 0. Note that (Y k , Z k , w k ) k≥0 is obtained from (Λ k ) k≥0 . Proposition 2. For any (t, x) ∈ [0, T ) × D, (16) |u(t, x) -E x [ u(t, x)]| ≤ δ E x [W ν ν exp(M θ ν )] ,
where δ is defined in (H5), ν is the number of steps that (Λ k ) k≥0 takes to reach the boundary Γ d ∩ {T } × D and

M = sup (s,y)∈[t,T )×D c(s, y).
Remark 4. Note that the weak error in [START_REF] Gobet | Efficient schemes for the weak approximation of reflected diffusions[END_REF] does not depend on the choice of the importance sampling technique, while the Monte Carlo error depends on this choice. If the coefficients a, b, f and c are constant on the domain, one can choose δ = 0 and the simulation becomes exact.

Proof. To the Markov chain (Λ k ) k≥0 is associated a random sequence of parallelepipeds (P (k) ) k=0,...,ν . Let us denote by τ (k) the successive times the diffusion process X reaches the boundary of the P (k) 's.

Since Z 0 = 0, Y 0 = 1 and u = φ on the boundary of Q, we get ( 17)

E x [ u(t, x)] = E x [W ν Y ν φ(θ ν , Ξ ν ) + W ν Z ν ] = u(t, x) + E x W ν ν-1 k=0 (Z k+1 -Z k + Y k+1 u(θ k+1 , Ξ k+1 ) -Y k u(θ k , Ξ k )) .
Let (G k ) k≥0 be the filtration generated by the Markov chain (Λ k ) k≥0 . We remark that Y k and Z k are measurable with respect to G k , while w k is measurable with respect to G k+1 (since it is obtained from θ k , Ξ k , θ k+1 and Ξ k+1 ).

By using the Markov property, after setting

W k+1,ν = E x [w k+1 • • • w ν | G k+1 ],
we get

E x [W ν (Z k+1 -Z k )] = E x [W k+1,ν E x [w k (Z k+1 -Z k ) | G k ]W k-1 ], E x [W ν (Y k+1 u(θ k+1 , Ξ k+1 ) -Y k u(θ k , Ξ k ))] = E x [W k+1,ν E x [w k (Y k+1 u(θ k+1 , Ξ k+1 ) -Y k u(θ k , Ξ k )) | G k ]W k-1 ].
Let us denote by (X (k) , P (k) t,x ) the process generated by the operator L (k) with constant coefficients a (k) and b (k) on P (k) . Define recursively (t (0) , x (0) ) = (t, x) and (t (k+1) , x (k+1) ) = (τ (k) , X (k) τ (k) ), where τ (k) is, as above, the first exit time from P (k) for the diffusion X (k) . Let also f (k) and c (k) be the values that approach f and c on P (k) , and define also recursively y (0) = 1 and y

(k) = y (k-1) exp(c (k) (t (k+1) -t (k) )).
By using the properties of P x and the Itô formula we obtain

E x [w k (Y k+1 u(θ k+1 , Ξ k+1 ) -Y k u(θ k , Ξ k )) | G k ] = y (k) E (k) t (k) ,x (k) [(e c (k) (t (k+1) -t (k) ) u(t (k+1) , X (k+1) t (k+1) ) -u(t (k) , x (k) ))] = y (k) E (k) t (k) ,x (k) t (k+1) t (k) e c (k) (s-t (k) ) ∂ ∂t + L (k) + c (k) u(s, X (k) s ) ds .
Also,

E x [w k (Z k+1 -Z k ) | G k ] = y (k) E (k) t (k) ,x (k) f (k) t (k+1) t (k)
e c (k) s ds .

Under the hypothesis on the coefficients and the parallelepiped P (k) we have

E x [w k (Y k+1 u(θ k+1 , Ξ k+1 ) -Y k u(θ k , Ξ k ) + Z k+1 -Z k ) | G k ] = y (k) E (k) t (k) ,x (k) t (k+1) t (k) e c (k) (s-t (k) ) ∂ ∂t + L (k) + c (k) u(s, X (k) s ) + f (k) ds ≤ y (k) δ ≤ E x [δw k Y k G k ] ,
since the Y k 's (and so the y (k) 's) are positive. Hence, from [START_REF] Heath | A variance reduction technique based on integral representations[END_REF] and the Jensen inequality applied to | • |, we obtain

| E x [Y ν φ(θ ν , Ξ ν ) + Z ν ] -E x [ u(t, x)]| ≤ δ E x W ν ν-1 k=0 Y k .
As 0 < Y k ≤ e M θ k for k = 0, . . . , ν, we deduce [START_REF] Gobet | Efficient schemes for the weak approximation of reflected diffusions[END_REF].

4.4. The Monte Carlo error. In order to compute the solution u(t, x) of ( 11), we have constructed the estimator u(t, x) given by ( 15), whose variance is

Var b Px u(t, x) = 1 N Var b Px (W ν φ(θ ν , Ξ ν )Y ν + W ν Z ν ).
The Monte Carlo error depends on this variance s 2 = Var b Px u(t, x), since asymptotically for N → ∞ the true mean E x [ u(t, x)] lies in the interval [ u(t, x)-2s, u(t, x)+ 2s] with a confidence of 95.4 %.

We denote by P n the distribution of (Λ k ) k≥0 with respect to the real distribution of the exit time and position of the rectangles. In this case the weights are equal to 1. Any event Φ measurable with respect to (Λ k ) k≥0 satisfies P n

[Φ] = P x [W Φ]. We get thus Var b Px (W ν φ(θ ν , Ξ ν )Y ν + W ν Z ν ) = Ψ + Var b P n ( u(t, x)) with Ψ = E n [(W ν -1)(φ(θ ν , Ξ ν )Y ν + Z ν ) 2 ].
This shows that a good choice for the density of the exit time and position from the parallelepipeds is such that Ψ ≤ 0 is as small as possible. By the way, reducing the variance is a difficult task, and requires some automatic selection/optimization techniques, as explained in the introduction part.

In addition, the numerical experiments we performed up to now highlight another difficulty. W ν may take large values, and this implies meaningless values for u(t, x). That's why we suggest to keep track also of the empirical distribution, or at least of the variance of W ν .

In order to illustrate this, let us assume that the diffusion process X has no drift and that for the simulation, the right parallelepipeds we use are squares centered on the particle, and consider the same density for the exit time and position. By a scaling argument, the distribution of the weight w k at the k-th step does not depend on the size of the squares, so that the w k 's are independent and identically distributed under P x .

Let us fix an integer n such that ν ≥ n a.s. (for example, the minimal number of steps needed to reach an absorbing boundary). We set ξ i = log(w i ), so that W n = exp ( n i=1 ξ i ). As the ξ i are independent and identically distributed, let us note S n = n i=1 ξ i , then S n / √ n converges to some normal random variable χ with mean m and variance s 2 . For n large enough, the distribution of W n is close to the distribution of exp( √ nχ). We obtain, with the expression of the Laplace transform for the normal distribution, for j ∈ {1, 2},

E x [(W n ) j ] ≈ E x [exp(j √ nχ)] = exp mj √ n + n j 2 s 2 2 .
This leads us to the following approximation

Var b Px (W n ) ≈ exp 2m √ n + 2ns 2 -exp m √ n + n 2 s 2 ≈ exp(2ns 2 ) exp m √ n + 1 -exp m 2 √ n - 3n 2 s 2 ∼ n→∞ exp(1 + 2ns 2 ).
So, for large n, the variance of W n explodes, while E x [W n ] = 1 for any n ≥ 1.

In [START_REF] Glynn | Importance sampling for stochastic simulations[END_REF] (see also [START_REF] Glasserman | Monte Carlo methods in financial engineering[END_REF]), P. Glynn and D. Iglehart exhibit another argument that shows that the simulation performs badly if too much steps are used. 4.5. Population Monte Carlo. In order to overcome the explosion of the variance due to the weights one can use a population Monte Carlo method. This kind of method, also known as quantum Monte Carlo, sequential Monte Carlo, Green Monte Carlo, ... has been used from a long time in physical simulations (see for example [START_REF] Iba | Population Monte Carlo algorithms Trans[END_REF] for a brief survey) but also in signal theory, statistics, ... A probabilistic point of view is developed in the book [START_REF] Del Moral | Feynman-Kac formulae: Genealogical and interacting particle systems with applications[END_REF] of P. Del Moral.

In our case, instead of simulating the particles one after another, the idea is to keep track of the whole population of N particles (y (i) ) i∈{1,...N } with time and space coordinates (t (i) , x (i) ) and a weight w (i) according to the algorithm given below. Each particle has two possible states : still running or stopped. A particle is stopped either at the first time it reaches an absorbing boundary, or if its time is equal to the finite final time T . Otherwise, the particle is still running.

Algorithm 3. This algorithm computes an approximation of the quantity E x [f (T ∧ τ, X T ∧τ )] when X 0 = 0 by using a population of N particles.

(1) Set n = 0, n is the number of steps.

(2) For i from 1 to N set (a) (w , where w

(i) 0 , t (i) 0 , x (i) 0 ) = (0, 0, x). (3) Set S = ∅ and R n = {(w (i) 0 , t (i) 0 , x (i) 0 )} i=1,...,N . ( 4
(j)
n is the weight of the particle after n iterations. (ii) The particle is moved in time and space according to the exit time and position from a time-space parallelepiped that contains (t

(j) n , x (j)
n ). Its new position is denoted (t

(j) n+1 , x (j)
n+1 ) and its associated weight w

(j) n+1 . (iii) If t (j) n+1 = T or if x (j)
n+1 belongs to an absorbing boundary, then (w

(j) n+1 , t (j) n+1 , x (j) n+1 ) is added to the set S of stopped particles. Oth- erwise, it is added to R n+1 . (c) Increment n by 1. (5) Return 1 N i=1 w (i) N i=1 w (i) f (t (i) , x (i) )
when S = {(w (i) , t (i) , x (i) )} i=1,...,N .

As we need to keep track of the positions of all the particles, this algorithm is memory consuming. On the other hand, it avoids the multiplication of the weights. In addition, this algorithm can be modified in the following way: instead of using #R n particles at step n, it is possible to use N particles, and in this case, one has to keep track of the number of still running particles, and to multiply the weights by the proportion of still running particles. The algorithm stops when the proportion of still running particles is smaller than a given threshold. This approach can be used for example for long time simulation, or to estimate rare events, as for example in [START_REF] Cérou | Genetic genealogical model in rare event analysis[END_REF][START_REF] Del Moral | Feynman-Kac formulae: Genealogical and interacting particle systems with applications[END_REF][START_REF] Del Moral | Genealogical particle analysis of rare events[END_REF][START_REF] Lejay | Computing the first eigenelements of some linear operators using a branching Monte Carlo method[END_REF]. 4.6. Estimation of the number of steps. Let us consider now the estimation of the number of steps. In order to do this we will use the techniques employed in [START_REF] Mil'shteȋn | Solution of the first boundary value problem for equations of parabolic type by means of the integration of stochastic differential equations[END_REF][START_REF] Milstein | Simulation of a space-time bounded diffusion[END_REF][START_REF] Milstein | Stochastic numerics for mathematical physics[END_REF].

In Algorithm 2, we have constructed the Markov chain (Λ k ) k≥0 which is absorbed when reaching

Γ k = Γ d ∩ {T } × D.
For a function u on D, we set

P u(t, x) = E n [u(Λ 1 ) | Λ 0 = (t, x)] and A = P u(t, x) -u(t, x).
The operator A is the generator of a Markov chain.

Lemma 1. If T < +∞ and

E n [θ 1 | (θ 0 , Ξ 0 ) = (t, x)] -t ≥ γ, then E n [ν | (θ 0 , Ξ 0 ) = (t, x)] ≤ 1 + T -t γ .
Proof. Consider the problem

Av(t, x) = -g(t, x) on Q, u(t, x) = 0 on [0, T ] × Γ whose solution is u(t, x) = E n ν-1 k=0 g(Λ k ) .
We remark that if u and g are well chosen this equality gives a good estimate of

E n [ν]. Let V (t, x) be the function V (t, x) = (T -t)1 (t,x)∈Q . For (t, x) in Q, we have AV (t, x) = E n [V (θ 1 , Ξ 1 ) | (θ 0 , Ξ 0 ) = (t, x)] -(T -t) ≤ -γ. Hence T -t ≥ E n [ ν-1 k=0 γ | (θ 0 , Ξ 0 ) = (t, x)
] and the result follows easily.

Lemma 2. With the previous notations, for every L > 0 fixed, we have

sup x∈Q P n [ν ≥ L (θ 0 , Ξ 0 ) = (t, x)] ≤ (1 + T -t) exp(-cγL/(1 + T -t))
where c is a constant depending on γ, more precisely c converges to 1 as γ decreases to 0.

Proof. The proof follows from the one of Theorem 7.2 in [START_REF] Milstein | Simulation of a space-time bounded diffusion[END_REF].

Lemma 3. If T = +∞, Q is bounded and E n [|x + Ξ 1 + c| 2 ] ≥ γ > 0, where c is such that min x∈Q |x + c| ≥ C > 0. Then E n [ν] ≤ B 2 -C 2 B 2 -γ with B > max{γ, sup x∈Q |x + c|}.
Proof. Let us proceed as in [START_REF] Mil'shteȋn | Solution of the first boundary value problem for equations of parabolic type by means of the integration of stochastic differential equations[END_REF]. Choose a vector c such that min x∈Q |x+c| ≥ C > 0 and set

V (t, x) = B 2 -|x + c| 2 if (t, x) ∈ R + × Q, 0 otherwise. Thus for B 2 > γ, AV (t, x) ≤ |x + c| 2 -E n [|x + Ξ 1 + c| 2 | (θ 0 , Ξ 0 ) = (t, x)] ≤ B 2 -γ
and the result follows.

Numerical examples

We present in this Section some numerical examples in order to test our algorithm.

5.1.

Speeding up the random walk on squares algorithm. In [START_REF] Milstein | Simulation of a space-time bounded diffusion[END_REF] (see also [START_REF] Milstein | Stochastic numerics for mathematical physics[END_REF]), G.N. Milstein and M.V. Tretyakov proposed a method to simulate Brownian motions and solutions of SDEs by using the first exit time and position from a hyper-cube or a time-space parallelepiped with cubic space basis. A similar method has been previously proposed by O. Faure in his PhD thesis [START_REF] Faure | Simulation du mouvement brownien et des diffusions[END_REF]. This method is a variation of the random walk on spheres method. Some authors already used random walk on squares and rectangles by using the explicit expression of the Green function but without simulating the exit time (see for example [START_REF] Simonov | Random walk algorithm for estimating effective properties of digitized porous media[END_REF]). One of the main feature of our approach is the simulation of the couple of non-independent random variables (exit time, exit position) by means of real valued random variables. We have explained in [START_REF] Deaconu | A random walk on rectangles algorithm[END_REF] how to extend this approach to rectangles and the starting point everywhere in the rectangle. This approach is still using only one-dimensional distributions. However, by using symmetry properties, we can notice that it is simpler to deal with squares centered on the current position of the particle than with a rectangle.

Nevertheless, the computation may be time consuming. We are looking now to speed up the computations by using a simple density for the exit position.

Let us consider here the d-dimensional hypercube C = [-1, 1] d , and a fixed time T > 0 (possibly T = +∞). Let B be a d-dimensional Brownian motion. We set

τ B = inf{t > 0 | B t ∈ C}. Let W be a one-dimensional Brownian motion. We set τ W [-1,1] = inf{t > 0 | W t ∈ [-1, 1]}, R(t) = P 0 [τ W [-1,1] < t], r the density of τ W [-1,1] , S(t, y) = P 0 [W t < y | t < τ W [-1,1]
] and s(t, y) = ∂ y S(t, y) the density of

W t given {t < τ W [-1,1] }.
Let us note that we can easily switch from C = [-1, 1] d to any hypercube [-L, L] d after a scaling argument in space and time. Thus, from a numerical point of view, we need only to implement the required functions r, s, R and S on [-1, 1]. Analytical expressions for these distribution functions are easily deduced from the series presented in Appendix A.

To simulate the exit time and position from [0, T ]×C, we proceed in the following steps:

• Compute the probability

β = 1 -(1 -R(T )) d that τ B < T . • With probability β, decide if {τ B < T } happens or not. • If {τ B < T } happened, -For a realization U of a uniform random variable U on [0, 1), set τ B = R -1 (1 -(1 -Uβ) 1/d ).
which is a realization of τ B given {τ B < T }. -Choose with probability 1/2d an exit side (J, ε) and set ξ J = ε.

-For each i = 1, . . . , d, i = J, set χ i = √ U i , where the U i 's are d -1 independent realizations of uniform random variables on [0, 1). With probability 1/2, set ξ i = χ i -1 and with probability 1/2, set ξ i = 1-χ i .

-Compute the weight

w = 1 1 -R(τ B ) i=1,...,d,i =J s(τ B , ξ i ) χ i . • If {τ B ≥ T } happened, then -Set τ B = T . -For i = 1, . . . , d, set χ i = √ U i
, where the U i 's are d -1 independent realizations of uniform random variables on [0, 1). With probability 1/2, set ξ i = χ i -1 and with probability 1/2, set

ξ i = 1 -χ i . -Compute the weight w = 1 1 -β i=1,...,d s(T, ξ i ) χ i .
(τ B , ξ 1 , . . . , ξ d ) represent the first exit time and position from [0, T ] × C, and w is the associated weight.

For the random walk on squares we can also use the idea proposed in [START_REF] Milstein | Simulation of a space-time bounded diffusion[END_REF] and in [START_REF] Faure | Simulation du mouvement brownien et des diffusions[END_REF]. This leads to the following algorithm :

• Compute the probability

β = 1 -(1 -R(T )) d that τ B < T .
• With probability β, decide if {τ B < T } happens or not.

• If {τ B < T } happened, -For a realization U of a uniform random variable U on [0, 1), set

τ B = R -1 (1 -(1 -Uβ) 1/d ).
which is a realization of τ B given {τ B < T }. -Choose with probability 1/2d an exit side (J, ε) and set ξ J = ε.

-For each i = 1, . . . , d, i = J, draw ξ i according to the distribution of

B i τ B given τ B i > τ B , where τ B i = inf{t > 0 | B i ∈ [-1, 1]}. • If {τ B ≥ T } happened, then -Set τ B = T . -For i = 1, . . . , d, draw ξ i according to the distribution of B i τ B given τ B i > τ B .
In both cases, we use tabulated values for R and R -1 . In order to simulate B i t given τ B i > t, we use the rejection method proposed by O. Faure in [START_REF] Faure | Simulation du mouvement brownien et des diffusions[END_REF] for t ∈ [0. [START_REF] Lions | Stochastic differential equations with reflecting boundary conditions[END_REF][START_REF] Arouna | Algorithmes stochastiques et méthodes de Monte Carlo[END_REF]. Otherwise, we draw B i t by using the fact that it is equal to S -1 (t, U ) for some random variable U with uniform distribution on [0, 1). This is the method proposed by G.N. Milstein and M.V. Tretyakov in [START_REF] Milstein | Simulation of a space-time bounded diffusion[END_REF]. For t > 2, the latter method is more efficient than the previous one. For t < 0.2, the rejection method may give wrong results. For t close to 0.2, the rejection method can be up to 6 times faster than the inversion method, while for t close to 2, they are comparable in the computation time.

If the Brownian motion reaches the side labelled by (1, -1) first at time τ B , then in order to simulate B i t for i = 2, . . . , d we use a random variable with density φ(x) = 1 + x if x ∈ (-1, 0] and φ(x) = 1 -x if x ∈ [0, -1). In this case, the weights w are close to 1 as we see in Table 1, and the execution time is usually divided by 10. For T = 0.1, the variance of w is too high and leads to some instabilities. In this case, it is preferable to simulate the exact distributions of B T given {T ≤ τ B }.

Solving a bi-harmonic problem.

To test the validity of our approach with respect to other algorithms, we consider first an example borrowed in [START_REF] Milstein | Simulation of a space-time bounded diffusion[END_REF] (see also [29, p. 332]). Let D = [-1, 1] 2 , and consider the bi-harmonic equation ( 18)

     1 2 2 u(x) = 1, x ∈ D, u(x) = φ(x) on ∂D, 1 2 u(x) = ψ(x) on ∂D,

Method

T Mean of w Variance of w Time (s) Walk on squares 0. 1. Speeding up the random walk on squares: experiments with 1,000,000 samples are used. with

φ(x 1 , ±1) = 1 + x 4 1 12 , φ(±1, x 2 ) = 1 + x 4 2 12 , (19)
ψ(x 1 , ±1) = 1 + x 2 1 2 , ψ(±1, x 2 ) = 1 + x 2 2 2 . ( 20
) After setting v(x) = 1 2 u(x), (18)
may be transformed into the system

1 2 v(x) = 1 on D with u(x) = ψ(x) on ∂D, 1 2 u(x) -v(x) = 0 on D with u(x) = φ(x) on ∂D, whose exact solution is u(x) = x 4 1 + x 4 2 12 , v(x) = x 2 1 + x 2 2 2 .
By Itô formula, it is easy to show that

u(x) = E[φ(x + B τ B )] -E[τ B ψ(x + B τ B )] + 1 2 E[(τ B ) 2 ], v(x) = E[ψ(x + B τ B )] -E[τ B ],
where B is a 2-dimensional Brownian motion, and τ B is, as above, its first exit time from D.

Here, in contrast with the values presented in [START_REF] Milstein | Simulation of a space-time bounded diffusion[END_REF], we only need to use one square, since we are not forced to start from its center. We compare the results given by our algorithm (first lines) with the one given by the random walk on rectangles (second line). Each side is chosen uniformly with probability 1/4. The time is drawn by using an exponential random variable of parameter 1/(1 -εx i ) if (i, ε) is the exit side. The position is drawn uniformly on the exit side. This strategy corresponds in some sense to a "naive" and simple way to choose the exit time and position.

As we evaluate quantities of the form E[f (τ B , B τ B )], we report the quantities µ n ± 2σ n / √ n, where µ n is the empirical mean of f (τ B , B τ B ) with n samples, and σ n is the corresponding empirical standard deviation. The interval [µ n -

2σ n / √ n, µ n + 2σ n / √ n] represents the 95.5 % confidence interval for E[f (τ B , B τ B )
]. The estimation u(x) and v(x) of u and v for three points are given in Table 2. Table 2. Solution of the bi-harmonic equation: the first line of each row contains the results for our algorithm, the second line contains the results for the Random Walk on rectangles.

x n u(x) u(x) v(x) v(x) Time (s) (0.
Although a small numerical bias seems to appear, our algorithm provides results comparable with the random walk on rectangles method. The execution time is much more smaller than the one given by this method (also the one given by the random walk on squares, for which the simulation of one step takes less time, but where more steps are needed).

(x) = 1 if x ∈ S, 0 if x ∈ ∂D \ S.
We illustrate our method on the simple two dimensional domain D drawn in Figure 2 and we compute the value of p at the five points marked respectively by (a), (b), (c), (d) and (e) on Figure 2. To set up our algorithm, we use two rectangles as in Figure 3. The numbers marked on each side are the probabilities to reach each one of these sides.

In order to obtain the simulated exit time we draw an exponential random variable with parameter α where α is given by α = 1/( L i /2). The L i notes the length of the rectangle in the direction perpendicular to the boundary that the particle hits. We perform 100,000 samples, each computation takes around 1s on our computer (a MacBook 12", 2GHz with a code written in C). The values for p are given in Table 3. We perform a comparison with the value given by MATLAB/PDEtool, where ( 21) is solved by using a finite element method, and with the method of random walk on rectangles [START_REF] Deaconu | A random walk on rectangles algorithm[END_REF] which is exact (up to the Monte Carlo error), for such a domain. In this case, with a sample of size n, the variance of the empirical mean is p(x)(1 -p(x))/n.

We notice that the results given by our method are close to the one given by the finite element method. As one can expect, the random walk on rectangles (and any other methods that do not rely on importance sampling or variance reduction Point Import. Sampling Finite Element Walk on rect.

(a) 3.32 techniques) is not efficient to estimate the values of p(x) when they are of the same order as the standard deviation of the empirical mean.

In order to test the validity of our method for the simulation of rare events, we use the domain D as in Figure 4. The numerical results are reported in Table 4. p n is the empirical mean with n = 100,000 samples, and s 50 (p n ) is the empirical standard deviation computed over 50 realizations of p n . We obtain really good results even while computing small probabilities of order of magnitude 10 -10 .

Point

p n s 50 (p n ) Finite Element (a) 1.00 • 10 -10 2.

X t = x + t 0 σ(X s ) dB s with σ(x) = 1 1 2 sin(x 1 + x 2 ) 0 1
which is driven by a two-dimensional Brownian motion B. The process X is killed when it exits from the domain D which is represented in Figure 5. In order to simulate X, we use either an Euler scheme with a time step of 0.0025 or a (possibly modified) random walk on squares. The squares sides lengths are smaller than 2L with L = 0.05 (note that the time step of the Euler scheme corresponds to 0.05 2 , which is close to the average exit time of the square [0.1, 0.1] 2 ). As the diffusion moves in a bounded domain, we use to deal with the boundary condition and apply the technique proposed in [START_REF] Campillo | A Monte Carlo method without grid for a fractured porous domain model[END_REF]: If the distance between the position of the particle and the boundary is smaller than 2L, we choose the square such that one of its sides is included in the boundary when it is possible to do so.

Unless the coefficients of the SDE are constant, one needs to simulate many couples of exit times and positions from small squares, and the computational time becomes very large and is not competitive with respect to the Euler scheme. In addition, when the random walk on squares is coupled with importance sampling, the weights are growing quickly (See Section 4.4).

When the Euler scheme is used, we simply stop the algorithm when the particle leaves the domain D. This is a crude way to proceed, and some refinements can be done (see [START_REF] Gobet | Weak approximation of killed diffusion using Euler schemes[END_REF] for example). Note that the exit time is then over-estimated.

The idea is to mix the two methods and to use the Euler scheme inside the domain, and a random walk on squares when the particle is closed to the boundary. We improve thus the simulation as in this case the behavior of the particle is taken into account. In addition, it is possible by making a change of measure, to increase or to decrease the probability that the particle hits the boundary.

Our aim is here to increase the number of particles which are not killed before a given time T . When one side of the square is set on the boundary, we use a probability p that the particle reaches the side of the square that is opposite to the boundary, and q = (1 -p)/3 for any other side. We have thus a "repulsing" effect.

We use P 1 = {p = 0.7, q = 0.1} and P 2 = {p = 0.91, q = 0.03}. In order to avoid the explosion of the variance of the weight, we have used a limitation N max for the number of times this procedure is used. The variance of the weight for each time this procedure is used is 0.04 for the set P 1 and 0.34 for the set P 2 .

All the simulations are done with 100, 000 particles. The results are summarized in table 5. For T = 1, the proportion of particles still alive is of order 0.19% (using the Euler scheme without specific treatment on the boundary, we get an estimation of 0.33%, yet for a quicker simulation of 7s). With a population Monte Carlo method, we obtain an estimate of 0.17%, using the set P 1 and a running time of 126s. We see that our scheme allows one to get much more alive particles.

Appendix A. Appendix : How to get densities for different situations?

We present in this section analytical expressions for the density in different cases. Except for the case of a drifted Brownian motion with Dirichlet boundary condition at one endpoint of [-1, 1] and a Neumann boundary condition at the other endpoint of [-1, 1], we obtain two expressions, one which follows from the images method and the other one from the spectral decomposition. From a numerical point of view, the spectral decomposition gives rise to series that converge very quickly for large times. It is worth using the expressions given by the method of images for small times.

A.1. Brownian motion without drift. We are interested in this section in writing down some useful formulas for the calculations. Let us consider first the case of the standard one dimensional Brownian motion starting from x ∈ [-1, 1] which is killed or reflected when hitting the boundaries -1 or 1. We shall write D for Dirichlet condition on the boundary and N for Neumann condition, which of course correspond to killing and respectively reflection. Furthermore we shall note, for example, p DN (t, x 1 , x 2 ) the density of the Brownian motion on [-1, 1] killed when hitting -1 and reflected on 1, more precisely the order in the indices indicates the boundary condition in -1 and 1 respectively. A.1.1. Reflected Brownian Motion on [-1, 1]. Let p N N (t, x 1 , x 2 ) denote the probability density function of a Brownian motion at time t, starting from x 1 and reflected at -1 and 1. By using the method of images we get the following formula for the transition density: Table 5. Simulations of the proportions (in %) of the particles reaching a given part of the boundary as well as the surviving particles at time T . We write "unstable" in the column for the variance of weights when the mean of the global weights is rather far from 1.

p N N (t, x 1 , x 2 ) = 1 √ 2πt ∞ n=-∞ e -(x 1 -x 2 +4n) 2 2t + e -(x
The spectral representation of this density writes: where τ DD = inf{t ≥ 0; B t / ∈ [-1, 1]}. Then, by the images' method we have:

p N N (t,
p DD (t, x 1 , x 2 ) = 1 √ 2πt ∞ n=-∞ e -(x 1 -x 2 +4n) 2 2t -e -(x 1 +x 2 +4n+2) 2 2t
.

For the law of the exit time we get:

P x 1 [τ DD ∈ dt] = 1 √ 2πt 3 ∞ n=-∞
(-1) n (x 1 + 2n + 1)e -(x 1 +2n+1) 2 2t

dt.

The spectral representation can be also written and yields:

p DD (t, x 1 , x 2 ) = ∞ n=1
e -n 2 π 2 8 t sin nπ 2 (x 1 + 1) sin nπ 2 (x 2 + 1) .

The law of the exit time is given by: These expressions may be found for example in [START_REF] Beck | Heat conduction using Green's functions[END_REF] or in [START_REF] Milstein | Simulation of a space-time bounded diffusion[END_REF].

P x 1 [τ DD ∈ dt] =
A.1.3. Mixed boundary conditions for the Brownian Motion on [-1, 1]. We give here explicit solutions for the Brownian motion killed on -1 and reflected on 1. Let p DN (t, x 1 , x 2 ) denote the probability density function of a Brownian motion at time t, starting from x 1 and killed when it hits -1 and reflected on 1. Then, by the images' method, one gets

p DN (t, x 1 , x 2 ) = 1 √ 2πt ∞ n=-∞ (-1) n e -(x 1 -x 2 +4n) 2 2t -e -(x 1 +x 2 +4n+2) 2 2t
.

Let us denote also by τ DN the killing time for the Brownian motion on [-1, 1) killed on -1 and reflected on 1. Hence,

P x 1 [τ DN ∈ dt] = 1 √ 2πt 3 ∞ n=-∞
(-1) n (x 1 + 4n + 1)e -(x 1 +4n+1) 2 2t dt.

The spectral representation can be also written and yields:

p DN (t, x 1 , x 2 ) = ∞ n=0
e -(2n+1) 2 π 2 32 t sin (2n + 1)π 4 (x 1 + 1) sin (2n + 1)π 4 (x 2 + 1) .

Then we get from the spectral representation the law of this exit time: The dual situation (reflection on -1 and absorption on 1) can be obtained easily by the transformation p N D (t, x 1 , x 2 ) = p DN (t, -x 1 , x 2 ).

P x 1 [τ DN ∈ dt] =
These expressions may be found for example in [START_REF] Beck | Heat conduction using Green's functions[END_REF]. This formula can be obtained also from the results in D. Veestraeten [START_REF] Veestraeten | The conditional probability density function for a reflected Brownian motion[END_REF].

By the spectral method (see for example [START_REF] Linetsky | On the transition densities for reflected diffusions[END_REF]), we have after some calculations: A.2.2. Brownian motion with drift µ on [-1, 1] killed at the boundary. We keep the same notations as before. By using classical properties of the Brownian motion and the results from Milstein and Tretyakov [START_REF] Milstein | Simulation of a space-time bounded diffusion[END_REF] we have the following transformation:

p µ N N (t,
p µ DD (t, x 1 , x 2 ) = e µ(x 2 -x 1 )-µ 2 t 2 p DD (t, x 1 , x 2 ).

Then, by the images' method, p µ DD (t, x 1 , x 2 ) = e µ(x 2 -x 1 )-µ 2 t .

We write down both distribution and density for the exit time. The distribution writes: .

P x 1 [τ µ DD < t] = 1 -
The spectral representation can be also written and yields:

p µ DD (t, x 1 , x 2 ) = e µ(x 2 -x 1 )-µ 2 t 2 ∞ n=1
e -n 2 π 2 8 t sin nπ 2 (x 1 + 1) sin nπ 2 (x 2 + 1) .

4. 2 .

 2 The hypotheses. Let us consider a domain Q in R + × R d . For the sake of simplicity, we assume that Q is the cylinder [0, T ] × D (with possibly T = +∞), where D is an open, bounded domain of R d with piecewise smooth boundary. Let us consider a function a with values in the space of d×d-symmetric matrices which is continuous on D and everywhere positive definite, together with some functions b

  and Γ n ∩Γ d = ∅. This means that the boundary of [0, T ]×∂D is split in two distinct parts, where either the Dirichlet or the Neumann boundary conditions hold. More precisely a side of a parallelepiped in P D contained in ∂Q is either from Γ n or from Γ d .

  belongs to the basis of P (s is possibly infinite if for example T = +∞ and the coefficients are time-inhomogeneous). On P (k) , consider the differential operator L (k) as well c (k) and f (k) which approximate L, c and f as in (H5). (b) Draw a realization of a random variable (θ k+1 , Ξ k+1 ) with values in ({s} × U ) ∪ ((θ k , s) × ∂U ) and compute its associated weight w k as shown in Sections 2 and 4.1 by considering the exit time and position from the parallelepiped P (k) . (c) Compute W k = W k-1 w k and

  ) While the set R n of still running particles at step n is non-empty do (a) Set R n+1 = ∅. (b) Do #R n times the following operations (i) Pick a still running particle of index j at random according to a family of discrete probability distribution p j = w (j) n k index of particles in Rn w (k) n

5. 3 .

 3 Estimation of rare events: computing hitting probabilities. Let us consider the following problem: what is the probability p(x) that starting from a point x in a domain D a Brownian motion reaches a part S of the boundary ∂D? It is well known that p is solution of the Dirichlet problem (21) 1 2 p(x) = 0 on D and p

Figure 2 .

 2 Figure 2. A simple domain D.

Figure 3 .

 3 Figure 3. Decomposition of D into rectangles.

Figure 4 .

 4 Figure 4. A simple domain D .

Figure 5 .

 5 Figure 5. Domain D with the label of the sides and the starting point.

A. 2 .

 2 Brownian motion with drift µ. As in the previous part of the appendix we consider here the case of the Brownian motion with drift on the interval [-1, 1] which is killed or reflected on -1 and 1. If we note by p L,µ .. (t, x 1 , x 2 ) the law of the process with drift µ and living on [-L, L] and p µ .. (t, x 1 , x 2 ) the corresponding law on [-1, 1] then, by the properties of the Brownian motion we have: p L,µ .. (t, x 1 , x 2 in the indices can take the value D for a Dirichlet condition or N for a Neumann condition, as previously noted. A.2.1. Brownian motion with drift µ reflected on [-1, 1]. We keep the same notations as before. The use of images' method gives the following representation of the density:

e 2 2t-µe 2µx 2 ∞

 22 -2µx 1 e -µ(4n+2) e -(x 1 +x 2 -µt+4n+2) n=-∞ e µ(4n+2) erfc x 1 +x 2 +µt+4n+2 √ 2t .

2 √ 2πt 3 ∞

 23 2µx 1 +µ(4n+2)) erfc x 1 -µt+4n+1 √ 2t -erfc x 1 -µt+4n+3 √ 2t ,while for the density we obtain:P x 1 [τ µ DD ∈ dt] = e -µx 1 -µ 2 t n=-∞ e -µ (x 1 + 4n + 1)e -(x 1 +4n+1) 2 2t--e µ (x 1 + 4n -1)e -(x 1 +4n-1) 2 2t

Table 3 .

 3 Computation of p(x) at given points of D.

		• 10 -6	3.39 • 10 -6	0.00
	(b)	2.31 • 10 -5	2.23 • 10 -5	1.00 • 10 -5
	(c)	1.70 • 10 -4	1.77 • 10 -4	1.90 • 10 -4
	(d)	4.43 • 10 -5	4.64 • 10 -5	3.00 • 10 -5
	(e)	2.79 • 10 -3	2.81 • 10 -3	2.36 • 10 -3

Table 4 .

 4 3 • 10 -11 1.15 • 10 -10 (b) 7.67 • 10 -10 1.6 • 10 -10 8.13 • 10 -10 (c) 5.19 • 10 -9 1.0 • 10 -9 6.61 • 10 -9 (d) 1.31 • 10 -9 2.8 • 10 -10 1.73 • 10 -9 (e) 2.27 • 10 -7 4.9 • 10 -8 2.29 • 10 -7 Computation of p(x) at given points of D . 5.4. Simulation of SDEs: approximation close to the boundary. Let us consider the two-dimensional SDE

  x 1 , x 2 ) = 1 2 +These expressions may be found for example in[START_REF] Beck | Heat conduction using Green's functions[END_REF].A.1.2. Killed Brownian Motion on [-1, 1]. Let p DD (t, x 1 , x 2 ) denote the probability density function of a Brownian motion at time t, starting from x 1 and killed when it exits from the interval [-1, 1]. That is p DD (t, x 1 , x 2 )dx 2 = P x 1 [B t ∈ dx 2 ; t < τ DD ]

	∞ n=1	e -n 2 π 2 8 t cos	nπ 2	(x 1 + 1) cos	nπ 2	(x 2 + 1) .

  x 1 , x 2 ) = 2µe 2µx 2 e 2µ -e -2µ + e µ(x 2 -x 1 )-µ 2

	× ×	πn 2 πn 2	cos cos	πn 2 πn 2	∞ n=1 (x 1 + 1) + µ sin e -n 2 π 2 8 t 2 t µ 2 + n 2 π 2 4 πn (x 1 + 1) 2 (x 2 + 1) + µ sin πn 2 (x 2 + 1) .

The authors are grateful to the referee for his helpful remarks and suggestions.

This author has been partially supported by the GdR MOMAS (funded by ANDRA, BRGM, CEA, CNRS, EDF and IRSN)..

The distribution of the exit time is given by:

(e -µ -(-1) n e µ ) 2nπ 4µ 2 +n 2 π 2 e -n 2 π 2 8 t sin nπ 2 (x

and

In a more detailed expression we can write this on the form:

These expressions may be found for example in [START_REF] Beck | Heat conduction using Green's functions[END_REF] or in [START_REF] Milstein | Simulation of a space-time bounded diffusion[END_REF].

A.2.3. Mixed boundary condition for the Brownian Motion on [-1, 1] with drift µ.

The aim is to express some explicit solutions for the Brownian motion killed on -1 and reflected on 1. We are solving now the following eigenvalue problem:

We can remark first that if ϕ λ is an eigenfunction for the eigenvalue λ for the preceding PDE then λ is negative. We associate to this problem the corresponding second degree equation and note ∆ = µ 2 + 2λ. After a detailed calculus with respect to the sign of ∆ we can express the countable set of eigenfunctions and eigenvalues with respect to the possible values of µ. There are three different situations, expressed in the Table 6 (see the book [START_REF] Pinsky | Partial Differential Equations and Boundary-Value Problems with Applications[END_REF] for example). The density p DN (t, x 1 , x 2 ) is obtained by using the spectral expansion p DN (t,

The density q DN (t, x 1 ) of the exit time is also expressed by


Stefan Manegold

Ioana Manolescu

Loredana Afanasiev

Jianling Feng

Gang Gou

Marios Hadjieleftheriou

Stavros Harizopoulos

Panos Kalnis

Konstantinos Karanasos

Dominique Laurent

M Lupu

N Onose

V Sans

P Senellart

T Wu

D Shasha

Repeatability & Workability Evaluation of SIGMOD 2009

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.

THE GOAL

On a voluntary basis, authors of accepted SIGMOD 2009 papers provided their code/binaries, experimental setups and data to be tested for: repeatability of the experiments described in each accepted paper;

workability of the software by running different/more experiments with different/more parameters than shown in the accepted paper; by a repeatability/workability committee (which we call the RWC ), under the responsibility of the repeatability/workability editors-in-chief (which we call the RWE ).

THE PEOPLE

The RWE were Ioana Manolescu and Stefan Manegold. The 2009 RWC consisted of the other authors of this paper, along with D. Shasha. * http://homepages.cwi.nl/ ~manegold/ SIGMOD-2009-RWE/

THE PLAN

Several lessons learned from the first repeatability evaluation with SIGMOD 2008 [START_REF] Manolescu | The repeatability experiment of SIGMOD[END_REF] led us to improve and extend the process. The following paragraphs describe the details.

Accepted papers, only

The SIGMOD 2009 repeatability & workability committee evaluated accepted papers only. The primary reason for this change was to reduce the workload for the evaluation by avoiding evaluation of papers that would eventually not be accepted. A second reason was that authors had commented that they wouldn't mind the extra work of preparing their repeatability & workability submission once their papers were accepted.

Adapted schedule

Focussing on accepted papers only required an adaptation of the general schedule for the repeatability & workability evaluation. After the SIGMOD 2009 program committee had announced the accepted papers, the contact authors of all accepted research papers were personally invited via email to prepare and submit their experiments including code, data sets and detailed instructions. This later start of the evaluation did not leave enough time to finish the evaluation before the camera ready deadline, thus preventing authors from mentioning the result of the evaluation in the final versions of their papers. In fact, the evaluation was completed just before the conference to give the authors the chance to mention the results in their presentations at SIG-MOD 2009.

Refined submission method

In contrast to the push-based submission in 2008 via upload to a FTP server, the submission in 2009 was pull-based. Authors were asked to make their submissions available for download by the RWE. This helped to avoid problems with uploading large (sometimes tens of gigabytes) submissions to a single FTP server. The RWE then made the submissions available for download to the assigned reviewers.

Refined submission instructions

To give the reviewers some information to better plan their evaluation, the authors were asked to include in their submission information the length of time their experiments were expected to run. In addition, in order to facilitate the workability evaluation, the authors were asked to extend their repeatability instructions with suggestions as to how to extend their experiments beyond the contents of their paper. Possibilities ranged from explanations of how to use different data sets, query work loads, tuning and/or configuration parameters to compilation, and installation instructions for alternative hardware/software environments.

Refined reviewing process

As in 2008, the assignment of papers to reviewers was mainly determined by the need to match the papers' hardware and software requirements with the reviewers resources. Of course, (potential) conflicts of interest were avoided. In contrast to 2008, each paper was assigned two reviewers: a primary reviewer to do the actual repeatability and workability evaluation, and a secondary reviewer as back-up and to double-check the primary reviewers report.

Author-reviewer-interaction

The 2008 experiences revealed that successful repeatability evaluation can be hindered or even prevented by minor problems in setting up and running the experiments due to missing details in the provided instructions. To solve this problem, the 2009 effort provided a web-based anonymous communication channel to allow interaction between authors and reviewers to resolve problems as early as possible. All communication has been archived. With standard WIKI or BLOG software either not providing convenient and effective means for anonymous peer-to-peer communication, we wrote a PHP script to efficiently provide the basic functionality required.

THE PROCESS

After the announcement of the accepted papers, the contact authors of all 64 accepted research papers were invited by email to prepare and submit their contribution. By the (extended) deadline of April 22 2009, 19 authors had provided their contribution. The remaining 45 authors chose not to reply at all. In contrast to 2008, authors were not asked to provide an explanation why they could not submit their code, data and experiments for evaluation.

Each RWC member was assigned three papers, either two for primary review and one for secondary review, or one for primary and two for secondary review. Assigned reviewers met the software and hardware requirements of the experiments though sometimes at significant effort. For example, some reviewers installed extra software or even (re-)install complete machines. In one case, the reviewer's group (re-)installed a 40-node Linux PC cluster to repeat a scaleddown version of experiments that were originally run on a 100-node cluster.

In nearly all cases, the anonymous web-based communication channel between authors and reviewers was successfully used to resolve problems ranging from missing gnuplot files to insufficient specification of the versions of required software. In only two cases was the discussion insufficient to solve all problems, resulting in only a partial repeatability evaluation for those papers.

The reviewing process stretched over the complete two month period, with the final reviews being finished only the day before SIGMOD 2009 started. The long time partly due to (i) the installation of extra hardware and software, as well as configuration work required; (ii) author-reviewer communication to solve initial problems; and (iii) experiments that took several days or even weeks to run. Not all authors provided hints how to modify and/or extend their experiments for workability evaluation. In all but 5 cases, the reviewers managed to find their own ways to modify/extend the respective experiment to assess their workability. Even when workability suggestions were provided, the reviewers volunteered to go beyond the authors' suggestions.

THE RESULTS

Overall, the results of the evaluation for the 19 submissions were rather positive:

• For 10 papers, the presented experiments could be fully repeated and workability was confirmed.

• For 1 paper, repeatability was fully confirmed and workability was mostly confirmed.

• For 4 papers, all original experiments were successfully repeated, but workability was not, mostly due to missing or insufficient instructions on how to modify the original setup conveniently.

• For 1 paper, the experiments were mostly repeated, but workability could not be evaluated.

• For 1 paper, the repeatability evaluation was successful, but the workability evaluation failed.

• For 2 papers, major technical problems could not be solved within the two months reviewing period, preventing most or all of the repeatability and workability evaluation.

The authors were informed before the conference about the results for their papers, and thus given the opportunity to mention the results during their presentation at SIGMOD 2009.

THE ASSESSMENT

With many of the lessons learned from the 2008 effort, the 2009 repeatability and workability evaluation went much smoother than the previous round. In particular focussing on accepted papers only (an idea suggested by Donald Kossmann), pull-based submission, and the web-mediated discussions between reviewers and authors to solve minor technical problems proved to be successful.

Though creating a higher workload for the reviewers, the newly introduced workability evaluation (when successful) gave even more credibility to the authors than a pure repeatability evaluation.

The unexpectedly low submission rate appears to be due to the fact that the authors were not aware of the SIG-MOD 2009 repeatability & workability evaluation by the paper submission deadline. Due to several delays and issues, the SIGMOD 2009 repeatability & workability evaluation was not announced in any call for papers, nor mentioned on the SIGMOD 2009 web site. Several personal communications with authors during the conference revealed that many authors were caught by surprise when invited to submit repeatability material for their accepted papers, or were simply not sure how "official" the evaluation was. In other words, there was probably insufficient publicity around the SIGMOD 2009 repeatability & workability evaluation.

While serving its primary purpose, the PHP script for the reviewer-author-communication could be improved. Not being a standard tool, the "look-and-feel" was considered "unusual" and the automatic email notification of new postings did not always work reliably.

Given the diversity of the papers and their experiments, the reviewers were not given a strict format for their reviews, but rather allowed to freely determine the format, structure and content of their reviews themselves, to accommodate the process they followed as well as their findings and final verdict.

RECOMMENDATIONS

Here are some lessons for 2010 and beyond:

• Publicize the effort well in advance of the submission deadline.

• Improve the author instructions, in particular to ask more explicitly for workability instructions. More generally, collecting, improving and disseminating guidelines for the preparation of repeatable experiments requires more work in the community; tutorials such as [START_REF] Manolescu | Performance Evaluation in Database Research: Principles and Experience[END_REF] are a step in this direction.

• Improve the reviewer guidelines to unify the results. Given the diversity of the experiments, this is not a trivial task, as any guidelines and/or format still need to leave sufficient room for all cases.

• Improve the visibility of the evaluation results. The SIGMOD PubZone server [START_REF]PubZone: scientific publication discussion forum[END_REF] is a promising tool for this purpose.

• Improve the software support for author-reviewer discussions. With respect to the last item above, we are currently considering the extension of the MyReview conference management tool [START_REF]The MyReview Conference Management System[END_REF] to accomodate the specific needs of our process. The main feature we need, and which is not yet supported by MyReview and other similar tools, is the possibility for reviewers and authors to exchange an unbounded number of messages, over the whole period of reviewing (as opposed to one single exchange, at a specific point in the process, as currently used for conferences such as ACM SIGMOD, and supported by the Microsoft Research Conference Management Tool [START_REF]The Microsoft Research Conference Management Tool[END_REF]).

SUMMARY

Our community can be justly proud of its practical impact over the years. This is due to a confluence of good ideas with good engineering. Repeatability and Workability help to ensure the validity of good ideas and provide paradigms and platforms for good engineering.

Acknowledgments We are thankful to a set of external referees: Tuyet Tram Dang Ngoc, Tao-Yuan Jen and Dan Vodislav from ETIS -Université de Cergy-Pontoise; Guangxishui Yang, Yi Zhou, and Junyu Liu from Sun Yat-Sen University; Jeffrey Xu Yu, Lijun Chang and Lu Qin from the Chinese University of Hong Kong; and YongChul Kwon from the University of Washington, Seattle.

Appendix

The following SIGMOD 2009 accepted papers passed all the repeatability tests, and also passed workability test: