Christophe Profeta
email: christophe.profeta@univ-evry.fr

Abass Sagna
email: abass.sagna@ensiie.fr.

Conditional hitting time estimation in a nonlinear filtering model by the Brownian bridge method

The model consists of a signal process X which is a general Brownian diffusion process and an observation process Y , also a diffusion process, which is supposed to be correlated to the signal process. We suppose that the process Y is observed from time 0 to s > 0 at discrete times and aim to estimate, conditionally on these observations, the probability that the non-observed process X crosses a fixed barrier after a given time t > s. We formulate this problem as a usual nonlinear filtering problem and use optimal quantization and Monte Carlo simulations techniques to estimate the involved quantities.

Introduction

We consider in this work a nonlinear filtering model where the signal process X and the observation process Y evolve following the stochastic differential equations: dX t = b(X t , t)dt + σ(X t , t)dW t , X 0 = x 0 , dY t = h(Y t , X t , t)dt + ν(Y t , t)dW t + δ(Y t , t)d W t , Y 0 = y 0 .

(

In these equations, W and W are two independent standard real valued Brownian motions. We suppose that the functions b, σ, h, ν, and δ are Lipschitz and that, for every (x, t) ∈]0, +∞) 2 , δ(x, t) > 0, ν(x, t) > 0 and σ(x, t) > 0.

Let a be a real number such that 0 < a < x 0 and let τ X a = inf{u ≥ 0, X u ≤ a} be the first hitting time of the barrier a by the signal process. As usually, we consider that inf ∅ = +∞. Our aim is to estimate the distribution of the conditional hitting time

P τ X a > t| F Y s ( 2
)
for t ≥ s > 0 and where (F Y t ) t≥0 is the filtration generated by the observation process Y :

F Y s = σ(Y u , u ≤ s).
More generally, we shall denote by (F Z t ) t≥0 the filtration generated by the process Z. Such a problem arises for example in credit risk when modeling a credit event in a structural model as the first hitting time of a barrier a by the firm value process X. Investors are supposed to have no access to the true value of the firm but only to the observation process Y , which is correlated to the value of the firm (see e.g. [START_REF] Coculescu | Valuation of default sensitive claims under imperfect information[END_REF][START_REF] Duffie | Term structures of credit spreads with incomplete accounting information[END_REF]). We will typically suppose that we observe the process Y at regular discrete times t 0 = 0 < t 1 < . . . < t m = s over the time interval [0, s] and intend for estimating the quantity P τ X a > t| Y t 0 , . . . , Y tm for every t ≥ s. Note that, if t < s, then, applying the Markov property to the diffusion Y , the computations boil down to the case s = t. In [START_REF] Callegaro | An application to credit risk of optimal quantization methods for nonlinear filtering[END_REF], the quantity

P inf u∈[s,t]
X u > a|Y t 0 , . . . , Y tm has been estimated by a hybrid Monte Carlo-Optimal quantization method in the case where the observation process dynamics is given by:

dY t = Y t h(X t )dt + ν(t)dW t + δ(t)d W t , Y 0 = y 0 ,
where ν and δ are deterministic functions. However, the approach used in the previous work does not apply to our framework because we want to compute the conditional distribution of a function of the whole trajectory of the signal process from 0 to t given the observations from 0 to s, with s ≤ t.

Example 1. A particular case of model ( 1) one may consider is the following "Black-Scholes" case:

dX t = X t (µdt + σdW t ), X 0 = x 0 , dY t = Y t (rdt + νdW t + δd W t ), Y 0 = y 0 , (3)
so that

dY t Y t = ν σ dX t X t + r -µ ν σ dt + δd W t , (4)
or

Y t = y 0 (x 0 ) ν/σ X ν/σ t exp δ W t + r - ν 2 + δ 2 2 - µν σ + νσ 2 t . (5)
Observe that setting r = µ and σ = ν yields

dY t Y t = dX t X t + δd W t ,
meaning that the return on Y is the return on X affected by a noise (see e.g. [START_REF] Coculescu | Valuation of default sensitive claims under imperfect information[END_REF]). Of course, in this case, we may compute theoretically the expression (2) by noticing that:

τ X a = inf{u ≥ 0, X u ≤ a} = inf u ≥ 0, δ W u - δ 2 2 u ≥ ln(Y u ) + ln x 0 y 0 -ln(a)
and that, conditionally to F Y s , the process W u , u ≤ s has the same law as

σ √ σ 2 + δ 2 B u + δ σ 2 + δ 2 ln(Y u ) + ln x 0 y 0 + σ 2 + δ 2 2 u , u ≤ s
where B is a standard Brownian motion independent from the process Y . This follows from the fact that, since W and W are two independent Brownian motion, so are

B = σW + δ W √ σ 2 + δ 2 and B = δW -σ W √ σ 2 + δ 2 ,
and from the relations :

B t = 1 √ σ 2 + δ 2 ln(Y t ) + ln x 0 y 0 + σ 2 + δ 2 2 t and W t = δB t -σ B t √ σ 2 + δ 2 .
Therefore, in this particular setting, the problem boils down to the computation of the first passage time of a Brownian motion to a curved boundary. We refer to [START_REF] Coculescu | Valuation of default sensitive claims under imperfect information[END_REF] for some similar, and far more general, considerations.

Example 2. Another particularly simple case is given when both the signal and the observation process evolve following the Ornstein-Uhlenbeck dynamics:

dX t = λ(θ -X t )dt + σdW t , X 0 = x 0 , dY t = λ(θ -Y t )dt + σdW t + δd W t , Y 0 = y 0 . (6)
In this case, we have :

X t = Y t -(y 0 -x 0 )e -λt -δe -λt t 0 e λu d W u , so that τ X a = inf u ≥ 0, Y u -a -(y 0 -x 0 )e -λu ≤ δe -λu u 0 e λv d W v
and, conditionally to F Y s , the process e -λu u 0 e λv d W v , u ≤ s has the same law as

σ √ σ 2 + δ 2 U u + δe λu σ 2 + δ 2 Y u -y 0 e -λu -θ(1 -e -λu ) , u ≤ s
where U is an Ornstein-Uhlenbeck process with parameters λ and 0 started from 0, and independent from Y . This follows from Knight's representation theorem combined with the same ideas as above.

The rest of the paper is organized as follows: in Section 2, we state and prove the main theorems, i.e. we give an approximation of the expectation (2) when X is replaced by its continuous Euler scheme X, see especially Subsection 2.3. Then, in Section 3, we introduce some numerical tools in order to compute the quantities involved. We finally conclude the paper by a few simulations.

2 Estimation of the conditional survival probability

Preliminaries results

To deal with the computation of the conditional hitting time (2), we introduce the first hitting time from time s:

τ (s) = inf{u ≥ s, X u ≤ a}.
Define furthermore the filtration (G t ) t≥0 by, for every t ≥ 0,

G t = σ{(W u , W u ), u ≤ t}. (7)
We have the following result.

Lemma 2.1. For every t ≥ s,

P τ X a > t| F Y s = E 1 1 {τ X a >s} F (s, t, X s )| F Y s (8) with F (s, t, X s ) = P inf u∈[s,t] X u > a|X s . Proof. If t ≥ s, one has {τ X a > t} = {τ X a > s} ∩ {τ (s) > t}, so that P τ X a > t| F Y s = P {τ X a > s} ∩ {τ (s) > t}| F Y s = E E 1 1 {τ X a >s} 1 1 {τ (s) >t} |G s |F Y s ,
where the σ-algebra G s is defined in [START_REF] Gobet | Schémas dÕEuler pour diffusion tuée. Application aux options barrière[END_REF]. Hence, using the fact that W is independent from W and the Markov property of the process X:

P τ X a > t| F Y s = E 1 1 {τ X a >s} E 1 1 {τ (s) >t} |G s |F Y s = E 1 1 {τ X a >s} E 1 1 {τ (s) >t} |F W s |F Y s = E 1 1 {τ X a >s} P τ (s) > t|F X s |F Y s = E 1 1 {τ X a >s} P inf u∈[s,t] X u > a|X s |F Y s .
This shows assertion [START_REF] Graf | Foundations of Quantization for Probability Distributions[END_REF].

Note that in general, there is no closed-form expression for computing P inf

u∈[s,t]
X u > a|X s , except in some very special cases, such as Brownian motion or Bessel processes... We shall therefore need an approximation of this expression, which is the purpose of the next subsection.

The continuous Euler scheme

Let us denote by X the continuous Euler scheme process associated to the signal process X and defined by

Xs = Xs + b( Xs , s)(s -s) + σ( Xs , s)(W s -W s ), X0 = x 0 , with s = t k if s ∈ [t k , t k+1 ).
To estimate the distribution of the conditional hitting time given in [START_REF] Callegaro | An application to credit risk of optimal quantization methods for nonlinear filtering[END_REF], consider that we observe the process Y at discrete and regularly spaced times: t 0 , t 1 , . . . , t m , with

0 = t 0 < • • • < t m = s < t m+1 < • • • < t n = t.
In our general model (see Equation ( 1)), the discrete time observation processes { Xt k , k = 0, . . . , n} and {Y t k , k = 0, . . . , m} are obtained from Euler scheme as:

Xt k+1 = Xt k + b( Xt k , t k )∆ k + σ( Xt k , t k )(W t k+1 -W t k ) Y t k+1 = Y t k + h(Y t k , Xt k , t k )∆ k + ν(Y t k , t k )(W t k+1 -W t k ) + δ(Y t k , t k )( W t k+1 -W t k ) (9)
where k = 0, . . . , m -1 for the observation process, k going till n -1 for the signal process and ∆ k = t k+1t k . Note that if t > s, the number of discretization steps over [0, s] may differ from the number of discretization steps over [s, t] so that we choose

   t k = ks m for k = 0, . . . , m t k = s + k(t-s) n for k = m + 1, . . . , n.
Supposing that we have observed the trajectory (Y t 0 , . . . , Y tm ) of the observation process Y we aim to estimate

P(τ X a > t n | Y t 0 , . . . , Y tm ) by P τ X a > t n | Y t 0 , . . . , Y tm , (10)
where τ X a = inf{u ≥ 0, Xu ≤ a}. To this end, we will use a useful result called the regular Brownian bridge method. This result recalled below allows us to compute the distribution of the minimum (or the maximum) of the continuous Euler scheme X of the process X over the time interval [0, t], given its values at discrete and regular time observation points 0 = t 0 < t 1 < • • • < t n = t (see, e.g. [START_REF] Glasserman | Asymptotically optimal importance sampling and stratification for pricing path-dependent options[END_REF][START_REF] Gobet | Schémas dÕEuler pour diffusion tuée. Application aux options barrière[END_REF]).

Lemma 2.2. Let (X t ) t≥0 be a diffusion process with dynamics given by

X t = x + t 0 b(X u , u)du + t 0 σ(X u , u)dW u
and let ( Xt ) t≥0 be its associated continuous Euler process. Then, the following equality in law holds:

L min u∈[0,t] Xu | Xt k = x k , k = 0, • • • , n = L min k=0,••• ,n-1 H x k ,x k+1 ∆ k σ 2 (x k ,t k ) -1 (Λ k ) , (11)
where (Λ k ) k=0,••• ,n-1 are i.i.d. random variables uniformly distributed over the unit interval and

H x k ,x k+1 ∆ k σ 2 (x k ,t k ) -1
is the inverse function of the conditional cumulative function H

x k ,x k+1 ∆ k σ 2 (x k ,t k ) , defined by H x k ,x k+1 ∆ k σ 2 (x k ,t k ) (u) := exp - 2 ∆ k σ 2 (x k ,t k ) (u -x k )(u -x k+1 ) if u ≤ min(x k , x k+1 ) 1 otherwise. ( 12
)
In the following, we shall replace the expression σ(x k , t k ) in the expression of H by σ k and rather write:

H x k ,x k+1 ∆ k σ 2 k . Short proof. Observe first that, conditionally to { Xt k = x k ; k = 0, • • • , n}, the random variables min u∈[t k ,t k+1 ]
Xu ; k = 0, . . . , n -1 are mutually independent, thanks to the independent increments property of Brownian motion. Then, it suffices to notice that computing the law of min

u∈[t k ,t k+1 ] Xu conditionally to { Xt k = x k , Xt k+1 = x k+1 }
amounts to computing the law of the hitting times of the bridge of a Brownian motion with drift. This law is well-known to be independent from the drift (i.e. from the function b here) and to be given by an expression such as [START_REF] Pagès | Optimal quadratic quantization for numerics: the gaussian case[END_REF].

In the rest of the paper, the function 1 -H

x k ,x k+1 ∆ k σ 2 k
(a) will be often used. We shall denote it by

G x k ,x k+1 ∆ k σ 2 k (a), so that G x k ,x k+1 ∆ k σ 2 k (a) = 1 -exp - 2(x k -a)(x k+1 -a) ∆ k σ 2 (x k , t k ) 1 1 {x k ≥a; x k+1 ≥a} . (13)
Now, in our set-up, we shall apply the Brownian bridge method to obtain the following lemma, which is taken from [START_REF] Callegaro | An application to credit risk of optimal quantization methods for nonlinear filtering[END_REF]:

Lemma 2.3. We have

P inf u∈[tm,tn] Xu > a ( Xt k ) k=0,...,m = n-1 k=m G Xt k , Xt k+1 ∆ k σ 2 k (a). (14)
Proof. It follows from Lemma 2.2 that

P inf u∈[tm,tn] Xu > a ( Xt k ) k=0,...,m = P min k=m,••• ,n-1 H Xt k , Xt k+1 ∆ k σ 2 k -1 (Λ k ) > a .
Since the function H

z k ,z k+1 ∆ k σ 2 k (•)
Xu > a ( Xt k ) k=0,...,m = n-1 k=m P Λ k ≥ H Xt k , Xt k+1 ∆ k σ 2 k (a) = n-1 k=m 1 -H Xt k , Xt k+1 ∆ k σ 2 k (a) .
This completes the proof.

Main theorem

We now state and prove the main theorems of this section. We first show that the conditional hitting time of the continuous Euler process X given the discrete path observations Y t 0 , . . . , Y tm may be written as an expectation of an explicit functional of the discrete path Xt 0 , . . . , Xtm of the signal given the observations Y t 0 , . . . , Y tm . Therefore, we reduce the initial problem to the characterization of the conditional distribution of Xt 0 , . . . , Xtm given Y t 0 , . . . , Y tm .

Theorem 2.4. Set Xm = ( Xt 0 , . . . , Xtm ) and Y m = (Y t 0 , . . . , Y tm ). We have:

P τ X a > t n | Y t 0 , . . . , Y tm = Ψ(Y t 0 , . . . , Y tm ), (15)
where for y m := (y 0 , . . . , y m ),

Ψ(y m ) = E Φ Xm Y m = y m . ( 16
)
The function Φ is defined for every z m = (z 0 , . . . , z m ) by

Φ(z m ) = F (t m , t n , z m ) m-1 k=0 G z k ,z k+1 ∆ k σ 2 k (a), with
F (t m , t n , x) = E n-1 k=m G Xt k , Xt k+1 ∆ k σ 2 k (a) Xtm = x .
Proof. We may show similarly to Equation [START_REF] Graf | Foundations of Quantization for Probability Distributions[END_REF]

Xu > a ( Xt k ) k=0,...,m Xtm = E n-1 k=m G Xt k , Xt k+1 ∆ k σ 2 k (a) Xtm .
On the other hand, combining a successive conditioning rule with the independence between W and W gives :

P τ X a > t n | Y t 0 , . . . , Y tm = E E 1 1 {τ X a >tm} F (t m , t n , Xtm )|(W t k , W t k ) k=0,...,m | Y m = E E 1 1 {τ X a >tm} F (t n , t m , Xtm )|(W t k ) k=0,...,m | Y m = E F (t m , t n , Xtm )P τ X a > t m |(W t k ) k=0,...,m | Y m = E F (t m , t n , Xtm )P τ X a > t m |( Xt k ) k=0,...,m | Y m . Now, using Lemma 2.3 once again, we have P τ X a > t m |( Xt k ) k=0,...,m = m-1 k=0 G Xt k , Xt k+1 ∆ k σ 2 k (a) so that P τ X a > t n | Y t 0 , . . . , Y tm = E Φ( Xm ) | Y t 0 , . . . , Y tm where for every x m := (x 0 , . . . , x m ), Φ(x m ) = F (t m , t n , x m ) m-1 k=0 G x k ,x k+1 ∆ k σ 2 k (a).
This completes the proof.

It remains now to characterize the law of Xm |Y m = y m appearing in (16). This is the purpose of the following theorem, in which we shall write the conditional survival probability in a usual form with respect to standard nonlinear filtering problems. We set from now on

σ k = σ(x k , t k ), b k = b(x k , t k ), δ k = δ(y k , t k ), h k = h(y k , x k , t k ), and ν k = ν(y k , t k ).
We next give the main result of the paper. Theorem 2.5. We have:

P τ X a > t n | Y t 0 , . . . , Y tm = Ψ(Y t 0 , . . . , Y tm ), (17)
where for y m := (y 0 , . . . , y m ),

Ψ(y m ) = E F (t m , t n , Xtm )K m L m y E[L m y ] , (18)
with

K m = m-1 k=0 G Xt k , Xt k+1 ∆ k σ 2 k (a) and L m y = m-1 k=0 g k ( Xt k , y k ; Xt k+1 , y k+1 ).
The function g k is defined by

g k (x k , y k ; x k+1 , y k+1 ) = 1 (2π∆ k ) 3/2 σ 2 k δ k exp - ν 2 k 2δ 2 k ∆ k x k+1 -m 1 k σ k - y k+1 -m 2 k ν k 2 (19) with m 1 k := x k + b k ∆ k and m 2 k := y k + h k ∆ k .
Proof. It follows from Theorem 2.4 that

P τ X a > t n | Y t 0 , . . . , Y tm = E Φ( Xm ) | Y t 0 , . . . , Y tm where for every x m := (x 0 , . . . , x m ), Φ(x m ) = F (t m , t n , x m ) m-1 k=0 G x k ,x k+1 ∆ k σ 2 k (a).
It remains to characterize the conditional distribution of Xm given Y m . Recall that the dynamics of the processes ( Xt k ) and (Y t k ) are given for k = 0, . . . , m by

Xt k+1 = Xt k + b( Xt k , t k )∆ k + σ( Xt k , t k )(W t k+1 -W t k ) Y t k+1 = Y t k + h(Y t k , Xt k , t k )∆ k + ν(Y t k , t k )(W t k+1 -W t k ) + δ(Y t k , t k )( W t k+1 -W t k ) (20) Let P x k ,y k ( Xt k+1 ∈ dx k+1 , Y t k+1 ∈ dy k+1 ) denote the density function of the couple ( Xt k+1 , Y t k+1 ) given ( Xt k , Y t k ) = (x k , y k ).
Using Bayes formula and the Markov property of the process ( Xt k , Y t k ) k we get:

P x 0 ,y 0 ( Xt 1 ∈ dx 1 , . . . , Xtm ∈ dx m |Y t 1 ∈ dy 1 , . . . , Y tm ∈ dy m ) = P x 0 ,y 0 ( Xt 1 ∈ dx 1 , Y t 1 ∈ dy 1 ); . . . ; ( Xtm ∈ dx m , Y tm ∈ dy m ) P x 0 ,y 0 (Y t 1 ∈ dy 1 , . . . , Y tm ∈ dy m ) = P x 0 ,y 0 Xt 1 ∈ dx 1 , Y t 1 ∈ dy 1 ) × • • • × P x m-1 ,y m-1 ( Xtm ∈ dx m , Y tm ∈ dy m ) P x 0 ,y 0 (Y t 1 ∈ dy 1 , . . . , Y tm ∈ dy m ) .
For every x m = (x 0 , . . . , x m ) and y m = (y 0 , . . . , y m ), set :

N k (x m ; y m ) := P x 0 ,y 0 Xt 1 ∈ dx 1 , Y t 1 ∈ dy 1 ) × • • • × P x m-1 ,y m-1 ( Xtm ∈ dx m , Y tm ∈ dy m ).
With this notation, the denominator reads:

P x 0 ,y 0 (Y t 1 ∈ dy 1 , . . . , Y tm ∈ dy m ) = N k (x m ; y m )dx 1 . . . dx m .
On the other hand, the random vector

( Xt k+1 , Y t k+1 ) ( Xt k , Y t k ) = (x k , y k
) has a Gaussian distribution with mean m k and covariance matrix Σ k given for every k = 0, . . . , m -1 by

m k = x k + b k ∆ k y k + h k ∆ k , Σ k = ∆ k σ 2 k σ k ν k σ k ν k ν 2 k + δ 2 k . (21)
Then, for every k = 0, . . . , m -1, the density

f k (x k , y k ; x k+1 , y k+1 ) of ( Xt k+1 , Y t k+1 ) ( Xt k , Y t k ) = (x k , y k ) reads f k (x k , y k ; x k+1 , y k+1 ) = 1 2πσ k δ k ∆ k exp - ν 2 k + δ 2 k 2δ 2 k ∆ k (x k+1 -m 1 k ) 2 σ 2 k + (y k+1 -m 2 k ) 2 ν 2 k + δ 2 k - 2ν k σ k (ν 2 k + δ 2 k ) (x k+1 -m 1 k )(y k+1 -m 2 k ) , (22)
where

m 1 k := x k + b k ∆ k and m 2 k := y k + h k ∆ k . As a consequence, P τ X a > t n | (Y t 0 , . . . , Y tm ) = y m = E Φ( Xm ) | (Y t 0 , . . . , Y tm ) = y m = Φ(x m )N k (x m ; y m )dx 1 . . . dx m N k (x m ; y m )dx 1 . . . dx m . (23)
Now, we know that the random variable Xt k+1 | Xt k = x k has a Gaussian distribution with mean m 1 k and variance σ 2 k ∆ k . Its density P k (x k , x k+1 )dx k+1 is therefore given by :

P k (x k , x k+1 )dx k+1 = 1 √ 2π∆ k σ k exp - (x k+1 -m 1 k ) 2 2σ 2 k ∆ k dx k+1 .
Then, by definition of N k , we may write:

N k (x m ; y m )dx 1 . . . dx m = m-1 k=0 f k (x k , y k ; x k+1 , y k+1 )dx k+1 = m-1 k=0 f k (x k , y k ; x k+1 , y k+1 ) P k (x k , x k+1 ) P k (x k , x k+1 )dx k+1 = m-1 k=0 f k (x k , y k ; x k+1 , y k+1 ) P k (x k , x k+1 ) m-1 k=0 P k (x k , x k+1 )dx k+1 = m-1 k=0 g k (x k , y k ; x k+1 , y k+1 ) P x 0 ,y 0 ( Xt 1 ∈ dx 1 , . . . , Xtm ∈ dx m ) = E[L m y ],
where g k is defined by ( 19) and the next-to-last equality follows from the Markov property of the process ( Xt k , k = 0, . . . , m). Now, looking at the numerator of ( 23), similar computations lead to:

Φ(x m )N k (x m ; y m )dx 1 . . . dx m = F (t m , t n , x m ) m-1 k=0 G x k ,x k+1 ∆ k σ 2 k (a)f k (x k , y k ; x k+1 , y k+1 )dx k+1 , = F (t m , t n , x m ) m-1 k=0 G x k ,x k+1 ∆ k σ 2 k (a) × m-1 k=0 g k (x k , y k ; x k+1 , y k+1 ) P x 0 ,y 0 ( Xt 1 ∈ dx 1 , . . . , Xtm ∈ dx m ) =E F (t m , t n , Xtm )K m L m y .
Finally, going back to (23), we obtain

P τ X a > t n | (Y t 0 , . . . , Y tm ) = y m = E F (t m , t n , Xtm )K m L m y E[L m y ] which is the announced result.
Remark 2.1. Remark that more generally we have for every real valued bounded function f ,

E[f ( Xtm )|Y t 0 = y 0 , . . . , Y tm = y m ] = E f ( Xtm )L m y E[L m y ] . (24)
In this case we can estimate the right hand side quantity of (24) using recursive algorithms. In our setting, we can adapt these algorithms by noting that the expression (18) may be read in the similar form of (24) as:

P τ X a > t n | (Y t 0 , . . . , Y tm ) = y m = E F (t m , t n , Xtm )L ′ m y E[L m y ]
where

L ′ m y = m-1 k=0 g k ( Xt k , y k ; Xt k+1 , y k+1 )G Xt k , Xt k+1 ∆ k σ 2 k (a)
and where the involved functions are defined in Theorem 2.5.

Considering the model given in Example 1, one may apply the result of Theorem 2.5 using the continuous Euler process X of the signal. However since in this framework the solutions of the stochastic differential equations are explicit, we refrain from using the Euler scheme to avoid adding additional error. In the next result we deduce a similar representation of the conditional survival probability using the explicit solutions of the stochastic differential equations of the model (3).

Corollary 2.6. Consider that the signal process X and the observation process Y evolve following the stochastic differential equations given in [START_REF] Coculescu | Valuation of default sensitive claims under imperfect information[END_REF]:

dX t = X t (µdt + σdW t ), X 0 = x 0 , dY t = Y t (rdt + νdW t + δd W t ), Y 0 = y 0 . ( 25
)
Then,

Ψ(y m ) = E F (t m , t n , X tm )K m L m y E[L m y ] , (26)
with

K m = m-1 k=0 G Xt k ,Xt k+1 ∆ k σ 2 k (a) and L m y = m-1 k=0 g k (X t k , y k ; X t k+1 , y k+1 ).
The function g k is defined by

g k (x k , y k ; x k+1 , y k+1 ) = 1 (2π∆ k ) 3/2 π k exp - ν 2 2δ 2 ∆ k log x k+1 -m 1 k σ - log y k+1 -m 2 k ν 2 (27) with π k := σ 2 δx 2 k+1 y k+1 , m 1 k := log x k + (µ - 1 2 σ 2 )∆ k and m 2 k := log y k + (r - 1 2 ν 2 - 1 2 δ 2 )∆ k .
Proof. It follows from Itô formula that for every s ≤ t,

X t = X s exp (µ -1 2 σ 2 )(t -s) + σ(W t -W s ) Y t = Y s exp (r -1 2 ν 2 -1 2 δ 2 )(t -s) + ν(W t -W s ) + δ( W t -W s ) . (28)
Then, for every k = 0, . . . , m -1 the random vector

(X t k+1 , Y t k+1 )|(X t k , Y t k ) = (x k , y k ) has a
bivariate lognormal distribution with mean m k and covariance matrix Σ given by

m k = log x k + (µ -1 2 σ 2 )∆ k log y k + (r -1 2 ν 2 -1 2 δ 2 )∆ k , Σ = ∆ k σ 2 σν σν ν 2 + δ 2 .
Hence, its density reads (setting

µ k = σδx k+1 y k+1 ) f k (x k , y k ; x k+1 , y k+1 ) = 1 2π∆ k µ k exp - ν 2 + δ 2 2δ 2 ∆ k (log x k+1 -m 1 k ) 2 σ 2 - 2ν σ(ν 2 + δ 2 ) (log x k+1 -m 1 k )(log y k+1 -m 2 k ) + (log y k+1 -m 2 k ) 2 ν 2 + δ 2 ,
where

m 1 k = log x k + (µ -1 2 σ 2 )∆ k and m 2 k := log y k + (r -1 2 ν 2 -1 2 δ 2 )∆ k .
Furthermore, for every k = 0, . . . , m -1, the random variable X t k+1 |X t k = x k has a lognormal distribution with mean m 1 k and variance σ 2 so that its density distribution P k (x k , x k+1 ) is given by

P k (x k , x k+1 ) = 1 √ 2π∆ k σx k+1 exp - (log x k+1 -m 1 k ) 2 2σ 2 ∆ k .
We then conclude the proof by using the same arguments than those of the proof of Theorem 2.5.

There exist several methods to estimate the above representation of the conditional survival probability. These methods involve, amount others, Monte Carlo simulations and optimal quantization methods. Owing to the numerical performance of the optimal quantization method due for example to its fast performability as soon as the optimal grids are obtained, we will use optimal quantization methods to estimate the conditional survival probability.

The use of optimal quantization methods to estimate the filter supposes to have numerical access to optimal (or stationary) quantizers of the marginals of the signal process. One may use the optimal vector quantization method (as done in the seminal work [START_REF] Pagès | Optimal quantization methods for nonlinear filtering with discrete time observations[END_REF] and used in [START_REF] Callegaro | An application to credit risk of optimal quantization methods for nonlinear filtering[END_REF]) to estimate these marginals. This method requires the use of some algorithms, like stochastic algorithms or Lloyd's algorithm, to obtain numerically the optimal (or stationary) quantizers of the marginals of the signal process. Given these optimal quantizers, the filter estimation is obtained quite instantaneously. However, the step of search of stationary quantizers is very time consuming.

We propose here an alternative method, the (quadratic) marginal functional quantization method (introduce in [START_REF] Sagna | Pricing of barrier options by marginal functional quantization method[END_REF] to price barrier options), to quantize the marginals of the signal process. The marginal functional quantization method consists first in considering the ordinary differential equation (ODE) resulting to the substitution of the Brownian motion appearing in the dynamics of the signal process by a quadratic quantization of the Brownian motion. Then, by constructing some "good" marginal quantization of the signal process based on the solution of the previous ODE's, we will show how to estimate the nonlinear filter. Since this procedure is based on the quantization of the Brownian motion and skips the use of algorithms to perform the stationary quantizers, the computation of the marginal quantizers is quite instantaneous. This reduce drastically the time computation of the procedure with respect to the vector quantization method since we skip the step of the use of algorithms search of marginal quantizers by using instead, the marginals of the functional quantization of the signal process.

In the rest of the paper we deal with the estimations methods of the conditional survival probability.

Numerical tools

Our aim in this section is to derive a way to compute numerically the conditional survival probability using Equation (18). To this end, we have to estimate three quantities: the quantity F (t m , t n , Xtm ), the expectation E[ F (t m , t n , Xtm )K m L m y ] as soon as F (t m , t n , Xtm ) is estimated, and finally, the expectation E[L m y ]. Both expectations will be estimated using marginal functional quantization method. The probability F (t m , t n , Xtm ) will be estimated by Monte Carlo simulations.

Before dealing with the estimation tools, we shall recall first some basic results about both optimal vector quantization and functional quantization. Then, we will show how to construct the marginal functional quantization process which will be used to estimate the quantities of interest.

Overview on optimal quantization methods

The optimal vector quantization on a grid Γ = {x 1 , • • • , x N } (which will be called a quantizer) of an R d -valued random vector X defined on a probability space (Ω, A, P) with finite r-th moment and probability distribution P X consists in finding the best approximation of X by a Borel function of X taking at most N values. This turns out to find the solution of the following minimization problem:

e N,r (X) = inf { X -X Γ r , Γ ⊂ R d , |Γ| ≤ N }, (29)
where X r := (E[|X| r ]) 1/r and where X Γ := N i=1 x i 1 1 {X∈C i (Γ)} is the quantization of X (we will write X N instead of X Γ ) on the grid Γ and (C i (Γ)) i=1,...,N corresponds to a Voronoi tessellation of R d (with respect to a norm | • | on R d ), that is, a Borel partition of R d satisfying for every i,

C i (Γ) ⊂ {x ∈ R d : |x -x i | = min j=1,...,N |x -x j |}.
We know that for every N ≥ 1, the infimum in (29) is reached at one grid Γ ⋆ at least, called a L r -optimal N -quantizer. It is also known that if card(supp(P X )) ≥ N then |Γ| = N (see e.g. [START_REF] Graf | Foundations of Quantization for Probability Distributions[END_REF] or [START_REF] Pagès | A space quantization method for numerical integration[END_REF]). Moreover, the L r -mean quantization error e N,r (X) decreases to zero at an N -1/d -rate as the size N of the grid Γ goes to infinity. This convergence rate has been investigated in [START_REF] Bucklew | Multidimensional asymptotic quantization theory with r-th power distribution measures[END_REF] and [START_REF] Zador | Asymptotic quantization error of continuous signals and the quantization dimension[END_REF] for absolutely continuous probability measures under the quadratic norm on R d , and studied in great details in [START_REF] Graf | Foundations of Quantization for Probability Distributions[END_REF] under an arbitrary norm on R d for absolutely continuous measures and some singular measures.

From the numerical integration viewpoint, finding an optimal quantization grid Γ ⋆ may be a challenging task. In practice (we will only consider the quadratic case, i.e. when r = 2) we are sometimes led to find some "good" quantizations Proj Γ (X) (with Γ = {x 1 , . . . , x N }) which are close to X in distribution, so that for every Borel function

F : R d → R, we can approximate E[F (X)] by E F X N = N i=1 F (x i ) p i , (30)
where p i = P Proj Γ (X) = x i . Amount "good" quantizations of X we have stationary quantizers. A grid Γ inducing the quantization X N of X is said stationary if

∀ i = j, x i = x j and P (X ∈ ∪ i ∂C i (Γ)) = 0 (31) and E X| X N = X N .
The stationary quantizers search is based on zero search recursive procedures like Newton algorithm in the one dimensional framework and some algorithms as Lloyd's I algorithms (see e.g. [START_REF] Gersho | Vector Quantization and Signal Compression[END_REF]), the Competitive Learning Vector Quantization (CLVQ) algorithm (see [START_REF] Gersho | Vector Quantization and Signal Compression[END_REF]) or stochastic algorithms (see [START_REF] Pagès | Optimal quadratic quantization for numerics: the gaussian case[END_REF]) in the multidimensional framework. Note that optimal quantizers estimates of the multivariate Gaussian random vector are available in the website www.quantize.math-fi.com.

We next recall some error bounds induced from approximating E[F (X)] by (30). Let Γ be a stationary quantizer and F be a Borel function on

R d . (i) If F is convex then E F ( X N ) ≤ E[F (X)].
(ii) Lipschitz functions:

-If F is Lipschitz continuous then (this error bound doesn't require the quantizer Γ to be stationary)

E[F (X)] -E[F ( X N )] ≤ [F ] Lip X -X N 2 , where [F ] Lip := sup x =y |F (x) -F (y)| |x -y| .
-Let θ : R d → R + be a nonnegative convex function such that θ(X) ∈ L 2 (P). If F is locally Lipschitz with at most θ-growth, i.e. |F (x) -

F (y)| ≤ [F ] Lip |x -y| (θ(x) + θ(y)) then F (X) ∈ L 1 (P) and E[F (X)] -E[F ( X N )] ≤ 2[F ] Lip X -X N 2 θ(X) 2 .
(iii) Differentiable functionals:

If F is differentiable on R d with an α-Hölder differential DF (α ∈ [0, 1]), then E[F (X)] -E[F ( X N )] ≤ [DF ] α X -X N 1+α 2 .
Other error bounds related to the regularity of F may be found in [START_REF] Pagès | Functional quantization for numerics with an application to option pricing[END_REF]).

The optimal vector quantization may be extended to random vectors with values in a set of infinite dimension, in particular to stochastic processes viewed as random variables with values in L 2 ([0, T ], dt) endowed with the norm

E |X| 2 L 2 T = T 0 E[X 2 s ]ds < +∞.
The functional quantization of the stochastic process (X t ) t∈[0,T ] with dynamics

dX t = b(X t , t)dt + σ(X t , t)dW t , X 0 = x ∈ R,
is based on the functional quantization of the Brownian motion W . One way to quantize the Brownian motion is to use the optimal product quantization using its Karhunen-Loève expansion which reads :

W L 2 T = n≥1 λ n ξ n e n ,
where ξ n := W, e n / √ λ n , n ≥ 1, is a sequence of i.i.d. random variables with standard normal distribution and

e n (t) = 2 T sin π n - 1 2 t T , λ n = T π n -1 2 2 , n ≥ 1.
In fact, from the previous expansion, a functional quantization W of the process W of size at most N is defined by

W N t = n≥1 λ n ξx Nn n e n (t) (32)
where ξx Nn n (with x Nn = {x Nn 1 , . . . , x Nn Nn }) is the optimal N n -quantization of ξ n and N 1 ו • •×N n ≤ N , with N i ≥ 2 for i ≤ n, and N k = 1 for k ≥ n + 1, so that the expansion defined in (32) is a finite sum. The product quantizer χ that produces the above Voronoi quantization W is defined by

χ i (t) = n≥1 λ n x Nn in e n (t), i = (i 1 , • • • , i n , • • • ) ∈ n≥1 {1, • • • , N n } .
and for every multi

-index i ∈ n≥1 {1, • • • , N n }, the associated Voronoi cell of χ is C i (χ) = n≥1 λ n C in (x Nn ).
The optimal product quantizer of size at most N , denoted W N , of the Brownian motion is defined as the solution of the following optimization problem:

min W -W N 2 , N 1 , • • • , N n ≥ 1, N 1 × • • • × N n ≤ N, N ≥ 1 . ( 33
)
Moreover this optimal product quantizer induces a rate-optimal sequence of quantizers (see e.g. [START_REF] Luschgy | Functional quantization of a class of Brownian diffusions: A constructive approach[END_REF] for more details), i.e.

W -W N 2 ≤ K W T (log N ) 1 2
for some real constant K W > 0. To define a functional quantization of the stochastic process (X t ) t≥0 consider (χ N ) N ≥1 a sequence of rate-optimal product quantizers of the Brownian motion and, for every multi

-index i ∈ n≥1 {1, • • • , N n }, with N 1 × • • • × N n ≤ N , consider
x i the solution of the following integral equation

dx i (t) = b(x i (t), t) - 1 2 σσ ′ (x i (t), t) dt + σ(x i (t), t)dχ i (t), (34)
where σ ′ is the derivative of σ. We define the functional (non-Voronoi) quantization X process of the stochastic process X of size at most N by

X N = N i=1 x i 1 C i (χ) (W ), N ≥ 1.
We have the following result.

Proposition 3.1 (See [START_REF] Luschgy | Functional quantization of a class of Brownian diffusions: A constructive approach[END_REF]). Under some suitable conditions on the coefficients of the diffusion, which in the homogeneous case are equivalent to: b is differentiable, σ is positive twice differentiable and b

′ -b σ ′ σ -1 2 σσ ′′ is bounded, we have X -X N 2 = O (log N ) -1 2 .
We observe that for any initial value x of the quantized process X, the marginals X t are of size d N = N 1 × • • • × N n where (N 1 , . . . , N n ) is the solution of (33). To define the marginal functional quantization process (still be denoted by ( X t )) of the process (X t ) we order the values of the marginals and define the marginal quantizer x d N (t) for every t > 0 by

x d N (t) = {x d N 1 (t), . . . , x d N d N (t)} and x d N (0) = x and define likewise the marginal quantizers of W t by χ d N (t) = {χ d N 1 (t), . . . , χ d N d N (t)}.
The marginal functional quantization process of size at most N is then defined by

X N t = d N i=1 x d N i (t)1 {Wt∈C i (χ d N (t))} .
Remark that the use of the marginal functional quantization method can not be justified from the theoretical point of view since we do not know yet the rate of convergence of the marginals of the quantized process to the marginals of the initial process. However, this method has proved its efficiency from the numerical viewpoint when used to estimate barrier option by optimal quantization, see [START_REF] Sagna | Pricing of barrier options by marginal functional quantization method[END_REF] (not that the considered marginal quantization in [START_REF] Sagna | Pricing of barrier options by marginal functional quantization method[END_REF] is a little bit different from the one considered in this paper, nevertheless the numerical results are the same, up to at least a 10 -3 absolute error order).

Let us come back to the problem of interest and let the functionals π y,m and ̟ y,m be defined for every bounded measurable function f by

π y,m f = E f ( Xtm )L m y K m and ̟ y,m f = E f ( Xtm )L m y .
Then,

Π y,m F (t m , t n , •) := P τ X a > t n | Y t 0 = y 0 , . . . , Y tm = y m = π y,m F (t m , t n , •) ̟ y,m 1 ,
where 1(x) = 1, for every real x. Then, it will be enough to show how to estimate π y,m F (t m , t n , •) since ̟ y,m 1 is estimated similarly. We discuss in the section below the estimation of π y,m .

Estimation of π y,m by marginal functional quantization

Our aim is to estimate

π y,m F (t m , t n , •) = E F (t m , t n , Xtm ) m-1 k=0 g a k ( Xt k , y k ; Xt k+1 , y k+1 )
by marginal functional quantization (with X ≡ X in the model ( 3)), where

g a k ( Xt k , y k ; Xt k+1 , y k+1 ) = g k ( Xt k , y k ; Xt k+1 , y k+1 ) × G Xt k , Xt k+1 ∆ k σ 2 k (a).
We deal with a general setting: the estimation of π y,m f for any bounded and measurable function f . Our main reference is [START_REF] Pagès | Optimal quantization methods for nonlinear filtering with discrete time observations[END_REF] where π y,m f been estimated using marginal vector quantization methods. Let us define for every k = 1, . . . , m the transition kernel H y,k by

H y,k f (x) = E f ( Xt k )g a k-1 (x, y k-1 ; Xt k , y k )| Xt k-1 = x = f (z)g a k-1 (x, y k-1 ; z, y k )P k (x, dz)
where

P k (x, dz) is the density of the random variable Xt k | Xt k-1 = x. We set H y,0 f = E[f ( X0 )] = f (z)µ(dz). (35)
Then, we have (setting for every k,

F X t k = σ({ Xt i , i = 0, . . . , k})) π y,k f = E E f ( Xt k ) k-1 i=0 g a i ( Xt i , y i ; Xt i+1 , y i+1 ) F X t k-1 = E E f ( Xt k )g a k-1 ( Xt k-1 , y k-1 ; Xt k , y k )|F X t k-1 k-2 i=0 g a i ( Xt i , y i ; Xt i+1 , y i+1 ) = E E f ( Xt k )g a k-1 ( Xt k-1 , y k-1 ; Xt k , y k )| Xt k-1 k-2 i=0 g a i ( Xt i , y i ; Xt i+1 , y i+1 ) ,
where the last equality is a consequence of the Markov property of the process ( Xt k ). Thus, we deduce that for every k = 1, . . . , m,

π y,k f = E H y,k f ( Xt k-1 ) k-2 i=0 g a i ( Xt i , y i ; Xt i+1 , y i+1 ) = π y,k-1 H y,k f.
It follows that π y,m f can be computed by the following recursive formula:

π y,m f = (H y,0 • H y,1 • • • • • H y,m )f. (36)
Therefore, to achieve the estimation of π y,m , it remains to estimate the kernels H y,k . This will be done by marginal functional quantization. Consider time discretization steps t k , k = 0, . . . , m and let χ N k := {χ 1 k , . . . , χ N k k } be a N kquantizer of W t k (we will consider the marginal functional quantization of Brownian motion so that N k = d N for every k). Suppose that we have also access to the marginal functional quantization process ( X t k ) k of the process ( Xt ) t≥0 the time steps t k , k = m, . . . , n:

{x 1 k , . . . , x N k k }, of sizes N k (keep in mind that N k = d N for every k).
It follows that the transition kernels H y,k may be estimated for every k = 1, . . . , m by

H y,k = N k j=1 H ij y,k δ x i k-1
, where

H ij y,k = g a k-1 (x i k-1 , y k-1 ; x j k , y k ) pij k , i = 1, . . . , N k-1 ; j = 1, . . . , N k (37)
and where the pij k 's correspond to the estimation of the transition probabilities from

X t k-1 = x i k-1 to X t k = x j k : pij k = P( X t k = x j k | X t k-1 = x i k-1 ), i = 1, . . . , N k-1 ; j = 1, . . . , N k . (38)
Finally, one will perform the estimation

π y,m = H y,0 • H y,1 •• • •• H y,m of π y,m
as soon as we will be able to compute pij k . In the proposition below we show how to compute these probabilities from the cumulative distribution function of the random variable W t k |W t k-1 = x, which has a Gaussian distribution with mean x and variance ∆ k (we refer to [START_REF] Sagna | Pricing of barrier options by marginal functional quantization method[END_REF] for a similar result).

Proposition 3.2. The transition probabilities can be estimated for every

k = 1, . . . , m by pij k ≈ N χ j+ k ; χ i k-1 -N χ j- k ; χ i k-1 , (39)
where N (• ; χ i k-1 ) is the cumulative distribution function of the normal distribution with mean χ i k-1

and variance ∆ k and where for every k = 1, . . . , m,

     χ j+ k = χ j+1 k +χ j k 2 ; χ j- k = χ j-1 k +χ j k 2 j = 1, . . . , N k χ 1- k = -∞; χ N + k k = +∞.
Proof. One has for every k = 1, . . . , m,

pij k = P W t k ∈ C j (χ N k ) W t k-1 ∈ C i (χ N k-1 ) = P W t k ≤ χ j+ k W t k-1 ∈ C i (χ N k-1 ) -P W t k ≤ χ j- k W t k-1 ∈ C i (χ N k-1
) . On the other hand, we have for every x ∈ R,

P W t k ≤ x W t k-1 ∈ C i (χ N k-1 ) = P W t k ≤ x; W t k-1 ∈ C i (χ N k-1 ) P W t k-1 ∈ C i (χ N k-1
) .

Set C i = C i (χ N k-1
). Then the numerator on the right hand side of the previous equation may be expressed as

P W t k ≤ x; W t k-1 ∈ C i = x -∞ C i P(W t k ∈ dx|W t k-1 = y)P(W t k-1 ∈ dy) dx = C i N (x; y)P(W t k-1 ∈ dy) (40) ≈ N (x; χ i k-1 )P W t k-1 ∈ C i ,
where N (•; y) is the cumulative distribution function of W t k |W t k-1 = y. The last quantity is the approximation of (40) by optimal quantization with one grid point, considering that {χ i k-1 } is the quantizer of size one of the random variable W t k-1 over the cell C i (χ N k-1 ).

Following the previous approach, the estimations π y,m and ̟ y,m of π y,m and ̟ y,m are computed from the following recursive formulae:

   π y,0 = H y,0 π y,k = π y,k-1 H y,k := N k-1 i=1 H i,j y,k π i y,k-1 j=1,...,N k , k = 1, . . . , m (41)
where

H ij y,k = g a k-1 (x i k-1 , y k-1 ; x j k , y k ) pij k , i = 1, . . . , d N ; j = 1, . . . , d N ; and    ̟ y,0 = Υ y,0 ̟ y,k = ̟ y,k-1 Υ y,k := N k-1 i=1 Υ i,j y,k ̟ i y,k-1 j=1,...,N k , k = 1, . . . , m (42)
where

Υ ij y,k = g k-1 (x i k-1 , y k-1 ; x j k , y k ) pij k , i = 1, . . . , N k-1 ; j = 1, . . . , N k .
Then we approximate Π y,m by Π y,m given by

Π y,m = Nm i=1 Π i y,m δ x i m with Π i y,m = π i y,m Nm j=1 ̟ i y,k-1 , i = 1, . . . , N m .
Recall that our aim is to estimate

Π y,m F (t m , t n , •) = π y,m F (t m , t n , •) ̟ y,m 1 ,
where the function F (t m , t n , •) is defined for every x ≥ 0 by

F (t m , t n , x) = E n-1 k=m G Xt k , Xt k+1 ∆ k σ 2 k (a) Xtm = x .
We shall take the estimation:

Π y,m F (t m , t n , •) = Nm i=1 Π i y,m F (t m , t n , x i m ). (43)
Remark that the initial function F (s, t, X s ) = P inf s≤u≤t X u > a|X s (which has been estimated by F (s, t, Xs )) has semi-closed expression in some specific models, like in the model [START_REF] Coculescu | Valuation of default sensitive claims under imperfect information[END_REF] in which case it is given by

P inf s≤u≤t X u > a|X s = N (h 1 (X s , t -s)) - a X s σ -2 (µ-σ 2 /2) N (h 2 (X s , t -s)) (44)
where

h 1 (x, u) = 1 σ √ u log(x/a) + µ - 1 2 σ 2 u , h 2 (x, u) = 1 σ √ u log(a/x) + µ - 1 2 σ 2 u
and where N (•) is the cumulative distribution function of the standard Gaussian distribution.

Except in these specific cases, the function F (t m , t n , •) has to be estimated, and we shall do it by Monte Carlo methods.

Estimation of F (t m , t n , •) by Monte Carlo

It remains to estimate the function F (t m , t n , x) defined by

F (t m , t n , x) = E n-1 k=m G Xt k , Xt k+1 ∆ k σ 2 k (a) Xtm = x
by Monte Carlo simulations. The steps of the Monte Carlo procedure are the following.

1. Let us consider regular time discretization steps {t m , t m+1 , . . . , t n } over [t m , t n ] and let M be the number of trials. We simulate for every j = 1, . . . , M the sample path ( Xj t k ) k=m,...,n , with Xj tm = x for every j.

Setting

p j m,n (x, a) := n-1 k=m G Xj t k , Xj t k+1 ∆ k σ 2 k (a),
we estimate F (t m , t n , x) by

F M (t m , t n , x) = 1 M M j=1 p j m,n (x, a). (45)
Consequently, integrating both formulae (43) and (45), the conditional survival probability Π y,m F (t m , t n , •) will be estimated (for a fixed trajectory (y 0 , . . . , y m ) of the observation process (Y t 0 , . . . , Y tm )) by

Π y,m F M (t m , t n , •) = 1 M Nm i=1 M j=1 Π i y,m p j m,n (x i m , a). ( 46
)
Remark 3.1. It follows from Monte Carlo error analysis that for every x ≥ 0,

F (t m , t n , x) -F M (t m , t n , x) 2 = O 1 √ M . (47)

Error analysis

We target to give in this section the error bound resulting from the estimation of P(τ X a > t|Y t 0 , . . . , Y tm ) by

Π y,m F M (t m , t n , •) = 1 M Nm i=1 M j=1 Π i y,m p j m,n (x i m , a).
We shall first give the error bound due to the fact that we compute the survival function of the random variable τ X a instead of τ X a , and then the error bound due to the simulations procedures. To this end, the following definitions and assumptions are needed. Remark that that in our framework, the transition operators P k (x, dy), k = 1, • • • , m are Lipschitz, so that we set

[P ] Lip := max k=1,••• ,m [P k ] Lip < +∞.
We furthermore define some useful quantities which appear in the error bound in the following.

(i) For every k = 1, • • • , m, we set

K m g := max k=1,••• ,m g k ∞ ,
where g k ∞ is the supremum norm of the functions g k defined by (19).

(ii) For every

k = 1, • • • , m, let [g 1 k ]
Lip and [g 2 k ] Lip be so that for every x, x ′ , x, x ′ ∈ R and y, y ′ ∈ R,

|g a k (x, y, x ′ , y ′ ) -g a k ( x, y, x ′ , y ′ )| ≤ [g 1 k ] Lip (y, y ′ ) |x -x| + [g 2 k ] Lip (y, y ′ ) |x ′ -x ′ |.
Let us make now some assumptions which will be used to compute (see [START_REF] Gobet | Schémas dÕEuler pour diffusion tuée. Application aux options barrière[END_REF]) the convergence rate of the quantity E 1 1 {τ X >t} -1 1 {τ >t} towards 0.

(H1) b is a C ∞ b (R, R) function and σ is in C ∞ b (R, R). (H2) there exists σ 0 > 0 such that ∀x ∈ R, σ(x) 2 ≥ σ 2
0 (uniform ellipticity). Before giving the error bound associated to our estimation we recall the following useful results. Consider in this scope that τ X = inf{u ≥ 0, X u ∈ D} where D = (a, +∞) and X is the signal process. Let ( Xt k ) k=0,...,m the continuous Euler process taken at discrete times t k , k = 0, . . . , m and τ X = inf{u ≥ 0, Xu ∈ D}.

We have the following result.

Proposition 3.3 (see [START_REF] Gobet | Schémas dÕEuler pour diffusion tuée. Application aux options barrière[END_REF]). Let t > 0. Suppose that Assumptions (H1) and (H2) are fulfilled. Then, for every η ∈ (0, 1 2 [ there exists an increasing function K(T ) such that for every t ∈ [0, T ] and for every

x ∈ R, E x 1 1 {τ X >t} -1 1 {τ X >t} ≤ 1 n 1 2 -η K(T ) √ t ,
where n is the number of time discretization steps over [0, t].

The convergence rate of the filter approximation is given by the following theorem. Since we do not know the convergence rate and some properties as the stationary property of the marginals of the functional quantization we consider here that for every k, X N k t k denotes the marginal quantization of (X t k ) of size N k , obtained from the marginal vector quantization method (see [START_REF] Pagès | Optimal quantization methods for nonlinear filtering with discrete time observations[END_REF]). Theorem 3.4. (see [START_REF] Pagès | Optimal quantization methods for nonlinear filtering with discrete time observations[END_REF] for a similar result). We have for every p ≥ 1,

|Π y,m F (s, t, •) -Π y,m F (s, t, •)| ≤ K m g φ m (y) ∨ φm (y) m k=0 B m k ( F (s, t, •), y, p) Xt k -X N k t k p where φ m (y) := π y,m 1, φ m (y) := π y,m 1
and where

B m k (f, y, p) := (2 -δ 2,p )[P ] m-k Lip [f ] Lip + 2 f ∞ K m g [g 1 k+1 ] Lip (y k , y k+1 ) + [g 2 k ] Lip (y k-1 , y k ) + (2 -δ 2,p ) f ∞ K m g m j=k+1 [P ] j-k-1 Lip [g 1 j ] Lip (y j-1 , y j ) + [P ] Lip [g 2 j ] Lip (y j-1 , y j )
(with the convention that g 0 = g m+1 ≡ 0 and δ n,p is the usual Kronecker symbol).

Proof. The proof is similar to the proof of Theorem 3.1 in [START_REF] Pagès | Optimal quantization methods for nonlinear filtering with discrete time observations[END_REF].

Remark 3.2. One may remark that the function B m k (f, y, p) involves the norm f ∞ which for f = F (s, t, •) is bounded by 1.

Let us give now the error bounds induced from the approximation of P(τ X a > t n |Y t 0 , . . . , Y tm ) by Π y,m F M (t m , t n , •).

Theorem 3.5. (See [START_REF] Callegaro | An application to credit risk of optimal quantization methods for nonlinear filtering[END_REF] for a similar result) Suppose that the coefficients b and σ of the continuous signal process X are such that Assumptions (H1) and (H2) are satisfied and let η ∈ (0, 1 2 ]. Then,

P(τ X a > t n |Y t 0 , . . . , Y tm ) - 1 M Nm i=1 M j=1 Π i y,m p j m,n (x i m , a) ≤ O n -1 2 +η + O M -1 2 + K m g φ m (y) ∨ φm (y) m k=0 B m k ( F (s, t, •), y, p) Xt k -X x N k t k p ,
where the functions φ m , φ m , B m k (•, y, p) are introduced in Theorem 3.4.

Proof. We have

P(τ X a > t n |Y t 0 , . . . , Y tm ) - 1 M Nm i=1 M j=1 Π i y,m p j m,n (x i m , a) ≤ P(τ X a > t n |Y t 0 , . . . , Y tm ) -P(τ X a > t n |Y t 0 , . . . , Y tm ) + Π y,m F (t m , t n , •) -Π y,m F (t m , t n , •) + Π y,m F (t m , t n , •) -Π y,m F M (t m , t n , •) . (48)
Owing to Remark 3.1 and to the fact that Nm i=1 Π i y,m ≤ 1 we get

Π y,m F (s, t, •) -Π y,m F M (t m , t n , •) = Nm i=1 Π i y,m F (t m , t n , x i m ) - Nm i=1 Π i y,m F M (t m , t n , x i m ) ≤ Nm i=1 Π i y,m sup x≥0 | F (t m , t n , x) -F M (t m , t n , x)| ≤ sup x≥0 | F (t m , t n , x) -F M (t m , t n , x)| = O 1 √ M .
On the other hand, the error bound of the term Π y,m F (t m , t n , •) -Π y,m F (t m , t n , •) is given by Theorem 3.4. Now, let us consider the first term of the right hand side of Equation ( 48). We have:

E P(τ X a > t n |Y t 0 , . . . , Y tm ) -P(τ X a > t n |Y t 0 , . . . , Y tm ) = E E 1 1 {τ X a >tn} -1 1 {τ X a >tn} |Y t 0 , . . . , Y tm ≤ E 1 1 {τ X a >tn} -1 1 {τ X a >tn} = O(n -1 2 +η
), the last statement following from Proposition 3.3.

Numerical examples

We deal with numerical simulations in this section by considering two example of models.

In both examples we fix t m = 1 and, given a (simulated) trajectory of the observation process Y from 0 to t m , we estimate the conditional cumulative function P(τ X a > t n |Y t 0 , . . . , Y tm ) using formula (46), for t n varying 0.1 by 0.1 from 1.1 to 11 (where the time unit is expressed in years). Furthermore, we set the number m of discretization points over [0, t m ] equal to 50 and for every k = 1, • • • , m, the quantization grid size N k is set to 966 (as a consequence of the numerical solution of the Problem 33 for N = 10000, with the optimal decomposition (23, 7, 3, 2), see [START_REF] Pagès | Functional quantization for numerics with an application to option pricing[END_REF] for more detail), with N 0 = 1. All the programs have been coded using the C language on a CPU 2.7 GHz and 4 Go memory computer.

Example 3 (The "Black-Scholes" example). The first model is the one considered in Example 1 and Corollary 2.6 where the dynamics of the signal process X and the observation process Y are given by dX t = X t (µdt + σdW t ), X 0 = x 0 , dY t = Y t (µdt + σdW t + δd W t ), Y 0 = y 0 (49)

or equivalently

dY t Y t = dX t X t + δd W t .
We choose the following parameters (as in [START_REF] Coculescu | Valuation of default sensitive claims under imperfect information[END_REF]) of the model: µ = 0.03, σ = 0.03, x 0 = y 0 = 86.3 and a = 76. The numerical results are depicted in Figure 1 and Figure 2. In Figure 1, we draw three trajectories of the observation process Y for the same δ = 0.1 (left side graphic) and the corresponding cumulative functions P(τ X a ≤ t n |Y t 0 , . . . , Y tm ) on the right hand side graphics, t m = 1 and t n ∈ [1. [START_REF] Bucklew | Multidimensional asymptotic quantization theory with r-th power distribution measures[END_REF][START_REF] Pagès | Optimal quantization methods for nonlinear filtering with discrete time observations[END_REF], in years. We remark that, for a fixed time t n , the lower the trajectory is, the higher its probability is to hit the barrier.

The left hand side graphics of Figure 2 corresponds to three trajectories of the observation process Y for δ ∈ {0.1, 0.3, 0.5} and the right hand side graphics, to (a zoom of) the corresponding cumulative functions P(τ X a ≤ t n |Y t 0 , . . . , Y tm ), with t m = 1 and t n ∈ [1. [START_REF] Bucklew | Multidimensional asymptotic quantization theory with r-th power distribution measures[END_REF][START_REF] Pagès | Optimal quantization methods for nonlinear filtering with discrete time observations[END_REF], in years. We observe in this example that the noisier the observations are, the higher the probability is to hit the barrier a before a fixed time t n .

Note that in both examples, the function F (t m , t n , •) has been computed using formula (44) and the computation time to get one cumulative function P(τ X a ≤ t n |Y t 0 , . . . , Y tm ) for a given t n is about of 24 seconds. Example 4 (The Ornstein-Uhlenbeck example).

In the second model, we suppose that both the signal and the observation process evolve following the Ornstein-Uhlenbeck dynamics:

dX t = λ(θ -X t )dt + σdW t , X 0 = x 0 , dY t = λ(θ -Y t )dt + σdW t + δd W t , Y 0 = y 0 ( 50
)
or setting Z t = Y t -X t dZ t = -λZ t dt + δd W t , meaning that Z is still an Ornstein-Uhlenbeck process with mean value θ = 0 and with volatility δ.

The parameters are chosen as follows: λ = 0.18, θ = 0.35, σ = 0.12, x 0 = y 0 = 0.35 (as in [START_REF] Coculescu | Valuation of default sensitive claims under imperfect information[END_REF]) and a = 0.2. The numerical results are represented in Figure 3 where we depict three trajectories of the observation process Y for δ = 0.16 (left hand side graphics of Figure 3) and the associated cumulative functions P(τ X a ≤ t n |Y t 0 , . . . , Y tm ), with t m = 1 and t n ∈ [1. [START_REF] Bucklew | Multidimensional asymptotic quantization theory with r-th power distribution measures[END_REF][START_REF] Glasserman | Asymptotically optimal importance sampling and stratification for pricing path-dependent options[END_REF] in years (right hand side graphics of Figure 3). Once again, we remark that, as in the "Black Scholes example", for a fixed time t n , the lower the trajectory is, the higher its probability is to hit the barrier.

In this example, the function F (t m , t n , •) has been computed using Monte Carlo simulations of size M = 10 5 (see the formula (45)) with 50 discretization steps over [t m , t n ]. The computation time to get one cumulative function P(τ X a ≤ t n |Y t 0 , . . . , Y tm ) for a given t n is about 6.5 minutes.

Definition 3 . 1 .

 31 A probability transitionP on E is C-Lipschitz (with C > 0) if for any Lipschitz function f on E with ratio [f ] Lip , P f is Lipschitz with ratio [P f ] Lip ≤ C[f ]Lip . Then, one may define the Lipschitz ratio [P ] Lip by[P ] Lip = sup [P f ] Lip [f ] Lip , f a nonzero Lipschitz function < +∞.

Figure 1 :

 1 Figure 1: ("Black Scholes example") Three trajectories for the observation process Y and for δ = 0.1 (on the left), and the corresponding cumulative functions P(τ X a ≤ tn|Yt 0 , . . . , Yt m ) with tm = 1 and tn ∈ [1.1, 11] years (on the right).

Figure 2 :

 2 Figure 2: ("Black Scholes example") Left graphics corresponds to three trajectories for the observation process Y , for δ = 0.1 (Y UP), δ = 0.3 (Y MID), δ = 0.5 (Y DOWN), and the right hand side graphics correspond to a zoom of the cumulative functions P(τ X a ≤ tn|Yt 0 , . . . , Yt m ) with tm = 1 and tn ∈ [1.1, 11] years.

Figure 3 :

 3 Figure 3: ("Ornstein-Uhlenbeck example") Left graphics correspond to three trajectories for the observation process Y , for δ = 0.16 and the right hand side graphics correspond to the cumulative functions P(τ X a ≤ tn|Yt 0 , . . . , Yt m ) with tm = 1 and tn ∈ [1.1, 6] years.

Both author's researches are supported by an AMaMeF exchange grant and the "Chaire Risque de Crédit" of the French Banking Federation.


Valerio Cosentino
email: valerio.cosentino@fr.ibm.com

Macros Didonet Del Fabro

Adil El Ghali
email: elghali@lutin-userlab.fr

A model driven approach for bridging ILOG Rule Language and RIF

Nowadays many companies run their own business using Business Rule Management Systems (BRMS), that offer: a clear separation between decision logic and procedural structure; the ability to modify a rulebase set rather than processes and the reusability of rules across applications. All these factors allow a company to quickly react and align its policies to the ever-changing market needs. Despite these advantages, different BRMSs can be found on the market, each of them implementing a proprietary business rule language (ex.: JBoss uses Drools, IBM uses Ilog Rule Language, etc.). Rule Interchange Format(RIF) is a W3C open standard aiming at reducing the heterogeneity among business rule languages, which makes rules less reusable and interchageable. Our work is focused on providing an implementation based on a Model Driven approach for bridging Ilog Rule Language to RIF.

Introduction

The field of BRMS is characterized by a number of normative, open source, or proprietary systems and languages (Ilog JRules, JBoss Drools, etc.), allowing the expression of various solutions to business problems at a high abstraction level, but with heterogeneous sets of capabilities and languages. Rule Interchange Format (RIF [1]) and Production Rules Representation (PRR [2]) are two standard proposals respectively developed at the W3C and at the OMG aiming at providing a level of standardization to the domain. An important difference between RIF and PRR is that RIF is a standard for distribution of shared rules at runtime, whereas PRR is a standard for interchanging design-time rules.

This work presents a case study of bridging a business rule language [START_REF] Gottesdiener | Business Rules Show Power, Promise[END_REF]: Ilog Rule Language (IRL) to a standard format: Rule Interchange Format (RIF) and viceversa. We show how to make it possible with the help of Model Driven Engineering (MDE [3]).

RIF is the W3C Rule Interchange Format and it is part of the infrastructure composing the semantic web. It is an XML language for expressing rules that computers can execute. One of the main goals of RIF is to promote a standard format to interchange rules between existing rule languages.

Because of the serious tradeoffs in the design of rule language, RIF provides multiple versions, called dialects:

-Core: it is the fundamental RIF language. It is designed to be the common subset of most rule engines -BLD: it adds to Core dialect logic functions, equality in the then-part, and named arguments -PRD: adds a notion of forward-chaining rules, where a rule fires and then performs some action (adding, updating or retracting some information)

In this work we have implemented transformations for the Core and PRD dialects. An example of a RIF example is shown in (1). The primary software artifacts of the MDE approach are models, which are considered as firstclass entities. Every model defines a concrete syntax and conforms to a metamodel (or grammar, or schema), which defines its abstract syntax. One of the most common operations applied on models is transformation, which consists in the creation of a set of target models from a set of source models according to specific rules.

In the work presented in this paper three main tools have been used:

-Ecore allows to handle and create metamodels in Eclipse Modeling Framework (EMF) [START_REF]Eclipse Modeling Framework[END_REF] -TCS (Textual Concrete Syntax) ([8], [START_REF] Jouault | TCS: a DSL for the Specification of Textual Concrete Syntaxes in Model Engineering[END_REF]) is a bidirectional mapping tool between metamodels and grammars. It is able to perform both text-to-model (injections) and model-to-text (extractions) translations from a single specification. TCS is used to map context-free concrete syntaxes to metamodels. -ATL (Atlas Transformation Language) [START_REF] Jouault | ATL: a Model Transformation Tool[END_REF] is a model transformation language specified as both a metamodel and a textual concrete syntax. It allows developers to produce a number of target models from a set of source models by writing rules that define how to create target models from source model elements.

Transformations

The RIF2IRL transformation takes as input a RIF file, a Business Object Model (BOM [6]) and a Dictionary model. The output generates an IRL model and the related IRL file, generated by a TCS parser.

From the first input, using a TCS parser, a model conforming to an implementation of the RIF metamodel is extracted. A BOM file is passed to the transformation to resolve the domain information declared inside the RIF file/model. The related BOM model is extracted out of the BOM file by a TCS parser and it conforms to the BOM metamodel implemented in [START_REF] Del Fabro | Industrial-strength Rule Interoperability using Model Driven Engineering[END_REF]. A Dictionary is added as input parameter in the transformation in order to provide a translation from the element names declared in the RIF file to the corresponding ones appearing in the IRL file. The transformation has been implemented in ATL [START_REF]Atlas Transformation Language[END_REF] and it defines how to convert RIF to IRL.

The IRL2RIF transformation takes as input a IRL file and two BOMs and returns a RIF file. From the first input, a model conforming to the IRL metamodel is extracted by means of a TCS parser. The first BOM contains the information of the domain model, while the second one is used to map IRL functions to the built-in RIF functions for handling dates, numbers, strings, etc. The output is a RIF model, that is translated into a RIF file by a TCS parser.

The demonstration

The IRL2RIF and RIF2IRL transformations will be first demonstrated using a set of simple examples (Sect. A) that illustrate tha main functionalities of the transformation. We will then show a real world application of the transformation using data from the Arcelor usecase studied in the ONTORULE project [START_REF] Gonzalez-Moriyon | Final steel industry public demonstrators[END_REF]. In the context of the Rule Xchanger application depicted in the (3), we will focus on the transformation of IRL rule edited using JRules into RIF rules that will be executed using the Tight engine.

Conclusion

The transformation we intend to demonstrate is a key component in the interoperability between two rule languages that has been used in some real world applications such us the Arcelor usecase. it shows that rule application developed using commercial BRMS such as JRules can be exported to other platforms that support the RIF standard.

A Annex

Fig. 1 .

 1 Fig. 1. Example of a rule in RIF

Fig. 2 .

 2 Fig. 2. Example of a rule in IRL

Fig. 3 .

 3 Fig. 3. Rule Xchanger scenario

Fig. 4 .Fig. 5 .

 45 Fig. 4. Example 1


Iadine Chadès
email: iadine.chades@csiro.au

Josie Carwardine
email: josie.carwardine@csiro.au

Tara G Martin
email: tara.martin@csiro.au

Samuel Nicol
email: scnicol@alaska.edu

Régis Sabbadin
email: sabbadin@toulouse.inra.fr

Olivier Buffet
email: olivier.buffet@loria.fr

MOMDPs: a Solution for Modelling Adaptive Management Problems

In conservation biology and natural resource management, adaptive management is an iterative process of improving management by reducing uncertainty via monitoring. Adaptive management is the principal tool for conserving endangered species under global change, yet adaptive management problems suffer from a poor suite of solution methods. The common approach used to solve an adaptive management problem is to assume the system state is known and the system dynamics can be one of a set of pre-defined models. The solution method used is unsatisfactory, employing value iteration on a discretized belief MDP which restricts the study to very small problems. We show how to overcome this limitation by modelling an adaptive management problem as a restricted Mixed Observability MDP called hidden model MDP (hmMDP). We demonstrate how to simplify the value function, the backup operator and the belief update computation. We show that, although a simplified case of POMDPs, hm-MDPs are PSPACE-complete in the finite-horizon case. We illustrate the use of this model to manage a population of the threatened Gouldian finch, a bird species endemic to Northern Australia. Our simple modelling approach is an important step towards efficient algorithms for solving adaptive management problems.

Introduction

In conservation biology and natural resource management, adaptive management or 'learning by doing' is an iterative process of improving management by reducing uncertainty via monitoring. Coined by [START_REF] Walters | Ecological optimization and adaptive management[END_REF], adaptive management has gained notoriety as an approach to manage ecosystems to conserve biodiversity. The key elements of adaptive management include clear definition of objectives, specification of alternative models to achieve these objectives, implementation of two or more models in a comparative experimental framework, monitoring to assess the relative merits and limitations of alternative models, and iterative modification of management to determine the true model, if any [START_REF] Keith | Uncertainty and adaptive management for biodiversity conservation[END_REF]. Despite its virtues, there are few examples of the application of adaptive management in practice. Several factors have been proposed to explain the widespread implementation difficulties in adaptive manage-Copyright c 2012, Association for the Advancement of Artificial Intelligence (www.aaai.org). All rights reserved. ment programs. Amongst them is the inefficiency of methods used to solve adaptive management problems and the inability to experimentally manipulate endangered populations.

To date, adaptive management problems have been tackled using discretized belief MDP methods [START_REF] Williams | Markov decision processes in natural resources management: Observability and uncertainty[END_REF]. The belief space is discretized using p subintervals for each belief variable. The updating rule does not guarantee that the updated belief falls on one of these grid points and therefore an interpolation rule must be used to define the transition probabilities for the belief states. This technique has also been studied to sidestep the intractability of exact POMDP value iteration, using either a fixed grid [START_REF] Lovejoy | Computationally feasible bounds for partially observed Markov decision processes[END_REF][START_REF] Bonet | An e-optimal grid-based algorithm for partially observable Markov decision processes[END_REF] or a variable grid [START_REF] Brafman | A heuristic variable grid solution method for POMDPs[END_REF][START_REF] Zhou | An improved grid-based approximation algorithm for POMDPs[END_REF]. The grid based methods differ mainly in how the grid points are selected and what shape the interpolation function takes. In general, regular grids do not scale well in problems with high dimensionality and non-regular grids suffer from expensive interpolation routines.

In this paper, we propose a transparent and formal way of representing an adaptive management problem as a special case of POMDP. We demonstrate for the first time how to model and solve an adaptive management problem as a restricted Mixed Observability MDP (MOMDP) called hidden model MDP (hmMDP). Our approach benefits from recent developments in the field of robotics and decisionmaking under uncertainty [START_REF] Ong | Planning under uncertainty for robotic tasks with mixed observability[END_REF][START_REF] Araya-López | A closer look at MOMDPs[END_REF]. Our framework is particularly relevant in the case of endangered species where it may not be possible to undertake a replicated experiment to learn the true model.

Case Study

We illustrate our method on managing a population of a threatened bird species, the Gouldian finch. The most pervasive threats to wild populations of Gouldian finches are habitat loss and degradation caused by inappropriate fire and grazing regimes and introduced predators such as feral cats. The response of the population to different management actions is uncertain. Each of four experts provided a possible model, which are probability distributions describing how the population might respond to four alternative threat management actions. Our objective is to implement the management action that is most likely to lead to a high persistence probability for the Gouldian finch population. An optimal adaptive strategy will provide the best decision by determining which model (1 to 4) is most likely the real model, and hence which action is optimal, over time.

The next section presents background on POMDPs and their solution. Then we present the hmMDP model and study its complexity. The following section shows how the specificities of hmMDPs allow for simplified computations. Our simple case study is finally described and used in experiments to show the benefit of advanced solution techniques.

POMDPs

Partially Observable Markov Decision Processes (POMDPs) are a convenient model for solving sequential decisionmaking optimization problems when the decision-maker does not have complete information about the current state of the system. Formally, a discrete POMDP is specified as a tuple S, A, O, T, Z, r, H, b 0 , γ , where:

• S is the set of states s that might be partially observed or imperfectly detected by a manager;

• A is the set of actions (or decisions) a the manager needs to choose from at every time step;

• O is the set of observations o the manager perceives;

• T is a probabilistic transition function describing the stochastic dynamics of the system; an element T (s, a, s ′ ) represents the probability of being in state

s ′ at time t + 1 given (s, a) at time t, T (s, a, s ′ ) = p(s t+1 = s ′ |s t = s, a t = a); • Z is the observation function, with Z(a, s ′ , o ′ ) = p(o t+1 = o ′ |a t = a, s t+1 = s ′
) representing the conditional probability of a manager observing o ′ given that action a led to state s ′ ;

• r : S × A → ℜ is the reward function identifying the benefits and costs of being in a particular state and performing an action;

• H is the -finite or infinite-horizon (this section focuses on the infinite case);

• b 0 is an initial belief, a probability distribution over states;

• γ ∈ [0, 1] is a discount factor (may be 1 if H is finite).

The optimal decision at time t may depend on the complete history of past actions and observations. Because it is neither practical nor tractable to use the history of the actionobservation trajectory to compute or represent an optimal solution, belief states, i.e., probability distributions over states, are used to summarize and overcome the difficulties of imperfect detection. A POMDP can be cast into a fully observable Markov decision process defined over the (continuous) belief state space. In our case, solving a POMDP means finding a strategy (or policy) π : B → A mapping the current belief state (b ∈ B) to an allocation of resources. An optimal strategy maximizes the expected sum of discounted rewards over an infinite time horizon

E[ t γ t R(b t , a t )]
where b t and a t denote the belief state and action at time t, R(b, a) = s b(s)r(s, a). For a given belief state b and a given policy π this expected sum is also referred to as the value function V π (b). A value function allows us to rank strategies by assigning a real value to each belief b. An optimal strategy π * is a strategy such that,

∀ b ∈ B, ∀ π, V π * (b) ≥ V π (b)
. Several strategies can be optimal and share a same optimal value function V * which can be computed using the dynamic programming operator for a POMDP represented as a belief MDP [START_REF] Bellman | Dynamic Programming[END_REF], i.e., ∀b ∈ B:

V * (b) = max a∈A s∈S r(s, a)b(s) + γ o ′ p(o ′ |b, a)V * (b ao ′ ) ,
where b ao ′ is the updated belief given that action a was performed and o ′ is observed. This function can be computed recursively using Bellman's principle of optimality (Bellman 1957):

V n+1 (b) = max a∈A s∈S r(s, a)b(s) + γ o ′ p(o ′ |b, a)Vn(b ao ′ ) , (1)
where b ao ′ is updated using Bayes' rule:

b ao ′ (s ′ ) = p(o ′ |a, s ′ ) p(o ′ |b, a) s∈S p(s ′ |s, a)b(s). (2)
Sondik has shown that the finite time horizon value function is piecewise linear and convex (PWLC) and that the infinite time horizon value function can be approximated arbitrary closely by a PWLC function [START_REF] Sondik | The Optimal Control of Partially Observable Markov Decision Processes[END_REF]). An alternative way of representing V is to use vectors:

V (b) = max α∈Γ α • b, (3)
where Γ is a finite set of vectors called α-vectors, b is the belief represented as a finite vector, and α • b denotes the inner product of an α-vector and b. The gradient of the value function at b is given by the vector α b = arg max α∈Γ b • α.

The policy can be executed by evaluating (3) at b to find the best α-vector: π(b) = a(α b ). Exact methods like Incremental Pruning (IP) rely on regularly pruning dominated hyperplanes to reduce their number [START_REF] Cassandra | Exact and Approximate Algorithms for Partially Observable Markov Decision Processes[END_REF].

While various algorithms from the operations research and artificial intelligence literatures have been developed over the past years, exact resolution of POMDPs is intractable: finite-horizon POMDPs are PSPACE-complete [START_REF] Papadimitriou | The complexity of Markov decision processes[END_REF] and infinite-horizon POMDPs are undecidable [START_REF] Madani | On the undecidability of probabilistic planning and related stochastic optimization problems[END_REF].

In recent years, approximate methods have been developed successfully to solve large POMDPs. Amongst them, point based approaches approximate the value function by updating it only for some selected belief states [START_REF] Pineau | Point-based value iteration: An anytime algorithm for POMDPs[END_REF][START_REF] Spaan | Perseus: Randomized point-based value iteration for POMDPs[END_REF]. Typical point-based methods sample belief states by simulating interactions with the environment and then update the value function and its gradient over a selection of those sampled belief states.

From MOMDPs to hmMDPs

In this section we demonstrate how to model adaptive management problems using a specific factored POMDP model [START_REF] Boutilier | Computing optimal policies for partially observable decision processes using compact representations[END_REF]. A factored POMDP is a POMDP that explicitly represents the independence relationships between variables of the system.

MOMDPs

A Mixed Observability MDP (MOMDP) [START_REF] Ong | Planning under uncertainty for robotic tasks with mixed observability[END_REF]) is specified as a tuple X, Y, A, O, T x , T y , Z, r, H, b 0 , γ where:

• S = X × Y is the factored set of states of the system with X a random variable representing the completely observable components and Y a random variable representing the partially observable components. A pair (x, y) specifies the complete system state;

• A is the finite set of actions; • Z is the observation function with

• O = O x × O y is set of observations with O x = X
Z(a, x ′ , y ′ , o ′ x , o ′ y ) = p(o ′ x , o ′ y |a, x ′ , y ′ ) the probability of observing o ′ x , o ′ y
given action a was performed, leading to system state (x ′ , y ′ ). In a MOMDP we assume the variable X ′ is perfectly observable so we have p

(o ′ x |a, x ′ , y ′ , o ′ y ) = 1 if o ′ x = x ′ ,

Hidden Model MDPs

Our adaptive management problem assumes managers can perfectly observe the state of the studied system but are uncertain about its dynamics. In ecology, this problem has traditionally been solved by assuming that the real but unknown MDP model is one of a finite set of known models. This can be treated as a hidden model MDP (hmMDP), i.e., a MOMDP where the partially observable state variable corresponds to the hidden model. In this setting, the following assumptions can be incorporated to simplify the solution:

1. The real model of the dynamics of the system y r is an element of a finite set Y of predefined models;

2. The finite set of actions A effects on the completely observable variable X and has no effect on Y ; Y is also independent from X, i.e., p(y ′ |x, y, a, x ′ ) = p(y ′ |y);

(A) x = x ′ and 0 otherwise. Figure 1 illustrates the comparison of the standard POMDP and hmMDP models.

o t o t-1 a t-1 s t-1 s t (B) x t-1 x t y t o t o t-1 y t-1 s t-1 s t a t-1

Definition of an hmMDP

An hmMDP problem is thus a simplified version of MOMDP. In fact, it can be seen as a finite set of (completely observed) MDP models, differing only in their transition and reward functions. Namely, an hmMDP is a tuple X, Y, A, {p y } y∈Y , {r y } y∈Y , H, γ, b 0 , where Y is the finite set of potential models, p y (x ′ |x, a) = p(x ′ |x, a, y) and r y (x, a) = r(x, a, y). Here, b 0 is our initial belief on which model is the true one.

When solving an hmMDP problem, we are looking for a policy π(x, b) → a maximizing the usual value function:

V π (x, b) = E H t=0 γ t r(X t , π(X t , B t ), B t )|b 0 , π , = x 1 ...x H ,y b 0 (y)P y (x 1 . . . x H |x 0 = x) H t=0 γ t r y (x t , a t ),
where the X t are random variables over X and the B t are random vectors over possible beliefs.

Complexity of hmMDPs

Although hmMDPs are simplified MOMDPs, they are still computationally complex. We show that solving an hmMDP in the finite-horizon case is a PSPACE-complete problem, thus as hard to solve as a classical finite-horizon POMDP. Proposition 1. Deciding whether a finite-horizon hmMDP problem admits a solution policy of value greater than a predefined threshold is a PSPACE-complete problem.

Finite-horizon hmMDP clearly belongs to PSPACE, as a particular case of POMDP. The hardness proof results from a reduction of the PSPACE-hard State Disambiguation problem (SD) [START_REF] Conitzer | Definition and complexity of some basic metareasoning problems[END_REF] to the hmMDP problem.

Definition 1 (STATE-DISAMBIGUATION). We are given:

• A set Θ = {θ 1 , . . . , θ n } of possible states of the world and a uniform probability distribution p over Θ. • A utility function u : Θ → [0; +∞[. u(θ i ) is the utility of knowing for sure that the state of the world is θ i . • A set Q = {q 1 , . . . , q m } of queries. q j = {q j1 , . . . , q jpj } is a set of subsets of Θ, such that 1≤k≤pj q jk = Θ. If the true state of the world is θ i and q j is asked, an answer is chosen (uniformly) randomly among the answers q jk containing θ i . • A maximum number N of queries that can be asked and a target real value G > 0.

The STATE DISAMBIGUATION problem consists in deciding whether there exists a policy asking at most N queries that gives expected utility at least G. If π δ (θ i ) denotes the probability of identifying θ i by using policy δ, the SD problem amounts to deciding whether there exists δ such that

v(δ) = 1≤i≤n p(θ i )π δ (θ i )u(θ i ) ≥ G,
where p is the uniform distribution over Θ.

Proof. The following transformation is a reduction from SD to hmMDP:

• X = {q jkj } j=1...m,kj =1...pj , Y = Θ. • A = Q ∪ {{θ 1 }, . . . , {θ n }}. • If a ∈ Q, p(x ′ |x, a, y) = p(q rk ′ |q jk , q r , θ i ) = p(q rk ′ |q r , θ i ) = 1 |qr| if θ i ∈ q rk ′ , and 0 else. If a = {θ i }, p(x|x, {θ i }, θ i ) = 1, ∀x, θ i .
• < N, r t (x, y, a) = r t (q jk , θ i , q l ) = 0, ∀q jk , θ i , q l and r t (x, y, a) = r t (q jk , θ i , {θ j }) < 0 (this, to prevent choosing a ∈ Θ for t < N ). • r N (x, y, a) = r N (q jk , θ i , {θ i }) = u(θ i ), ∀q jk , θ i and r N (q jk , i , a) = -max θj u(θ j ), ∀a = {θ i } (this ensures that, if a policy does not disambiguate θ i , it has expected value less than or equal to 0).

• H = N, b 0 is a uniform distribution over Y .
• value threshold G (identical to SD threshold).

One can easily check that the policy spaces of the SD and hmMDP problems are isomorphic and there exists a policy δ such that v(δ) ≥ G in the SD problem if and only if there is a policy π which has value at least G in the corresponding hmMDP.

Solution Methods

We have provided a new framework to model adaptive management problems in ecology. This framework is simpler than the MOMDP framework and sufficient, as we will show in the application section. However, we have also shown that the computational complexity of this new problem is the same as general POMDPs. In this section we demonstrate that solution algorithms, such as the MO-SARSOP algorithm, can benefit from the model simplification to gain computational efficiency, even though exact resolution remains time-exponential.

Building on [START_REF] Ong | Planning under uncertainty for robotic tasks with mixed observability[END_REF] and Araya-López et al. ( 2010), we show how to simplify the calculation of the value function (1) and belief update (2) for existing POMDP algorithms. First, [START_REF] Ong | Planning under uncertainty for robotic tasks with mixed observability[END_REF] have shown that V (b) = V (x, b y ) for any MOMDP, and Eq. ( 1) can be rewritten:

V n+1 (x, b y ) = max a∈A [ y∈Y r y (x, a)b y (y) + γ y,x ′ ,y ′ ,o ′ b y (y) p(x ′ |x, y, a, y ′ )p(y ′ |y, a)p(o ′ |a, x ′ , y ′ )V n (x ′ , b ao ′ y )]. (4)
Accounting for the stationarity assumptions (2-3) and p(o ′ |x ′ , y ′ ) = 1, o ′ = x ′ , (assumption 4) leads to:

V n+1 (x, b y ) = max a∈A [ y∈Y r y (x, a)b y (y)+ γ y,x ′ b y (y)p(x ′ |x, y, a)V n (x ′ , b a,o ′ =x ′ y )],
which can be rewritten substituting

V n (x ′ , b a,o ′ =x ′ y ) = max α∈Γn b a,o ′ =x ′ y • α: V n+1 (x, b y ) = max a∈A [ y∈Y r y (x, a)b y (y)+ (5) γ y,x ′ b y (y)p(x ′ |x, y, a) max α∈Γn b a,o ′ =x ′ y • α] with b a,o ′ =x ′ y
the updated belief given a, o ′ and b y simplified as follows:

b a,o ′ =x ′ y (y) = p(x ′ |x, y, a)b y (y) y ′′ b y (y ′′ )p(x ′ |x, y ′′ , a) . (6)
Substituting ( 6) in max α∈Γn b a,o ′ =x ′ y

• α we can rearrange (5):

V n+1 (x, b y ) = max a∈A [ y∈Y r y (x, a)b y (y) + γ y,x ′ b y (y)× p(x ′ |x, y, a) max α∈Γn y ′ p(x ′ |x, y ′ , a)b y (y ′ )α(y ′ ) y ′′ b y (y ′′ )p(x ′ |x, y ′′ , a) ] = max a [ y∈Y r y (x, a)b y (y)+ γ x ′ max α∈Γn y ′ p(x ′ |x, y ′ , a)b y (y ′ )α(y ′ )] = max a [b y • α a 0 + γ x ′ max α∈Γn b y • G ax ′ α (x)], (7)
where

G ax ′ α (x) = y ′ p(x ′ |x, y ′ , a)α(y ′
) and the reward function r(s, a) is represented as a set of |A| vectors α a 0 = (α a 0 (1), ..., α a 0 (|S|)), one for each action a, α a 0 (s) = r(s, a). Using the identity max yj x • y j = x • arg max yj x • y j twice we get:

V n+1 (x, b y ) = b y • arg max a b y • [α a 0 + γ x ′ arg max α∈Γn b y • G ax ′ α (x)].
Algorithm 1: MO-SARSOP α-vector backup at a node

(x, b y ) of T R BACKUP(T R , Γ, (x, b y )); 1 forall a ∈ A, x ′ ∈ X, o ′ ∈ O do 2 α a,x ′ ,o ′ ← arg max α∈Γy(x ′ ) (α • τ (x, b y , a, x ′ , o ′ )); 3 forall a ∈ A, y ∈ Y, do 4 α a (y) ← r(x, y, a) + 5 γ x ′ ,o ′ ,y ′ T x (x, y, a, x ′ )T y (x, y, a, x ′ , y ′ ) × Z(x ′ , y ′ , a, o ′ )α a,x ′ ,o ′ (y ′ ); α ′ ← arg max a∈A (α a • b y ); 6 Insert α ′ into Γ y (x); 7 Algorithm 2: hm-SARSOP α-vector backup at a node (x, b y ) of T R BACKUP(T R , Γ, (x, b y )); 1 forall a ∈ A, x ′ ∈ X do 2 α a,x ′ ← arg max α∈Γy(x ′ ) (α • τ (x, b y , a, x ′ )); 3 forall a ∈ A, y ∈ Y, do 4 α a (y) ← r y (x, a) + γ x ′ T x (x, y, a, x ′ )α a,x ′ (y); 5 α ′ ← arg max a∈A (α a • b y ); 6 Insert α ′ into Γ y (x); 7
Finally we can define the vector backup β(b y ) as the vector whose inner product with b y yields V n+1 (b y ):

V n+1 (x, b y ) = b y • β(b y ), where (8)
β(b y ) = arg max a b y • [α a 0 + γ x ′ arg max α∈Γn b y • G ax ′ α (x)].
By simplifying the backup operator calculation (8) and the belief update operation (6) we have shown how existing POMDP solving algorithms based on a PWLC approximation can be adapted when looking at problems with hidden models. These two procedures are common to most exact algorithms such as Witness or Incremental Pruning (Cassandra 1998), and most point-based POMDP algorithms, such as PBVI [START_REF] Pineau | Point-based value iteration: An anytime algorithm for POMDPs[END_REF], Perseus [START_REF] Spaan | Perseus: Randomized point-based value iteration for POMDPs[END_REF], symbolic Perseus [START_REF] Poupart | Exploiting Structure to Efficiently Solve Large Scale Partially Observable Markov Decision Processes[END_REF] and SARSOP [START_REF] Ong | Planning under uncertainty for robotic tasks with mixed observability[END_REF].

Algorithm 1 presents the backup procedure used at each belief point by MO-SARSOP, and Algorithm 2 shows how-under the adaptive management assumptions-it can be simplified at two levels: i) we do not need to consider the set of observations O, and the belief update calculation τ is simplified (line 2); ii) the observation function Z and the model dynamics T y are not required (line 5). The sampling procedure also benefits from the simplified belief update (6).

Management of a Threatened Bird

We apply our adaptive management model to a threatened Gouldian finch population in the Kimberley, Australia. Our management objective is to maximize the likelihood of a high persistence probability of this population over the medium term, in an area where we can measure the population state but are uncertain about the state response to the management actions. We asked four experts to assess the likelihood of a high (and conversely low) probability of persistence under four plausible management actions. Following our hidden model MDP model, we define the set of states S = X × Y where X represents the local probability of persistence X = {Low, High}. We pose Y = {Expert 1, Expert 2, Expert 3, Expert 4} the set of possible true models that predict the state dynamics of our studied species. The set of actions A = {DN,FG,C,N} represents the management actions to choose from at every time step: do nothing (DN), improve fire and grazing management (FG), control feral cats (C) and provide nesting boxes (N). We assume that the same amount of funds are spent on each action. We elicited from the experts the transition probabilities T x for each model and state transition. Finally we define the reward function so that r(Low,DN) = 0, r(High,DN) = 20, r(Low,{FG,C,N}) = -5 and r(High,{FG,C,N}) = 15. We assume that each expert is as likely to be correct at the beginning of our adaptive management program. We first attempted to solve our adaptive management problem using IP1 . Due to a lack of memory space we were not able to provide an exact solution for a time horizon greater than 5 time steps with 3753 α-vectors computed (Problem 1, Table 1). We solved the same problem using MO-IP (Araya-López et al. 2010) and reached the 5th horizon 3 times faster than IP. We then computed approximate solutions using the algorithms Grid-bMDP 1 , SARSOP and MO-SARSOP [START_REF] Ong | Planning under uncertainty for robotic tasks with mixed observability[END_REF]. We interrupted Grid-bMDP after ∼30 mins (horizon 26, 4402 α-vectors and an estimated error of 1.28). Using the same computational time, SARSOP and MO-SARSOP provided the best performances by far with a higher precision for MO-SARSOP. Experiments conducted with a version of MO-IP implementing the modifications specific to adaptive management showed no speed up at horizon 5 (106.69s vs. 106.80s for the original MO-IP). We also ran these algorithms on a version of our Gouldian finch problem with a larger set of states (|X| = 81) and a smaller set of experts (Experts 1 and 2). In this case we explicitly distinguish how the interactions between 4 species (dingo, cats, long-tailed finch and Gouldian finch) could drive the management success (Problem 2, Table 1).

The solutions of Problem 1 provide guidance on which management action should be performed in absence of precise knowledge of what the real model is. Expert 1 and 4 provided similar transition probabilities; both believe that the appropriate management of fire and grazing (FG) will generate a higher probability of the population reaching and remaining in the 'High' persistence state. Consequently, in absence of prior information on what the real model is, the first action to perform is FG. Expert 2 believes the management of feral cats (C) and the provision of nesting boxes (N) have more benefits than fire and grazing management. Expert 3 believes that habitat loss and intra-species competition is the limiting factor of persistence and therefore attributes the highest probabilities of persistence when action 'provide nesting boxes' (N) is implemented. Note that Expert 3 values the management of cats and fire and grazing as 'good' actions. We assessed the quality of our adaptive strategy by simulations. When we assume that cats are the main threat and Expert 2 holds the real model, our adaptive strategy quickly finds the best action to perform. However when we assume that fire and grazing are the main threats, our adaptive strategy provides the best action to perform but has trouble identifying the real model, because Experts 1 and 4 are similarly supportive of the combined FG action. When we assume Expert Model 3 holds the real model (see Figure 2), the simulations first lead to believing that the management of cats is most beneficial and Expert 2 holds the real model. Again this ambiguity is due to the high success rate given by Expert 3 when feral cats are managed. After several time steps, our adaptive strategy favors Expert 3.

Discussion

Adaptive management problems as currently solved in the ecology literature suffer from a lack of efficient solution methods. We demonstrate and explain how to model a classic adaptive management problem as a POMDP with mixed observability focusing on the special case of hidden model MDPs (hmMDPs). We show that hmMDPs are PSPACEcomplete in the finite-horizon case. However, if data were available, MO-SARSOP could approximately solve AM problems with up to ∼250,000 states in 20-30 minutes [START_REF] Ong | Planning under uncertainty for robotic tasks with mixed observability[END_REF]. The assumption that the real model is contained within the model set is simplistic (Assumption 1), however it is currently the way AM is solved. Point based methods help us account for a large set of models so that we do not risk being too far from the real model, but, when managing threatened species, having many models makes it difficult to be confident on the real model as observations are few. In this situation, we would recommend incorporating a model for each of the management options so that the best action can be chosen in the face of uncertainty.

While our examples may seem simplistic in AI they are complex and reasonably realistic in ecology (Chadès et al. 2008). Practitioners are reluctant to adopt new methods. Our aim is to bridge the gap between AI and ecological management by solving real, tractable problems at a scale defined by practitioners. Conservation practitioners need to understand the model and the management rules need to be simple enough to be implemented (Chadès et al. 2011). Moreover, we need to be able to populate these models with data, which in ecology is often lacking. Both examples are based on a real application with real experts, with only 4 management actions currently possible for the Gouldian finch. Example 2 is complex due to interacting species and requires eliciting more information to populate the corresponding DBN.

POMDPs are a powerful model that does not require the restrictive assumptions of the traditional adaptive management literature and offers new and exciting adaptive management applications. First, there is no need to assume a static model (Assumption 3). If we allow the real model to change over time we can tackle systems that are influenced by climate change or seasonal fluctuations. Second, there is no need to assume perfect detection of the state (X) (Assumption 4) and detection probabilities can be included. Not assuming that the real model is one of a predefined set of models (Assumption 1) remains an unsolved challenge. To our knowledge this cannot be tackled with current POMDP solutions but model-based Bayesian Reinforcement Learning [START_REF] Poupart | An analytic solution to discrete Bayesian reinforcement learning[END_REF]) could be a relevant alternative.

  and 0 otherwise; • r, H, b 0 and γ are defined as for POMDPs.The belief space B now represents our beliefs on y only since the state variable X is fully observable. Any belief b in B on the complete system state s = (x, y) is represented as (x, b y ). In algorithms based on PWLC approximations, accounting for visible state variables allows for substantial speed-ups by reasoning on multiple low-dimensional belief spaces instead of the original high-dimensional one, as demonstrated with MO-SARSOP[START_REF] Ong | Planning under uncertainty for robotic tasks with mixed observability[END_REF]-based on SARSOP, a state of the art point-based solver-and MO-IP (Araya-López et al. 2010)-based on IP.

Figure 1 :

 1 Figure 1: Standard POMDP (A) vs. hidden model MDP (B)

  Figure 2: Simulation of our adaptive management strategy when Expert 3 holds the real model. The graphs represent the state dynamics (top), the belief in each expert's predictions over time (middle), and the action performed (bottom).

  the completely observable component, and O y the set of observations of the hidden variables;• T x (x, y, a, x ′ , y ′ ) = p(x ′ |x, y, a) gives the probability that the fully observable state variable takes the value x ′ at time t + 1 if action a is performed in state (x, y) at time t and has already led to y ′ ; T y (x, y, a, x ′ , y

′ ) = p(y ′ |x, y, a, x ′ ) gives the probability that the value of the partially observable state variable changes from y to y ′ given that action a is performed in state (x, y) and the fully observable state variable has value x ′ ;

Table 1 :

 1 Performance for the Gouldian finch problems. Experiments conducted on a 2.40 GHz Core 2 computer with 3.45 GB of memory. ( + ) err represents the Bellman residual error between 2 successive value functions (precision).

		Problem 1: |S| = 8 Problem 2: |S| = 162
		|X| = 2, |Y | = 4	|X| = 81, |Y | = 2
	IP	h= 5 |α| = 3753	h=4 |α| = 1181
		t=349.2s	t=703.7s
	MO-IP	h=5 |α| = 2052	h=4 |α| = 218
		t=106.8s	t=0.59s
	Grid	h=26 |α| = 4402	h=23 |α| = 426
	bMDP +	t=1831s err=1.28	t=1849s err= 3.84
	SARSOP +	|α| = 25153	|α| = 6995
		t=1831s err=0.066	t=422.85s err<0.001
	MO-	|α| = 38137	|α| = 3861
	SARSOP + t=1831s err=0.055	t=19.36s err<0.001

Using Cassandra's pomdp-solver toolbox.

Acknowledgements

This research was funded through an Australian Government's NERP, an Australian Research Council's Centre of Excellence (I.C., J.C., T.G.M.) and a CSIRO Julius Career Award (T.G.M.).