S Herrmann
email: samuel.herrmann@u-bourgogne.fr

C Zucca

Exact simulation of first exit times for one-dimensional diffusion processes

Keywords: Exit time, Brownian motion, diffusion processes, Girsanov's transformation, rejection sampling, exact simulation, randomized algorithm, conditioned Brownian motion. 2010 AMS subject classifications: primary 65C05; secondary: 65N75, 60G40

The simulation of exit times for diffusion processes is a challenging task since it concerns many applications in different fields like mathematical finance, neuroscience, reliability... The usual procedure is to use discretization schemes which unfortunately introduce some error in the target distribution. Our aim is to present a new algorithm which simulates exactly the exit time for one-dimensional diffusions. This acceptance-rejection algorithm requires to simulate exactly the exit time of the Brownian motion on one side and the Brownian position at a given time, constrained not to have exit before, on the other side. Crucial tools in this study are the Girsanov transformation, the convergent series method for the simulation of random variables and the classical rejection sampling. The efficiency of the method is described through theoretical results and numerical examples.

Introduction

First exit time distributions for general stochastic processes are of prime importance in many contexts. In mathematical finance, they permit to appreciate the risk of default for a given path-dependent option; in neuroscience, they characterize for instance the interspike time distribution... Since diffusion processes (solutions of stochastic differential equations) form an important family of stochastic processes, we aim to describe quite precisely their exit times. Unfortunately a simple and explicit expression of their distribution is not generally available which leads us to consider numerical approximations. Our task is to point out an algorithm which permits to simulate the exit time of the diffusion. This objective was already concerned by several studies introducing a discretization scheme for the corresponding stochastic differential equation. Most of them are based on improvements of the classical Euler scheme (see for instance [START_REF] Broadie | A continuity correction for discrete barrier options[END_REF], [START_REF] Gobet | Stopped diffusion processes: boundary corrections and overshoot[END_REF], [START_REF] Gobet | Weak approximation of killed diffusion using Euler schemes[END_REF]) which essentially consist in reducing the error stem from the approximation procedure. Let us also note another point of view which consists in approximating the probability density function of the exit time and therefore to approximate the solution of an integral equation [START_REF] Sacerdote | Joint densities of first hitting times of a diffusion process through two time-dependent boundaries[END_REF].

Our approach is completely different since we emphasize an exact simulation procedure: the distribution of the random outcome of the algorithm is identical to the distribution of the first exit time of the diffusion process. For such simulation methods based on an acceptance-rejection sampling, the challenge is to describe and reduce if possible the time consumption of the simulation. Beskos & Roberts, in their founding work [START_REF] Beskos | Exact simulation of diffusions[END_REF], already introduced the exact simulation for diffusion paths on some fixed time interval. Meanwhile several modifications of this algorithm have been proposed [START_REF] Beskos | Retrospective exact simulation of diffusion sample paths with applications[END_REF][START_REF] Beskos | A factorisation of diffusion measure and finite sample path constructions[END_REF][START_REF] Jenkins | Exact simulation of the sample paths of a diffusion with a finite entrance boundary[END_REF]. The basic idea of the rejection sampling is to sequentially observe independent random objects generated according to a proposal distribution until a condition is satisfied. That means that each object is accepted or rejected according to a certain probability weight.

For diffusion processes on a fixed interval [0, T ], the proposed object pointed out by Beskos & Roberts is the Brownian paths (easy to simulate) and the rejection weight can be computed using the Girsanov transformation which essentially requires the simulation of Brownian bridges. In a previous work [START_REF] Herrmann | Exact simulation of the first-passage time of diffusions[END_REF], the authors proposed a similar procedure in order to exactly simulate the first passage time of a one-dimensional diffusion through a given threshold: the proposal object is then the Brownian first passage time (inverse gaussian random variable) and the rejection weight requires the simulation of Brownian bridges conditioned to stay under a given threshold (Bessel processes). In order to adapt such a procedure to exit times, it suffices to generate Brownian exit times and to use the rejection weight suggested by the Girsanov transformation. Unfortunately this random weight requires the simulation of a Brownian bridge conditioned to stay in a given interval, which corresponds to a SDE with singular coefficients: there is no way to exactly simulate such a stochastic paths. Such an intuitive generalization leads therefore to a deadlock.

That's why we aim to present a quite different rejection sampling based on a similar concept (Girsanov's transformation) and avoiding the simulation of the whole conditioned Brownian paths. Let us consider the stochastic process (X t , t ≥ 0), solution of the SDE:

dX t = µ(X t )dt + σ(X t )dB t , X 0 = x ∈ (a, b),
where (B t , t ≥ 0) stands for the standard one-dimensional Brownian motion. We denote by τ a,b the first time the diffusion exits from the interval [a, b]: τ a,b (X) := inf{t > 0 :

X t / ∈ [a, b]}. (0.1)
The aim of the study is to propose an efficient algorithm in order to simulate the first exit time τ a,b (X). Under suitable conditions, the Lamperti transform permits to reduce the area of investigation to the constant diffusion case: σ(x) ≡ 1 providing to change the interval [a, b]. That's why we shall only focus our attention on the diffusion process:

dX t = µ(X t )dt + dB t , X 0 = 0, (0.2)
for t ≤ τ a,b (X) with a < 0 < b. The exact algorithm presented here (Section 3) is essentially based on the Girsanov transformation which permits to relate the diffusion (0.2) to a standard one-dimensional Brownian motion. Let us roughly describe the crucial tools needed to execute such an algorithm. The proposal random variable is the Brownian exit time of the interval [a, b], denoted T prop .

The weight used in order to accept or reject this proposal depends both on the conditional distribution of the Brownian paths (B t ) given that τ a,b (B) = T prop and on a Poisson process defined on the time axis R + which is independent of the Brownian motion (B t ) and whose events occur at time T 1 , . . . , T n , . . . Then conditionally to τ a,b (B) = T prop , only the values of B T k , for any T k ≤ T prop , and B Tprop (dots in Figure 1) are involved in the specific conditions associated to the acceptance of T prop . The exact simulation of τ a,b (X) therefore requires to describe:

1. The conditional distribution of B t for a given t > 0, given that the first exit of the interval τ a,b (B) is larger than t (Section 1). Using the classical convergent series method, we propose an algorithm denoted by CON-DITIONAL_DISTR (x,[a,b],t) for the simulation of the constrained Brownian motion B t with initial condition x.

The distribution of the

E[F (B t , t ≤ τ a,b (B))|B 0 = y] = E F b -a 2 B t + a + b 2 , t ≤ τ -1,1 (B) B 0 = x ,
where x = 2y-a-b b-a . So if we denote B x the Brownian motion starting in x, we deduce easily the following distribution identity

τ a,b (B y ), B y τ a,b (B y ) (d) = (b -a) 2 4 τ -1,1 (B x ), b -a 2 B x τ-1,1(B x ) + a + b 2 . (1.1)
In other words, it suffices to study precisely the first Brownian exit problem from the normalized interval [-1, 1] with the initial condition

B 0 = x ∈ [-1, 1].
For notational convenience, we denote by

τ x = inf{t > 0 : B x t / ∈ I} with I = [-1, 1]. (1.2)
Let us first emphasize a classical result on the exit position: since (B t , t ≥ 0) is a martingale, the identity function s(x) = x corresponds to the scale function of the Brownian motion. The optimal stopping theorem permits therefore to describe the probability to exit from the interval on one particular side:

P(B x τ = 1) = s(-1) -s(x) s(-1) -s(1) = x + 1 2 , ∀x ∈ [-1, 1]. (1.3)

Distribution of B t given τ > t

Let us now focus our attention to the distribution of the constrained Brownian motion. We introduce the killed Brownian motion: as soon as the process hits the boundary of the interval, it jumps to a cemetery point x * / ∈ [-1, 1]. Indeed let us denote by p(t, x, dy) := P(B x t ∈ dy, τ > t) the transition probability of the Brownian motion joint with τ > t. It can be written as p(t, x, y) = q(t, x, y)P(τ x > t) where q(t, x, y) is the transition probability of the Brownian motion conditioned to stay in the interval [-1, 1] till time t. For any non negative or bounded measurable function ψ, we get

E[ψ(B x t ), τ x > t] = 1 -1 ψ(y)p(t, x, y)dy, x ∈ [-1, 1].
It is well known (see, for instance [START_REF] Itô | Diffusion processes and their sample paths[END_REF] Section 4.11 or [8] Section 5.7) that (t, y) → p(t, x, y) satisfies an Initial-Boundary value problem associated to the heat equation:

       ∂p ∂t (t, x, y) = 1 2 ∂ 2 p ∂y 2 (t, x, y) p(t, x, y) → δ x (y) as t → 0, p(t, x, y) = 0 for y ∈ {-1, 1}, t > 0.
(1.4) Milstein and Tretyakov [START_REF] Milstein | Simulation of a space-time bounded diffusion[END_REF] recall that the solution of (1.4) takes two different expressions. The first one is obtained by the method of images and the second one is based on a spectral decomposition of the heat equation. We introduce the standard gaussian pdf and cdf : φ

(x) = 1 √ 2π e -x 2 /2 and Φ(x) = x -∞ φ(y)dy. Then p(t, x, y) = 1 √ t +∞ n=-∞ φ((x -y -4n)/ √ t) -φ((x + y -2 -4n)/ √ t) (1.5)
which is a convenient formula for small values of t. For large times t, we would prefer the second formula:

p(t, x, y) = n≥1 exp - n 2 π 2 t 8 sin nπ 2 (x + 1) sin nπ 2 (y + 1) . (1.6)
Let us fix (t, x): since we get an explicit expression of p(t, x, y), we can point out exact simulation algorithms based on the acceptance/rejection method even if p(t, x, y) is not a probability density function. Indeed it suffices to divide by P(τ x > t) in order to obtain such a density (cf. Section 2.1). The method applied here is a particular acceptance/rejection method presented in [START_REF] Devroye | Nonuniform random variate generation[END_REF] as the convergent series method.

Inequalities related to the series expansion

Let us consider f a non negative function satisfying a series expansion:

f (y) = n≥0 f n (y).
We assume that I(f ) := R f (y) dy < ∞. The convergent series method permits to simulate a random variable with the corresponding probability distribution function f /I(f ). It requires two important features:

• the existence of both a pdf h and a constant κ > 0 such that f (y) ≤ κh(y), y ∈ R.

(1.7)

• the existence of a positive sequence (r n ) n≥0 which converges toward 0 and satisfies the following reminder upper-bound

|R n (y)| ≤ r n , ∀y ∈ R where R n (y) := k≥n+1 f k (y). (1.8)
CONVERGENT SERIES METHOD

Initialization. N c = 0.
Step 1. Generate a random variable Y with density h.

Step 2. Generate a random variable U uniformly distributed on [0, 1] and define W = κU h(Y ).

Step 3. Initialization. N l = 0, S = 0, Test = 0.

Step 4. While (Test = 0) do:

• S ← S + f N l (Y ) • Test = 1 {|S-W |>r N l } • N l ← N l + 1 and N c ← N c + 1.
Step 5. If W ≤ S then X = Y otherwise go to Step 1.

Outcome: the random variable X with density f /I(f ) and the global number of terms of the series expansions N c used.

In order to simulate q(t, x, y), i.e. the conditional distribution of B t given t < τ , we need to choose the density h, to point out some constant κ and the sequence (r n ) n≥0 for each explicit expression of f given by (1.5) and (1.6). Let us just note that it is challenging to find the smallest constant κ as possible since the number of iterations of Step 1 in order to get an outcome is geometrically distributed with average κ/I(f ).

Some comments on the series expansion (1.5)

We note that (1.5) is not an alternating series and can be rewritten as p(t, x, y) = a 0 (t, x, y) + ∞ n=1 (a n (t, x, y) + a -n (t, x, y)), where a n (t, x, y)

= 1 √ t (φ((x -y -4n)/ √ t) -φ((x + y -2 -4n)/ √ t)).
We first observe that a -n (t, x, y) < 0 and a n (t, x, y) > 0, ∀n ≥ 1, ∀(x, y) ∈ [-1, 1] 2 .

(1.9)

Moreover both increments of the function φ are computed on intervals whose length does not depend on the variables x and n:

(x + y -2 -4n) -(x -y -4n) = (x + y -2 + 4n) -(x -y + 4n) = 2y -2.
For n large enough, the increments are decreasing and therefore a n (t, x, y) + a -n (t, x, y) becomes negative. Indeed by considering the four terms:

(x + y -2 -4n) < (x -y -4n) < 0 < (x + y -2 + 4n) < (x -y + 4n),
it is straightforward to see that the smallest one in absolute value is |x + y -2 + 4n|. Since the change of convexity of the curve x → φ(x) appears for x = 1, the sum a n (t, x, y) + a -n (t, x, y) is negative as soon as x + y -2 + 4n ≥ √ t, which is satisfied in particular if n ≥ n 0 with

n 0 = √ t/4 + 2.
(1.10)

Since the terms of the series are negative for n ≥ n 0 , the series (1.5) is obviously not an alternating series.

Bound of the series reminder in (1.5)

Let us fix (t, x) ∈ R * + × [-1, 1].
We introduce the reminder of the series (1.5):

R n (y) := k≥n+1 (a k (t, x, y) + a -k (t, x, y)), n ≥ 0.

In order to apply the convergent series method, we need to bound this reminder. Since φ is a decreasing function, we obtain the following bound:

1 √ t k≥n+1 φ((α + 4k)/ √ t) ≤ 1 √ t ∞ n φ((α + 4z)/ √ t)dz,
as soon as α + 4n ≥ 0. Moreover, by symmetry,

1 √ t k≥n+1 φ((α -4k)/ √ t) ≤ 1 √ t ∞ n φ((α -4z)/ √ t)dz,
as soon as α -4n ≤ 0. If α = x -y or α = x + y -2, then the inequalities are satisfied for any n ≥ 1. Finally, for n ≥ 1, we get

|R n (y)| ≤ 1 √ t k≥n+1 φ((x -y -4k)/ √ t) + φ((x + y -2 + 4k)/ √ t) ≤ 1 √ t ∞ n φ((x -y -4z)/ √ t)dz + 1 √ t ∞ n φ((x + y -2 + 4z)/ √ t)dz.
Let us now observe that we can obtain a bound which is uniform with respect to both the variable x and y.

|R n (y)| ≤ 1 8 1 -erf 4n -2 √ 2t + 1 8 1 -erf 4n + 4 √ 2t ≤ r n (t) := 1 4 1 -erf 4n -2 √ 2t .
Let us note that this upper-bound is efficient for small values of t.

Proposal distribution for the rejection method using (1.5)

The aim is to use the acceptance/rejection algorithm in order to simulate a random variable with the target density q(t, x, y) = p(t, x, y)/P(τ x > t) where p(t, x, y) is given by (1.5). That's why we are looking for a proposal pdf h and a constant κ(t, x) independent of y satisfying

p(t, x, y) ≤ κ(t, x)h t,x (y), ∀y ∈ [-1, 1],
as explained in (1.7). We suggest here the following choice of the proposal distribution:

h t,x (y) = 1 √ t φ((x -y)/ √ t),
that means that we choose a gaussian distribution centered in x with variance t: it corresponds to the distribution of B t without any conditioning. The rejection method shall permit to go from this initial distribution to the conditional distribution with respect to the event {τ > t}. Let us note that this bound is uniform with respect to the variable x (starting value of the Brownian paths).

Some comments on the series expansion (1.6)

Since q(t, x, y) can be characterized by two different series expansions, it is useful to understand which kind of acceptance/rejection algorithm each series produces. Let us now focus on the second series (1.6):

p(t, x, y) = n≥1 β n (t)a n (x, y), t > 0, (x, y) ∈ [-1, 1] 2 ,
where β n (t) = e -n 2 π 2 t/8 and a n (x, y) = sin nπ 2 (x+1) sin nπ 2 (y +1) . Let us first notice that, for fixed t > 0, (β n (t)) n≥1 is a decreasing sequence of positive real numbers which converges towards 0. Since a n is defined as the product of two sine functions, one idea to bound the reminder of the series is to use the Abel transform and to prove that N n=1 a n (x, y) is bounded as N tends to infinity. It is quite easy to obtain some bound when x = y but for the particular case x = y, the sequence of partial sums is increasing (sum of squares) and tends to infinity. In conclusion, we cannot apply technics based on the Abel transform for the series (1.6). We shall therefore present an other approach.

Bound of the series reminder in (1.6) Let us define the reminder of the series (1.6):

Rn (y) := k≥n+1 β k (t)a k (x, y), t > 0, (x, y) ∈ [-1, 1] 2 .
Let us find a bound of this reminder. Since |a k (x, y)| ≤ 1 for any k ∈ N and any (x, y) ∈] -1, 1[ 2 , we obtain

| Rn | ≤ k≥n+1 β k (t), n ≥ 0.
Moreover, u → e -u 2 π 2 t/8 is a decreasing function on R + , and therefore we get the following bound

| Rn | ≤ ∞ n e -u 2 π 2 t/8 du = 2 πt 1 -erf nπ √ t 2 √ 2 := rn (t), n ≥ 0. (1.13)
Note that this bound is sharp for large value of t.

Proposal distribution for the rejection method using (1.6)

The aim is to use an acceptance /rejection method in order to simulate a random variable with the target density p(t, x, y)/P(τ > t) based on the equation (1.6), as already done for (1.5). We are looking for a probability distribution function h(y) and a constant κ(t, x) > 0 independent of y such that:

p(t, x, y) ≤ κ(t, x)h(y), ∀y ∈ [-1, 1],
as explained in (1.7). Our particular choice for the function h is

h(y) = π 4 sin π 2 (y + 1) , y ∈ [-1, 1],
which corresponds to the invariant probability measure of the Brownian motion conditioned to stay in the interval [-1, 1]. Let us now determine a constant κ(t, x). We have

p(t, x, y) = 4 π h(y) sin π 2 (x + 1) n≥1 β n (t)s n (x)s n (y),
where

s n (x) = sin( nπ 2 (x + 1)) sin( π 2 (x + 1)) . (1.14) Using the formula a n -b n = (a -b)(a n-1 + a n-2 b + . . . + b n-1 ) applied to the sine function, we observe that |s n (x)| ≤ n for all x ∈ [-1, 1]. Hence p(t, x, y) ≤ 4 π h(y) sin π 2 (x + 1) n≥1 n 2 β n (t), t > 0.
Moreover the function x → x 2 e -x 2 π 2 t/8 is decreasing as soon as x ≥ 2 √ 2 π √ t , we obtain the following upper bound:

n≥n0+1 n 2 β n (t) ≤ ∞ n0 z 2 e -z 2 π 2 t/8 dz, (1.15)
where

n 0 = 2 √ 2 π √ t + 1.
The integration by parts formula and a change of variables lead to

∞ n0 z 2 e -z 2 π 2 t/8 dz = 4n 0 π 2 t e -n 2 0 π 2 t/8 + 4 π 2 t ∞ n0 e -z 2 π 2 t/8 dz = 4 π 2 t e -n 2 0 π 2 t/8 n 0 + ∞ 0 e -u 2 π 2 t/8-n0uπ 2 t/4 du ≤ 4 π 2 t e -n 2 0 π 2 t/8 n 0 + 2 πt ≤ C 1 (t) := 8n 0 π 2 t e -n 2 0 π 2 t/8 , since 2 πt ≤ √ πn 0 /2 ≤ n 0 . We just note that C 1 (t) = 8 π 2 t 2 √ 2 π √ t + 1 e -n 2 0 π 2 t/8 ∼ 8 π 2 t e -π 2 t/8 , as t becomes large. Moreover n0 n=1 n 2 β n (t) ≤ C 2 (t) := n 3 0 β 1 (t) = 2 √ 2 π √ t + 1 3 e -π 2 t/8 . (1.16)
Finally we can choose the constant κ(t, x)

:= 4 π (C 1 (t) + C 2 (t)) sin( π 2 (x + 1)). Let us note that κ(t, x) ∼ 4 π sin( π 2 (x+1
))e -π 2 t/8 as t tends to infinity. Since the averaged number of iterations in the acceptance/rejection method corresponds to κ(t, x)/P x (τ > t) and since both κ(t, x) and P x (τ > t) are of the same order in the large time limit, the efficiency of the algorithm using (1.6) still remains strong when the time variable enlarges. Indeed let us decompose the following probability:

P x (τ > t) = n≥0 R n (t, x) where R n (t, x) := 4 π 1 2n + 1 exp - (2n + 1) 2 π 2 8 t sin (2n + 1)(x + 1) π 2 .
Using the definition (1.14), we obtain

P x (τ > t) R 0 (t, x) = 1 + n≥1 s 2n+1 (x) 2n + 1 exp - n(n + 1) 2 π 2 t .
The inequality |s n (x)| ≤ n for any x ∈ [-1, 1] and n ≥ 1 leads to

P x (τ > t) R 0 (t, x) ≥ 1 - n≥1 exp - n(n + 1) 2 π 2 t .
We observe that each term of the series converges in a monotonous way towards 0 as t tends to ∞ which implies that the ratio tends towards 1 by the Lebesgue theorem. We therefore deduce that for any θ > 1, there exists t θ > 0 such that the average number of iterations in the acceptance/rejection algorithm is smaller than θ as soon as t ≥ t θ .

Algorithm and numerics

Let us now describe the algorithms used in the following in order to simulate q(t, x, dy), the conditional distribution of the Brownian motion at time t given τ > t. We know that p(t, x, dy) admits two different series expansion presented in (1.5) and (1.6). As introduced in Section 1.2, the classical convergent series method requires both a proposal distribution h, satisfying the inequality (1.7) that is p(t, x, y) ≤ κ(t, x)h t,x (y) (also necessary for classical acceptance/rejection methods), and a precise description of the series reminder characterized by the sequence (r n (t)) n and (r n (t)) n , see (1.8).

We just recall the results obtained in the previous section:

Table 1: Series expansion (1.5)
Proposal distribution: h t,x (y)

1 √ t φ((x -y)/ √ t) Constant: κ(t, x) 3 + √ t/4
Reminder bounds: r n (t)

1 4 1 -erf 4n-2 √ 2t Table 2: Series expansion (1.6) h t,x (y) π 4 sin π 2 (y + 1) κ(t, x) 4 π sin π 2 (x + 1) 8n0 π 2 t e -n 2 0 π 2 t/8 + n 3 0 e -π 2 t/8 where n 0 = 2 √ 2 π √ t + 1 rn (t) 2 πt 1 -erf( nπ √ t 2 √ 2 )
It is straightforward that the convergent series method using (1.5) is convenient for small values of t while (1.6) is rather convenient for large t. That's why we choose a threshold t c > 0 (threshold for the conditional distribution) such that (1.5) is used for t ≤ t c and (1.6) otherwise. For practical purposes, we fix t c = 0.7, this choice is motivated by the curves in Fig. 2 and will be held for all numerical illustrations presented in this study. Applying the convergence series algorithm either for small times or large times permits to obtain the simulations presented in Fig. 3. We observe that, even if x = 0, the condition t < τ leads to a distribution which looks like symmetric as t becomes large and which converges towards h t,x (y) the invariant measure of the diffusion conditioned to stay in the interval [-1, 1].

Efficiency of the algorithms

Let us just recall classical results concerning the efficiency of the convergent series method introduced in Section 1.2. We introduce N c the random number of computations of terms f n used in order to simulate just one random variable X. In fact N c also corresponds to the number of random variable S generated before the algorithm halts. We also define N l the local counter depending on Y which represents the number of terms f n used till the decision of acceptance or rejection of the variable Y can be taken. Theorem IV.5.2 in Devroye [START_REF] Devroye | Nonuniform random variate generation[END_REF] emphasizes the following upper-bound:

E[N l |Y ] ≤ 2 κ h(Y ) ∞ n=0 R n (Y ) (1.17)
where R n (y) is the reminder of the series expansion defined in the general framework (1.8). Let us also notice that the starting idea behind this bound is quite classical: it suffices to use the classical expansion:

E[N l |Y ] ≤ ∞ n=0 P(N l > n|Y ).
In the study developed in the previous section, we obtained precise bounds for the reminder terms only for n ≥ 1. So we shall modify the bound isolating the term n = 0 -instead of (1.17) -which plays a crucial role for the description of the algorithm efficiency:

E[N l |Y ] ≤ 1 + 2 κ h(Y ) ∞ n=1 R n (Y ).
Since the uniform bound |R n (y)| ≤ r n holds for n ≥ 1, and since the average number of iterations is κ/I(f ), Wald's inequality immediately implies that the total number N c satisfies for general positive integrable functions f :

E[N c ] ≤ κ I(f ) 1 -1 E[N l |Y = y]h(y) dy ≤ 1 I(f ) κ + 4 ∞ n=1 r n . (1.18)
Let us describe the consequences of this general statement to our particular algorithms.

For small values of the time variable t.

For small t, it is useful to use the series expansion (1.5) and the associated reminder bounds presented in Table 1.

Proposition 1.1. For t > 0 and x ∈ [-1, 1], we observe the following bound for the number of computations in the convergent series algorithm associated to the expansion (1.5):

E[N c ] P x (τ > t) ≤ U 1 (t) := 3 + √ t/4 + √ t 2 √ 2π ϑ 2 (0, e -8/t ), (1.19)
where ϑ 2 stands for the Jacobi theta function given by

ϑ 2 (z, q) := 2q 1/4 ∞ n=0 q n(n+1) cos((2n + 1)z),
see [START_REF] Whittaker | A course of modern analysis[END_REF] for instance. In particular the r.h.s of the previous inequality tends to 3 as t tends to 0.

Proof. The arguments are based on the upper-bound (1.18). Let us recall that

r n (t, x) = 1 4 1 -erf 4n -2 √ 2t , for n ≥ 1.
Using the bound

1 -erf(x) ≤ 1 x √ π e -x 2 , (1.20)
we obtain

I(f )E[N c ] ≤ κ + √ 2t √ π ∞ n=1 e -(4n-2) 2 2t 4n -2 ≤ κ + √ t √ 2π e -2/t ∞ n=1 e -(4n-2) 2 -4 2t ≤ κ + √ t √ 2π e -2/t ∞ n=0 e -8 t n(n+1) ≤ κ + √ t 2 √ 2π ϑ 2 (0, e -8/t ).
In order to conclude, it suffices to replace κ by (1.12) and to notice that I(f

) = P x (τ > t).
For large values of the time variable t.

For large values of t it is more convenient to use the series expansion (1.6). A similar approach to Proposition 1.1 easily leads to the upper bound:

Proposition 1.2. Let t ≥ 0 and x ∈ [-1, 1].
The number of computations N c for the convergent series algorithm associated to the expansion (1.6) satisfies

E[N c ] P x (τ > t) ≤ U 2 (t) := κ(t, x) + 8 π 2 t ϑ 3 (0, e -π 2 t/8 ) -1 , (1.21)
where the constant κ(t, x) associated to (1.6) is described in Section 1.3 and ϑ 3 stands for the Jacobi theta function given by

ϑ 3 (z, q) := 1 + 2 ∞ n=1 q n 2 cos(2nz).
Since the proofs of Proposition 1.2 and Proposition 1.1 are similar, we let the details of the proof to the reader.

Distribution of the Brownian first exit time

In this section, we focus our attention to the distribution of the first Brownian exit (τ, B x τ ), where B x stands for the one-dimensional Brownian motion starting in x and τ the exit time of the interval [-1, 1] as defined in (1.2).

Distribution of the first exit time τ

Since p(t, x, y) satisfies two different series expansions (1.5) and (1.6), we can also obtain two series for the cumulative distribution of the exit time τ :

P x (τ ≤ t) = 1 -P x (τ > t) = 1 - 1 -1 p(t, x, y)dy.
Therefore (1.5) becomes

P x (τ ≤ t) = 1 - +∞ n=-∞ -Φ((x -1 -4n)/ √ t) + Φ((x + 1 -4n)/ √ t) -Φ((x + 1 -2 -4n)/ √ t) + Φ((x -1 -2 -4n)/ √ t) = 1 -2 +∞ n=-∞ -Φ((x -1 -4n)/ √ t) + Φ((x + 1 -4n)/ √ t) = 1 + 2 +∞ n=-∞ (-1) n Φ((x -(2n + 1))/ √ t).
We deduce easily the expression of the pdf :

p τ (t) = +∞ n=-∞ (-1) n+1 (x -(2n + 1)) t 3/2 φ((x -(2n + 1))/ √ t). (2.1)
Of course the particular case x = 0 ensures simplifications: the symmetry of the function φ implies φ(-

(2n + 1)/ √ t) = φ(-(2 × (-n -1) + 1)/ √ t) and therefore p τ (t) = +∞ n=0 (-1) n R 1 (2n + 1, t) with R 1 (n, t) := 2n t 3/2 φ n √ t . (2.2)
Such an expression for the pdf of τ is of prime interest for simulation purposes. Let us note that the reminder of the series is small for small values of t. Indeed the sequence (R 1 (2n + 1, t)) n≥1 is a decreasing sequence under the assumption t ≤ 9. For large values t, it is more convenient to consider (1.6) which leads to

P x (τ > t) = 4 π n≥0 1 2n + 1 exp - (2n + 1) 2 π 2 8 t sin (2n + 1)(x + 1)π/2 ,
(see also (3.3) in Milstein and Tretyakov [START_REF] Milstein | Simulation of a space-time bounded diffusion[END_REF]). We deduce the expression

p τ (t) = π 2 +∞ n=0 (2n + 1) exp - (2n + 1) 2 π 2 8 t sin (2n + 1)(x + 1)π/2 . (2.3)
Here also the case x = 0 plays a crucial role due to its simple expression:

p τ (t) = +∞ n=0 (-1) n R 2 (2n + 1, t) with R 2 (n, t) := πn 2 exp - n 2 π 2 8 t . (2.4)
This expression is of prime interest as soon as the variable t is large. Indeed the series becomes alternating with a decreasing sequence (R 2 (2n + 1)) n≥1 as soon as t ≥ 4/(9π 2 ).

Algorithms and numerics

The aim is to describe an algorithm which permits to simulate both the Brownian exit time τ and the exit position B x τ of the interval [-1, 1]. In the previous section, we have pointed out two different series (2.1) and (2.3) corresponding to the pdf of τ . We could therefore apply the classical convergent series method for simulation purposes. But here, we prefer to introduce another method based on an iterative procedure, the advantage of our approach is to deal with alternative series rather than general convergent series and therefore the reminder of the series is easier to bound.

First of all, we shall focus our attention to the case x = 0. In this case the interval is symmetric and the expression of the pdf of τ is simplified, see (2.2) and (2.4). Moreover these series are alternating for suitable conditions on the time variable t. So we can apply the alternating series method for the simulation of the exit time τ , the exit location being just uniformly distributed in {-1, 1}. The alternating series method is an acceptance/rejection method. We need therefore a proposal distribution and an acceptance procedure. For the proposal distribution, let us first fix t e > 0 (threshold for the exit distribution). Let us consider a standard gaussian random variable G. We define a new variable Y as follows: if 1/G 2 ≤ t e then Y = 1/G 2 otherwise Y = t e -8 π 2 log(U ) where U is uniformly distributed on [0, 1] and independent of G. The density function of Y which corresponds to the proposal distribution of the algorithm satisfies:

ĥ(t) = 1 t 3/2 φ 1 √ t 1 {t≤te} + π 4κ e -π 2 8 t 1 {t>te} , (2.5)
where κ -1 = π erf( 1/(2t e )) e π 2 te/8 /2. Using (2.2) and (2.4), we deduce the following link between the proposal distribution and the target distribution:

p τ (t) =    2 ĥ(t) 1 -R1(3,t) R1(1,t) + R1(5,t) R1(1,t) -. . . for t ≤ t e , 2κ ĥ(t) 1 -R2(3,t) R2(1,t) + R2(5,t) R2(1,t) -. . . for t > t e . (2.6)
Let us present now the acceptance/rejection method applied to this particular situation which permits to simulate the exit time τ in the symmetric case.

BROWNIAN EXIT TIME FOR SYMM. INTERVALS (with parameter t e )

BROWNIAN_EXIT_SYMMETRIC First initialization: N s = 0.

Step 0: Second initialization. n = 0, Test = 0, L 0 = 0, U 0 = 1.

Step 1. Generate a random variable Y with pdf ĥ given by (2.5). If Y ≤ t e then set i = 1 and C = 1 else set i = 2 and C = κ.

Step 2. Generate a random variable V uniformly distributed on [0, 1].

Step 3. While (V < CU n ) & (Test = 0) do:

• n ← n + 1 and N s ← N s + 1, • L n = U n-1 -Ri(4n-1,Y ) Ri(1,Y ) , U n = L n + Ri(4n+1,Y ) Ri(1,Y ) • Test = 1 {V ≤CLn}
Step 4. If Test = 1 then Z = Y otherwise go to Step 0.

Outcome: the random variable Z with density p τ and the number of incrementations needed N s .

The outcome variable distribution p τ corresponds to the target one as soon as the sequences (R i (2n + 1, t)/R i (1, t)) n≥1 appearing in (2.6) are decreasing. This is an easy adaptation of the classical alternating series method (proof left to the reader). Such a property leads to the condition:

4 9π 2 ≤ t e ≤ 1. (2.7)
Let us just note that the probability of acceptance A = {Accept Y } in the acceptance/rejection algorithm satisfies:

P(A) = P(A, Y ≤ t e ) + P(A, Y > t e ) = P V ≤ C p τ (Y ) 2C ĥ(Y ) , Y ≤ t e + P V ≤ C p τ (Y ) 2C ĥ(Y ) , Y > t e = 1 2 . (2.8)
We deduce that the number of random variables Y simulated in order to obtain Z is geometrically distributed with average 2 and does not depend on the choice of t e . Nevertheless the parameter t e has an influence on the efficiency of the algorithm as illustrated in Fig. 5. This figure represents the averaged value of the number of iterations of Step 3 needed by the algorithm in order to simulate one r.v. with the p.d.f. p τ . This figure suggest to choose a parameter t e of the order of 1/2. For this parameter, we obtain: 1 -e -5π 2 te , (2.9)

κ -1 ≈ 2.
where x → erfc(x) is the complementary error function and t e is the parameter appearing in the algorithm and satisfying (2.7). In particular, for t e = 1/2, the value of the upper bound is approximatively 1.027 which emphasizes that the inequality is quite sharp in comparison with the estimated number of iterations in Figure 5 (solid line/stars). Moreover the bound is reasonably small which permits to confirm the efficiency of the proposed algorithm.

Proof. The arguments are quite similar to those developed in Theorem 5.1 in Devroye [START_REF] Devroye | Nonuniform random variate generation[END_REF] for the efficiency of the alternating series method. Let us denote by N loop s the number of steps of type 3 used in order to go from Step 1 to Step 4, in other words: the number of increments n ← n + 1 during one loop. By Wald's equation and since the number of Steps 0 in this algorithm is geometrically distributed with parameter 1/2 (see (2.8)), we get

E[N s ] = 2E[N loop s ].
Let us first compute the probability of the following event {N loop s > 0}. This event corresponds to {V < C}. In the following, we shall just recall that the constant C so as L 1 , U 1 ... depend on Y , that's why we use from now on the notation C = C Y . We obtain:

P(N loop s > 0|Y ) = P(V < C Y |Y ) = C Y , since 0 ≤ C Y ≤ 1 a.s.
Let us now consider the event {N loop s > 1}. This event corresponds to {CL 1 < V < CU 1 } and therefore

P(N loop s > 1|Y ) = C Y (U 1 -max(L 1 , 0)) ≤ C Y (U 1 -L 1 ).
The definition of L 1 and U 1 leads to

P(N loop s > 1|Y ) ≤ C Y • R i (5, Y ) R i (1, Y ) ,
where i = 1 for Y ≤ t e and i = 2 otherwise. The same arguments lead to

P(N loop s > k|Y ) = P(C Y L k < V < C Y U k |Y ) = C Y (U k -L k ) = C Y • R i (4k + 1, Y ) R i (1, Y ) .
We deduce

E[N loop s |Y ] = k≥0 P(N loop s > k|Y ) = C Y k≥0 R i (4k + 1, Y ) R i (1, Y ) .
Using the distribution of the random variable Y with density ĥ in (2.5), we obtain By the change of variable z = 1/w, we get

E[N loop s ] = k≥0 te 0 R 1 (4k + 1, y) R 1 (1, y) ĥ(y)dy + κ k≥0 ∞ te R 2 (4k + 1, y) R 2 (1, y) ĥ(y)dy = 1 2 k≥0 te 0 R 1 (4k + 1, y) dy + 1 2 k≥0 ∞ te R 2 (4k + 1, y) dy =: 1 2 A 1 (t e ) + 1 2 A 2 (t e ). ( 2
t e 8π te 0 1 √ w e -1 2w w 3/2 dw = t e 8π ∞ 1/te e -z/2 dz = t e 2π e -1 2te . (2.11)
Moreover, using the change of variable w = r -2 /2, we have Let us now come to the exit simulation in the asymmetric case. Since we know how to handle with symmetric intervals, we shall use this first algorithm BROWNIAN_EXIT_SYMMETRIC in an iteration procedure in order to solve the asymmetric case. The main idea consists in the following :

3 √ 8π te 0 e -1 2w w 3/2 dw = 3 √ π ∞ 1/ √ 2te e -r 2 dr = 3 2 erfc 1 √ 2t e . ( 2
• starting in x, we consider the largest interval centered in x of the type [x -δ, x + δ] and included in [-1, 1].

• We simulate the exit time of this interval denoted by T 1 and the exit position will be uniformly distributed in {x -δ, x + δ}. Step 0: Initialization. X = x, T = 0, test = 0, L = a, U = b and N as = 0.

• If B x T1 ∈ {-
While (test=0) do:

Step 1. Set D = min(X-L, U -X). Generate a random variable τ (and the number of iterations N s ) using the algorithm BROWNIAN_EXIT_SYMMETRIC and define S = D 2 τ . Set T ← T + S and N as ← N as + N s .

Step 2. Generate a random variable V uniformly distributed on {-D, D}. Set X ← X + V Figure 6: Histogram of the algorithm counter corresponding to the Brownian exit from the interval [-1.5, 2] (left), p.d.f. and histograms of the exit time when exiting at the top or at the bottom of this interval (right). A sample of 100 000 simulations has been used for these figures and te = 0.5.

Step 3. If X ∈ {a, b} then test = 1 else set L ← X -L and U ← U -X.

End While

Outcome: the exit time T and the exit location X of the interval [a, b] and the total number of iterations N as .

The number of iterations is stochastically upper-bounded by a geometrically distribution with parameter 1/2 since we only deal with symmetric intervals. Moreover this random number is independent of the generation cost of any Brownian symmetric exit time and position. Let us finally note that (T, X) and (τ, B x τ ) are obviously identically distributed. This algorithm can be illustrated by Fig. 6 Corollary 2.2. The algorithm BROWNIAN_EXIT_ASYMM involves a random number of calls to BROWNIAN_EXIT_SYMMETRIC whose efficiency is characterized by the number of iterations N s . That's why the total number of iterations N as is directly linked to the efficiency of BROWNIAN_EXIT_ASYMM. We have

E[N as ] ≤ 2E[N s ] ≤ 2t e π e -1 2te + 3 erfc 1 √ 2t e + 8 π e -π 2 te 8 + 8 5π
e -25 π 2 te 8

1 -e -5π 2 te , where x → erfc(x) is the complementary error function and t e is the parameter appearing in BROWNIAN_EXIT_SYMMETRIC.

The statement is a direct consequence of the geometrical distributed upperbound of the number of calls to the symmetric case algorithm on one hand (Wald's identity therefore leads to E[N as ] ≤ 2E[N s ]) and of Proposition 2.1 on the other hand.

First exit time for one-dimensional diffusion

This section is concerned with the exit problem for a one-dimensional diffusion. Let us first consider (X t , t ≥ 0) the solution of the following stochastic differential equation:

dX t = µ(X t )dt + σ(X t )dB t , X 0 = x ∈ [a, b], (3.1)
where (B t , t ≥ 0) stands for the standard one-dimensional Brownian motion.

As already introduced in (0.1), we denote by τ a,b (X) the first exit time of the interval [a, b]. It suffices to assume the existence of a unique weak solution to the equation, see for instance [START_REF] Karatzas | Brownian motion and stochastic calculus[END_REF] for the corresponding conditions. In order to simplify the presentation of all algorithms, we restrict our study to the constant diffusion case: σ(x) ≡ 1. Using the classical Lamperti transform, we observe that this restriction is not sharp at all. That is why, from now on, X stands for the unique solution of

dX t = µ(X t )dt + dB t , t ≥ 0, X 0 = x ∈ [a, b]. (3.2)
Let us note that, for the particular Brownian case: µ(x) ≡ 0, the first exit time has already been presented in Section 2. The aim is to use the results developed in the previous sections when considering the general diffusion case and the important tool for such a strategy is Girsanov's formula.

First Exit Time Algorithm (DET)

Let us now consider the exact simulation algorithm which permits to handle with the diffusion exit problem.

DIFFUSION EXIT TIME (DET) Parameter: γ 0 , input functions γ(•) and β(•)

First initialization: N tot = 0.

Step 0: Initialization. Z = x, T = 0, test = 0. Here x stands for the initial value of the diffusion.

While (test=0) do:

Step 1. Generate an expon. distr. random variable E with parameter γ 0 and U and V two random variables uniformly distributed on [0, 1], the variables E, U and V being independent.

Step 2. Simulate the Brownian exit time and location

(S, Y, N as ) = BROWNIAN_EXIT_ASYMM (Z,[a,b])
and set N tot ← N tot + N as .

Step 3. If S < E then

• if U ≤ β(Y ) then set test = 1, Z ← Y
E[f (T )g(Z)] = E x f (τ a,b (X))e -ρ τ a,b (X) g(X τ a,b (X) ) E x [e -ρ τ a,b (X) ] , (3.4)
where X stands for the diffusion defined by (3.2). In particular, if µ + µ 2 is a non-negative function on the interval [a, b] (we set ρ = 0), then the outcome (Z, T ) has the same distribution as (X τ a,b (X) , τ a,b (X)).

Proof. The proof of Theorem 3.1 is organized as follows: first we introduce a stochastic theoretical model which is based on the one-dimensional Brownian motion and on an independent Poisson process. This model is closely related to the exact simulation introduced by Beskos and Roberts in [START_REF] Beskos | Exact simulation of diffusions[END_REF]. Secondly we consider in details the outcome of Algorithm (DET) and point out the strong link between the outcome and random variables derived from the theoretical model of the first step. Finally we shall use the Girsanov transformation in order to conclude the proof.

Step 1. A theoretical model. Let us consider for any non-negative functions f and g: where λ(•) is the Lebesgue measure. In particular, if A is defined as follows:

I(x, f, g) := E x f (τ B )g(B τ B )e -
A := (t, y) ∈ [0, τ B ] × [0, γ + ] : y ≤ γ(B t )
the independence of both the Brownian motion and the Poisson process implies

I(x, f, g) = E x f (τ B )g(B τ B )1 {N (A)=0} = n≥0 I n (x, f, g)
where

I n (x, f, g) :=      E x f (τ B )g(B τ B )1 {N (A)=0} 1 {ξn<τ B ≤ξn+1} , for n ≥ 1, E x f (τ B )g(B τ B )1 {N (A)=0} 1 {τ B ≤ξ1} , n = 0. Let us first observe that {τ B ≤ ξ 1 } ⊂ {N (A) = 0} which implies that I 0 (x, f, g) = E x f (τ B )g(B τ B )1 {τ B ≤ξ1} . (3.6)
For n ≥ 1, due to the Markov property of the Brownian path and since the Poisson process is independent of the Brownian motion, we have

I n (x, f, g) = E x f (τ B )g(B τ B )1 {U1>γ(B ξ 1 ),...,Un>γ(B ξn )} 1 {ξn<τ B ≤ξn+1} = E x I n-1 (B ξ1 , f (ξ 1 + •), g)1 {U1>γ(B ξ 1 ), ξ1<τ B } . (3.7)
Since ξ 1 is exponentially distributed and independent of the Brownian motion, we get

I n (x, f, g) = ∞ 0 E x [I n-1 (B t , f (t + •), g)1 {U1>γ(Bt),t<τ B } ]γ + e -γ+t dt = ∞ 0 E x [I n-1 (B t , f (t + •), g)1 {U1>γ(Bt)} |t < τ B ]p(t, x)γ + e -γ+t dt = R+×[a,b] I n-1 (y, f (t + •), g)Γ(y)p(t, x)M(x, dt, dy), (3.8)
where Γ(y) := P(U 1 > γ(y)) and p(t, x) = P x (τ B > t). The probability measure M(x, dt, dy) represents the distribution of the couple (T, Y ) where T is an exponentially distributed r.v. with parameter γ + and Y represents the distribution of B T the value of the Brownian motion at time T starting in x and given {τ B > T }, T and (B t ) being independent. The recurrence relations (3.7) and (3.8) are satisfied for any n ≥ 1.

Step 2. Relation between Algorithm (DET) and the Brownian-Poisson theoretical model presented in Step 1.

If we denote by N 0 the number of Step 0 used during the procedure which leads to the computation of the outcome (Z, T ) and N 1 the number of exponentially distributed random variables generated (Step 1 ), then we obviously obtain

E x [f (T )g(Z)1 {N0=1, N1=1} ] = E[f (S)g(Y )1 {U ≤β(Y ), S<E} ] = E[f (S)g(Y )β(Y )1 {S<E} ] = E x f (τ B )g(B τ B )β(B τ B )1 {τ B ≤ξ1} = I 0 (x, f, g × β) (3.9) using (3.7).
Let us now consider that N 1 > 1. On the event {N 0 = 1} ∩ {N 1 > 1}, we observe that

• the first exponentially distr. r.v. E satisfies E < S where S is given by

Step 2 and corresponds to the Brownian exit time of the interval [a, b].

• Y c given that E < S satisfies V > γ(Y c ), where V is uniformly distributed and Y c becomes the new starting point for the next use of the exit problem.

In other words, if we denote by

I n := E x [f (T )g(Z)1 {N0=1, N1=n+1} ] then (3.6) leads to I 0 (x, f, g) = I 0 (x, f, g • β)
. Moreover by (3.7)

I 1 (x, f, g) = E x [f (T )g(Z)1 {N0=1, N1=2} ] = E x [I 0 (Y c , f (E + •), g • β)1 {V >γ(Yc),E<S} ] = E x [I 0 (B ξ1 , f (ξ 1 + •), g × β)1 {U1>γ(B ξ 1 ),ξ1<τ B } ] = I 1 (x, f, g × β),
and using the same arguments, we prove easily that I n satisfies the recurrence relations (3.7) and (3.8). We deduce that

I n (x, f, g) = I n (x, f, g × β), ∀n ≥ 0. Therefore I(x, f, g × β) = n≥0 I n (x, f, g × β) = n≥0 I n (x, f, g) = E x [f (T )g(Z)1 {N0=1} ].
(3.10) Since the Algorithm (DET) is an acceptance/rejection algorithm, (3.10) leads to

E x [f (T )g(Z)] = E x [f (T )g(Z)1 {N0=1} ] P(N 0 = 1) = I(x, f, g × β) I(x, 1, β)
.

By (3.5) E x [f (T )g(Z)] = E x f (τ B )g(B τ B )β(B τ B )e -τ B 0 γ(Bs) ds E x β(B τ B )e -τ B 0 γ(Bs) ds . (3.11)
Step 3. The Girsanov transformation. In this last part of the proof, the aim is to link the distribution of (T, Z) described in (3.11) to the distribution of

(X τ a,b (X) , τ a,b (X))
where X is the diffusion defined by (3.2). Using the definition of the function β and Itô's formula, we obtain that

M t := β(B t ) β(B 0 ) e ρ t-t 0 γ(Bs) ds = exp Bt B0 µ(y) dy - 1 2 t 0 µ (B s ) + µ 2 (B s ) ds = exp t 0 µ(B s ) dB s - 1 2 t 0 µ 2 (B s ) ds , (3.12)
is an exponential martingale. Moreover the stopped martingale (M t∧τ B ) t≥0 is bounded. The stopping theorem therefore implies E x [M τ B ] = 1 which means that

E x β(B τ B )e ρ τ B -τ B 0 γ(Bs) ds = β(x).
The expression (3.11) and Girsanov's transformation permits to obtain

E x [f (T )g(Z)] = E x f (τ B )g(B τ B )M τ B e -ρτ B E x M τ B e -ρτ B = E x f (τ a,b (X))e -ρ τ a,b (X) g(X τ a,b (X) ) E x [e -ρ τ a,b (X) ]
,

where X stands for the diffusion defined by (3.2). In particular, if ρ = 0 that is µ 2 + µ is a non-negative function, we get

E x [f (T )g(Z)] = E x f (τ a,b (X))g(X τ a,b (X) ) .

Efficiency of the algorithm

Let us now focus our attention to the analysis of Algorithm (DET). Of course since the algorithm permits to simulate exactly the random variables desired, there is no error terms to deal with, it suffices therefore to describe the time needed by the algorithm. We introduced in Algorithm (DET) the random variable N tot which permits to have a precise idea of the efficiency. Let us just add the information concerning the starting position of the Brownian motion B 0 = x with the following notation N tot = N x tot . Let us also note that the efficiency of Algorithm (DET) in particular depends on two different parameters: t e which appears in the use of the algorithm BROWNIAN_EXIT_SYM and therefore also in BROWNIAN_EXIT_ASYMM and t c ∈ R + which appears in CONDITIONAL_DISTR (if t ≤ t c we consider the algorithm associated to the small values of t and for t > t c the algorithm associated to the large values).

Using informations concerning the cost of each part of the algorithm, namely Proposition 1.1, Proposition 1.2 and Proposition 2.1, we obtain a bound for

E[N x tot ].
Theorem 3.2. The random variable N x tot which is one of the outcomes of Algorithm (DET) and represents its cost satisfies the following bound: there exists a constant C(t c , t e ) > 0 independent of both the interval [a, b] and the starting position x, such that Each term appearing in the denominator, that is E x [e -ρτ a,b (X) ] on one hand and E (a+b)/2 [e -γ+τ a,b (B) ] on the other hand, tends to 0 when the interval size b -a tends to infinity. It is therefore important to choose ρ and γ + as small as possible in order to obtain a sharper bound. The exit problem of the Brownian motion can be precisely described using classical results on Laplace transforms and differential equations, see for instance [START_REF] Darling | The first passage problem for a continuous Markov process[END_REF] or [START_REF] Borodin | Handbook of Brownian motion-facts and formulae[END_REF]. We obtain

E[N x tot ] ≤ C(t c , t e ) E x [e -
E (a+b)/2 [e -γ+τ a,b (B) ] = cosh γ + 2 (b -a) -1
.

Let us also note that E x [e -ρτ a,b (X) ] can be linked to the two linearly independent solutions of the differential equation (see, for instance [START_REF] Darling | The first passage problem for a continuous Markov process[END_REF])

1 2 d 2 u dx 2 + µ(x) du dx -λu = 0,
which permits to describe the asymptotic behaviour as b-a tends to infinity. For the link between the Laplace transform and the speed measure, see for instance [START_REF] Bass | Diffusions and elliptic operators. Probability and its Applications[END_REF].

Proof. The counter N tot introduced in the algorithm can be decomposed as follows:

N x tot = k≥1 N x,k tot ,
where N x,k tot represents the number of counter increases observed in-between the k-th and (k+1)-th passage through the item Step 0. We recall that we defined N 0 in the proof of Theorem 3.1: it corresponds to the number of Step 0 necessary to obtain the desired outcome. Since Algorithm (DET) is an acceptance-rejection algorithm, the random variable N 0 is geometrically distributed. Let us also note that N x,k tot = 0 a.s. on the event {N 0 < k} and conditionally to {N 0 ≥ k}, N x,k tot has the same distribution as N x,1 tot . Hence

E[N x tot ] = k≥1 E[N x,k tot ] = k≥1 E[N x,k tot 1 {N0≥k} ] = E[N x,1 tot ] k≥1 P(N 0 ≥ k) = E[N x,1 tot ] P(N 0 = 1) = E[N x,1 tot ] E x [e -ρτ a,b (X) ]
.

(3.13)

The last equality is related to the third step in the proof of Theorem 3. • if S < E then N x,1 tot = N as .

• if S > E then N x,1 tot = N as + N c + N Yc,1 tot 1 {γ0V >γ(Yc)} where N x,1 tot is an independent copy of N x,1 tot . Such a property is essentially based on the Markov property of the Brownian paths.

We deduce ) is a concave function whose derivative vanishes for y = 0 (see the expression (2.3)), we get

E[N x,1 tot ] ≤ E[N as ] + E[N c 1 {S>E} ] + E[ N Yc,1 tot 1 {S>E} 1 {γ0V >γ(Yc)} ] ≤ E[N as ] + E[N c 1 {S>E} ] + E[ N Yc,1 tot 1 {S>E} ] ≤ E[N as ] + E[E[N c |E]1 {S>E} ] + P(S > E) sup x∈[a,b] E[N x,1 tot ]. ( 3
P x (τ a,b (B) > ξ) ≤ P 0 (τ > 4ξ(b -a) -2 ).
• Since N as is the cost of the function BROWNIAN_EXIT_ASYMM (x, [a, b]) and using the scaling property (1.1), we obtain that it is equal to the cost of 1 -e -5π 2 te , where x → erfc(x) is the complementary error function and t e is the parameter appearing in BROWNIAN_EXIT_SYMMETRIC. This parameter t e satisfying (2.7) can be chosen in order to minimize this average.

• Moreover we need some information on E x

E x [N c |E = t] ≤ Θ(4t(b -a) -2 ) P y (τ > 4t(b -a) -2 ) with Θ(t) = U 1 (t)1 t≤tc + U 2 (t)1 t>tc . Let us introduce C 1 (t c ) := sup t≥0 Θ(t) < ∞, then E[E[N c |E]1 {S>E} ] = ∞ 0 γ + E[E[N c |E = t]1 {S>t} ]e -γ+t dt ≤ ∞ 0 γ + P(S > t) Θ(4t(b -a) -2 ) P y (τ > 4t(b -a) -2 ) e -γ+t dt ≤ C 1 (t c ) ∞ 0 γ + e -γ+t dt = C 1 (t c ).
Using these three items and (3.14), we obtain

E[N x,1 tot ] ≤ E[N as ] + E[E[N c |E]1 {S>E} ] P 0 (τ ≤ 4ξ(b -a) -2 ) ≤ C 0 (t e ) + C 1 (t c ) E 0 [e -γ+(b-a) 2 τ /4 ]
.

The identity (3.13) permits to conclude the proof.

Modifications of Algorithm (DET)

Under the drift condition µ + µ 2 ≥ 0, Algorithm (DET) permits to simulate in an exact way the exit time of the interval [a, b] for diffusion processes. Let us now improve this algorithm in order to deal with any one dimensional diffusion. This generalization first requires a modified algorithm with outcome Step 0: Initialization. K = κ, Z = x, T = 0, test = 0. Here x stands for the initial value of the diffusion and κ the time upper-bound.

(X τ a,b (X)∧κ , τ a,b (

While (test=0) do:

Step 1 Generate an expon. distr. random variable E with parameter γ 0 and U , V and W three random variables uniformly distributed on [0, 1], the variables E, U , V and W being independent.

Step 2. Simulate the Brownian exit time and location

(S, Y, N as ) = BROWNIAN_EXIT_ASYMM (Z,[a,b])
and set N tot ← N tot + N as .

Step 3. If S = min(K, E, S) then

• if U ≤ β m (Y ) and W ≤ e -ρ (K-S) then set test = 1, Z ← Y and T ← T + K else go to Step 0 end if, elseif K = min(K, E, S)
• simulate the conditional Brownian location at time K:

(Y c , N c ) = CONDITIONAL_DISTR (Z,[a,b],K) and set N tot ← N tot + N c . • if U ≤ β m (Y c ) then set test = 1, Z ← Y and T ← T + S else go to Step 0 end if, elseif E = min(K, E, S)
Let us introduce the following series expansion:

I(x, f, g, κ) = n≥1 I n (x, f, g, κ),
the definition of I n (x, f, g, κ) is similar to the definition appearing in the proof of Theorem 3.1, it suffices to replace in the definition τ B by τ B ∧ κ. Therefore

I 0 (x, f, g, κ) = E x [f (τ B ∧ κ)g(B τ B ∧κ )1 {τ B ∧κ≤ξ1} ],
where (ξ 1 , U 1 ) stands for the coordinates of the Poisson process point with the smallest abscissa. Moreover we obtain the following step by step property:

I n (x, f, g, κ) = E x I n-1 (B ξ1 , f (ξ 1 + •), g, κ -ξ 1 )1 {U1>γ(B ξ 1 ),ξ1<τ B ∧κ} . (3.16)
Step 2. Relation between Algorithm (κ-DET) and the Brownian-Poisson theoretical model presented in Step 1.

If we denote by N 0 the number of Step 0 used during the procedure which leads to the computation of the outcome (Z, T ) and N 1 the number of exponentially distributed random variables generated (Step 1 ), then we obviously obtain for

I n (x, f, g, κ) := E x [f (T )g(Z)1 {N0=1, N1=n+1} ], I 0 (x, f, g, κ) = E f (S)g(Y )1 {S<E∧K, U ≤βm(Y ), W ≤e -ρ(κ-S) } + E f (κ)g(Y c )1 {κ≤E∧S, U ≤βm(Yc)} = E f (τ B ∧ κ)e -ρ(κ-τ B ∧κ) g(B τ B ∧κ )β m (B τ B ∧κ )1 {ξ1≥τ B ∧κ} = I 0 (x, f × e -ρ(κ-•) , g × β m , κ).
Using similar arguments, we can prove that I n satisfies the same step by step relation as (3.16). So by identification, we get I n (x, f, g, κ) = I n (x, f ×e -ρ(κ-•) , g× β m , κ) for all n ≥ 0 and therefore

E[f (T )g(Z)1 {N0=1} ] = I(x, f, g, κ) = n≥0 I n (x, f, g, κ) = I(x, f × e -ρ(κ-•) , g × β m , κ).
Since Algorithm (κ-DET) is an acceptance rejection algorithm, we have

E[f (T )g(Z)] = E[f (T )g(Z)1 {N0=1} ] P(N 0 = 1) = I(x, f × e -ρ(κ-•) , g × β m , κ) I(x, e -ρ(κ-•) , β m , κ) .
The link between β m and β leads to

E[f (T )g(Z)] = I(x, f × e -ρ(κ-•) , g × β, κ) I(x, e -ρ(κ-•) , β, κ) = E x f (τ B ∧ κ)g(B τ B ∧κ )M τ B ∧κ E x [M τ B ∧κ ]
,

where M t is the exponential martingale defined in (3.12). It is actually important to note that since the time interval is bounded (upper-bounded by κ)

we can use the time reversal expression e -ρ(κ-•) and the ratio in the previous equation verifies a cancellation of the terms e -ρκ . Finally it suffices therefore to use the Girsanov transformation in order to obtain the announced result:

E[f (T )g(Z)] = E x [f (τ X ∧ κ)g(X τ X ∧κ )].
Let us now focus our attention on the exact simulation of (X τ a,b , τ a,b ) for any diffusion process with regular drift term and constant diffusion coefficient. We have already seen that Algorithm (DET) permits to reach such objective however it is restricted to drift terms satisfying µ + µ 2 ≥ 0 on the interval [a, b]. If such a condition is not verified, we propose the following procedure: first we choose some time parameter κ > 0, then we apply: GENERAL DIFFUSION EXIT TIME (GDET) parameters: κ, a and b.

Initialization. Z = x, T = 0, N it = 0

While Z / ∈ {a, b} do:

• Simulate (S, Y ) = κ-DET(Z, κ) • T ← T + S, Z ← Y and N it ← N it + 1.

End While

Outcome: the couple of random variables (Z, T ) and the number of iterations N it .

Due to the Markov property, it is obvious that the outcome of such an algorithm and (X τ a,b (X) , τ a,b (X)) are identically distributed. Moreover the number of iterations of GDET has the same distribution as

τ a,b (X) κ + 1. Hence E x [N it ] ≤ 1 + E x [τ a,b (X)] κ .
It is evident that the number of iterations decreases as κ becomes large but we should be careful for a clever choice of κ since the number of rejections in Algorithm (κ-DET) grows exponentially fast when κ enlarges, see Proposition 3.4. A reasonable choice is therefore κ ≈ ρ -1 .

Examples and numerics

The aim of this section is to emphasize the efficiency of Algorithm (DET) and Algorithm (κ-DET) through the analysis of two examples. The first situation concerns a diffusion whose drift term b satisfies the condition µ + µ 2 ≥ 0 and consequently only requires the basic Algorithm (DET). The second situation concerns the Ornstein-Uhlenbeck process which plays an essential role in several applications namely in neuroscience. For both examples, we set the parameters appearing in the algorithms: t c = 0.7 and t e = 0.5.

Example with the Algorithm (DET)

We first consider a stochastic differential equation which was already presented in [START_REF] Herrmann | Exact simulation of the first-passage time of diffusions[END_REF] for the simulation of the first passage time. Here the objective is clearly different since we focus our attention to the exit time and exit position of the diffusion and the algorithm is different too. Let us also note that a similar diffusion process was also introduced in [START_REF] Beskos | Exact simulation of diffusions[END_REF].

We consider the following stochastic differential equation: dX t = (2 + sin(X t )) dt + dB t , t ≥ 0, X 0 = 0.

(3.17)

We first observe that γ(x) = (µ 2 (x) + µ (x))/2 = ((2 + sin(x)) 2 + cos(x))/2 satisfies 0 ≤ γ ≤ 5. We deduce that we can apply Theorem 3.1 with the particular choice ρ = 0: the outcome (Z, T ) of Algorithm (DET) has therefore the same distribution as (X τ a,b (X) , τ a,b (X)). The algorithm permits to obtain the histograms of respectively the exit time (Fig. 7 When considering larger intervals like for instance [-1, 2] (Fig. 8), the counter becomes large (average: 1205) and the algorithm DET is rather time consuming (C++ progaming: CPU 1,77 sec for 100 000 samples of the exit from the interval [-0.5, 0.5] and CPU 230,8 sec for the interval [-1, 2]). Figure 9 is an illustration of the high level of rejection for the algorithm (DET) as the interval size increases. Let us now consider the Ornstein-Uhlenbeck process given by the following SDE dX t = -µ 0 X t dt + dB t , X 0 = 0.

(3.18)

We shall determine the distribution of the exit time and location from the interval [-a, a] with a > 0. The drift term is defined by µ(x) = -µ 0 x, we can therefore find ρ such that Since the coefficient ρ is strictly positive (for positive µ 0 > 0), we cannot use the (DET)-algorithm in order to simulate the exit time. We shall therefore use the (κ-DET) algorithm in order to simulate exactly the couple (X τ-a,a(X)∧κ , τ -a,a (X)∧ κ) for any given constant time κ > 0. Fig. 10 represents the distribution of X τ-a,a(X)∧κ and Fig. 11 illustrates the time distribution and describes the counter values N tot . In order to simulate exactly the exit time for the Ornstein- Uhlenbeck process, we use the (GDET)-algorithm. The histogram in Figure 12 emphasizes the distribution of the exit time for the particular case: µ 0 = 2 and a = 1 and the efficiency of the algorithm for such a simulation. We can easily observe that the exit time increases as the parameter µ 0 increases, see Fig. 13. This result can be generalized for any K ≤ κ and therefore we obtain .

Figure 1 :

 1 Figure 1: Decomposition of the Brownian paths

4 1

 4 Brownian exit time of the interval [a, b] as the Brownian trajectory starts in x. The corresponding algorithm, presented in Section 2, is denoted BROWNIAN_EXIT_ASYMM (x,[a,b]) . The two main algorithms (DET) and (GDET) presented in Section 3 use these two previous algorithms as basic components. They permit the exact simulation of diffusion exit times under weak conditions on the coefficient µ. Efficiency results are pointed out (Theorem 3.2 and Proposition 3.4) and numerical illustrations complete the study in Section 3.Brownian motion constrained to stay in an interval [a, b] Let us first consider the Brownian motion in the space interval [a, b]. We need to describe the distribution of B t for a given t > 0, given that the first exit of the interval τ a,b (B) is larger than t. We fix the starting position of the Brownian motion B 0 = y ∈ (a, b). Due to both the translation invariance and the scaling property of the Brownian motion, we get

Figure 2 :

 2 Figure 2: Comparison between the constants κ(t, x) obtained either from the series expansion (1.5) or (1.6) as time elapses ad for different initial values x ∈ [0, 1).

Figure 3 :

 3 Figure 3: Conditional distribution (probability density function and normalized histogram) for x = 0.5, t = 0.2 (left) or t = 1 (right) and 100 000 simulations.

Figure 4 :

 4 Figure 4: Average value of the algorithm counter Nc versus the time variable, for the starting position x = 0.5 (solid line) and x = 0.2 (dashed line) with tc = 0.7 and 100 000 simulations.

Figure 5 :

 5 Figure 5: The estimated average value of the counter Ns versus the barrier te with 100 000 simulations for the estimation (stars) and the upper-bound function given by (2.9) (solid line).

. 10 )

 10 Combining the definition of R 1 in (2.2) and the change of variable w =

1 ,

 1 1} then we set T = T 1 and X = B x T1 else we start a new simulation for the exit time and position of a Brownian motion Bx starting in x = B T1 from the largest symmetric interval centered in x and included in [-1, 1]. This exit time is denoted by T 2 . • As above, if Bx T2 ∈ {-1, 1} then we set X = Bx T2 and T = T 1 + T 2 else we start again with a new initial position in the interval [-1, 1]. BROWNIAN EXIT TIME FOR ASYMMETRIC INTERVALS [a, b] BROWNIAN_EXIT_ASYMM Input: x (initial value of the Brownian paths) and [a, b].

  for the standard Brownian exit time of the asymmetric interval [-1.5, 2]. Of course the algorithm is not restricted to the standard Brownian case, it takes into account any initial position x belonging to the interval [a, b].

2 + 1 . 3 . 1 .

 2131 and T ← T + S else go to Step 0 end if, else • simulate the Brownian location at time E given that the exit time is larger than E: (Y c , N c ) = CONDITIONAL_DISTR (Z,[a,b],E), and set N tot ← N tot + N c . • if γ 0 V ≤ γ(Y c ) then go to Step 0 else set Z ← Y c and T ← T + E. end if. End While Outcome: the exit location Z, the exit time T of the interval [a, b] and the efficiency index N tot . Let ρ ≥ 0 such that γ(x) := µ 2 (x)+µ (x) ρ is a non negative function on the interval [a, b]. We introduce γ + = sup x∈[a,b] γ(x) and define ∆ = b a µ(y)dy and β(x) = min 1, e -∆ exp x a µ(y) dy. (3.3) Let us just notice that 0 ≤ β(a) ≤ 1 and 0 ≤ β(b) ≤ Theorem The outcome (Z, T ) of the Algorithm (DET) with parameter γ 0 = γ + and input functions γ(•) and β(•) satisfies: for any non-negative measurable functions f and g,

  τ B means τ a,b (B) for notational simplicity and B represents a Brownian motion starting in x ∈ [a, b] and γ a non negative continuous function. Let us now introduce a Poisson point process N on the space R + × [0, γ + ] with Lebesgue intensity measure and independent of the Brownian paths. If we reorganize the points of the Poisson process with respect to the abscissa, we obtain a sequence of points defined by (ξ n , U n ) n≥1 where (U n ) are independent uniformly distributed r.v. on [0, γ + ] and ξ n := n j=1 e j with (e j ) j≥1 a sequence of independent exponentially distributed r.v. with parameter γ + . Both sequences (U n ) n and (e n ) n are independent. Therefore, for any subset A of R + × [0, γ + ], we get N (A) = n≥1 1 {(ξn,Un)∈A} and P(N (A) = 0) = e -λ(A) ,

  ρτ a,b (X) ]E (a+b)/2 [e -γ+τ a,b (B) ] , ∀x ∈]a, b[. We recall that τ a,b stands for the first exit time of the interval [a, b]. The parameters ρ and γ + are defined in the introduction of Theorem 3.1.

1 .

 1 Let us now describe E[N x,1 tot ]. Let (S, Y, N as ) stands for the result of the first use of the function BROWNIAN_EXIT_ASYMM and (Y c , N c ) of the first use of CONDITIONAL_DISTR , we can therefore distinguish two different cases.

. 14 )•

 14 Let us note that P(S > E) = P x (τ a,b (B) > ξ) where τ a,b (B) is the exit time of the Brownian motion from the interval [a, b] and ξ is exponentially distributed with parameter γ + . Using the scaling property (1.1), we have P x (τ a,b (B) > ξ) = P y (τ > 4ξ(b -a) -2 ), where τ is the first Brownian exit time of the normalized interval [-1, 1] and y = 2x-a-b b-a . Since y → P y (τ > 4ξ(b -a) -2

  [N c |E = t] where N c is the cost of the function CONDITIONAL_DISTR (x, [a, b], t). Due to the scaling property (1.1), we know that the cost of this function is identical as the cost of CONDITIONAL_DISTR (y, [-1, 1], 4t(b-a) -2 ). We choose a parameter t c > 0 such that 4t(b -a) -2 ≤ t c corresponds to the algorithm for small time values and 4t(b -a) -2 > t c corresponds to large times. Combining Proposition 1.1 and Proposition 1.2 permits to have the following bound: for any x ∈ [a, b],

  X) ∧ κ) for any κ > 0. The simulation of (X τ a,b (X) , τ a,b (X)) can then be obtained by iteration due to the diffusion Markov property.Let us present Algorithm (κ-DET).MODIFIED DIFFUSION EXIT TIME (κ-DET) Parameters: ρ and γ 0 , input functions γ(•) and β m (•) First initialization: N tot = 0.

  left) and the counter N tot illustrating the efficiency of the algorithm DET (right). These histograms use a sample of size 100 000 and concerns the interval [a, b] = [-0.5, 0.5]. The average value of the counter is 8.5 and its estimated standard deviation is 8.88. Let us just compare the approximated results obtained by Algorithm (DET) for the interval [a, b] = [-0.5, 0.5] with a classical Euler method with step size 0.0001 and 100 000 samples. Algo. E[τ a,b ] σ(τ a,b ) E[τ a,b 1 {Xτ a,b =a} ] P(X τ a,b = a)

Figure 7 :

 7 Figure 7: Simulation of the first exit time from the interval [a, b] = [-1/2, 1/2] for the diffusion (3.17) (100 000 samples). Histograms of the exit time variable (left) and of the counter Ntot (right).

Figure 8 :

 8 Figure 8: Simulation of the first exit time from the interval [a, b] = [-1, 2] for the diffusion (3.17) (100 000 samples). Histograms of the exit time variable (left) and of the counter Ntot (right).

Figure 9 :

 9 Figure 9: FET from the interval [a, b] = [-1, 2] for the diffusion (3.17) (100 000 samples). Empirical c.d.f. of the exit times when the exit occurs at the top -lower curve -or at the bottom of the interval -upper curve -(left). Simulation of the FET from the interval [-a, a] for the diffusion (3.17) (10 000 samples). Average counter E[Ntot] in logarithmic scale versus a (right).

3. 4 . 2

 42 The Ornstein-Uhlenbeck case.

It suffices to choose ρ = µ0 2 .

 2 Such a choice implies sup x∈[-a,a] γ(x) ≤ µ 2 0 a 2 2 .

Figure 10 :

 10 Figure 10: Distribution of X τ -a,a (X)∧κ with a sample of size 10 000, for the drift parameter µ0 = 2, the interval size a = 1 and the constant time upper-bound κ = 0.5: histogram (left) and cumulative distribution (right).

Figure 11 :

 11 Figure 11: Cumulative distribution of τ-a,a(X) ∧ κ with a sample of size 10 000, for the drift parameter µ0 = 2, the interval size a = 1 and the constant time upper-bound κ = 0.5 (left, blue curve) or κ = 1 (left, black curve) and the corresponding histogram of the counter Ntot for κ = 0.5 (right).

Figure 12 :

 12 Figure 12: Histogram of the exit time distribution for the Ornstein-Uhlenbeck process with parameter µ0 = 2 and interval [-1, 1] (left), histogram of the associated counter Ntot (right).

Figure 13 :

 13 Figure 13: Cumulative distribution on the time interval [0, 15] for the O.-U.-exit time from the interval [-1, 1]

Figure 14 :

 14 Figure 14: Average counter for the (GDET)-algorithm versus the intensity of the drift (µ0 varies between 0.3 and 6) in logarithmic scale (time constant: κ = 0.5, sample size: 10 000) (left), average counter versus the time parameter κ used in the (GDET)-algorithm (right).

P

  (a+b)/2 (S ∧ κ ≤ E) sup x∈[a,b], K≤κ E[N x,1 K,tot ] ≤ C 0 (t e ) + C 1 (t c ). (3.21) Let us note that P (a+b)/2 (S ∧ κ ≤ E) = P (a+b)/2 (τ a,b (B) ∧ κ ≤ ξ) where τ a,b (B)is the exit time of the Brownian motion from the interval [a, b] and ξ is exponentially distributed with parameter γ + . Using the scaling property (1.1), we haveP (a+b)/2 (τ a,b (B) ∧ κ ≤ ξ) = E (a+b)/2 [e -γ+(τ a,b (B)∧κ) ].(3.22) Combining (3.21) (3.22) and (3.19), we obtain E[N x κ,tot ] ≤ e ρκ E[N x,1 κ,tot ] β m (x) ≤ e ρκ C 0 (t e ) + C 1 (t c ) β m (x) E (a+b)/2 [e -γ+(τ a,b (B)∧κ) ]

4529458

  Proposition 2.1. Let us note that N s is the random number of iterations of Step 3 used in the previous algorithm called BROWNIAN_EXIT_SYMMETRIC in order to simulate just one random variable Z with density p τ . Then

	E[N s ] ≤	t e 2π	e -1 2te +	3 2	erfc	1 √ 2t e	+	4 π	e -π 2 te 8	+	5π 4	e -25 π 2 te 8

  Let us note that A 1 (t e ) becomes small as t e becomes small. Let us now focus our attention to A 2 defined in (2.10).

	A 2 (t e ) = 4	k≥0	exp(-(4k + 1) 2 π 2 te 8 ) (4k + 1)π
	≤	4 π	e -π 2 te 8	+	4 5π	e -25 π 2 te 8	k≥1	exp -((4k + 1) 2 -5 2 )	π 2 t e 8
	≤	4 π	e -π 2 te 8	+	4 5π	e -25 π 2 te 8	k≥0	e -5kπ 2 te
	=	4 π	e -π 2 te 8	+	4 5π	e -25 π 2 te 8 1 -e -5π 2 te .	(2.14)
	Combining (2.13) and (2.14) permits to obtain the announced bound.
										.12)
	Combining (2.11) and (2.12) we obtain the following upper-bound:
					A 1 (t e ) ≤	t e 2π	e -1 2te +	3 2	erfc	1 √ 2t e	.	(2.13)

  BROWNIAN_EXIT_ASYMM (y, [-1, +1]) which satisfies due to Corollary 2.2:E[N as ] ≤ C 0 (t e ) with

	C 0 (t e ) :=	2t e π	e -1 2te + 3 erfc	1 √ 2t e	+	8 π	e -π 2 te 8	+	5π 8	e -25 π 2 te 8

Acknowledgment: We thank Antoine Lejay for helpful discussions about the Brownian exit times and their simulation. We particularly appreciated his working paper [17].

• simulate the conditional Brownian location at time E:

(Y c , N c ) = CONDITIONAL_DISTR (Z,[a,b],E) and set N tot ← N tot + N c .

• if V ≤ γ(Y c ) then go to Step 0 else set Z ← Y c , T ← T + E and K ← K -E.

end if.

End While

Outcome: the stopping time T corresponding to the minimum between the diffusion exit time of the interval [a, b] and the constant time κ, the location Z of the stopped diffusion at time T and the efficiency index N tot .

+ ρ is a non negative function on the interval [a, b]. We recall that γ + and β are defined in (3.3). We introduce Proposition 3.4. The random variable N x tot which is one of the outcomes of Algorithm (κ-DET) satisfies the following bound: there exists a constant C(t c , t e ) > 0 independent of both the interval [a, b] and the starting position x, such that

Appendix

Proof of Proposition 3.4. The counter N tot depends on the parameter κ. So we set N tot = N κ,tot . This number can be decomposed as follows:

where N x,k κ,tot represents the number of counter increases observed in-between the k-th and (k + 1)-th passage through the item Step 0. Since Algorithm (κ-DET) is an acceptance-rejection algorithm, the random variable N 0 is geometrically distributed. Let us also note that N x,k κ,tot = 0 a.s. on the event {N 0 < k} and conditionally to {N 0 ≥ k}, N x,k κ,tot has the same distribution as N x,1 κ,tot . Hence

since

the martingale (M t ) being defined in (3.12).

Let us now describe E[N x,1 κ,tot ]. Let (S, Y, N as ) stands for the result of the first use of the function BROWNIAN_EXIT_ASYMM and (Y t c , N t c ) of the first use of CONDITIONAL_DISTR at time t, we can therefore distinguish three different cases.

κ,tot = N as .

is an independent copy of N x,1 κ,tot . We deduce


Pierre David

Mikaël Le Pendu

Christine Guillemot

SPARSE TO DENSE SCENE FLOW ESTIMATION FROM LIGHT FIELDS

Keywords: Scene flow, optical flow, light field, sparse to dense

The paper addresses the problem of scene flow estimation from sparsely sampled video light fields. The scene flow estimation method is based on an affine model in the 4D ray space that allows us to estimate a dense flow from sparse estimates in 4D clusters. A dataset of synthetic video light fields created for assessing scene flow estimation techniques is also described. Experiments show that the proposed method gives error rates on the optical flow components that are comparable to those obtained with state of the art optical flow estimation methods, while computing a more accurate disparity variation when compared with prior scene flow estimation techniques.

INTRODUCTION

Light fields have proven to be very useful for scene analysis, and in particular for scene depth estimation [START_REF] Wanner | Variational light field analysis for disparity estimation and super-resolution[END_REF][START_REF] Zhang | Robust depth estimation for light field via spinning parallelogram operator[END_REF][START_REF] Jeon | Accurate depth map estimation from a lenslet light field camera[END_REF][START_REF] Huang | Empirical bayesian light-field stereo matching by robust pseudo random field modeling[END_REF][START_REF] Jiang | Depth estimation with occlusion handling from a sparse set of light field views[END_REF]. However, while a lot of work has been dedicated to scene analysis from static light fields, very little effort has been invested on the problem of dynamic scene analysis from video light fields. The concept of scene flow has first been defined in [START_REF] Vedula | Three-dimensional scene flow[END_REF] as describing the 3D geometry as well as the 3D motion of each scene point. Considering a multi-view set-up, the scene flow is estimated using an optical flow estimator for each view. The 3D scene flow is then computed by fitting its projection on each view to the estimated optical flows. However, in the recent literature (e.g. [START_REF] Menze | Object scene flow for autonomous vehicles[END_REF][START_REF] Vogel | 3d scene flow estimation with a piecewise rigid scene model[END_REF][START_REF] Behl | Bounding boxes, segmentations and object coordinates: How important is recognition for 3d scene flow estimation in autonomous driving scenarios?[END_REF][START_REF] Sun | Layered rgbd scene flow estimation[END_REF]), the scene flow is instead defined as a direct extension of the optical flow, where the depth (or disparity) variation ∆d of objects along time is represented in addition to the apparent 2D motion (∆x, ∆y).

Considering this definition, the scene flow analysis problem has first been addressed for stereo video sequences: in [START_REF] Menze | Object scene flow for autonomous vehicles[END_REF][START_REF] Vogel | 3d scene flow estimation with a piecewise rigid scene model[END_REF][START_REF] Behl | Bounding boxes, segmentations and object coordinates: How important is recognition for 3d scene flow estimation in autonomous driving scenarios?[END_REF] a scene flow is estimated assuming that the scene can be decomposed into rigidly moving objects. Several methods based on RGB-D videos have also been developed [START_REF] Sun | Layered rgbd scene flow estimation[END_REF], [START_REF] Quiroga | Local/global scene flow estimation[END_REF], [START_REF] Jaimez | Fast odometry and scene flow from rgb-d cameras based on geometric clustering[END_REF]. The first methods for scene flow analysis from dense light fields have been proposed in [START_REF] Heber | Scene flow estimation from light fields via the preconditioned primal-dual algorithm[END_REF] and [START_REF] Navarro | Variational scene flow and occlusion detection from a light field sequence[END_REF]. Both are based on variational models. The authors in [START_REF] Srinivasan | Oriented light-field windows for scene flow[END_REF] proposed oriented light field windows to estimate the scene flow. All these methods rely on epipolar images which only provide sufficient information for scene flow estimation of densely sampled light fields (as those captured with plenoptic cameras). However, they do not address the case of sparse light fields (i.e. with large baselines).

In this paper, we focus on the problem of scene flow analysis from sparse video light fields. This problem is made difficult due to the large temporal and angular occlusions. To cope with this difficulty, we propose a method for interpolating, in the 4D ray space, a sparse scene flow into a dense one. The sparse-to-dense approach naturally handles occlusions: while the sparse estimation can only be obtained on non-occluded parts of the image, the dense interpolation extends these estimations to every pixel (occluded or not). The proposed method is based on an affine model in the 4D ray space that allows us to estimate the dense flow from sparse estimates in nearby clusters. Note that an advantage of our sparse-to-dense approach is the possibility to use it as a post-processing step in other scene flow estimation methods, in the same way as EpicFlow [START_REF] Revaud | Epicflow: Edge-preserving interpolation of correspondences for optical flow[END_REF] has been used in other methods (e.g. [START_REF] Bailer | Flow fields: Dense correspondence fields for highly accurate large displacement optical flow estimation[END_REF]) to interpolate the optical flow of pixels detected as outliers.

In order to validate the proposed algorithm, we have generated a synthetic light field video dataset based on the Sintel movie (used in the optical flow benchmark [START_REF] Butler | A naturalistic open source movie for optical flow evaluation[END_REF][START_REF] Wulff | Lessons and insights from creating a synthetic optical flow benchmark[END_REF]). This dataset is composed of two synthetic light fields ('Bamboo2' and 'Temple1') of 3 × 3 views of 1024 × 536 pixels and 50 frames. The light field views are provided with the corresponding ground truth scene flow (optical flow and disparity variation).

The proposed method has been assessed using this synthetic dataset in comparison with the oriented window method in [START_REF] Srinivasan | Oriented light-field windows for scene flow[END_REF], and with two methods proposed for optical flow estimation: the sparse to dense estimation method EpicFlow [START_REF] Revaud | Epicflow: Edge-preserving interpolation of correspondences for optical flow[END_REF] and a state-of-the-art technique based on a deep learning architecture called PWC-Net [START_REF] Sun | PWC-Net: CNNs for optical flow using pyramid, warping, and cost volume[END_REF]. Note that, since the depth estimation component of the method in [START_REF] Srinivasan | Oriented light-field windows for scene flow[END_REF] has been designed for dense light fields and is hardly applicable when the baseline is large, when using this method, we have coupled the optical flow estimation performed by the method with the ground truth disparity. The two optical flow estimation techniques [START_REF] Sun | PWC-Net: CNNs for optical flow using pyramid, warping, and cost volume[END_REF] and [START_REF] Revaud | Epicflow: Edge-preserving interpolation of correspondences for optical flow[END_REF] are used for separately estimating the optical flow in each view as well as the disparity between views. Experimental results show that the proposed method gives optical flow endpoint error rates that are comparable to those obtained with state of the art optical flow estimation methods, while computing at the same time the disparity variation with lower mean absolute errors. Although the accuracy of our optical flow is slightly lower than that of PWC-Net, the depth variation maps estimated by our method have considerably lower errors.

NOTATIONS AND METHOD OVERVIEW

Let us consider the 4D representation of Light Fields proposed in [START_REF] Levoy | Light field rendering[END_REF] and [START_REF] Gortler | The lumigraph[END_REF] to describe the radiance along the different light rays. This 4D function, at each time instant, is denoted LF t (x, y, u, v), where the pairs x = (x, y) and u = (u, v) respectively denote the spatial and angular coordinates of light rays. A view u of a light field at t is written L t u . We denote by C : (I, I ) → (S, S , f ) the function that determines two sets S and S of matching points in I and I respectively. The associated bijection f is such that f (S) = S .

f2 f 0 f4 f5 f3 f1 L t u L t c L t+1 u L t+1 c u c x u x LF t LF t+1
Fig. 2. Matching between a view u and the central view c at t and t + 1, shown in epipolar images at t and t+1.

The inverse bijection is noted f -1 . The proposed method takes as inputs two consecutive frames of a video light field LF t (x, y, u, v), LF t+1 (x, y, u, v), a matching function C and estimates a consistent dense scene flow on each view. It produces disparity maps at t and t+1 and a disparity variation map between t and t+1 (in addition to the optical flow components). The algorithm proceeds as follows (see Fig. 1). The first step of the method is the robust sparse scene flow estimation. All light field rays are then grouped in 4D clusters to guide the sparse to dense interpolation of the scene flow. An affine model is computed from the sparse flow estimates in neighboring clusters, and then used to compute a flow value in each light field pixel.

SPARSE SCENE FLOW ESTIMATION

The method pivots on the central view L t c . We first compute (S0, S 0 , f0) = C(L t c , L t+1 c ). Then, for every non central view u, we perform three kinds of matching (see Fig. 2), angular matching to estimate the sets of points (S1,

S 1 , f1) = C(L t+1 c , L t+1 u ) and (S3, S 3 , f3) = C(L t u , L t c ), temporal matching to estimate the sets of points (S2, S 2 , f2) = C(L t+1 u , L t u ), and temporo- angular matching to estimate (S4, S 4 , f4) = C(L t c , L t+1 u ) and (S5, S 5 , f5) = C(L t+1 c , L t u ).
In the experiments, we used [START_REF] Li | Coarse-to-fine patchmatch for dense correspondence[END_REF] as function C.

The next step is to link the different sets together based on the distance between the points computed as in [START_REF] Achanta | Slic superpixels compared to state-of-theart superpixel methods[END_REF]. We define a distance D based on color, angular and spatial proximity in the 7D space [labuvxy] as :

D (Pi, Pj) = dc 2 + ms 2 S 2 ds 2 + ma 2 A 2 da 2 (1)
with :

S = h × w/|S0| A = √ Nu × Nv (2)
where dc, ds and da are the color, spatial and angular distances defined as euclidean distances in the CIELAB colorspace, the [xy] space and the [uv] space respectively. The variables w and h are the width and height of a view, Nu and Nv are the angular size of the light field and ms is a parameter controlling the balance between color and spatial distances that we fix to 10 in the experiments. For this step, we have da = 0 because we are aiming at linking points in the same view, but this distance will be used later in the clustering step. Based on the above distances, for each point P0 ∈ S0, we build a chain of points (P0, P 0 , P1, P 1 , P2, P 2 , P3, P 3 ) by searching for the nearest neighbors successively in S1, S2 and S3 (see Eq. ( 3) and Fig. 3a). We call this step a forward nearest neighbor search and denote NN the simple nearest neighbor search function for a point Pi ∈ Si in a set Sj, according to the distance D.

P0 f 0 -→ P 0 NN -→ P1 f 1 -→ P 1 NN -→ P2 f 2 -→ P 2 NN -→ P3 f 3 -→ P 3 (3)
To ensure the robustness of the chain of points, we then perform a backward nearest neighbor check as in Eq.( 4). Given a point P i ∈ S i whose nearest neighbor in Sj is Pj, we check if P i is reciprocally the nearest neighbor of Pj in S i.

P 3 f -1 3 -→ P3 NN? -→ P 2 f -1 2 -→ P2 NN? -→ P 1 f -1 1 -→ P1 NN? -→ P 0 f -1 0 -→ P0 (4)
Let A1 be the set of points P0 ∈ S0 that passed the backward nearest neighbor test. We have A1 ⊂ S0. We then perform a similar forward nearest neighbor search and backward nearest neighbor test for each point P0 ∈ S0 but this time with the two temporo-angular matching sets S4 and S5, building the chain of points (P0, P 0 , P4) and (P0, P 0 , P5) (see Fig. 3b).

Let A2 be the set of points P0 ∈ S0 that passed the backward nearest neighbor test in S4 and S5. We have A2 ⊂ S0. For each point P0 ∈ A = A1 ∩ A2, we build a complete chain of points g = (P0, P1, P2, P3, P4, P5).

Given a chain of points, we check the consistency of each matching point P i = fi(Pi) with the point Pj in the same view and instant. More precisely, we discard the chain if either D(P0, P 3 ) > τ , D(P2, P 4 )) > τ , or D(P3, P 5 ) > τ , where τ is a distance threshold that was fixed to 10 in our experiments.

Any remaining chain g is used to estimate a sparse scene flow. From g, we can compute the scene flow on two points of the light field LF t : on the central view L t c at P0 and on the view L t u at P 2 :

P 0 -P0 = ∆xc ∆yc and P 2 -P2 = - ∆xu ∆yu (5)
Let dt and dt+1 be the disparities at t and t + 1. We have:

P 1 -P1 = (u -c) • dt+1 and P 3 -P3 = -(u -c) • dt (6)
With the temporo-angular matching, we also have:

P 4 -P4 = -(P 2 -P2) -(P 3 -P3) = (P 0 -P0) + (P 1 -P1) (7)
P 5 -P5 = -(P 0 -P0) -(P 3 -P3) = (P 1 -P1) + (P 2 -P2)
We solve this linear overdetermined system with a least square method in order to determine the values ∆xc, ∆yc, ∆xu, ∆yu, dt and dt+1 for each chain g.

L t c L t+1 c L t u L t+1 u f0 S0 S 0 f1 S1 S 1 f2 S2 S 2 f3 S3 S 3 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • NN NN?
(a) Chain of angular and temporal correspondences

L t c L t+1 c L t u L t+1 u f0 S0 S 0 f5 S5 S 5 f4 S4 S 4 • • • • • • • • • • • • • • • • • • • • • • • • NN NN?
(b) Chain of temporo-angular correspondences Fig. 3. Forward nearest neighbor search and backward nearest neighbor check. The green dots represent the points that passed the backward nearest neighbor check.

CLUSTERING THE LIGHT FIELD

In order to estimate a dense scene flow from sparse estimates, we need to make an assumption about the scene. We assume that rays which have similar colors and are located near each other in the light field should have similar motions. As a consequence, similarly to [START_REF] Hu | Robust interpolation of correspondences for large displacement optical flow[END_REF] but extended to 4D light fields, we propose to decompose the light field into clusters based on color, spatial and angular proximity. Then, the clustering will be used to estimate one scene flow model per cluster. We generalize in 4D the SLIC superpixels developed in [START_REF] Achanta | Slic superpixels compared to state-of-theart superpixel methods[END_REF] to cluster the light field LF t . We use as distance D , already defined in Eq.( 1). We only modify the definition of S :

S = h × w/K (8)
K denotes the number of desired clusters. We use the parameters ms and ma to control the balance between spatial, angular and color proximity. S (in Eq. ( 8)) and A (in Eq. ( 2)) are the maximum spatial and angular distances expected in a cluster. The spatial positions of the centroids C k of the clusters are initialized on a regular grid of step S, while their initial angular positions are randomly sampled. Then, each pixel i is assigned to the nearest cluster center C k . We fix the search region to Nu × Nv × 2S × 2S because the expected approximate spatial size of each cluster is S × S and the approximate angular size is Nu × Nv. Indeed, each cluster is excpected to contain pixels in all the Nu × Nv views.

Once the assignment step is over, we update the cluster centers by computing the new centroid of each cluster. We only need to iterate this second step Niteration = 10 in order to have stabilized cluster centroids.

SPARSE TO DENSE INTERPOLATION

Once we have divided the light field into clusters, we estimate a model for each one of them. For each cluster, we look for the points estimated in Section 3 that are inside the cluster, we compute the mean value of the sparse estimates and associate it with the centroids of the cluster. For clusters that have no seeds inside them, we do not do anything. We end up with initialized clusters and unitialized ones.

As in [START_REF] Revaud | Epicflow: Edge-preserving interpolation of correspondences for optical flow[END_REF], we build a weighted graph where each node is a cluster (initialized or not) and where an edge between two nodes means that the two clusters are adjacent. The associated weight is defined as the distance D (as in Eq.( 1)) between the centroids of the two clusters. We then look for the N nearest initialized neighbors using Dijkstra's algorithm on the graph [START_REF] Dijkstra | A note on two problems in connexion with graphs[END_REF], discarding every empty cluster.

For each cluster, we estimate the parameters of the affine model in Eq.( 9) in the 4D space by fitting the model on the initialized centroids using the RANSAC algorithm [START_REF] Fischler | Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography[END_REF].

a1u + b1v + c1x + d1y + e1 = x + ∆x a2u + b2v + c2x + d2y + e2 = y + ∆y (9) a3u + b3v + c3x + d3y + e3 = ∆d
We then apply the model to every pixel belonging to the cluster to compute its (∆x, ∆y, ∆d) scene flow value.

As a final post-processing step, we perform an energy minimization like in [START_REF] Revaud | Epicflow: Edge-preserving interpolation of correspondences for optical flow[END_REF] for the optical flow of each subaperture image independently. To regularize the disparity variation map, we then perform joint bilateral filtering using the optical flow as a guide.

EVALUATION

Scene Flow Dataset

For our experiments, we have prepared a synthetic video Light Field dataset 1 . For that purpose, we have used the production files of the open source movie Sintel [START_REF]Sintel[END_REF] and have modified them in the Blender 3D software [START_REF]Blender (3d software)[END_REF] in order to render an array of 3x3 views. Similarly to the MPI Sintel flow dataset [START_REF] Butler | A naturalistic open source movie for optical flow evaluation[END_REF][START_REF] Wulff | Lessons and insights from creating a synthetic optical flow benchmark[END_REF], we have modified the scenes to generate not only the 'final' render, but also a 'clean' render without lighting effects, motion blur, or semi-transparent objects. Ground truth optical flow and disparity maps were also generated for each view. Since disparity variation maps could not be rendered within Blender, we have computed them using the disparity and optical flow. However, this process requires projecting the disparity map of a frame to the next frame using the optical flow, which results in unavailable disparity variation information in areas of temporal occlusion. We have processed two scenes of 3 × 3 views of 1024 × 536 pixels and 50 frames corresponding to the scenes 'Bamboo2' and Epic [START_REF] Revaud | Epicflow: Edge-preserving interpolation of correspondences for optical flow[END_REF] PWC [START_REF] Sun | PWC-Net: CNNs for optical flow using pyramid, warping, and cost volume[END_REF] Fig. 4. Visual comparison of our method with [START_REF] Srinivasan | Oriented light-field windows for scene flow[END_REF][START_REF] Revaud | Epicflow: Edge-preserving interpolation of correspondences for optical flow[END_REF][START_REF] Sun | PWC-Net: CNNs for optical flow using pyramid, warping, and cost volume[END_REF]. The optical flows are visualized with the Middlebury color code, and the disparity variations are visualized using a gray-scale representation. The red pixels are the occlusion mask where there is no ground truth disparity variation available.

'Temple1' in [START_REF] Butler | A naturalistic open source movie for optical flow evaluation[END_REF]. The disparities (in pixels) between neighboring views are in the range [-8, +52] for 'Bamboo2' and [-22, +9] for 'Temple1'.

Results

The proposed method is first assessed in comparison with the method in [START_REF] Srinivasan | Oriented light-field windows for scene flow[END_REF], referred to here as OW (Oriented Window). This latter was designed for dense light fields captured with plenoptic cameras. However, the optical flow searched via the oriented window can be combined with disparity maps estimated by methods suitable for sparse light fields. In the test reported here, we used for this method ground truth disparity maps, thus showing the best results it can give for the estimated scene flow. We also compare the method with a naive approach that would consist in separately estimating the disparity maps at t and t+1, the optical flows between t and t+1 and finally computing the disparity variation as the difference of disparities along the optical flows. We consider two methods proposed for optical flow estimation: the sparse to dense estimation method EpicFlow [START_REF] Revaud | Epicflow: Edge-preserving interpolation of correspondences for optical flow[END_REF] and a state-of-the-art technique based on a deep learning architecture called PWC-Net [START_REF] Sun | PWC-Net: CNNs for optical flow using pyramid, warping, and cost volume[END_REF] Note that this separate disparity and optical flow estimation does not handle occlusions. So, the disparity variation in occluded or disoccluded areas will never be consistent. The results are summarized in Table 1. For each successive light field frame of the four sequences (Bamboo2 and Temple1, both rendered as clean and final), we compute the endpoint error for the optical flow (EPE OF) and separately the mean absolute error for the disparity variation (MAE ∆d). The latter is only computed for disoccluded pixels because there is no ground truth on occluded pixels. We compute these two errors on every ray of the light field (Global) and because the method proposed in [START_REF] Srinivasan | Oriented light-field windows for scene flow[END_REF] only gives the scene flow for the central view, we also compute the errors on the central view only. We can observe that although [START_REF] Sun | PWC-Net: CNNs for optical flow using pyramid, warping, and cost volume[END_REF] yields the most accurate optical flows, our method has a lower optical flow end point error than the two other methods and that it gets close to [START_REF] Sun | PWC-Net: CNNs for optical flow using pyramid, warping, and cost volume[END_REF] for some sequences like Temple1 clean and Bamboo2 clean. Regarding the disparity variation, our method provides much lower errors than any other tested method. In the end, we propose a method that is slightly less precise than state of the art optical flow methods but more accurate in terms of disparity variation.

CONCLUSION

In this paper, we have presented a new method to estimate scene flows from sparsely sampled video light fields. The method is based on three steps: first a sparse scene flow estimation, then a 4D clustering of the light field, and finally a sparse to dense scene flow interpolation for each cluster. For the performance evaluation, we have generated a synthetic dataset from the open source movie Sintel in order to extend the popular MPI Sintel dataset to light fields and scene flow. Our method gives comparable performances with stateof-the art approaches only considering the horizontal and vertical displacements (i.e. the optical flow). However, significant improvements were obtained in the estimation of the disparity variation.

Fig. 1 .

 1 Fig. 1. Overview of the sparse to dense scene flow estimation

Table 1 .

 1 EPE OF MAE ∆d EPE OF MAE ∆d EPE OF MAE ∆d Results on our light field scene flow dataset for the optical flow endpoint error EPE OF and for the disparity variation mean absolute error MAE ∆d. This latter is only computed on pixels that remain visible between two adjacent frames because of the lack of ground truth in occluded and disoccluded areas. The lowest errors are in red, the second lowest errors in orange.

				Bamboo2			Temple1
			clean		final		clean		final
	Ours EPE OF MAE ∆d Central 1.007 0.136 OW [15] 1.421 0.356 view Epic [16] 1.078 0.685 only PWC [20] 0.945 0.576	1.102 1.462 1.280 1.017	0.140 0.345 0.663 0.596	1.042 2.061 5.405 1.042	0.109 0.152 0.664 0.291	1.383 2.374 2.472 1.322	0.128 0.162 1.009 0.325
		Ours	1.090	0.140	1.169	0.142	1.109	0.111	1.453	0.131
	Global	OW [15] Epic [16]	/ 1.071	/ 1.020	/ 1.264	/ 1.015	/ 5.676	/ 0.730	/ 2.609	/ 1.208
		PWC [20]	0.947	0.873	1.017	0.884	1.038	0.385	1.323	0.429
	Central view	Ground Truth	Ours		OW[15]
	t	(∆x, ∆y)
	t + 1	∆d

http://clim.inria.fr/Datasets/SyntheticVideoLF/index.html

This work was supported by the EU H2020 Research and Innovation Program under grant agreement N • 694122 (ERC advanced grant CLIM).