Benjamin Jourdain
email: jourdain@cermics.enpc.fr

Mohamed Sbai
email: sbai@cermics.enpc.fr

Exact retrospective Monte Carlo computation of arithmetic average Asian options

Taking advantage of the recent literature on exact simulation algorithms (Beskos et al. [1]) and unbiased estimation of the expectation of certain functional integrals (Wagner [23], Beskos et al. [2] and Fearnhead et al. [6]), we apply an exact simulation based technique for pricing continuous arithmetic average Asian options in the Black & Scholes framework. Unlike existing Monte Carlo methods, we are no longer prone to the discretization bias resulting from the approximation of continuous time processes through discrete sampling. Numerical results of simulation studies are presented and variance reduction problems are considered.

Introduction

Although the Black & Scholes framework is very simple, it is still a challenging task to efficiently price Asian options. Since we do not know explicitly the distribution of the arithmetic sum of log-normal variables, there is no closed form solution for the price of an Asian option. By the early nineties, many researchers attempted to address this problem and hence different approaches were studied including analytic approximations (see Turnball and Wakeman [START_REF] Turnball | A quick algorithm for pricing european average options[END_REF], Vorst [START_REF] Vorst | Prices and hedge ratios of average exchange rate options[END_REF], Levy [START_REF] Levy | Pricing european average rate currency options[END_REF] and more recently Lord [START_REF] Lord | Partially exact and bounded approximations for arithmetic Asian options[END_REF]), PDE methods (see Vecer [START_REF] Vecer | A new pde approach for pricing arithmetic asian options[END_REF], Rogers and Shi [START_REF] Rogers | The value of an Asian option[END_REF], Ingersoll [START_REF] Ingersoll | Theory of Financial Decision Making[END_REF], Dubois and Lelievre [START_REF] Dubois | Efficient pricing of asian options by the pde approach[END_REF]), Laplace transform inversion methods (see Geman and Yor [START_REF] Geman | Bessel processes, asian option and perpetuities[END_REF], Geman and Eydeland [START_REF] Geman | Domino effect[END_REF]) and, of course, Monte Carlo simulation methods (see Kemna and Vorst [START_REF] Kemna | A pricing method for options based on average asset values[END_REF], Broadie and Glasserman [START_REF] Broadie | Estimating security price derivatives using simulation[END_REF], Fu et al. [START_REF] Fu | Pricing continuous asian options : a comparison of monte carlo and laplace transform inversion methods[END_REF]).

Monte Carlo simulation can be computationally expensive because of the usual statistical error. Variance reduction techniques are then essential to accelerate the convergence (one of the most efficient techniques is the Kemna&Vorst control variate based on the geometric average). One must also account for the inherent discretization bias resulting from approximating the continuous average of the stock price with a discrete one. It is crucial to choose with care the discretization scheme in order to have an accurate solution (see Lapeyre and Temam [START_REF] Lapeyre | Competitive Monte Carlo methods for pricing asian options[END_REF]). The main contribution of our work is to fully address this last feature by the use, after a suitable change of variables, of an exact simulation method inspired from the recent work of Beskos et al. [START_REF] Beskos | Retrospective exact simulation of diffusion sample paths[END_REF][START_REF] Beskos | Exact and computationally efficient likelihood-based estimation for discretely observed diffusion processes[END_REF] and Fearnhead et al. [START_REF] Paul Fearnhead | Particle filters for partially observed diffusions[END_REF].

In the first part of the paper, we recall the algorithm introduced by Beskos et al. [START_REF] Beskos | Retrospective exact simulation of diffusion sample paths[END_REF] in order to simulate sample-paths of processes solving one-dimensional stochastic differential equations. By a suitable change of variables, one may suppose that the diffusion coefficient is equal to one. Then, according to the Girsanov theorem, one may deal with the drift coefficient by introducing an exponential martingale weight. Because of the one-dimensional setting, the stochastic integral in this exponential weight is equal to a standard The constants α and β are two given non-negative parameters.

When α > 0, we remark that, by a change of variables inspired by Rogers and Shi [START_REF] Rogers | The value of an Asian option[END_REF], αS T + β T 0 S t dt has the same law as the solution at time T of a well-chosen one-dimensional stochastic differential equation. Then it is easy to implement the exact methods previously presented. The case α = 0 of standard Asian options is more intricate. The previous approach does not work and we propose a new change of variables which is singular at initial time. It is not possible to implement neither the exact simulation algorithm nor the method based on the unbiased estimator of Wagner [START_REF] Wagner | Unbiased Monte Carlo evaluation of certain functional integrals[END_REF] and we propose a pseudo-exact hybrid method which appears as an extension of the exact simulation algorithm. In both cases, one first replaces the integral with respect to the time variable in the function f by an integral with respect to time in the exponential function. Because of the nice properties of this last function, exact computation is possible.

1 Exact Simulation techniques 1.1 The exact simulation method of Beskos et al. [START_REF] Beskos | Retrospective exact simulation of diffusion sample paths[END_REF] In a recent paper, Beskos et al. [START_REF] Beskos | Retrospective exact simulation of diffusion sample paths[END_REF] proposed an algorithm which allows to simulate exactly the solution of a 1-dimensional stochastic differential equation. Under some hypotheses, they manage to implement an acceptance-rejection algorithm over the whole path of the solution, based on recursive simulation of a biased Brownian motion. Let us briefly recall their methodology. We refer to [START_REF] Beskos | Retrospective exact simulation of diffusion sample paths[END_REF] for the demonstrations and a detailed presentation.

Consider the stochastic process (ξ t ) 0≤t≤T determined as the solution of a general stochastic differential equation of the form :

dξ t = b(ξ t )dt + σ(ξ t )dW t ξ 0 = ξ ∈ Ê (1)
where b and σ are scalar functions satisfying the usual Lipschitz and growth conditions with σ non vanishing.

To simplify this equation, Beskos et al. [START_REF] Beskos | Retrospective exact simulation of diffusion sample paths[END_REF] suggest to use the following change of variables :

X t = η(ξ t )
where η is a primitive of 1 σ (η(x) =

x .

1 σ(u) du).

Under the additional assumption that 1 σ is continuously differentiable, one can apply Itô's lemma to get

dX t = η ′ (ξ t )dξ t + 1 2 η ′′ (ξ t ) d < ξ, ξ > t = b(ξ t ) σ(ξ t ) dt + dW t - σ ′ (ξ t ) 2 dt = b(η -1 (X t )) σ(η -1 (X t )) - σ ′ (η -1 (X t )) 2 a(Xt) dt + dW t So ξ t = η -1 (X t )
where (X t ) t is a solution of the stochastic differential equation

dX t = a(X t )dt + dW t X 0 = x. (2)
Thus, without loss of generality, one can start from equation (2) instead of (1).

Let us denote by (W x t ) t∈[0,T ] the process (W t + x) t∈[0,T ] , by É W x its law and by É X the law of the process

(X t ) t∈[0,T ]
. From now on, we will denote by (Y t ) t∈[0,T ] the canonical process, that is the coordinate mapping on the set C([0, T ], Ê) of real continuous maps on [0, T ] (see Revuz and Yor [17] or Karatzas and Shreve [START_REF] Karatzas | Brownian motion and stochastic calculus[END_REF]).

One needs the following assumption to be true Assumption 1 : Under É W x , the process

L t = exp t 0 a(Y u )dY u - 1 2 t 0 a 2 (Y u )du is a martingale.
According to Rydberg [START_REF] Rydberg | A note on the existence of unique equivalent martingale measures in a markovian setting[END_REF] (see the proof of Proposition 4 where we give his argument on a specific example), a sufficient condition for this assumption to hold is -Existence and uniqueness in law of a solution to the SDE (2).

-∀t ∈ [0, T ],

t 0 a 2 (Y u )du < ∞, É X and É W x almost surely on C([0, T ], Ê).
Thanks to this assumption, one can apply the Girsanov theorem to get that É X is absolutely continuous with respect to É W x and its Radon-Nikodym derivative is equal to

dÉ X dÉ W x = exp T 0 a(Y t )dY t - 1 2 T 0 a 2 (Y t )dt .
Consider A the primitive of the drift a, and assume that Assumption 2 : a is continuously differentiable.

Since, by Itô's lemma,

A(W x T ) = A(x) + T 0 a(W x t )dW x t + 1 2 T 0 a ′ (W x t )dt, we have dÉ X dÉ W x = exp A(Y T ) -A(x) - 1 2 T 0 a 2 (Y t ) + a ′ (Y t )dt .
Before setting up an acceptance-rejection algorithm using this Radon-Nikodym derivative, a last step is needed. To ensure the existence of a density h(u) proportional to exp(A(u)

-(u-x) 2
2T ), it is necessary and sufficient that the following assumption holds

Assumption 3 : The function u → exp(A(u) -(u-x) 2 2T ) is integrable.
Finally, let us define a process Z t distributed according to the following law É Z

É Z = Ê L (W x t ) t∈[0,T ] |W x T = y h(y)dy
where the notation L(.|.) stands for the conditional law. One has

dÉ X dÉ Z = dÉ X dÉ W x dÉ W x dÉ Z = C exp - 1 2 T 0 a 2 (Y t ) + a ′ (Y t )dt
where C is a normalizing constant. At this level, Beskos et al. [START_REF] Beskos | Retrospective exact simulation of diffusion sample paths[END_REF] need another assumption

Assumption 4 : The function φ : x → a 2 (x)+a ′ (x)
2 is bounded from below.

Therefore, one can find a lower bound k of this function and eventually the Radon-Nikodym derivative of the change of measure between X and Z takes the form

dÉ X dÉ Z = Ce -kT exp - T 0 φ(Y t ) -k dt .
The idea behind the exact algorithm is the following : suppose that one is able to simulate a continuous path Z t (ω) distributed according to É Z and let M (ω) be an upper bound of the mapping t → φ(Z t (ω)) -k.

Let N be an independent random variable which follows the Poisson distribution with parameter T M (ω) and let (U i , V i ) i=1...N be a sequence of independent random variables uniformly distributed on

[0, T ] × [0, M (ω)].
Then, the number of points (U i , V i ) which fall below the graph

{(t, φ(Z t (ω)) -k); t ∈ [0, T ]} is equal to zero with probability exp - T 0 φ(Z t (ω)) -k dt .
Actually, simulating the whole path (Z t ) t∈[0,T ] is not necessary. It is sufficient to determine an upper bound for φ(Z t )k since, as pointed out by the authors, it is possible to simulate recursively a Brownian motion on a bounded time interval by first simulating its endpoint, then simulating its minimum or its maximum and finally simulating the other points2 . For this reason, one needs the following assumption for the algorithm to be feasible :

Assumption 5 : Either lim sup u→+∞ φ(u) < +∞ or lim sup u→-∞ φ(u) < +∞.
Suppose for example that lim sup u→+∞ φ(u) < +∞. The exact algorithm of Bekos et al. [START_REF] Beskos | Retrospective exact simulation of diffusion sample paths[END_REF] then takes the following form : Algorithm 1 1. Draw the ending point Z T of the process Z with respect to the density h.

2.

Simulate the minimum m of the process Z given Z T .

Fix an upper bound

M (m) = sup{φ(u) -k; u ≥ m} for the mapping t → φ(Z t ) -k.

Draw N according to the Poisson distribution with parameter T M (m) and draw

(U i , V i ) i=1...N , a sequence of independent variables uniformly distributed on [0, T ] × [0, M (m)].
5. Fill in the path of Z at the remaining times (U i ) i=1...N .

Evaluate the number of points

(V i ) i=1...N such that V i ≤ φ(Z Ui ) -k.
If it is equal to zero, then return the simulated path Z.

Else, return to step 1.

This algorithm gives exact skeletons of the process X, solution of the SDE (2). Once accepted, a path can be further recursively simulated at additional times without any other acceptance/rejection criteria. We also point out that the same technique can be generalized by replacing the Brownian motion in the law of the proposal Z by any process that one is able to simulate recursively by first simulating its ending point, its minimum/maximum and then the other points. Also, the extension of the algorithm to the inhomogeneous case, where the drift coefficient a in (2), and therefore the function φ, depend on the time variable t, is straightforward given that the assumptions presented above are appropriately modified.

The unbiased estimator (U.E)

In finance, the pricing of contingent claims often comes down to the problem of computing an expectation of the form

C 0 = (f (X T )) ( 3
)
where X is a solution of the SDE (2) and f is a scalar function such that f (X T ) is square integrable. In a simulation based approach, one is usually unable to exhibit an explicit solution of this SDE and will therefore resort to numerical discretization schemes, such as the Euler or Milstein schemes, which introduce a bias. Of course, the exact algorithm presented above avoids this bias. Here, we are going to present a technique which permits to compute exactly the expectation (3) while assumptions 4 and 5 on the function a 2 +a ′ 2 which appears in the Radon-Nikodym derivative are relaxed.

Using the previous results and notations, we get, under the assumptions 1 and 2, that

C 0 = f (W x T ) exp A(W x T ) -A(x) - 1 2 T 0 a 2 (W x t ) + a ′ (W x t )dt . (4)
In order to implement an importance sampling method, let us introduce a positive density ρ on the real line and a process (Z t ) t∈[0,T ] distributed according to the following law É Z

É Z = Ê L (W x t ) t∈[0,T ] |W x T = y ρ(y)dy.
By (4), one has

C 0 = ψ(Z T ) exp - T 0 φ(Z t )dt (5)
where

ψ : z → f (z) e A(z)-A(x)-(z-x) 2 2T √ 2πρ(z)
and φ : z → a 2 (z)+a ′ (z)

2

. We do not impose ρ to be equal to the density h of the previous section. It is a free parameter chosen in such a way that it reduces the variance of the simulation.

In his first paper, Wagner [START_REF] Wagner | Unbiased Monte Carlo evaluation of certain functional integrals[END_REF] constructs an unbiased estimator of the expectation (5) when ψ is a constant, (Z t ) t∈[0,T ] is an Ê d -valued Markov process with known transition function and φ is a measurable function such that e T 0 |φ(Zt)|dt < +∞. His main idea is to expand the exponential term in a power series, then, using the transition function of the underlying Markov process and symmetry arguments, he constructs a signed measure ν on the space Y = +∞ n=0 ([0, T ] × Ê d ) n+1 such that the expectation at hand is equal to ν(Y). Consequently, any probability measure µ on Y that is absolutely continuous with respect to ν gives rise to an unbiased estimator ζ defined on (Y, µ) via ζ(y) = dν dµ (y). In practice, a suitable way to construct such an estimator is to use a Markov chain with an absorbing state. Wagner also discusses variance reduction techniques, specially importance sampling and a shift procedure consisting on adding a constant c to the integrand φ and then multiplying by the factor e -cT in order to get the right expectation. Wagner [START_REF] Wagner | Unbiased multi-step estimators for the Monte Carlo evaluation of certain functional integrals[END_REF] extends the class of unbiased estimators by perturbing the integrand φ by a suitably chosen function φ 0 and then using mixed integration formulas representation. Very recently, Beskos et al. [START_REF] Beskos | Exact and computationally efficient likelihood-based estimation for discretely observed diffusion processes[END_REF] obtained a simplified unbiased estimator for [START_REF] Dubois | Efficient pricing of asian options by the pde approach[END_REF], termed Poisson estimator, using Wagner's idea of expanding the exponential in a power series and his shift procedure. To be specific, the Poisson estimator writes

ψ(Z T )e cpT -cT N i=1 c -φ(Z Vi ) c P ( 6
)
where N is a Poisson random variable with parameter c P and (V i ) i is a sequence of independent random variables uniformly distributed on [0, T ]. Fearnhead et al. [START_REF] Paul Fearnhead | Particle filters for partially observed diffusions[END_REF] generalized this estimator allowing c and c P to depend on Z and N to be distributed according to any positive probability distribution on AE. They termed the new estimator the generalized Poisson estimator. We introduce a new degree of freedom by allowing the sequence (V i ) i to be distributed according to any positive density on [0, T ]. This gives rise to the following unbiased estimator for ( 5) :

Lemma 1 -Let p Z and q Z denote respectively a positive probability measure on AE and a positive probability density on [0, T ]. Let N be distributed according to p Z and (V i ) i∈AE * be a sequence of independent random variables identically distributed according to the density q Z , both independent from each other conditionally on the process (Z t ) t∈[0,T ] . Let c Z be a real number which may depend on Z. Assume that

|ψ(Z T )|e -cZ T exp T 0 |c Z -φ(Z t )|dt < ∞.
Then

ψ(Z T )e -cZ T 1 p Z (N ) N ! N i=1 c Z -φ(Z Vi ) q Z (V i ) (7)
is an unbiased estimator of C 0 .

Proof : The result follows from

ψ(Z T )e -cZ T 1 p Z (N ) N ! N i=1 c Z -φ(Z Vi ) q Z (V i ) (Z t ) t∈[0,T ] = ψ(Z T )e -cZ T +∞ n=0 T 0 c Z -φ(Z t )dt n p Z (n) n! p Z (n) = ψ(Z T ) exp - T 0 φ(Z t )dt .
2 Using (7), one is now able to compute the expectation at hand by a simple Monte Carlo simulation. The practical choice of p Z and q Z conditionally on Z is studied in the appendix 4.1.

As pointed out in Fearnhead et al. [START_REF] Paul Fearnhead | Particle filters for partially observed diffusions[END_REF], this method is an extension of the exact algorithm method since, under assumptions 3, 4 and 5, the reinforced integrability assumption of Lemma 1 is always satisfied. Indeed, suppose for example that lim sup

|ψ(Z T )|e -(MZ +k)T e T 0 |MZ +k-φ(Zt)|dt = |ψ(Z T )|e -(MZ +k)T e T 0 MZ +k-φ(Zt)dt = |ψ(Z T )|e -T 0 φ(Zt)dt < ∞
and hence, one is allowed to write that

C 0 = ψ(Z T )e -(MZ +k)T 1 p Z (N )N ! N i=1 M Z + k -φ(Z Vi ) q Z (V i ) .
Better still, the random variable

ψ(Z T )e -(MZ +k)T 1 pZ (N )N ! N i=1 MZ +k-φ(ZV i ) qZ (Vi)
is square integrable when p Z is the Poisson distribution with parameter M Z T + k and q Z is the uniform distribution on [0, T ] since we have then

  ψ(Z T )e -(MZ +k)T 1 p Z (N )N ! N i=1 M Z + k -φ(Z Vi ) q Z (V i ) 2   = ψ 2 (Z T ) N i=1 1 - φ(Z Vi ) M Z + k 2 ≤ ψ 2 (Z T ) < ∞.
The last inequality follows from the square integrability of f : whenever one is able to simulate from the density h, introduced in the exact algorithm, by doing rejection sampling, there exists a density ρ such that ψ, which is equal to f (Z T ) h(ZT ) ρ(ZT ) up to a constant factor, is dominated by f and so is square integrable. The square integrability property is very important in that we use a Monte Carlo method. We see that, whenever the exact algorithm is feasible, the unbiased estimator of lemma 1 is a simulable square integrable random variable, at least for the previous choice of p Z and q Z . Remark 2 -One can derive two estimators of C 0 from the result of Lemma 1 :

δ 1 = 1 n n i=1 f (Z i T ) e A(Z i T )-A(x)- (Z i T -x) 2 2T √ 2πρ(Z i T ) e -cZ T 1 p Z (N i ) N i ! N i j=1 c Z -φ(Z i V i j ) q Z (V i j ) δ 2 = n i=1 f (Z i T ) e A(Z i T )-A(x)- (Z i T -x) 2 2T √ 2πρ(Z i T ) 1 p Z (N i ) N i ! N i j=1 c Z -φ(Z i V i j ) q Z (V i j ) n i=1 e A(Z i T )-A(x)- (Z i T -x) 2 2T √ 2πρ(Z i T ) 1 p Z (N i ) N i ! N i j=1 c Z -φ(Z i V i j ) q Z (V i j )
.

Application : the pricing of continuous Asian options

In the Black & Scholes model, the stock price is the solution of the following SDE under the risk-neutral measure È

dS t S t = (r -δ)dt + σdW t (8)
where all the parameters are constant : r is the short interest rate, δ is the dividend rate and σ is the volatility. Throughout, we denote γ = rδ -σ 2 2 . The path-wise unique solution of ( 8) is

S t = S 0 exp(σW t + γt) .
We consider an option with pay-off of the form

f αS T + β T 0 S t dt (9)
where f is a given function such that

f 2 αS T + β T 0 S t dt < ∞,
T is the maturity of the option and α, β are two given non negative parameters3 . Note that for α = 0, this is the pay-off of a standard continuous Asian option.

The fundamental theorem of arbitrage-free pricing ensures that the price of the option under consideration is

C 0 = e -rT f αS T + β T 0 S u du .
At first sight, the problem seems to involve two variables : the stock price and the integral of the stock price with respect to time. Dealing with the PDE associated with Asian option pricing, Rogers and Rogers and Shi [START_REF] Rogers | The value of an Asian option[END_REF] used a suitable change of variables to reduce the spatial dimension of the problem to one. We are going to use a similar idea.

Let

ξ t = αS 0 + βS 0 t 0 e -σWu-γu du e σWt+γt .
We have that As a consequence C 0 = e -rT f (ξ T ) .

ξ t = αS 0 e σWt+γt + βS 0 t 0 e σ(Wt-Wu)+γ(t-u) du = αS 0 e σBt+γt +
By applying Itô's lemma, we verify that the process (ξ t ) t≥0 is a positive solution of the following 1-dimensional stochastic differential equation for which path-wise uniqueness holds

dξ t = βS 0 dt + ξ t (σdW t + (γ + σ 2 2 )dt) ξ 0 = αS 0 . (10)
We are thus able to value C 0 by Monte Carlo simulation without resorting to discretization schemes using one of the exact simulation techniques described in the previous section. In the case α = 0, one has to deal with the fact that ξ t starts from zero which is the reason why we distinguish two cases.

The case α = 0

We are going to apply both the exact algorithm of Beskos et al. [START_REF] Beskos | Retrospective exact simulation of diffusion sample paths[END_REF] and the method based on the unbiased estimator of lemma 1.

We make the following change of variables to have a diffusion coefficient equal to 1 :

X t = log(ξ t ) σ ⇒ dX t = ( γ σ + βS0 σ e -σXt )dt + dW t X 0 = x with x = log(αS0) σ . (11)

Thus

C 0 = e -rT f (e σXT ) .

The following proposition ensures that assumption 1 is satisfied.

Proposition 4 -The process (L t ) t∈[0,T ] defined by

L t = exp T 0 ( γ σ + βS 0 σ e -σYt ) dY t - 1 2 T 0 ( γ σ + βS 0 σ e -σYt ) 2 dt is a martingale under É W x .
Proof : Under É W x , (L t ) t∈[0,T ] is clearly a non-negative local martingale and hence a super-martingale.

Then, it is a true martingale if and only if

É W x (L T ) = 1.
Checking the classical Novikov's or Kamazaki's criteria is not straightforward. Instead, we are going to use the approach developed by Rydberg [START_REF] Rydberg | A note on the existence of unique equivalent martingale measures in a markovian setting[END_REF] (see also Wong and Heyde [START_REF] Wong | On the martingale property of stochastic exponentials[END_REF]) who takes advantage of the link between explosions of SDEs and the martingale property of stochastic exponentials.

Let us define the following stopping times :

τ n (Y ) = inf t ∈ Ê + such that t 0 γ σ + βS 0 σ e -σYu 2 du ≥ n ,
with the convention inf{∅} = +∞.

The stopped process (L t∧τn(Y ) ) t∈[0,T ] is a true martingale under É W x since Novikov's condition is fulfilled. According to the Girsanov theorem, one can define a new probability measure É n X , which is absolutely continuous with respect to É W x , by its Radon-Nikodym derivative

dÉ n X dÉ W x = L T ∧τn(Y ) . Hence É n X ½ {τn(Y )>T } = É W x ½ {τn(Y )>T } L T ∧τn(Y ) .
Since (τ n (Y )) n∈AE is a non decreasing sequence, we can pass to the limit in the right hand side We get

lim n→+∞ É n X (τ n (Y ) > T ) = É W x ½ {τ∞(Y )>T } L T ∧τ∞(Y )
where τ ∞ (Y ) denotes the limit of the non decreasing sequence (τ n (Y )) n∈AE .

Under É W x , (Y t ) t∈[0,T ] has the same law as a Brownian motion starting from x so τ ∞ (Y ) = +∞ , É W x almost surely, and consequently

É W x L T = lim n→+∞ É n X (τ n (Y ) > T ) .
On the other hand, the Girsanov theorem implies that, under É n X , (Y t ) t∈[0,T ∧τn(Y )] solves a SDE of the form [START_REF] Ingersoll | Theory of Financial Decision Making[END_REF]. To conclude the proof, it is sufficient to check that trajectorial uniqueness holds for this SDE. Indeed, the law of (Y t ) t∈[0,T ∧τn(Y )] under É n X is the same as the law of (

Y t ) t∈[0,T ∧τn(Y )] under É X . Hence É n X (τ n (Y ) > T ) = É X (τ n (Y ) > T ) -→ n→+∞ É X (τ ∞ (Y ) > T ) .
Clearly,

t 0 γ σ + βS0 σ e -σYu 2 du < ∞, É X almost surely, so É W x L T = É X (τ ∞ (Y ) > T ) = 1
as required.

In order to check trajectorial uniqueness for the SDE [START_REF] Ingersoll | Theory of Financial Decision Making[END_REF], we consider two solutions X 1 and X 2 . We have that

d(X 1 t -X 2 t ) = βS 0 σ e -σX 1 t -e -σX 2 t dt ⇒ d|X 1 t -X 2 t | = βS 0 σ sign(X 1 t -X 2 t ) e -σX 1 t -e -σX 2 t dt. So |X 1 t -X 2 t | = βS 0 σ t 0 sign(X 1 s -X 2 s ) e -σX 1 t -e -σX 2 t ds ≤ 0.
The last inequality follows from the fact that x → e -σx is a decreasing function. Finally, almost surely, ∀t ≥ 0, X 1 t = X 2 t which leads to strong uniqueness. 2

Consequently, thanks to the Girsanov theorem, we have

dÉ X dÉ W x = exp     T 0 ( γ σ + βS 0 σ e -σYt ) a(Yt) dY t - 1 2 T 0 ( γ σ + βS 0 σ e -σYt ) 2 dt     . (12)
Set

A(u) = u 0 a(x)dx = γ σ u + βS0 σ 2 (1 -e -σu ). Then dÉ X dÉ W x = exp A(Y T ) -A(x) - 1 2 T 0 a 2 (Y t ) + a ′ (Y t )dt . The function u → exp A(u) -(u-Y0) 2 2T = exp γ σ u + βS0 σ 2 (1 -e -σu ) -(u-Y0) 2

2T

is clearly integrable so we can define a new process (Z t ) t∈[0,T ] distributed according to the following law É Z

É Z = Ê L (W t ) t∈[0,T ] |W T = y h(y)dy
where the probability density h is of the form

h(u) = C exp A(u) - (u -Y 0 ) 2 2T with C a normalizing constant. ( 13
)
Remark 5 -Simulating from this probability distribution is not difficult (see the appendix 4.2 for an appropriate method of acceptance/rejection sampling).

We have

dÉ X dÉ Z = C exp - T 0 1 2 (a 2 (Y t ) + a ′ (Y t ))dt . Set φ(x) = a 2 (x)+a ′ (x) 2 = ( γ σ + βS 0 σ e -σx ) 2 -βS0e -σx 2 . A direct calculation gives inf x∈Ê φ(x) =    γ 2 2σ 2 if 2γ ≥ σ 2 φ 1 σ log( 2βS0 σ 2 -2γ ) otherwise.
Set k = inf x∈Ê φ(x). Finally, we get

dÉ X dÉ Z = Ce -kT exp - T 0 φ(Y t ) -k dt .
We check that

lim x→+∞ φ(x) = γ 2 2σ 2 < ∞ lim x→-∞ φ(x) = +∞.
Hence we can apply the algorithm 1 to simulate exactly X T and compute C 0 = e -rT f (e σXT ) by Monte Carlo. On the other hand, using [START_REF] Karatzas | Brownian motion and stochastic calculus[END_REF] we get

C 0 = e -rT f (e σW x T ) exp A(W x T ) -A(x) - 1 2 T 0 a 2 (W x t ) + a ′ (W x t )dt
and we can also use the unbiased estimator presented in the previous section to compute this expectation.

Remark 6 -We also applied the exact algorithm based on a geometric Brownian motion instead of the standard Brownian motion which seems more intuitive given the form of the SDE [START_REF] Geman | Bessel processes, asian option and perpetuities[END_REF]. The algorithm is feasible because we can simulate recursively a drifted Brownian motion and therefore a geometric Brownian motion by an exponential change of variables. The results we obtained were not different from the first method.

Numerical computation

For numerical tests, we consider the case

f (x) = (x -K) +
which corresponds to the European call option with strike K. Using the exact simulation algorithm presented above, we can simulate the underlying αS T + β T 0 S t dt at maturity (see Figure 1). Then, all we have to do is a simple Monte Carlo method to get the price of the option under consideration. Using the unbiased estimator, we get

C 0 =   e -rT (e σZT -K) + e A(ZT )-A(x)-(Z T -x) 2 2T √ 2πρ(Z T ) e -(MZ +k)T 1 p Z (N )N ! N i=1 M Z + k -φ(Z Vi ) q(V i )  
where (Z t ) t∈[0,T ] , ρ, M Z , k, p Z and q Z are defined as in section 1.2. In order to ensure square integrability, we choose p Z to be a Poisson distribution with parameter M Z T + k and q Z to be the uniform distribution on [0, T ]. For the density ρ, a good choice is to consider the density that we use to simulate from the distribution h by rejection sampling. We test these exact methods against a standard discretization scheme with the variance reduction technique of Kemna and Vorst [START_REF] Kemna | A pricing method for options based on average asset values[END_REF]. As pointed out by Lapeyre and Temam [START_REF] Lapeyre | Competitive Monte Carlo methods for pricing asian options[END_REF], the discretization of the integral by a simple Riemannian sum is not efficient. Instead, we use the trapezoidal discretization. In the sequel, we will denote this method by Trap+KV. The table 1 gives the results we obtained for the following arbitrary set of parameters : S 0 = 100, K = 100, r = 0.05, σ = 0.3, δ = 0, T = 1, α = 0.6 and β = 0.4. The computation has been made on a computer with a 2.8 Ghz Intel Penthium 4 processor. We intentionally choose a large number of simulations in order to show the influence of the number of time steps when using a discretization scheme. Empirical evidence shows that the exact simulation method is quite slow. This is mainly due to the fact that the rejection algorithm has a little acceptance rate (24% according to table 1). Using a geometric Brownian motion instead of a standard Brownian motion did not improve the results. Also, simulating recursively a Brownian path conditionally on its terminal value and its minimum is time consuming.

The unbiased estimator is more efficient, especially when we can avoid the recursive simulation of the Brownian path. To do so, we choose for p Z a Poisson distribution with mean c P T where c P is a free parameter. If we assume that the integrability condition in lemma 1 holds, then we can write that

C 0 =   e -rT (e σZT -K) + e A(ZT )-A(x)-(Z T -x) 2 2T √ 2πρ(Z T ) e cP T -cZ T N i=1 c Z -φ(Z Vi ) c P   .
Regarding the dependence of the exact simulation method with respect to the parameters α and β, it is intuitive that whenever α >> β, the method performs well since the logarithm of the underlying is not far from the logarithm of the geometric Brownian motion on which we do rejection-sampling. The table 2 confirms this intuition. We see that we cannot apply the algorithm for small values of α and then let α → 0 to treat the case α = 0.

α α + β 0.
Obviously, the two variables ξ T and 1 T T 0 S u du have the same law. Hence, the price of the Asian option becomes

C 0 = e -rT f 1 T T 0 S u du = e -rT f (ξ T ) .
Remark 7 -The pricing of floating strike Asian options is also straightforward using this method. It is even more natural to consider these options since it unveils the appropriate change of variables as we shall see below.

Let us consider a floating strike Asian call for example. We have to compute

C 0 = e -rT 1 T T 0 S u du -S T + .
Using S t = S t e δt as a numéraire (see the seminal paper of Geman et al. [START_REF] Geman | Changes of numéraires, changes of probability measure and option pricing[END_REF]), we immediately obtain that

C 0 = È S S 0 e -δT 1 T T 0 S u S T du -1 +
where È S is the probability measure associated to the numéraire S t . It is defined by its Radon-Nikodym

derivative dÈ S dÈ = e σWT -σ 2 2 T .
Under È S , the process B t = W t -σt is a Brownian motion and we can write that

C 0 = È S S 0 e -δT 1 T T 0 e σ(Bu-BT )+(r-δ+ σ 2 2 )(u-T ) du -1 + = S 0 e -δT 1 T T 0 e σ(Wu-WT )+(r-δ+ σ 2 2 )(u-T ) du -1 + = e -δT ξ T -S 0 +
where ξ t is the defined by ( 14) but with γ = rδ + σ 2 2 . We see therefore that the problem simplifies to the fixed strike Asian pricing problem.

Let us write down the stochastic differential equation that rules the process (ξ t ) t∈[0,T ] . Using Itô's lemma, we get

dξ t = ξ0-ξt t dt + ξ t σdW t + (γ + σ 2 2 )dt ξ 0 = S 0 .
Note that we are faced with a singularity problem near 0 because of the term ξ0-ξt t . We are going to reduce its effect using another change of variables.

Using Itô's lemma, we show that

C 0 = e -rT f S 0 e XT (15)
where X t = log(ξ t /ξ 0 ) solves the following SDE

dX t = σdW t + γdt + e -X t -1 t dt X 0 = 0. ( 16
)
Lemma 8 -Existence and strong uniqueness hold for the stochastic differential equation [START_REF] Lord | Partially exact and bounded approximations for arithmetic Asian options[END_REF].

Proof : Existence is obvious since we have a particular solution X t . The diffusion coefficient being constant and the drift coefficient being a decreasing function in the spatial variable, we have also strong uniqueness for the SDE (see the proof of Proposition 4).

2

Because of the singularity of the term e -X t -1 t in the drift coefficient, the law of (X t ) t≥0 is not absolutely continuous with respect to the law of (σW t ) t≥0 . That is why we now define (Z t ) t≥0 by the following SDE with an affine inhomogeneous drift coefficient :

   dZ t = σdW t + γdt - Z t t dt Z 0 = X 0 = 0. ( 17
)
The drift coefficient exhibits the same behavior as the one in ( 16) in the limit t → 0 in order to ensure the desired absolute continuity property. It is affine in the spatial variable so that (Z t ) t≥0 is a Gaussian process and as such is easy to simulate recursively.

Lemma 9 -The process

Z t = σ t t 0 sdW s + γ 2 t ( 18
)
is the unique solution of the stochastic differential equation [START_REF] Revuz | Continuous martingales and Brownian motion[END_REF].

Proof : Using Itô's Lemma, we easily check that Z t given by ( 18) is a solution of [START_REF] Revuz | Continuous martingales and Brownian motion[END_REF]. Again, constant diffusion coefficient and decreasing drift coefficient ensures strong uniqueness.

2

Remark 10 -For the computation of the price C 0 = e -rT (S 0 e XT -K) + of standard Asian call option, the random variable e -rT (S 0 e ZT -K) + provides a natural control variate. Indeed, since Z T is a Gaussian random variable with mean γ 2 T and variance σ 2 T 3 , one has e -rT (S 0 e ZT -K)

+ = S 0 e ( γ 2 + σ 2 6 -r)T N d + σ 1 3 T -Ke -rT N (d)
where N is the cumulative standard normal distribution function and

d = log(S0/K)+ γ 2 T σ √ 1 3 T .
Notice that in Kemna and Vorst [START_REF] Kemna | A pricing method for options based on average asset values[END_REF], the authors suggest the use of the control variate e -rT S 0 exp 1 T T 0 σW t + γt dt -K + which has the same law than e -rT S 0 e ZT -K + as 1 T T 0 σW t + γt dt is also a Gaussian variable with mean γ 2 T and variance σ 2 T 3 .

In order to define a new probability measure under which (Z t ) t≥0 solves the SDE ( 16), one introduces

L t = exp t 0 e -Zs -1 + Z s σs dW s - 1 2 t 0 e -Zs -1 + Z s σs 2 ds .
Because of the singularity of the coefficients in the neighborhood of s = 0, one has to check that the integrals in L t are well defined. This relies on the following lemma Lemma 11 -Let ǫ > 0. In a random neighborhood of s = 0, we have

|Z s | ≤ cs 1 2 -ǫ and |X s | ≤ cs 1 2 -ǫ
where c is a constant depending on σ,γ and ǫ.

Since ∀ǫ > 0, ∀z ≤ cs 1 2 -ǫ , e -z -1 + z σs 2 ≤ Cs -4ǫ ,
we can choose ǫ < 1 4 to deduce that L t is well defined.

Proof : We easily check that the Gaussian process (B t ) t∈[0,T ] defined by B t = (3t)

1 3 0 sdW s is a standard
Brownian motion. Thanks to the law of iterated logarithm for the Brownian motion (see for example Karatzas and Shreve [START_REF] Karatzas | Brownian motion and stochastic calculus[END_REF] p. 112), there exists t 1 (ω) such that 4 ,

∀t ≤ t 1 (ω), |B t (ω)| ≤ t 1 2 -ǫ 3 .
4 ω is an element of the underlying probability space Ω.

Therefore,

∀t ≤ (3t 1 (ω)) 1 3 , |Z t (ω)| = σ t B t 3 3 (ω) + γ 2 t ≤ σ 3 1 2 -ǫ 3 t 1 2 -ǫ + γ 2 t.
Taking c = max( σ

3 1 2 -ǫ 3 , γ 2 ) yields ∀t ≤ (3t 1 (ω)) 1 3 ∧ 1, |Z t (ω)| ≤ ct 1 2 -ǫ .
On the other hand, recall that X t = log(ξ /ξ 0 ) = log 1 t e σWt+γt t 0 e -σWu-γu du . So, using the law of iterated logarithm for the Brownian motion, we deduce that there exists t 2 (ω) such that

∀t ≤ t 2 (ω), 0 ≤ 1 t e σWt(ω)+γt t 0 e -σWu(ω)-γu du ≤ 1 t e σt 1 2 -ǫ +γt t 0 e σu 1 2
-ǫ -γu du.

Denote g(t) = 1 t e σt -ǫ -γu du and let us investigate the order in time near zero of this function. We have that e σt 1 2

-ǫ +γt = 1 + σt 1 2 -ǫ + O(t 1-2ǫ ) t 0 e σu 1 2 -ǫ -γu du = t + σ 3 2 -ǫ t 3 2 -ǫ + O(t 2-2ǫ ) hence g(t) = 1 + (σ + σ 3 2 -ǫ )t 1 2 -ǫ + O(t 1-2ǫ ), so X t (ω) ≤ log (g(t)) ∼ t→0 (σ + σ 3 2 -ǫ )t 1 2 -ǫ , which ends the proof for X t . 2
Proposition 12 -(L t ) t∈[0,T ] is a martingale and, consequently, for all g : C([0, T ]) → Ê measurable, the random variables g((X t ) 0≤t≤T ) and g((Z t ) 0≤t≤T )L T are simultaneously integrable and then g((X t ) 0≤t≤T ) = g((Z t ) 0≤t≤T )L T .

Proof : The proof is similar to the proof of Proposition 4.

We have already shown existence and strong uniqueness for both SDE ( 16) and [START_REF] Revuz | Continuous martingales and Brownian motion[END_REF]. Showing that the stopping time

τ n (Y ) = inf t ∈ Ê + such that t 0 e -Ys -1 + Y s σs 2 ds ≥ n , with the convention inf{∅} = +∞,
have infinite limits when n tends to +∞, É X and É Z almost surely, follows from the previous lemma.

2

One has

L T = exp T 0 e -Zt -1 + Z t σ 2 t dZ t - T 0 e -Zt -1 + Z t σ 2 t e -Zt -1 + Z t 2t + γ - Z t t dt . Set A(t, z) = 1 -z + z 2 2 -e -z σ 2 t
. The function A : ]0, T ] × Ê → Ê is continuously differentiable in time and twice continuously differentiable in space. So, we can apply Itô's Lemma on the interval [ǫ, T ] for ǫ > 0 :

A(T, Z T ) = A(ǫ, Z ǫ ) + T ǫ e -Zt -1 + Z t σ 2 t dZ t - ǫ 1 -Z t + Z 2 t 2 -e -Zt σ 2 t 2 dt + T ǫ 1 -e -Zt 2t dt
Using the lemma 9, we let ǫ → 0 to obtain

A(T, Z T ) = T 0 e -Zt -1 + Z t σ 2 t dZ t - T 0 1 -Z t + Z 2 t 2 -e -Zt σ 2 t 2 dt + T 0 1 -e -Zt 2t dt.
Then

L T = exp A(T, Z T ) - T 0 φ(t, Z t )dt
where φ is the mapping

φ(t, z) = e -z -1 + z -z 2 2 σ 2 t 2 + 1 -e -z 2t + e -z -1 + z σ 2 t e -z -1 + z 2t + γ - z t . (19)
By [START_REF] Levy | Pricing european average rate currency options[END_REF] and Proposition 12, we get

C 0 = e -rT f (S 0 e ZT ) exp A(T, Z T ) - T 0 φ(t, Z t )dt . (20)
Since for each t > 0, lim z→-∞ φ(t, z) = +∞ and lim z→+∞ φ(t, z) = -∞, it is not possible to apply the exact algorithm. One can use the unbiased estimator, at least theoretically, if there exists a random variable c Z measurable with respect to Z such that e A(T,ZT )-(r+cZ )T |f (S 0 e ZT )|e

T 0 |cZ -φ(t,Zt)|dt < ∞.
Unfortunately, this reinforced integrability condition is never satisfied :

Lemma 13 -Assume that f is a non identically zero function. Let p Z and q Z denote respectively a positive probability measure on AE and a positive probability density on [0, T ]. Let N be distributed according to p Z and (U i ) i∈AE * be a sequence of independent random variables identically distributed according to the density q Z , both independent conditionally on the process (Z t ) t∈[0,T ] . Then the random variable

e A(T,ZT )-rT f (S 0 e ZT ) 1 p Z (N ) N ! N i=1 -φ(U i , Z Ui ) q Z (U i ) (21)
is non integrable.

Proof : By conditioning on Z, one has

∆ := e A(T ,Z T )-rT |f (S0e Z T )| pZ (N ) N ! N i=1 |φ(Ui,ZU i )| qZ (Ui) = e A(T,ZT )-rT |f (S 0 e ZT )|e T 0 |φ(t,Zt)|dt ≥ e A(T,ZT )-rT |f (S 0 e ZT )|e T T 2 |φ(t,Zt)|dt
One can easily show that, ∀z < 0 and ∀t ∈ [ T 2 , T ], φ(t, z) ≥ φ(z) where

φ(z) = e -z -1 + z -z 2 2 σ 2 ( T 2 ) 2 + e -z -1 + z σ 2 T 2 e -z -1 + z T + γ + -2 z T Since φ(z) ∼ -∞ 2 e -2z
σ 2 T 2 , there exists c < 0 such for all z < c, φ(z) ≥ e -2z σ 2 T 2 . Hence, ∆ ≥ e A(T,ZT )-rT |f (S 0 e ZT )|e

1 σ 2 T 2 T T 2 e -2Z t ½ {Z t <c} dt ≥ e A(T,ZT )-rT |f (S 0 e ZT )|e -e -2c 2σ 2 T e 1 σ 2 T 2 T T 2 e -2Z t dt Using Jensen's inequality we get ∆ ≥ e A(T,ZT )-rT |f (S 0 e ZT )|e -e -2c 2σ 2 T exp 1 2σ 2 T e -4 T T T 2

Ztdt

We have seen in the proof of lemma 11 that

Z t = σ t B t 3 3 + γ 2 t where (B t ) t≥0 is a standard Brownian motion. So, conditionally on Z T , T T 2
Z t dt is a gaussian random variable and hence ∆ = +∞.

2

We are in a situation where e A(T,ZT )-rT |f (S 0 e ZT )|

1 pZ (N ) N ! N i=1 -φ(Ui,ZU i ) qZ (Ui) (Z t ) t∈[0,T ] is non inte- grable while e A(T,ZT )-rT |f (S 0 e ZT )| 1 pZ (N ) N ! N i=1 -φ(Ui,ZU i ) qZ (Ui) (Z t ) t∈[0,T ] is integrable since e -rT |f (S 0 e ZT )| exp A(T, Z T ) - T 0 φ(t, Z t )dt < ∞.
Then, a natural idea would consist in considering, for a given n ∈ AE * , the random variable

e A(T,ZT )-rT |f (S 0 e ZT )|   1 n n j=1 1 p Z (N j ) N j ! Nj i=1 -φ(U j i , Z U j i ) q Z (U j i ) (Z t ) t∈[0,T ]  
where (N j ) 1≤j≤n are independent variables having the same law as N and (U j i ) i∈AE * 1≤j≤n are independent sequences having the same law as (U i ) i∈AE * , both independent conditionally on the process (Z t ) t∈[0,T ] . The following general result tells us that this is not sufficient to circumvent integrability problems.

Lemma 14 -Let Y and Z be two real random variables and g : Ê → Ê a given measurable function. On the other hand, ∀z ∈ Ê,

Assume that g(Z) (Y |Z) is integrable while g(Z) (|Y | |Z) is non integrable. Then, when (Y i )
e n (z) ≥ 1 n n i=1 Y i ½ {∀2≤j≤n,Yj ≥0} |Z = z + n i=1 Y i ½ {∀2≤j≤n,Yj <0} |Z = z ≥ 1 n Y + 1 |Z = z È (Y 1 ≥ 0|Z = z) n-1 + Y - 1 |Z = z È (Y 1 < 0|Z = z) n-1 = 1 n e1(z)+e(z) 2 È (Y 1 ≥ 0|Z = z) n-1 + e1(z)-e(z) 2 (Y 1 < 0|Z = z) n-1 ≥ 1 n e1(z) 4 ½ {e1(z)≥2|e(z)|} È (Y 1 ≥ 0|Z = z) n-1 + e1(z) 4 ½ {e1(z)≥2|e(z)|} È (Y 1 < 0|Z = z) n-1 ≥ e1(z) n2 n ½ {e1(z)≥2|e(z)|}
Hence, g(Z)

1 n n i=1 Y i |Z = Ê |g(z)|e n (z)È Z (dz) = +∞. 2
There is still hope yet. In the proof of Lemma 13, we saw that integrability problems appear when Z t takes large negative values so that φ(t, Z t ) tends rapidly towards +∞. Since lim z→+∞ φ(t, z) = -∞, one possible issue is to split the function φ(t, Z t ) into a positive part and a negative part. The first term can be handled by the exact simulation technique whereas the second term, which as we shall see in the following section presents no integrability problems, can be handled by the unbiased estimator technique.

An hybrid pseudo-exact method

We rewrite [START_REF] Turnball | A quick algorithm for pricing european average options[END_REF] in the following form

C 0 = e A(T,ZT )-rT f (S 0 e ZT )e T 0 φ -(t,Zt)dt e -T 0 φ + (t,Zt)dt . (22)
Let p Z and q Z denote respectively a positive probability measure on AE and a positive probability density on [0, T ]. Let N be distributed according to p Z and (U i ) i∈AE * be a sequence of independent random variables identically distributed according to the density q Z , both independent conditionally on the process (Z t ) t∈[0,T ] . Note that, since e A(T,ZT )-rT f (S 0 e ZT )e T 0 |φ -(t,Zt)|dt e -T 0 φ + (t,Zt)dt = e A(T,ZT )-rT f (S 0 e ZT )e -T 0 φ(t,Zt)dt is integrable, one has

C 0 = e A(T,ZT )-rT f (S 0 e ZT ) 1 p Z (N )N ! N i=1 φ -(U i , Z Ui ) q Z (U i ) e -T 0 φ + (t,Zt)dt . (23)
Remark 15 -There is no hope that this estimator is square integrable. Indeed, one can show as in Lemma 13 that e T 0 (φ -(t,Zt)) 2 dt = +∞ since (φ -(t, z)) 2 is of order z 4 for large positive z.

The idea then is to apply the exact simulation technique to simulate an event with probability e -T 0 φ + (t,Zt)dt . Since for each t > 0, lim z→-∞ φ + (t, z) = +∞, one needs to bound from above φ + (t, z), uniformly with respect to t ∈ [0, T ], for z > c where c < 0 is a given constant. Thanks to the following lemma, it is possible to do so but only uniformly with respect to t ∈ [ǫ, T ] for all ǫ > 0 :

Lemma 16 -For all 0 < t ≤ T , sup z≥0 φ + (t, z) ≤ γ 2 σ 2 + γ σ 2 t + 1 t 1 2 - γ σ 2 + and ∀c < 0, sup z∈[c,0] φ + (t, z) ≤ e -c -1 + c σ 2 t 2 (1 + γ + t) + (e -c -1) 2 2σ 2 t 2 - c 2 σ 2 t 2 .
Proof : Let z > 0. It is useful to distinguish two cases according to the sign of γ :

1. γ ≥ 0

We rewrite φ in the following form

φ(t, z) = e -z -1 + z z 2 2 σ 2 t 2 + 1 -e -z t 1 2 - γ σ 2 + γz σ 2 t - z 2 -(z ∧ 1) 2 2σ 2 t 2 + (e -z -1) 2 -(z ∧ 1) 2 2σ 2 t 2
First note that

e -z -1+z-z 2 2 σ 2 t 2 ≤ 0, 1-e -z t 1 2 -γ σ 2 ≤ 1 t 1 2 -γ σ 2
+ and (e -z -1) 2 -(z∧1) 2

2σ 2 t 2 ≤ 0. Moreover, γz σ 2 t - z 2 -(z ∧ 1) 2 2σ 2 t 2 = 1 σ 2 γ z t - 1 2 z t 2 + ( z t ∧ 1 t ) 2 2 ≤ γ σ 2 t if γt ≤ 1 γ 2 σ 2 otherwise Consequently, φ + (t, z) ≤ γ 2 σ 2 + γ σ 2 t + 1 t 1 2 -γ σ 2 + .
2. γ ≤ 0

Now we rewrite φ in the following form

φ(t, z) = e -z -1 + z -z 2 2 σ 2 t 2 + γ e -z -1 + z σ 2 t + (e -z -1) 2 -z 2 2σ 2 t 2 + 1 -e -z 2t
It is then easy to show that φ + (t, z) ≤ 1 2t .

Note that

1 2t ≤ γ 2 σ 2 + γ σ 2 t + 1 t 1 2 -γ σ 2
+ . Hence, gathering the two cases yields the first part of the lemma.

Let now z ∈ [c, 0] for a given negative constant c. We rewrite φ in the following form

φ(t, z) = e -z -1 + z σ 2 t 2 (1 + γ + t) + (e -z -1) 2 2σ 2 t 2 - z 2 σ 2 t 2 ≥0 for z<0 + 1 -e -z 2t -γ -e -z -1 + z σ 2 t ≤0 for z<0 . Since ∂ z e -z -1 + z σ 2 t 2 (1 + γ + t) + (e -z -1) 2 2σ 2 t 2 - z 2 σ 2 t 2 = 1 -e -2z -2z + tγ + (1 -e -z ) t 2 σ 2 is negative for all z < 0, one has that sup z∈[c,0] φ + (t, z) ≤ e -c -1 + c σ 2 t 2 (1 + γ + t) + (e -c -1) 2 2σ 2 t 2 - c 2 σ 2 t 2 .

2

This lemma suggests to apply the exact algorithm on [ǫ, T ] for a fixed positive threshold ǫ. It remains to handle the time interval [0, ǫ[. Thanks to the following lemma, we that φ + (t, Z t ) can be approximately bounded from above for small t, almost surely, by a function of t. The idea is then to extend the exact simulation algorithm by simulating an inhomogeneous Poisson process. Of course, this hybrid method is no longer exact since the positive threshold for which the upper bound holds is random.

Lemma 17 -For all η > 0, there exists a random neighborhood of t = 0 such that

φ + (t, Z t ) ≤ 2c 3 3σ 2 + c 2 t -1 2 -η (24)
where c = max( σ

3 1 2 - η 3
, γ 2 ).

Proof : We rewrite [START_REF] Rydberg | A note on the existence of unique equivalent martingale measures in a markovian setting[END_REF] this way

φ(t, z) = 1 -e -z 2 + γ e -z -1 + z σ 2 1 t - -z + z 2 2 -e -z -1 2 (e -z -1 + z)(e -z -1 -z) σ 2 1 t 2
and make the following Taylor expansions

1 -z + z 2 2 -e -z -1 2 (e -z -1 + z)(e -z -1 -z) σ 2 = 2 3σ 2 z 3 + O(z 4 ) and 1 -e -z 2 + γ e -z -1 + z σ 2 = 1 2 z + O(z 2 ).
On the other hand, we have seen in the proof of lemma 11 that there exists a random neighborhood of zero such that Z t ≤ ct 1 2 -η where c = max( σ

3 1 2 - η 3 , γ 2
). We conclude that, in a random neighborhood of zero,

φ + (t, Z t ) ≤ 2c 3 3σ 2 + c 2 t -1 2 -η . 2

Numerical computation

For numerical computation, we are going to use the following set of parameters : S 0 = 100, K = 100, σ = 0.2, r = 0.1, δ = 0 and T = 1. To fix the ideas, let us consider a call option. The price C 0 writes as follows

C 0 = e A(T,ZT )-rT S 0 e ZT -K + e cp N i=1 φ -(U i , Z Ui ) c p e -T 0 φ + (t,Zt)dt .
where N ∼ P(c p ) and (U i ) i≥1 is an independent sequence of independent random variables uniformly distributed in [0, T ]. The parameter c p > 0 is set to one in the following. We give a description of the hybrid method we implement :

Algorithm 2

On the time interval In table 3, we give the price obtained by our method for different values of the positive threshold ǫ = T 2 J+1 . The number M of Monte Carlo simulations is equal to 10 5 and the true price is equal to 7.042 (computed using a Monte Carlo method with a trapezoidal scheme and a Kemna-Vorst control variate technique).

I j := [ T 2 j+1 , T 2 j ], 1. Simulate Z T 2 j+1 , Z T 2
Price CPU ǫ = T

21s

Table 3: Price of the Asian call using the hybrid-pseudo exact method.

Clearly, the method is not yet competitive regarding computation time. Nevertheless, unlike the usual discretization methods, it is not prone to discretization errors.

Conclusion

In this article, we have applied two original Monte Carlo methods for pricing Asian like options which have the following pay-off : (αS T + β T 0 S t dt -K) + . In the case α = 0, we applied both the algorithm of Beskos et al. [START_REF] Beskos | Retrospective exact simulation of diffusion sample paths[END_REF] and a method based on the unbiased estimator of Wagner [START_REF] Wagner | Unbiased Monte Carlo evaluation of certain functional integrals[END_REF] and more recently the Poisson estimator of Beskos et al. [START_REF] Beskos | Exact and computationally efficient likelihood-based estimation for discretely observed diffusion processes[END_REF] and the generalized Poisson estimator of Fearnhead et al. [START_REF] Paul Fearnhead | Particle filters for partially observed diffusions[END_REF]. The numerical results show that the latter performs the best. The more interesting case α = 0, which corresponds to usual continuously monitored Asian options, can not be treated using neither the exact algorithm, nor the method of exact computation of expectation but we investigate an hybrid pseudo-exact method which combines the two techniques. More generally, this hybrid method is an extension of the two exact methods and can be applied in other situations.

From a practical point of view, the main contribution of these techniques is to allow Monte Carlo pricing without resorting to discretization schemes. Hence, we are no longer prone to the discretization bias that we encounter in standard Monte Carlo methods for pricing Asian like options. Even though these exact methods are time consuming, they provide a good and reliable benchmark.

Appendix

4.1 The practical choice of p and q in the U.E method

The best choice for the probability law p of N and the common density q of the variables (V i ) i≥1 is obviously the one for which the variance of the simulation is minimum. In a very general setting, it is difficult to tackle this issue. In order to have a first idea, we are going to restrict ourselves to the computation of 1 p(N ) N ! N i=1 g(V i ) q(V i ) where g : [0, T ] → Ê. Proof : Minimizing the variance in ( 7) comes down to minimizing the expectation of the square of 1 p(N ) N ! N i=1 g(V i ) q(V i ) .

Set

F (p, q) = 1 (p(N ) N !) 2 N i=1 g 2 (V i ) q 2 (V i ) = +∞ n=0 T 0 g 2 (t) q(t) dt n p(n) (n!) 2 .
Using Cauchy-Schwartz inequality we obtain a lower bound for F (p, q) F (p, q) = We easily check that this lower bound is attained for q opt and p opt . 2

The optimal probability distribution p opt is the Poisson law with parameter 4.2 Simulation from the distribution h given by [START_REF] Kemna | A pricing method for options based on average asset values[END_REF] Recall that

h(u) = C exp A(u) - (u -X 0 ) 2 2T = C exp γ σ u + βS 0 σ 2 (1 -e -σu ) - (u -X 0 ) 2 2T
where C is a normalizing constant. The expansion of the exponential e -σu at the first order yields

h(u) ≈ C exp γ σ u + βS 0 σ u - (u -X 0 ) 2 2T = C exp - (u -(X 0 + T (γ+βS0) σ )) 2 2T
.

This suggests to do rejection sampling using the normal distribution with mean X 0 + T (γ+βS0) σ and variance T as prior. Unfortunately, for a standard set of parameters, this method gives bad results. Even a second order expansion of e -σu which also modifies the variance does not work.

In order to get round this problem, we evaluate the mode u * of h. We have

h ′ (u * ) = C γ σ + βS 0 σ e -σu * - u * -X 0 T exp γ σ u * + βS 0 σ 2 (1 -e -σu * ) - (u * -X 0 ) 2 2T .
So, h ′ (u * ) = 0 if and only if γ σ + βS 0 σ e -σu * -u * -X 0 T = 0 which writes σ(u * -X 0 -γ σ T )e σ(u * -X0-γ σ T ) = T βS 0 e -σX0-γT .

The function x → xe x is continuous and increasing on [0, +∞[ and so is its inverse which we denote by W . Since T βS 0 e -σX0-γT ≥ 0, we deduce that h is unimodal and that its mode satisfies u * = γT + W βS 0 T e -γT -σX0 + σX 0 σ .

The function W is the well-known Lambert function, also called the Omega function. It is uniquely valued on [0, +∞[ and there are robust and fast numerical methods based on series expansion for approximating this function (see for example Corless et al. [START_REF] Corless | On the lambert W function[END_REF]).

Numerical tests showed that performing rejection sampling using a Gaussian distribution with variance T and mean u * instead of X 0 + T (γ+βS0) σ gives plain satisfaction. In table 4, we see that for arbitrary choice of the parameter α α+β , the acceptance rate of the algorithm is always high (of order 70%) and that the computation time is low.

  u→+∞ φ(u) < +∞ and let k be a lower bound of φ, m Z be the minimum of the process Z and M Z an upper bound of {φ(u)k, u ≥ m Z }. Then, taking c Z = M Z + k in Lemma 1 ensures the integrability condition :

2. 2

 2 Standard Asian options : the case α = 0 and β > 0 A standard Asian option is a European option on the average of the stock price over a determined period until maturity. An Asian call, for example, has a pay-off of the form ( 1 T T 0 S u du -K) + . With our previous notations, it corresponds to the case α = 0, β = 1 T and f (x) = (x -K) + . The change of variables we used above is no longer suitable because it starts from zero when α = 0. Instead, we consider the following new definition of the process ξ Wt-Wu)+γ(t-u) du 0 = S 0 .

1 2 -ǫ +γt t 0 e σu 1 2

 122

  1≤i≤n is a sequence of independent random variables having the same law as Y , ∀n ∈ AE * , the random variableg(Z) | 1 n n i=1 Y i | |Z is non integrable.Proof : Denote by e, e 1 and e n three functions satisfying ∀z ∈ Ê, e(z) = (Y |Z = z) , e 1 (z) = (|Y 1 | |Z = z) and e n (zOn the one hand, since Ê |g(z)| |e(z)|È Z (dz) < ∞ and Ê |g(z)| e 1 (z)È Z (dz) = +∞ , where È Z is the law of Z, we have that Ê |g(z)| e 1 (z)½ {e1(z)≥2|e(z)|} È Z (dz) = +∞.

Lemma 18 -

 18 When g is a measurable function on [0, T ] such that 0 < T 0 |g(t)|dt < +∞, the variance of1 p(N ) N ! N i=1 g(V i ) q(V i )is minimal forq opt (t) = |g(t)| T 0 |g(t)|dt ½ [0,T ] (t) and p opt (n) =

  |dt. This justifies our use of a Poisson distribution for p.

Table 4 :

 4 Acceptance rate of the rejection algorithm of simulating from the distribution h in[START_REF] Kemna | A pricing method for options based on average asset values[END_REF] with S 0 = 100, σ = 0.3, T = 2 and r = 0.1.

Table 1 :

 1 Price of the option (9) using a standard discretization technique and exact simulation methods.

		Exact Simulation of the underlying : αS T + β
	0.018							S T	T
	0.016							αS T + β	0	S t dt
	0.014
	0.012
	0.010
	0.008
	0.006
	0.004
	0.002
	0.000
	10	53	97	140	183	227	270	313	357	400
	Method		M	N	Acceptance rate Price	C.I at 95%	CPU
				10				11.46 [11.43, 11.48]	5 s
	Trap+KV		20 10 6	-		11.46 [11.43, 11.49]	9 s
				50				11.47 [11.44, 11.5]	21 s
	Exact Simulation		-	10 6	24%		11.46 [11.43, 11.5]	81 s
	U.E (c P = M Z , c Z = M Z + k)	-	10 6	-		11.46 [11.43, 11.49] 17 s
	U.E (c P = c Z = 1/T )	-	10 6	-		11.46 [11.43, 11.49]	6 s

T 0 S t dt Figure 1: Histogram of 10 5 independent realizations of αS T + β T 0 S t dt for α = 0.6 and β = 0.4 compared with the lognormal distribution of S T .

Table 2 :

 2 Influence of the parameter α α+β on the acceptance rate of the exact algorithm.

	3	0.4	0.5	0.6	0.7
	Acceptance Rate 0.003% 0.47% 5.66% 24.43% 53.85%

  j and a lower bound m j for the minimum of (Z t ) t∈Ij (use the fact that Z t = Find M j > 0 such that ∀t ∈ I j , φ + (t, Z t ) ≤ M j (useLemma 16).3. Simulate an homogeneous spatial Poisson process on the rectangleI j × [0, M j ]and accept (respectively reject) the trajectory simulated if the number of points falling below the graph (φ + (t, Z t )) t∈Ij is equal to (respectively different from) zero. Carry on this acceptance rejection algorithm until reaching a time interval I J for a chosen J ∈ AE * . On the remaining time interval [0, T 2 J+1 ], use the same acceptance/rejection algorithm but with an inhomogeneous spatial Poisson process this time (use Lemma 17).

	σ t B t 3 3	+ γ 2 t where (B t ) t≥0 is a standard Brownian motion).
	2.

In their paper, the authors explain how to do such a decomposition of the Brownian path.

The underlying of this option is a weighted average of the stock price at maturity and the running average of the stock price until maturity with respective weights α and βT .


Keywords: DHT, RFID, Smart memory placement, EPCGlobal DHT, RFID, Placement mémoire efficace, EPCGlobal

This paper introduces SENSATION, a novel self-organizing stock management structure. SENSATION is based on a double DHT mechanism and is inspired from existing works such as SOLIST and Tribe. It proposes an efficient unique structure which can be used for different purposes such as data replication, distributed storing and request managements in stock management. These features help in the scalability and reliability of new storing warehouse management which tends to increase in scale and to be more and more interconnected. Results show that SENSATION provides interesting and promising results in terms of reliability and scalability.

Introduction

Nowadays, trade industries are growing up and need more and more the support of technology for managing their goods, for both logistic and traceability concerns. Ordinary applications such as inventories need the support of new technologies to fit the scalability of companies. More and more servers and databases are deployed to store information and connected in networks to become accessible from every place. These databases are now supplied through the reading of bar-codes or RFID tags hold by objects. Shelves equipped by readers even consist in database themselves and can be connected to database network. Even mobile RFID readers can be directly connected to the network and transfers useful data.

Although the networking architecture might be general, it is desirable that the communication infrastructure takes application requirements into account. Here, we narrow our focus to applications of stock management with a large amount of products to be controlled in different geographic sites communicating, and a need for limited communication overhead and reliable data storing. The network should be standard compliant1 to be connected to an EPC2 network (e.g. which respects the EPC standards).

In order to identify products to be controlled, they are equipped with identifiers (RFID, bar code, etc.). To perform the identification, traceability and inventory actions, the network is composed of database servers, routers, and two kinds of readers: either fixed (as a smart shelf) or mobile. Mobile readers can read and store data without being connected to the network, and so transfer data when connected. In this case, routers or fixed readers act as gateways towards databases.

The first motivation is the scalability and reliability of such a network. Indeed, the number of entities can increase heavily and sets the need of large scale communications without overloading the network. Data have to be duplicated for ensuring reliability but in a smart way in order to avoid memory overhead and facilitate at the same time the query routing. Secondly, the different characteristics of stock management requests, coupled with the need of scalability, require a specific network infrastructure management. Different granularities for the different diffusion primitives are needed, such as broadcast ("Count number of this specific product"), k-cast ("Is there a quantity x of this product") or anycast ("Give me the price of this product"). At last, due to the mobility of some readers, the structure of the network has to organize itself and adapt at every modification in a transparent and local way.

In this paper, we propose SENSATION a reliable k-level structure for distributed redundant memory placement. This structure allows to smartly store data and fast request to be processed by limiting bandwidth overhead. This unique structure allow performing different operations: (i) distributed storing to allow scalability by reducing memory overhead, (ii) data replication to ensure data reliability and (iii) efficient request management for answering requests in a smart scalable way which limit traffic and flooding overhead. This is achieved through the use of two distributed hash tables. The first one allows to reach the correct level, i.e. the level corresponding to the specific object targeted by the request while the second one allows to reach the proper correspondent in this layer. For instance, if the request is "What is the price of a blue trouser ?", the first DHT directs the request towards the layer responsible for trousers while the second one then directs the request towards a node aware of the price of blue trousers. Results show that for a low replication level (3), only 25% of data are lost when more than 70% of servers fail. The remaining of this paper is organized as follows. Section 2 describes the case studies we are focusing on. It is worth noting that SENSATION may find several applications in stock management but we describe a specific case study in order to illustrate more easyly the behavior of SENSATION. Section 3 sets background works useful for the understanding of SENSATION. Section 4 describes the relative works useful in the construction of SENSATION while Section 5 explains in details how we use them in SENSATION. Section 6 presents SENSATION performance results. At last Section 7 concludes that work by prospecting some future works.

Case study

Our case study is represented by a wide area stock management application. It can be, for instance, a distributor with several warehouses spread all over one or more countries; in each warehouse the servers and the databases are supposed to be connected, thus forming a network. Such an organization may need several tools at different levels.

In order to have a glimpse on these different tools, let us firstly consider an user located in one warehouse: he may need to draw an inventory (i) of the whole set of all warehouses, or (ii) of a particular warehouse located in a specific area, or evenly (iii) of products laying in a small area of a specific warehouse. Moreover, the inventory may target either (i) a kind of product (e.g. "inventory of every trouser"), or (ii) products in a specific place (e.g. "inventory of everything located in Paris"), or (iii) even both (e.g. "inventory of every trouser located in Paris"). On its side, the application has to handle different granularity requests. Finally, requests can be of different natures, possibly requiring different mechanisms for their diffusion: consider, for instance, (i) broadcast queries ("how many items of the specific product do I have?"), (ii) k-cast queries ("do I have at least k items of this product?"), and (iii) anycast queries ("give me the price of this product").

With a plain database network infrastructure, in all the abovementioned situations, every reader and/or database has to be queried (i.e. the network is flooded). Upon answers, some filters may be applied and/or data are aggregated. This normally introduces a non-negligible latency and causes a network overload, which actually is not necessary. Indeed, (i) for broadcast queries, only databases storing information about the articles being enumerated should be contacted; (ii) a k-cast query should be delivered as soon as k products have been found; finally, (iii) an unicast query may be answered after having contacted only one server. Essentially, these are the features we are to give SENSATION, in order to ensure scalability and high quality of service.

A second, equally important, concern in this case study is enabling a tunable fault-tolerance level. In order to add some reliability to data, it is necessary to include a degree of redundancy (i.e. provide one or more alternate servers in which to store the same data). The main issues here are to decide where, how and how many times to replicate data. Indeed, in most cases, data are simply replicated in one or more backup servers, which are eventually used in the case that the primary server is experimenting a failure. On the other hand, SENSATION replicates efficiently data and takes advantage of this replication by routing requests towards the closest server storing the information thus decreasing the routing overhead.

Preliminaries

Indirect routing

Providing a scalable and efficient location service in the context of self-organizing systems is a non-trivial problem, due to the spontaneity of networks. This requires a dynamic association between identification and location of a node, and the specification of a mechanism to manage this association. Furthermore, there is the need for minimizing the control message overhead for routing or location discovery. An efficient solution is to perform an indirect routing [START_REF] Hubaux | Towards self-organized mobile ad hoc networks: the terminodes project[END_REF][START_REF] Li | A scalable location service for geographic ad hoc routing[END_REF][START_REF] Xue | A scalable location management scheme in mobile ad-hoc networks[END_REF].

A routing operation is referred as indirect when it is performed in two steps: (i) first locate the target and then (ii) communicate with the target. This allows the network to decouple the location of a node from the location itself. With this approach, the information can be totally distributed, which is important for achieving scalability in large scale networks.

Distributed Hash Tables (DHT) represent the basis of indirect routing. Basically, they provide a general mapping between any information and a location establishing then a location-independent routing layer. They use a virtual addressing space V. Partitions of this virtual space are assigned to nodes in the network. As a basic example, consider Fig. 1 where the addressing space V = [0, 30[ is shared among 3 nodes i, j, k. In the figure, node i is assigned partition q(i) = [20, 30[, node j gets partition q(j) = [10, 20[ and node k is responsible for q(k) = [0, 10[. The idea is to use a hash function to first distribute node location information among rendezvous points. This same hash function is known by every node. Each information is hashed into a key (hash(v) = key v ∈ V) of this virtual addressing space V and then stored in the node(s) responsible for the partition of the virtual space. Yet, in Fig. 1(a), node k has a content C to register. By applying hash(C) = 25 it gets the key of C, which is within the partition of node i. So, node k registers C in node i. Usually, node k also stores the location of C (i.e. it registers its own address). Later, when a node j is looking for content C, it applies the same hash function, obtaining the same result (i.e. hash(C) = 25) thus knowing that it has to contact node i (cf. Fig. 1(b)). Upon the request, node i answers node j with the information previously registered (Fig. 1(c)). If this information is the location of C, then node j gets the address of node k and can then directly communicate with k.

DHT Mechanisms

Node Arrival: When a node u enters the network, it has to join the DHT overlay. To do so, u retrieves the way to contact a DHT node v that is used as an entry point to the DHT. Then, u is assigned a partition in the logical address space3 . Routing information in the system is then updated to reflect the presence of u. Finally, u retrieves all (key, value) pairs under its responsibility from the node that stored them previously.

INRIA (a) (b) (c)
Node Failure: When a node fails, the application data that it used to store is lost unless the DHT uses replication to keep multiple copies on different nodes. Some DHTs follow the simpler soft-state approach which does not guarantee persistence of data. Data items are pruned from the DHT unless the application refreshes them periodically. Therefore, a node failure leads to a temporary loss of application data until the data is refreshed.

Node Departure: DHT implementations often require departing nodes to notify the system before leaving. This allows other nodes to copy application data from the leaving node and to immediately update their routing information leading to improved routing efficiency.

Related works

Because of page restrictions, we only focus on literature works useful to SEN-SATION.

SOLIST

SOLIST [START_REF] Busnel | Solist: A lightweight multi-overlay structure for wireless sensor networks[END_REF] presents a structured overlay network and provides an efficient * -cast suite for wireless sensor networks. SOLIST identifies a set of basics functionalities widely used in distributed applications and propose an efficient implementation of this suite in a multi-layer structured network. The aim here is to give a type to each sensor and to group them depending on this type either static (heterogeneous sensors) or dynamic (battery state, sensor data, etc.). A layer is composed by all sensors of the same type. Fig. 2 shows a SOLIST multi-layer projection. To be included in the overlay of the proper type t it belongs to, or to perform a query about type t, a given sensor has to contact at least one node from the corresponding layer. To find out such a close node, SOLIST proceeds in two steps. The first step consists in contacting an entry point for this layer. The entry point will not necessary belong to the correct layer but is aware of the closest sensor that belongs to that layer. The second step thus consists in contacting that contact node. Proceeding in two steps allows always reaching the closest node thus limiting energy, bandwidth and time spendings.

To identify entry points and contact nodes, SOLIST uses hash functions.

Content Addressable Network (CAN )

CAN [START_REF] Ratnasamy | A scalable contentaddressable network[END_REF] introduces the notion of multi-dimensional identifier spaces by which routing efficiency is greatly improved compared to linear neighbor traversal in a single dimension. The average path length in a system with n nodes and d dimensions scales as O(d(n

1 d ))
. Partition assignment: In CAN, each data item is assigned an identifier of the form < x, y, z > for d = 3. Each node is said to own a zone. CAN ensures that the entire space is divided into non-overlapping zones.Because a key represents a point P in the identifier space, (key, value) pairs are stored on the node owning the zone which covers P . A new node n joining a CAN system sends a Join message to node d which is the current node responsible for the zone where n lays. Then, d splits its zone in half and assigns one half to n (cf. Fig. 3(b)). Finally, d transfers the keys to n for which it has become responsible. Routing over the CAN overlay: For routing purposes, a CAN node stores information only about its immediate neighbors. Two nodes in a d-dimensional space are considered neighbors if their coordinates overlap in one dimension and are adjacent to each other in d -1 dimension. For instance, in Fig. 3(a), nodes N 1 and N 6 are neighbors unlike nodes N 5 and N 6 which are not neighbors. If the local node does not own the zone of the destination, it forwards the message to its neighbor with the coordinates closest to the destination (Fig. 3(a)). In a d-dimensional space equally partitioned into n zones, this procedure results in an average of O((d/4)(n 1 d )) routing steps. This expresses that increasing the number of dimensions significantly reduces the average route length. Node Failure: The zones of failing or leaving nodes must be taken over by alive nodes to maintain a valid partitioning of the identifier space. A node detects the failure of a neighbor when it ceases to send update messages. Through an efficient message advertising, the neighboring node with the smallest zone volume merges the deserted zone with its own zone if possible. Alternatively, it temporarily manages both zones.

Node Departure: When a node l deliberately leaves a CAN system, it notifies a neighbor n whose zone can be merged with l ′ s zone. If no such neighbor exists, l chooses the neighbor with the smallest zone volume. It then copies the contents of its hash table to the selected node so this data remains available.

Tribe

Tribe [START_REF] Viana | Indirect routing using distributed location information[END_REF] has been designed for large scale self-organizing networks. Without any central control entity or positioning mechanism, Tribe creates a topology that is a logical network representation, and describes the relative location of the nodes according to their neighborhood in the physical network.

Partition assignment: When a node arrives in the network, it receives a control region which will serve for two purposes: node identification and routing.

A new node u in the network receives the control region from its neighbors which control region is the largest. This latter one gives half part of its region to the new node v and becomes its parent node. Tribe thus builds a tree as illustrated by links between nodes on Fig. 4.

Routing in the Tribe structure: The routing operation needs to be performed to reach a node responsible for a given key returned by the hash function, i.e. either when a node needs to register an information or to find it. Let p(u) = [p ⊖ u , p ⊕ u [ be the control region of node u. When node u entered the network, it has been assigned the virtual space p init (u) = [p ⊖ u , p ⊕ uinit [ where p ⊕ u ≤ p ⊕ uinit since from then, u may have shared its initial space with its children. When node u needs to reach the node responsible for a key key ∈ V, it proceeds as follows.

• If key ∈ p(u), u can answer the request.

• If key /

∈ p init (u), then forwards to its parent. • If p ⊕ u < key ≤ p ⊕ uinit , then forwards to its child v such that key ∈ p init (v). Every node reiterates this process till reaching the node responsible for the key.

RR n°7192 (a) (b) (c) (d)

SENSATION

Solution overview

SENSATION aims at proposing a self-organizing structure to manage data dynamism in a stock oriented network deployed over several geographic sites. It aims at answering stock management requests which may concern either a product, or a family of products or a geographical site.

To do so, in SENSATION, data are assigned a type like in SOLIST [START_REF] Busnel | Solist: A lightweight multi-overlay structure for wireless sensor networks[END_REF]. This type may be a product kind (e.g. type = trousers) or a location (e.g. type = Warehouse A) all addressed with no distinctions. Each type is associated to a layer. To address requests concerning a family of products or a larger geographical area, requests will be sent to a collection of layers. Indeed, if a request aims at drawing an inventory of sport articles and that sport articles family includes bikes and rackets, independent requests are sent to the racket layer and to the bike layer and aggregated afterwards. Similarly, if a request concerns every warehouse in Lille and that there are 2 warehouses in Lille, independent requests are sent to warehouse A layer and warehouse B layer.

Each layer t is composed of an overlay network composed of nodes and links. Nodes are databases and servers detaining information relative to type t, e.g every server or database storing data about trousers. When joining the overlay network of a layer t, a database A will be assigned a virtual space partition together with "neighbors", e.g. other database from the layer to which A is connected to. Note that since this is a logical network, this link may not be direct in the physical underlying network. SENSATION uses Tribe [START_REF] Viana | Indirect routing using distributed location information[END_REF] to assign the space partition to databases and to create links. Although Tribe may be directly applied within a layer, the setting up of the SOLIST-like architecture needs some adaptation. In SENSATION, this first basis space needs to be defined and databases need to detain coordinates which allow locating them in this space. We use CAN [START_REF] Balakrishnan | Looking up data in P2P systems[END_REF] in this purpose. The structure is then used for three purposes: (i) Distributed storing: every database which has fresh data of a given type t, firstly retrieves the corresponding responsible node for the type t. Through the responsible node it can enter INRIA the Tribe t-layer, in order to reach the database responsible for storing its data, by using Tribe routing; (ii) Replication: any level of replication can be imposed to SENSATION. Indeed, once the Tribe layer is reached, just applying different hash functions to the item ID to store would change the server responsible for the storing. (iii) Efficient request management: thanks to the use of several layers, a request will not be flooded in the whole network but only on nodes concerned by the request. Resources and latency are saved. The first step which consists in reaching the proper layer allows this node selection. Then, the use of Tribe within an overlay facilitates the request management. Indeed, the tree-like structure of Tribe is useful to aggregate counting requests in each level of the tree and stop flooding in enough nodes have been reached.

Mapping of the virtual coordinate system

The virtual structure has to include all nodes of all sites of the stock oriented network. Considering a national scale, sites represent the cities where a store is located. Thus, the lookup structure will be divided in NC cells, being each cell corresponding to a site.

In the following, each cell (or city) is assigned a coordinate (X C , Y C ) representing its relative position in the virtual system, along the x and y dimensions. More in details, the cell located in the bottom left corner of the grid has a coordinate equal to (0, 0), the next ones have coordinates equal to (1, 0) and (0, 1) along the x and y dimensions respectively, and so on. Finally, each cell is logically divided into T square areas, where T is maximum number of categories (which as mentioned could be types of goods or sites) supported by the application, and should be a perfect square number. Each square represents a virtual location in the cell, with relative coordinate (X T , Y T ) within a cell; each location acts as an entry point for a different type in the logical cell (square of coordinates (X T , Y T ) in cell i and square of coordinates (X T , Y T ) in cell j are entry points for the same layer(s)). Fig. 5 clarifies this concept; here the geographical area is divided into 3 × 3 virtual cells and each cell is divided into 16 sites: each site in the cell represents a different entry point for the 16 managed types: this means that the application can manage up to 16 types, being each type represented by an integer value within 0 and 15. For the sake of simplicity in the following we will consider virtual grids divided into the same number of cells along the x and y dimensions. This assumption does not limit the applicability of the approach, since one can build a grid with some unused cells.

By construction, we have exactly T entry-points within each cell, one for each type managed by the application. 4 This means that whenever an action has to be performed (insertion, lookup, . . . ) on an item of type p, a requesting user u has firstly to reach the p entry point within the cell he belongs to. By contacting the server database which is responsible of the requested type p (which actually is the server responsible of p entry point), u can reach the layer p (in the following we will say that "u is in the layer p"), so he can reach the server which has actually stored the item, by using Tribe. The critic point is how to share the virtual space represented by the entry-point set. We simply propose to share the entry-points set among the servers in a cell, by using CAN. For this reason, at a steady state the servers in a cell share the entry-points so as each of them is responsible of an non-overlapping partition of the DHT, being the DHT the entry-point set 5 . For this, we use a bi-dimensional congruential mapping function in order to determine the entry point of type p in the cell (X C , Y C ) (in the following it will be denoted as ep

(XC ,YC ) p
), which is:

ep (XC ,YC ) p = Xep (XC ,YC ) p , Y ep (XC ,YC ) p = mod(p, d) + X C • d, p d + Y C • d (1)
where d = √ T , p d is the integer division between p and d and mod(p, d) is the remainder after this division. Finally, once ep (XC ,YC ) p is computed, the request can be routed (using CAN) towards the requested entry point p. Once reached the right layer, Tribe is used for finding the node storing the item. Virtual coordinate assignment: To be inserted in the structure, databases needs coordinates which are assigned by the user who arbitrarily assigns to each new node (database) a unique geographic position in the grid. On the basis on the assigned site, a unique ID in the virtual system is generated as a bit sequence containing information about the virtual cell coordinates (X C , Y C ) and the relative position within the cell (X T , Y T ):

ID = |bin(XC , nbC)|bin(YC , nbC)|bin(XT , nbT )|bin(YT , nbT )|
where bin(x, n) is the binary integer representation of x over n bits, nbC = ⌈log 2 (N C)⌉ 2 , nbT = ⌈log 2 (T )⌉ 2 and | represents a string concatenation. To better understand this translation from a virtual site into a unique address, let us consider an example in which a user plugs a new database to the network in the geographical coordinate (7, 11) of a grid of N C = 9 cells. If T = 16, the resulting virtual ID is 8-bits long: 2 bits being used for codifying each field. From Fig. 5 it is clear that coordinate [START_REF] Xue | A scalable location management scheme in mobile ad-hoc networks[END_REF]11) belongs to cell (1, 2) and to square 3, 3 in that cell. Thus, ID = 01|10|11|11 (01, being the binary representation of XC , 10 the binary representation of YC and 11 the binary representation of XT and YT , on 2 bits each). Adding a server: When a server u wants to join the infrastructure it has to (i) obtain a portion of the entry-point set (i.e. joining CAN) and (ii) obtain a portion of the Tribe addressing space in all the layers obtained by CAN and corresponding to the objects it has to register (i.e. join Tribe). To better clarify this, let us consider a server u wanting to join the infrastructure. The nearest CAN node v already in the structure has to give u a portion of the entry-points it currently manages. Let consider, for instance that u receives an entry-points set P = 1, 2, 7, 8. This means that u will be the entry point for any request of access layers 1, 2, 7, 8 within the cell, so that u has also to join the same layers in Tribe. In other words, it has to obtain a portion of the Tribe addressing space in each of the mentioned layers. Adding items: When an item of type i has to be stored in the infrastructure, the server which is currently managing the entry-point i is contacted by using CAN (note that this operation is purely local, in the sense that the responsible node of type i is local to the cell). Since this server is also registered in the Tribe i-layer, the request can be forwarded using Tribe. Notice that, when the right Tribe layer is reached, any level of replication of the data can be imposed, in order to increase the faulttolerance. Just applying different hash functions to same item ID, in fact, would result in a different server managing the hashed ID (and so the storing server). The only constrain to impose to have a level of replication of h is that the hash functions set h has to be known to all servers. With this technique a tunable fault tolerance level can be introduced. Moreover, this technique let the application further diminish the overhead cost introduced by Tribe, when querying. This aspect will be clearer in the next paragraph. Lookup phase: When a query has to be answered, firstly the responsible node of the requested type is contacted, by using CAN. This phase is common to the adding phase. Once in the right layer, in dependence on the type of query a different scheme is used. In particular, for unicast, the request is routed using Tribe until the server managing the requested hashed ID is reached. This server is the responsible of directly providing the answer to the asking user. When a replication level h is imposed, we know that actually the item (if it was added) was stored h times in the layer. For this reason, it is less expensive (from a latency and communication overhead point of view) computing the hash function h times and reach the tribe node at the minimum distance from the contact node. This nice feature of SENSATION is analyzed in the simulation experiments.

Example

Let us assume that a database B wishes to draw an inventory of trousers. To respond to it, by running SENSATION, B first contacts the closest entry point of the type 'trousers'. To identify this latter one, B applies Eq. 1 and CAN (which returns the relative location of the proper entry point corresponding to trousers layer in every cell). It determines the closest one ept. Then, B contacts ept and gets the location of the nearest node C managing trousers, called here contact node (Fig. 6). Note that ept knows node C since this latter one has registered at node ept by also using Eq. 1. The counting request is then sent to C. Since C has previously joined the layer corresponding to trousers, it can forward the message to all nodes managing trousers that have previously registered in this layer. The counting request will be sent in broadcast mode in the overlay and the answer will contain the number of trousers in the entire network. Another request is to know if a product is existent (here, trousers). As previously, node B contacts the closest entry point for the type 'trousers'. If the entry point returns a node id, the existence is proved. Otherwise, the entry point returns Not Found and there are no trousers.

If a user wants to know if there are at least 10 trousers to ensure a command, he still contacts the closest entry point, which in turn sends the nearest contact node C A query may also concern a geographical location. In such a case, the same mechanism is applied on the type 'Lille' if the query wishes to set an inventory of Lille warehouse, for instance.

Simulation

An exhaustive performance analysis is performed to analyze the routing cost of requests, the robustness at various replication levels and the distribution of the information in SENSATION. Testbed: In each cell 10 servers are uniformly randomly placed (so a total of 640 servers) is a square universe divided into 64 cells of equal size over the Matlab6 software. After having inserted the servers, items are added in the following way: 100 items of type i are inserted from randomly extracted servers, where i ∈ [START_REF] Balakrishnan | Looking up data in P2P systems[END_REF]400]. Each experiment was performed 30 times and results were averaged. Overhead: After the registration phase, each registered item is queried by a server. This action (which is actually what we have denoted as unicast query) in a classical database network should not have any cost, since the information is stored in the local server. On the other side, using SENSATION, the entry point of the requested type has to be contacted (routing with CAN); once the responsible node is reached, Tribe is used to locate the server which is actually storing the requested item. Fig. 8 shows the overhead introduced by CAN for reaching the requested entry points. For the sake of readability a histogram is shown: it can be easily seen that more than 35% of the entry points can be contacted with a routing cost of 2 hops, which is a very low overhead. Note that this overhead is independent on the replication level introduced in Tribe layers. Dealing with node failures: To evaluate SENSATION robustness, a percentage of failing servers varying within 0 and 100% is considered. Note that, when a server is experimenting a failure, the items registered there are unavailable (it is a database failure), while the routing in the layers is still possible. It means that during the failure the information is not available, but CAN and Tribe can react to the failure and auto-organize again. Fig. 7(a) shows, for different replication levels and different loss regimes the percentage of information lost (this is equivalent of showing the failure rate of unicast queries). It can be noticed that, with a level of replication equal to 1 the loss of information is linear. Indeed, the information is uniformly distributed as the number of failures. Moreover, as expected, increasing the level of replication in the Tribe layers increases the reliability of SENSATION. With a low replication level equal to 2, when 50% of servers are down, only 25% of data are lost. This same amount of losses is even achieved when 70% of servers are down with a replication level equal to 3. In addition, this is worth noting that this reflects the worst scenario since in these simulations, failures appear randomly while in realistic scenario, when a server fails, its closest servers are more likely to fail rather than other ones in the network. Since data are spread randomly, in realistic scenario, data losses will indeed be less numerous. Request performance: Fig. 7(b) plots the percentage of successful unicast queries, at different replication levels and loss regimes. Results show that, unicast queries are pretty well managed by SENSATION with a sufficiently low order of replication and even at relatively high loss regime. This has a simple explication: let us assume a replication level of k. Information about an item is lost if and only if the k servers in which the item was stored are simultaneously experimenting a failure; otherwise, after the re-organization the item is registered again. Lost are only temporary. Unfortunately, this behavior is not to be expected in broadcast queries. Indeed, a broadcast query, in fact, can be correctly served (i.e. produce the right answer) if and only if all the items of the requested type are present in any of the k servers in which they were registered. As a consequence, just the disappearing of an item of the given type (the item on which the unicast query is failing) produces a failure of the broadcast query. Nevertheless, as shown by Fig. 7(c), SENSATION well resists to these failures.

Conclusion

We have proposed SENSATION, a novel self-organizing stock management structure. SENSATION is based on existing works proposed for wireless sensor networks and INRIA adapts them to another context. It proposes an efficient unique structure which can be used for different purposes such as data replication, distributed storing and request managements. These features help in the scalability and reliability of new storing warehouse management which tends to increase in scale and to be more and more interconnected. Results show that SENSATION provides interesting and promising results in terms of reliability and scalability. Next steps would be to study in more details the cost of replication towards the gain in reliability. In addition, some real implementations and tests should be performed to complete the full analysis.

Figure 1 :

 1 Figure 1: (a) Node k stores content C, and registers information about C on node i which is its rendezvous node. (b) Lookup phase of node j to contact the rendezvous node of C. (c) Lookup answer with information about C.

Figure 2 :

 2 Figure 2: Projection of 3 groups of nodes ( , and •) into three layers.

Figure 3 :

 3 Figure 3: CAN mechanisms illustration.

Figure 4 :

 4 Figure 4: Tribe virtual space sharing. (a) node A is alone and responsible for the whole virtual space V = [0, 30[. (b) Node B pops up (b) and gets half of A partition. (c) Node C appears, A gives again half of its partition. (d) Node D arrives, B shares its partition since it has a larger partition than node A.

Figure 5 :

 5 Figure 5: Adding a node. Cell coordinates are reported in the bottom and left sides of the x and y dimensions, resp.

  node C, which is the nearest node in the ligh-t-layer Node B computes the coordinates of the nearest entry point (ep ) which returns the nearest node (C) identifier and coordinates belonging to the layer

Figure 6 :

 6 Figure 6: Node B retrieving the nearest node of Ligh-trousers-Layer.

  (a) Percentage of lost information. (b) successful unicast queries. (c) Percentage of failing broadcast queries.

Figure 7 :

 7 Figure 7: Performance of SENSATION varying the replication level in Tribe and the percentage of failing servers.

Figure 8 :

 8 Figure 8: Overhead in CAN.

RR n°7192

RFID relative standards are established by EPC Global.

http://www.epcglobalinc.org/ INRIA

Depending on the DHT implementation, u may choose arbitrarily partitions on its own or it determines one based on the current state of the system.RR n°7192

Note that there may be less than T entry points. If a position is empty, the closest server handles the role of the entry point for this empty position.RR n°7192

This is one of the main differences between SOLIST and SENSATION.INRIA

http://www.mathworks.fr/ INRIA

This work has been partially supported by the FP7 European project ASPIRE (http://fp7-aspire.eu) and the regional project ICOM.