Keywords: AMS Classification 60 J 25, 60 Bessel processes, Pseudo-inverses, Generalized Black-Scholes formulae

It is shown that the tail probability of a Bessel process is the distributio function of a random time which is related to first and last passage times of Bessel processes.

1 Introduction. Definition of a pseudo-inverse 1.0 The aim of this paper is to give some mathematical meaning to the assertion : "the stochastic process (R t , t ≥ 0), with values in R + , has a tendency to increase" (1.0) and, more precisely, to "measure this tendency" when (R t , t ≥ 0) is a Bessel process, by showing the existence of a "pseudo-inverse" of R and by studying its properties. First, here are several possible interpretations of (1.0) : (1.a) (R t , t ≥ 0) is a.s. increasing ; then, it admits an increasing inverse process.

(1.b) (R t , t ≥ 0) is a submartingale ; under some adequate conditions, it admits a Doob-Meyer decomposition :

R t = M t + A t (t ≥ 0)
with (A t , t ≥ 0) an increasing process and (M t , t ≥ 0) a martingale. Since a martingale is a "well balanced" process, (R t , t ≥ 0) has a tendency to increase, which is measured by (A t , t ≥ 0) .

(1.c) E(R t ), t ≥ 0 is an increasing function.

(1.d) (R t , t ≥ 0) is "stochastically increasing", i.e. : for every y ≥ 0, P (R t ≥ y) is an increasing function of t (in the case R 0 = x > 0, we need to modify the previous assertion by taking y ≥ x, to get a meaningful assertion). We note that, trivially : (1.a) =⇒ (1.b) =⇒ (1.c) ⇐= (1.d) (when x = 0). It is the assertion (1.d) which we have in mind throughout this study. It leads us to the definition of a pseudo-inverse of (R t , t ≥ 0), which we now present.

1.1 Let Ω, (R t , t ≥ 0), P x (x ∈ R + ) denote a process taking values on R + , which is a.s.

continuous and such that P x (R 0 = x) = 1. Definition 1.1 (R t , t ≥ 0) admits an increasing pseudo-inverse (resp. decreasing pseudoinverse) if : i) for every y > x,

lim t→∞ P x (R t ≥ y) = 1 (1.1)
ii) for every y > x, the application from R + into [0, 1] : t → P x (R t ≥ y) is increasing (1.2) resp. if : i') for every y < x,

lim t→∞ P x (R t ≤ y) = 1 (1.1 ′ )
ii') for every y < x, the application from R + into [0, 1] :

t → P x (R t ≤ y) is increasing (1.2 ′ )
Definition 1.2 Assume that (1.1) and (1.2) resp. (1.1') and (1.2') are satisfied. Then, there exists a family of positive r.v.'s (Y x,y , y > x) (resp. Y x,y , x > y) such that : P x (R t ≥ y) = P (Y x,y ≤ t) (t ≥ 0) (1.3) resp. P x (R t ≤ y) = P (Y x,y ≤ t) (t ≥ 0) (1.3 ′ ) We call the family of positive r.v.'s (Y x,y , y > x) resp. (Y x,y , y < x) the increasing (resp. decreasing) pseudo-inverse of the process (R t , t ≥ 0). We note that, a priori, the family (Y x,y , y > x) resp. (Y x,y , y < x) is only a family of r.v.'s, and does not constitute a process. Note that formulae (1.3) and (1.3') may be read, either considering x as fixed and y varying or by considering the two parameters x and y as varying.

Remark 1.3 i) If (R t , t ≥ 0) admits an increasing (resp. decreasing) pseudo-inverse, then R t -→ t→∞ +∞ in probability resp. R t -→ t→∞ 0 in probability .

ii) If (R t , t ≥ 0) admits an increasing pseudo-inverse, then, for every α > 0 (resp. for every α < 0), the process (R α t , t ≥ 0) admits an increasing (resp. decreasing) pseudo-inverse. If (R t , t ≥ 0) admits a decreasing pseudo-inverse, for every α > 0 (resp. for every α < 0), the process (R α t , t ≥ 0) admits a decreasing (resp. increasing) pseudo-inverse. iii) Justifying the term "pseudo-inverse" Condition (1.2) indicates that the process (R t , t ≥ 0) "has a tendency to increase". In case this process is indeed increasing, we introduce (τ l , l ≥ 0) its right-continuous inverse :

τ l := inf{t ≥ 0 ; R t > l}
Then, from (1.3), and for y > x ; P x (R t ≥ y) = P x (τ y ≤ t) = P (Y x,y ≤ t).

(1.4)

Thus Y x,y (law)
= τ y under P x . Hence, in this case, we may choose for the family (Y x,y , y > x) the process (τ y , y > x), and this justifies our terminology of pseudo-inverse. iv) A simple example Let (B t , t ≥ 0) be a Brownian motion started at 0, and a > 0. Define : R t := exp(B t + at). Then, the submartingale (R t , t ≥ 0) (which satisfies R 0 = 1) admits an increasing pseudoinverse. Indeed, for every y > 1 :

P (R t ≥ y) = P (B t + at ≥ log y) = P ( √ t B 1 ≥ log y -at) (by scaling) = P B 1 ≥ log y √ t -a √ t = 1 √ 2π ∞ log y √ t -a √ t e -x 2 2 dx
Obviously, P (R 0 ≥ y) = 0 and lim t→∞ P (R t ≥ y) = 1. Moreover :

∂ ∂t P (R t ≥ y) = 1 2 log y √ 2πt 3 + a √ 2πt e -1 2t (log y-at) 2 ≥ 0
Thus, with (Y z , z > 0) defined as or, rather, satisfying :

P (Y log y ≤ t) = P (R t ≥ y)
we have :

P (Y log y ∈ dt) = 1 2 log y √ 2πt 3 + a √ 2πt exp - 1 2t (log y -at) 2 dt = 1 2 P (T (
log y = inf{t ≥ 0 ; B t + at = log y} since as is well-known see e.g. [MRY] P (G 1.2 Pseudo-inverses and generalised Black-Scholes formulae

Our main reason to introduce this notion of a pseudo-inverse is, independently from the fact that we think it is of interest by itself, because it plays an important role in the set-up of some generalised Black-Scholes formulae. More precisely, let us denote by (M t , t ≥ 0 ; F t , t ≥ 0 ; P x , x ≥ 0) a Markov process taking values in R + , such that (M t , t ≥ 0) is a (F t , t ≥ 0, P x ) positive local martingale, hence a supermartingale, which satisfies P x (M 0 = x) = 1 and lim t→∞ M t = 0 P x a.s.

(1.5)

Let Π M : R + × R + -→ R + denote the put quantity associated with (M t , t ≥ 0), that is :

Π M (K, t) := E x (K -M t ) + (K, t ≥ 0) (1.6)
Then, it follows from our assumptions that, for every x > 0 i)

1 x Π M (y, ∞), 0 ≤ y ≤ x = y x , 0 ≤ y ≤ x is the distribution function of a uniform r.v. U x on [0, x] ii) 1 x Π M (x, t), t ≥ 0 is the distribution function of a positive r.v., since 1 x Π M (x, 0) = 0, 1 x Π M (x, ∞) = 1 and the application t → 1 x Π M (x, t) = 1 x E x (x -M t ) + is an increasing function of t, (x -M t ) + being a submartingale since (x -y) + is a convex function of y .
More precisely, it is proven in [MRY] that :

1 x Π M (x, t) = P x (G x ≤ t) (1.6)
with :

G x := sup{t ≥ 0 ; M t = x} (1.7) (= 0 if this set is empty) A natural question is then : x being fixed, is the function of K and t 1 x Π M (K, t), K ≤ x ; t < ∞ , the distribution function of a probability on [0, x] × [0, ∞], which
we denote, if it exists, by γ M (x) ? If so, one has

1 x Π M (K, t) = γ M (x) [0, K] × [0, t] (K ≤ x, t ≥ 0) (1.8)
On the other hand, note that, by Fubini, one has :

1 x Π M (K, t) = 1 x K 0 P x (M t ≤ y)dy (1.9)
Hence :

∂ ∂K 1 x Π M (K, t) = 1 x P x (M t ≤ K) (K ≤ x) (1.10)
Thus, the existence of the probability γ M (x) is equivalent to :

t → P x (M t ≤ K) is an increasing function of t (for K < x) (for K < x).
Clearly, this condition is equivalent to the existence of a decreasing pseudoinverse (Y x,y , y < x) for the process (M t , t ≥ 0 ; P x , x ≥ 0). At this point, we have obtained :

Proposition 1.4 For every x > 0, there exists a probability γ

M (x) on [0, x] × [0, ∞[ such that : 1 x Π M (K, t) = γ M (x) [0, K] × [0, t] (K < x, t ≥ 0) (1.11)
if and only if (M t , t ≥ 0 ; P x , x ≥ 0) admits a decreasing pseudo-inverse (Y x,y , y < x), i.e. if, for every 0 < y < x and t ≥ 0, one has :

P x (M t ≤ y) = P (Y x,y ≤ t) (1.12) Of course, if (1.11) is satisfied, then γ M (x)
is the law of a couple (U x , Λ x ) where U x is uniform on [0, x] and Λ x (law)

= G x , but U x and Λ x are dependent. Moreover, if Π M (x) is a regular function of K and t, then, denoting by f γ M (x) the density of γ M (x) :

f γ M (x) (K, t) = ∂ 2 ∂K∂t 1 x Π M (K, t) = 1 x ∂ ∂t P x (M t ≤ y) (1.13) from (1.10) = 1 x ∂ ∂t P (Y x,y ≤ t) = 1 x f Yx,y (t) (K ≤ x) (1.14)
where f Yx,y denotes the density of Y x,y .

Our aim in this paper

In the set-up of the classical Black-Scholes formulae, i.e. when

M t = E t := exp B t - t 2 ,
where under P x , (B t , t ≥ 0) is a Brownian motion starting from (log x), we have proven in [MRY] the existence of a decreasing pseudo-inverse of (M t , t ≥ 0) and thus established the existence of a probability γ E (x) characterized by :

1 x E x (K -E t ) + = γ E (x) [0, K] × [0, t] (K ≤ x, t ≥ 0) (1.15)
In fact, in [MRY], the study was done for x = 1, but reducing the study to the case x = 1 is easy and we have also described this probability γ E (x) in detail. Now, let Ω, (R t , F t ), t ≥ 0 ; P (ν)

x , x ≥ 0 denote a Bessel process of index ν (ν ≥ -1). We shall prove in Section 2 (see Theorem 2.1) the following results :

• when ν ≥ - 1 2
, the existence of an increasing pseudo-inverse (Y

(ν)
x,y , y > x) for (R t , t ≥ 0, P (ν)

x , x ≥ 0)

• when ν = -1, the existence of a decreasing pseudo-inverse (Y

(-1) x,y , y < x) for (R t , t ≥ 0, P (-1) x , x ≥ 0) • when ν ∈ -1, - 1 2 :
i) if x = 0, the existence of an increasing pseudo-inverse ii) if x > 0, the non-existence of a pseudo-inverse.

In Section 3, we describe the laws of the r.v.'s (Y

x,y , y > x). In particular, we obtain explicitely the Laplace transform of (Y (ν) x,y , y > x) :

E(e -λY (ν) x,y ) = I ν (x √ 2λ) (x √ 2λ) ν (y √ 2λ) ν+1 K ν+1 (y √ 2λ) ν ≥ - 1 2 (1.16)
where I ν and K ν are the classical modified Bessel functions see [Leb], p. 108 . We also answer (partially) the question : are the r.v.'s Y (ν)

x,y infinitely divisible ? For ν > 0, we apply (end of Section 2) the existence of a pseudo-inverse for the process (R t , t ≥ 0) to the generalised Black-Scholes formula relative to the local martingale (R -2ν t , t ≥ 0), when (R t , t ≥ 0) denotes a δ = 2(ν + 1) dimensional Bessel process.

In Section 4, we present two other families of Markovian positive submartingales which admit an increasing pseudo-inverse ; they are : i) the Bessel processes with drift, studied by S. Watanabe [Wata] with infinitesimal generator :

L (ν,a) = 1 2 d 2 dx 2 + 2ν + 1 2x + a I ν+1 I ν (ax) d dx for ν ≥ - 1 2
, and a ≥ 0 ;

ii) the generalized squares of Ornstein-Uhlenbeck processes see, e.g., [START_REF] Pitman | Yor A decomposition of Bessel bridges[END_REF] , also called Cox-Ingersoll-Ross ( : CIR) processes, with infinitesimal generator :

L (ν,β) = 2x d 2 dx 2 + 2βx + 2(ν + 1) d dx ,
again, with ν ≥ -1 2 , and β ≥ 0.

Of course, letting a → 0, resp. : β → 0, in i), resp. ii), we recover the result for the Bessel processes.

In Section 5, we exhibit a two parameter family (Y (ν,α) , ν ≥ 0, 0 ≤ α ≤ 1) of variables (indexed by (x, y), x < y) which extends the family Y (ν) , corresponding to α = 1.

The following paragraph aims at relating the results of Section 3, which presents the properties of the r.v.'s Y (ν)

x,y , and those of Section 5, where the r.v.'s Y

(ν,α)

x,y are defined and studied. Consider (R

(ν) t , t ≥ 0) and (R (ν ′ )
t , t ≥ 0) two independent Bessel processes starting from 0, and with respective indices ν and ν ′ , greater than -1, or dimensions : d = 2(ν + 1) and d ′ = 2(ν ′ + 1). From the additivity property of squares of Bessel processes, the process :

R (ν+ν ′ +1) t := (R (ν) t ) 2 + (R (ν ′ ) t ) 2 1 2 , t ≥ 0,
is a Bessel process with index ν + ν ′ + 1 (or dimension : d + d ′ ), starting from 0. Let, for y > 0 :

G (α) y := sup{t ≥ 0 : R (α) t = y} T (α) y := inf{t ≥ 0 : R (α) t = y} (α = ν, ν + ν ′ + 1)
It is then clear that :

T (ν+ν ′ +1) y ≤ T (ν) y ≤ G (ν) y , and
G (ν+ν ′ +1) y ≤ G (ν) y .
These inequalities invite to look for some r.v.'s Z (i)

y (i = 1, 2, 3), such that :

G (ν) y (law) = G (ν+ν ′ +1) y + Z (1)

y

(1.17)

G (ν) y (law) = T (ν) y + Z (2) y (1.18) T (ν) y (law) = T (ν+ν ′ +1) y + Z (3) y (1.19)
where the r.v.'s featured on the right-hand side of (1.17), (1.18) and (1.19) are independent. The existence of these r.v.'s

Z (i)
y , (i = 1, 2, 3) is obtained in Sections 3 and 5 as a sub product of the existence of pseudo-inverses for Bessel processes. Precisely, the identities (1.17), (1.18) and (1.19) are shown respectively as (5.12), in Proposition 5.1, resp. as (5.8), with y = z and z = 0, resp. as (5.11) in Proposition 5.1. In Theorem 3.1 point iv) and v) and Theorem 5.2, other equalities in law of the type of (1.17), (1.18) and (1.19), are established, when the starting points of the Bessel processes differ from 0.

Section 6 concludes with a small list of the ingredients we have used throughout the paper to obtain these pseudo-inverses.

2 Existence of pseudo-inverses for Bessel processes 2.1 Statement of our main result Let ν ≥ -1 and (R t , t ≥ 0 ; P (ν)

x , x ≥ 0) a Bessel process with index ν, i.e. with dimension δ = 2ν + 2, with δ ≥ 0.

Theorem 2.1 i) If ν ≥ - 1 2
, (R t , t ≥ 0) admits an increasing pseudo-inverse ; that is for every y > x > 0, there exists a positive r.v. Y

(ν)

x,y such that, for every t ≥ 0 : P (ν) x (R t ≥ y) = P (Y (ν) x,y ≤ t).

(2.1)

ii) Let ν ∈ -1, - 1 2 ; if x > 0, (R t , t ≥ 0)
does not admit a pseudo-inverse. If x = 0, (R t , t ≥ 0) admits an increasing pseudo-inverse (see point 2.3.2 of the proof of Theorem 2.1).

iii) If ν = -1, (R t , t ≥ 0) admits a decreasing pseudo-inverse, that is there exists, for every x > y > 0, a positive r.v. Y (-1)

x,y such that, for every t ≥ 0 :

P (-1) x (R t ≤ y) = P (Y (-1) x,y ≤ t) (2.2)
Thus, in Theorem 2.1, the value ν = -1 2 appears as a critical value. This is, in fact, rather natural. Indeed, for ν ≥ -1 2 , (R t , t ≥ 0) is a submartingale -therefore it is, in a way, increasing in t -and it may be written as :

R t = x + B t + 2ν + 1 2 t 0 ds R s if ν > - 1 2
where (B t , t ≥ 0) denotes a Brownian motion starting from 0, and for ν = -1 2 ,

R t = |x + B t | = |x| + β t + L -x t
, where (L -x t , t ≥ 0) denotes the local time of B at level -x and (β t , t ≥ 0) is a Brownian motion.

For ν ∈ -1, -1 2 , such a representation is no longer true. In fact :

R t = x + β t + 2ν + 1 2 k t with k t := p.v. t 0 ds R s := ∞ 0 a 2ν (L a t -L 0 t )da
where (L a t , t ≥ 0, a ≥ 0) denotes the jointly continuous family of diffusion local times associated with (R t , t ≥ 0) see [RY], Exerc. 1.26,Chap. XI,p. 451 ; in this case, (R t , t ≥ 0) is no longer a semimartingale, but it is a Dirichlet process. See Bertoin [Ber] for a deep study of excursion theory for the process (R t , k t ), t ≥ 0 . In order to prove Theorem 2.1, we gather a few results about Bessel processes.

2.2 A summary of some results about Bessel processes 2.2.a The density of the Bessel semi-group see [RY], Chap. XI We denote by p (ν) (t, x, y) the density of R t under P (ν)

x (ν > -1) ; one has :

p (ν) (t, x, y) = y t y x ν exp - x 2 + y 2 2t I ν xy t (ν > -1) (2.3) if x > 0 whereas for x = 0 p (ν) (t, 0, y) = 2 -ν t -(ν+1) Γ(ν + 1) y 2ν+1 exp - y 2 2t (ν > -1) (2.4)
where I ν denotes the modified Bessel function with index ν see [Leb], p. 108 .

2.2.b Density of last passage times See [PY1]

Let

G y := sup{t ≥ 0 ; R t = y} (2.5) (= 0 if this set is empty)
Then, for x ≤ y :

P (ν) x (G y ∈ dt) = ν y p (ν) (t, x, y)dt (ν > 0) (2.6)
It follows easily from (2.6) and (2.3) that, under P (ν) 0 , the law of G y is that of

y 2 2γ ν
, where γ ν is a gamma r.v. with parameter ν (ν > 0) :

G y (law) = y 2 2γ ν (G y being considered under P (ν) 0 ) (2.7)
and :

E (ν)
x (e -λGy ) = 2ν

y x ν I ν (x √ 2λ)K ν (y √ 2λ) (0 < x < y ; ν > 0) E (ν) 0 (e -λGy ) = 1 2 ν-1 Γ(ν) (y √ 2λ) ν K ν (y √ 2λ) (y > 0 ; ν > 0)
2.2.c Laplace transform of first hitting times See [K] and [START_REF] Pitman | Yor Bessel processes and infinitely divisible laws[END_REF] Let

T y := inf{t ≥ 0 ; R t = y} (2.8) (= +∞ if this set is empty)
Then :

E (ν)
x (e -λTy ) =

y x ν I ν (x √ 2λ) I ν (y √ 2λ) (x < y, λ ≥ 0) (2.9)
In particular, for x = 0

E (ν) 0 (e -λTy ) = 1 2 ν Γ(ν + 1) (y √ 2λ) ν I ν (y 2λ) (2.10) 2.2.d Resolvent kernel see [B.S] Let u (ν) λ (x, y) = ∞ 0 e -λt p (ν) (t, x, y)dt (λ ≥ 0)
Then, for every positive Borel function f :

∞ x u (ν) λ (x, y)f (y)dy = 2 x ν I ν (x √ 2λ) ∞ x y ν+1 K ν (y √ 2λ)f (y)dy (2.11) and u (ν) λ (x, y) = 2y y x ν I ν (x √ 2λ)K ν (y √ 2λ) (x < y)
where K ν denotes the Bessel-Mc Donald function with index ν see [Leb], p. 108 .

2.2.e Bessel realisations of the Hartman-Watson laws See [Y] and [PY1]

It was established by Hartman-Watson [H.W] that, for any ν ≥ 0 and r > 0, the ratio

λ -→ I √ λ+ν 2 I ν (r)
is the Laplace transform of a probability, say θ (ν) r on R + , which may then be called Hartman-Watson distribution. In [Y] and [START_REF] Pitman | Yor Bessel processes and infinitely divisible laws[END_REF] these laws were shown to be those of the Bessel clocks t 0 ds R 2 s for Bessel bridges (R s , s ≤ t) conditioned at both ends. Precisely : For ν > -1 and µ = 0, one has :

E (ν) x exp - µ 2 2 t 0 ds R 2 s R t = y = I √ µ 2 +ν 2 I ν xy t (x, y ≥ 0) (2.12)
This formula may be obtained from the particular case ν = 0 :

E (0) x exp - µ 2 2 t 0 ds R 2 s R t = y = I |µ| I 0
xy t together with the absolute continuity relationship :

P (ν) x Rt∩(t<T 0 ) = R t x ν exp - ν 2 2 t 0 ds R 2 s • P (0) x | Rt (2.13)
where (R t , t ≥ 0) denotes the natural filtration of (R t , t ≥ 0). If ν > 0, R t ∩ (t < T 0 ) may be replaced by R t in formula (2.13) since then T 0 = ∞ a.s., whereas if ν ∈] -1, 0[, then, for µ 2 > 0 :

E (ν) x exp - µ 2 2 t 0 ds R 2 s R t = y = E (ν) x exp - µ 2 2 t 0 ds R 2 s 1 T 0 >t |R t = y since, if T 0 < t, t 0 ds R 2 s = +∞.
We make several remarks : i) For ν < 0, letting µ → 0 in (2.12), we get :

P (ν) x (T 0 > t|R t = y) = I -ν I ν xy t (2.14) ii) For ν > - 1 2
, formula (2.12) becomes, taking µ 2 = 2ν + 1 :

E (ν) x exp -ν + 1 2 t 0 ds R 2 s R t = y = I ν+1 I ν xy t ν > - 1 2 (2.15)
From (2.15), we deduce that, for ν ≥ -1 2 and z ≥ 0 :

I ν+1 (z) ≤ I ν (z) (2.16) (2.16
) is an immediate consequence of (2.15) for ν > -1 2 and for ν = -1 2 :

I1 2 (z) = 2 πz 1 2 sinh(z) ≤ I -1 2 (z) = 2 πz 1 2
cosh(z) [Leb], p. 112 .

2.3.f The Hirsch-Song formula see [H.S]

Proposition 2.2 Let ν > - 1 2 • i) ∂ ∂x p (ν) (t, x, y) = - ∂ ∂y x y p (ν+1) (t, x, y) (2.17)
ii) Let H denote the operator defined on the space of C 1 functions, with derivative equal to 0 at 0, and which are bounded, as well as their derivative :

Hf (x) = 1 x f ′ (x) (2.18) Let (P (ν)
t , t ≥ 0) denotes the Bessel semi-group with index ν. Then :

HP (ν) t = P (ν+1) t H (2.19) iii) Let (Q (ν)
t , t ≥ 0) denotes the semi-group of the Bessel squared process, with index ν. Then :

DQ (ν) t = Q (ν+1) t D (2.20)
where D denotes the differentiation operator : Df (x) = f ′ (x), with domain the space of C 1 functions, bounded as well as their derivative.

It is easily verified that (2.19) and (2.20) are equivalent since

Q (ν) t f (z) = P (ν) t ( f )( √ z), with f (z) = f (z 2 ) .
On the other hand, denoting by L (ν) the infinitesimal generator of the semi-group

Q (ν) t : (see point 2.2.g below) L (ν) f (x) = 2x f ′′ (x) + 2(ν + 1)f ′ (x) (2.21)
an easy computation allows to obtain

DL (ν) = L (ν+1) D, hence (2.20).
We show (2.17), following [HS] For f regular with compact support, one has ;

P (ν) t f (x) = ∞ 0 p (ν) (t, x, y)f (y)dy, hence : ∂ ∂x P (ν) t f (x) = ∞ 0 ∂ ∂x p (ν) (t, x, y)f (y)dy (2.22)
On the other hand, with obvious notations :

∂ ∂x P (ν) t f (x) = ∂ ∂x E (ν) f (R x t ) = E (ν) f ′ (R x t ) ∂R x t ∂x
However, since :

R x t = x + B t + 2ν + 1 2 1 0 ds R x s ν > - 1 2
we obtain, by differentiation with respect to x :

∂R x t ∂x = 1 - 2ν + 1 2 t 0 ds (R x s ) 2 ∂R x s ∂x (2.23)
The linear equation (2.23) may be integrated, and we obtain :

∂R x t ∂x = exp - 2ν + 1 2 t 0 ds (R x s ) 2
(2.24)

Thus :

∞ 0 ∂ ∂x p (ν) (t, x, y)f (y)dy = E (ν) x f ′ (R t ) exp - 2ν + 1 2 t 0 ds R 2 s = ∞ 0 E (ν) x exp - 2ν + 1 2 t 0 ds R 2 s R t = y f ′ (y)p (ν) (t, x, y)dy = ∞ 0 p (ν) (t, x, y)f ′ (y) I ν+1 I ν xy t dy from (2.15) = - ∞ 0 f (y) ∂ ∂y I ν+1 I ν xy t p (ν) (t, x, y) dy = - ∞ 0 f (y) ∂ ∂y x y p (ν+1) (t, x, y) dy from (2.3) (2.25)
The comparison of (2.25) and (2.22) then implies point i) of Proposition 2.2.

2.2.g The Fokker-Planck formula

The infinitesimal generator L (ν) of the Bessel semi-group (P (ν) t , t ≥ 0) with index ν is given by :

L (ν) f (x) = 1 2 f ′′ (x) + 2ν + 1 2x f ′ (x) (2.26)
Its domain is the space of functions f such that L (ν) f is bounded continuous and satisfies : lim

x→0 x 2ν+1 f ′ (x) = 0. One has ∂ ∂t p (ν) (t, x, y) = L (ν) p (ν) (t, •, y) (x) = L (ν) * p (ν) (t, x, •) (y) (F okker -P lanck) (2.27)
where the operator L (ν) * , the adjoint of L (ν) , is defined by :

L (ν) * f (x) = 1 2 f ′′ (x) - ∂ ∂x 2ν + 1 2x f (x) (2.28)
The infinitesimal generator L (ν) of the semi-group (Q

(ν)
t , t ≥ 0) of the squared Bessel process with index ν, is given by :

L (ν) f (x) = 2xf ′′ (x) + (2ν + 2)f ′ (x) (2.29)
and its adjoint L (ν) * is defined by :

L (ν) * f (x) = ∂ 2 ∂x 2 2xf (x) -(2ν + 2)f ′ (x) (2.30)

2.2.h Sojourn time below level y see [PY1]

, also [MY], Th. 11.6, p. 180 Let

A - y = ∞ 0 1 (Rt≤y) dt = Gy 0 1 (Rt≤y) dt
Then :

E (ν) y (e -λ A - y ) = 2ν y √ 2λ I ν (y √ 2λ) I ν-1 (y √ 2λ) (ν > 0 ; λ > 0) (2.31) 2.3 Proof of Theorem 2.1 2.3.1 A useful Lemma Lemma 2.3 Let I ν denote the modified Bessel function with index ν see [Leb], p. 108 i) If ν ≥ - 1 2
, then, for every z ≥ 0,

I ν (z) ≥ I ν+1 (z) (2.32)
ii) If ν = -1, then, for every z ≥ 0,

I 1 (z) = I -1 (z) ≤ I 0 (z) (2.33) iii) If ν ∈] -1, -1 2 [ for z small enough I ν (z) > I ν+1 (z) whereas for z large enough, I ν (z) < I ν+1 (z)
Proof of Lemma 2.3 Point i) has been proven : it is relation (2.16). We now prove point ii) It follows from the previous point since I -1 = I 1 . We prove point iii) In the neighborhood of 0, one has :

I ν (z) I ν+1 (z) ∼ 1 Γ(ν + 1) z 2 ν Γ(ν + 2) ( z 2 ) ν+1 = 2 z (ν + 1) -→ z→0 +∞ (ν > -1) (2.34)
In the neighborhood of +∞, one has see [Leb], p. 122 and 123 :

I µ (z) = e z √ 2πz 1 - 4µ 2 -1 4 1 2z + o 1 z 2
Thus, with µ = -ν and ν ∈ 1 2 , 1 :

I -ν (z) = e z √ 2πz 1 - 4ν 2 -1 4 1 2z + o 1 z 2 I -ν+1 (z) = e z √ 2πz 1 - 4(1 -ν) 2 -1 4 1 2z + o 1 z 2 Now, for ν ∈] 1 2 , 1[, one has 4ν 2 -1 ≥ 4(1 -ν) 2 -1 ; indeed, this is equivalent to ν 2 ≥ 1 + ν 2 -2ν, i.e. ν ≥ 1 2 • Hence, for ν ∈] 1 2 , 1[, the ratio I -ν (z) I -ν+1 (z) -→ z→∞ 1 -as z → ∞.
As a conclusion, from this point and (2.34), for ν ∈ -1, -1 2 :

• for z large enough,

I ν (z) I ν+1 (z)
< 1 ;

• for z small enough

I ν (z) I ν+1 (z) > 1• 2.3.2
The case x = 0 When x = 0 and ν > -1, we shall show the existence of an increasing process (J t , t ≥ 0) such that J t -→ t→∞ +∞ a.s. and

P (ν) 0 (R t ≥ y) = P (J t ≥ y) (2.35)
This formula obviously implies the existence of a process which is a pseudo-inverse of (R t , t ≥ 0). Let us prove (2.35), which hinges in fact simply on the scaling property of (R t , t ≥ 0). Indeed, one has :

P (ν) 0 (R t ≥ y) = P (ν) 0 √ t ≥ y R 1 (by scaling) = P (ν) 0 t ≥ y 2 R 2 1 = P t ≥ y 2 2γ ν+1 since R 2 1 (law) = 2γ ν+1 under P (ν)
0 , where γ ν+1 is a gamma variable with parameter ν + 1 .

= P (ν+1) 0 (G y ≤ t) from (2.7) = P (ν+1) 0 inf u≥t R u ≥ y = P (ν+1) 0 (J t ≥ y) with J t := inf u≥t R u . Clearly (J t , t ≥ 0) is an increasing process and J t -→ t→∞ +∞ P (ν+1) 0 a.s. since ν + 1 > 0. 2.3.3 We now prove point i) of Theorem 2.1
Since, from the comparison theorem, (applied to Bessel processes) for ν ≥ -1 2 , one has :

P (ν) x (R t ≥ y) ≥ P (-1 2 ) 0 (R t ≥ y) = P |B t | ≥ y = P 0 |B 1 | > y √ t -→ t→∞ 1
where (B t , t ≥ 0) is a Brownian motion started from 0. It suffices, in order to prove point i) of Theorem 2.1 to show that : t → P (ν)

x (R t ≥ y) is an increasing function of t. Now, one has :

∂ ∂t P (ν) x (R t ≥ y) = ∞ y ∂ ∂t p (ν) (t, x, z)dz = ∞ y 1 2 ∂ 2 ∂z 2 p (ν) (t, x, z) - ∂ ∂z 2ν + 1 2z p (ν) (t, x, z) dz (from Fokker-Planck) = - 1 2 ∂ ∂y p (ν) (t, x, y) + 2ν + 1 2y p (ν) (t, x, y)
Now, from (2.3) : ν) x (R t ≥ y) is positive, since y > x and I ν (z) ≥ I ν+1 (z).

∂ ∂y p (ν) (t, x, y) = p (ν) (t, x, y) ν + 1 y - y t + x t I ′ ν ( xy t ) I ν ( xy t ) Hence : ∂ ∂t P (ν) x (R t ≥ y) = p (ν) (t, x, y) ν 2y + y 2t - x 2t
I ′ ν ( xy t ) I ν ( xy t ) = p (ν) (t, x, y) ν 2y + y 2t - x 2t
I ν+1 I ν xy t + νt xy since I ′ ν (z) I ν (z) = I ν+1 I ν (z) + ν z [Leb], p. 110 = p (ν) (t,

2.3.4

We now prove point ii) of Theorem 2.1 This follows immediately from (2.36) and from point iii) of Lemma 2.3 : for x = 0, ∂ ∂t P (ν) x (R t ≥ y) does not have a fixed sign, hence t → P (ν)

x (R t ≥ y) cannot be a monotone function.

2.3.5

We now prove point iii) of Theorem 2.1 We need to show that, for every x > 0 and 0 < y < x, P (-1) x (R t ≤ y) is an increasing function of t. We may of course replace R t by R 2 t . Now, the law of R 2 t , when R 2 0 = x, is given by :

P (R 2 t ∈ dy) = exp - x 2t
• δ 0 (dy) + q(t, x, y)dy with q(t, x, y) := 1 2t

x y exp - x + y 2t I 1 xy t (2.37)
Then, we need to show that, for every y < x, one has :

∂ ∂t exp - x 2t + ∂ ∂t y 0 q(t, x, z)dz ≥ 0 But, since ∂ ∂t exp - x 2t
≥ 0, it suffices to show that :

∂ ∂t y 0 q(t, x, z)dz ≥ 0.
Now, one has :

∂ ∂t y 0 q(t, x, z)dz = y 0 ∂ ∂t q(t, x, z)dz = y 0 ∂ 2 ∂z 2 2zq(t, x, z) dz (from Fokker-Planck) = ∂ ∂y 2y q(t, x, y) = 2q(t, x, y) + 2y q(t, x, y) - 1 2y - 1 2t + 1 2t
x y

I ′ 1 I 1 √ xy t from (2.37) = q(t, x, y) 1 - y t + √ xy t I 0 I 1 √ xy t - t √ xy since I ′ 1 I 1 (z) = I 0 I 1 (z) - 1 z = √ y t q(t, x, y) √ x I 0 I 1 √ xy t - √ y ≥ 0
from point ii) of Lemma 3.2 and since y < x. This ends the proof of Theorem 2.1.

Remark 2.4 i) The property : t → P (ν)

x (R t ≥ y) is an increasing function of t (y > x) is equivalent to : for every Borel increasing function ϕ with support in [x, ∞[, t → E (ν) x ϕ(R t ) is an increasing function of t. Indeed, we have, for such a ϕ :

E (ν) x ϕ(R t ) = E (ν) x Rt x dϕ(y) = E (ν) x ∞ x dϕ(y)1 (Rt≥y) = ∞ x dϕ(y)P (ν)
x (R t ≥ y)

ii) Around formula (2.36) Formula (2.36) :

∂ ∂t P (ν) x (R t ≥ y) = p (ν) (t, x, y) 2t I ν ( xy t ) y I ν xy t -x I ν+1 xy t
may still be written, with the help of (2.3) :

∂ ∂t P (Y (ν) x,y ≤ t) = ∂ ∂t P (ν) x (R t ≥ y) = y 2t p (ν) (t, x, y) - x 2 2ty p (ν+1) (t, x, y) (2.38)
We denote by θ (ν) (t, x, y) the RHS of (2.38). We now check directly that this positive function of t integrates to 1 i.e. : it is a density of probability :

∞ 0 θ (ν) (t, x, y)dt = 1 (2.39)
For this purpose, we use the classical Lipschitz-Hankel formula see, e.g. [Wat], p. 384 , for ν > 0 :

ν ∞ 0 du u exp(-au)I ν (u) = 1 (a + √ a 2 -1) ν (a > 1) (2.40)
Differentiating both sides with respect to a, we obtain :

∞ 0 exp(-au) • I ν (u)du = 1 (a + √ a 2 -1) ν • 1 √ a 2 -1 (a > 1, ν > 0) (2.41)
In order to check (2.39), we first prove that :

Λ (ν) x,y := ∞ 0 dt t p (ν) (t, x, y) = 2y y 2 -x 2 (y > x) (2.42)
Note that (2.42) (which we assume for a moment) shows that Λ (ν)

x,y does not depend on ν ! We deduce from (2.42) that :

∞ 0 y 2t p (ν) (t, x, y) - x 2 2ty p (ν+1) (t, x, y) dt = y 2 Λ (ν)
x,y -

x 2 2y Λ (ν+1) x,y = y 2 - x 2 2y 2y y 2 -x 2 = 1
which is the desired result (2.39). We now prove (2.42) From (2.3), we get :

Λ (ν) x,y = y x ν y ∞ 0 dt t 2 exp - x 2 + y 2 2t I ν xy t = y x ν y ∞ 0 du exp - x 2 + y 2 2 u I ν (xyu) (after the change of variable 1 t = u) = y x ν y xy ∞ 0 dv exp - x 2 + y 2 2xy v I ν (v) (after the change of variable xyu = v) = y x ν 1 x 1 (a + √ a 2 -1) ν 1 √ a 2 -1
with the help of (2.41), where a := x 2 + y 2 2xy (> 1 since (x < y). On the other hand, there are the elementary identities :

√ a 2 -1 = y 2 -x 2 2xy , a + √ a 2 -1 = y x
, which implies :

Λ (ν) x,y = 1 x 1 y 2 -x 2 2xy = 2y y 2 -x 2
which is the desired result (2.42). We shall come back to this computation, in Section 3, in order to establish certain properties of the r.v.'s Y (ν)

x,y .

2.4 Interpretation in terms of the local martingales (R -2ν t , t ≥ 0) We assume here ν > 0. We know that (R t , t ≥ 0) solves the S.D.E. :

R t = r + B t + 2ν + 1 2 t 0 ds R s (2.43)
where (B t , t ≥ 0) is a Brownian motion starting from 0. Let M (ν) t := R -2ν t and denote M (ν) 0 = x (= r -2ν ). Itô formula yields :

M (ν) t = x -2ν t 0 (M (ν) s ) 2ν+1 2ν dB s .
(2.44)

We denote by

Q (ν)
x the law of the solution of (2.44). Thus (M

(ν) t , t ≥ 0 ; Q (ν)
x , x ≥ 0) is the family of laws of a Markov process and (M (ν) t , t ≥ 0) is a positive local martingale, hence, it is a supermartingale, which converges a.s. to 0 as t → ∞. Gathering the results of Theorem 2.1, point ii) of Remark 1.3 and Proposition 1.4, we have obtained :

Theorem 2.4 i) For every ν > 0 and x > 0, there exists a probability measure γ M (ν,x) on [0, x] × [0, ∞[ such that :

1 x Q (ν) x (K -M (ν) t ) + = γ M (ν,x) [0, K] × [0, t] (K ≤ x, t ≥ 0) (2.

45)

ii) There exists a family of positive r.v.'s (Y

x,y , y < x) such that :

Q (ν) x (M t ≤ y) = P (Y (ν) x,y ≤ t) (y < x, t ≥ 0) (2.46)
Remark 2.5 i) Since :

Q (ν) x (M t ≤ y) = P (ν) x -1 2ν 1 R 2ν t ≤ y = P (ν) x -1 2ν (R t ≥ y -1 2ν ) (2.47)
we have :

Y (ν) x,y (law) = Y (ν) x -1 2ν ,y -1 2ν (x > y) (2.48)
In Section 3, we shall study in details the laws of the r.v.'s (Y

x,y , y > x) which, thanks to (2.48), allows to obtain easily the corresponding properties of the r.v.'s (Y

(ν) x,y , y < x). ii) When ν ∈]-1, 0[, (R -2ν t , t ≥ 0) is no longer a local martingale, but it is a submartingale see [D.M, RVY] such that : R -2ν t = N (ν) t + L (ν) t (2.49)
where (N (ν) t , t ≥ 0) is a martingale and (L (ν) t , t ≥ 0) an increasing process such that dL

(ν) t = 1 (Rt=0) dL (ν) t .
3 Some properties of the r.v.'s (Y (ν)

x,y , y > x) ν ≥ -1 2 We recall that these r.v.'s Y (ν)

x,y (or, more exactly, their laws) are defined (see Theorem 2.1) via :

P (ν) x (R t ≥ y) = P (Y (ν) x,y ≤ t) (3.1)
Theorem 3.1 Here are a few properties of the laws of

Y (ν) x,y v ≥ - 1 2 , y > x i) Laplace transform of Y (ν) x,y • E(e -λY (ν) x,y ) = I ν (x √ 2λ) (x √ 2λ) ν (y √ 2λ) ν+1 K ν+1 (y √ 2λ) (x > 0, λ ≥ 0) (3.2) = (y √ 2λ) y x ν I ν (x √ 2λ)K ν+1 (y √ 2λ)
For x = 0, we get :

• E(e -λY (ν) 0,y = 1 2 ν Γ(ν + 1) (y √ 2λ) ν+1 K ν+1 (y √ 2λ) (3.3)
ii) Scaling property

• The law of Y (ν)
x,y xy depends only on the ratio y x (> 1). In particular :

• Y (ν) x,y = x 2 Y 1, y x (y > x) (3.4)
and

• E e -λY (ν)
x,y xy

= √ 2λ b ν+ 1 2 I ν 2λ b K ν+1 ( √ 2λb) (3.5) where b = y x > 1.
iii) Further results about the laws of

Y (ν) x,y • Y (ν) 0,y (law) = y 2 2γ ν+1 (law) = G (ν+1) y (3.6)
where G (ν+1) y denotes the r.v. G y under P (ν+1) 0 . The following relations hold :

• P (Y (ν) x,y ≤ t) = P γ ν+1 ≥

y 2 2t + x 0 z y p (ν+1) (t, z, y)dz (3.7) • ∂ ∂x P (Y (ν) x,y ≤ t) = x y p (ν+1) (t, x, y) = x ν + 1 f G (ν+1) x,y (t) (3.8)
where f G (ν+1)

x,y denotes the density of the r.v. G y under P

(ν+1) x (x < y) • The density f Y (ν) x,y of Y (ν)
x,y is given by :

f Y (ν) x,y (t) = 1 2t y p (ν) (t, x, y) - x 2 y p (ν+1) (t, x, y) (3.9) • E(Y (ν) x,y ) = y 2 2ν - x 2 2(ν + 1) (ν > 0) (3.10) iv) An equation satisfied by Y (ν) x,y Y (ν)
x,z

(law) = T (ν) x,y + Y (ν) y,z (x < y < z) (3.11)
In particular :

Y (ν) 0,z (law) = T (ν) 0,y + Y (ν) y,z (0 < y < z) (3.12)
The r.v.'s which occur on the RHS of (3.11) and (3.12) are independent and T (ν)

x,y is the first hitting time of level y by the process (R t , t ≥ 0) starting from x.

v) T

(ν)

0,x + Y (ν)
x,y

(law) = G (ν+1) y (0 < x < y) (3.13)
The r.v.'s which occur on the LHS of (3.13) are independent see (2.5), (2.6) and (2.7) for the definition of G From (3.1), we deduce :

∞ 0 e -λt P (ν) x (R t ≥ y)dt = ∞ 0 dt ∞ y e -λt p (ν) (t, x, z)dz = ∞ y dz ∞ 0 e -λt p (ν) (t, x, z)dt = ∞ y u (ν) λ (x, z)dz see (2.2.d) = E ∞ 0 e -λt 1 (t≥Y (ν) x,y ) dt = 1 λ E(e -λY (ν)
x,y )

Hence :

E(e -λY (ν) x,y ) = λ ∞ y u (ν) λ (x, z)dz = 2λ x ν I ν (x √ 2λ)(2λ) -ν 2 -1 ∞ y √ 2λ K ν (h)h ν+1 dh from (2.11)
Now, since : [Leb], p. 110 we obtain :

-z ν+1 K ν (z) = ∂ ∂z z ν+1 K ν+1 (z)
E(e -λY (ν) x,y ) = (2λ) -1 2 I ν (x √ 2λ) x ν (y √ 2λ) ν+1 K ν+1 (y √ 2λ) = I ν (x √ 2λ) (x √ 2λ) ν (y √ 2λ) ν+1 K ν+1 (y √ 2λ)
Formula (3.3) may be obtained by letting x tend to 0 in (3.2) and (3.10) follows from (3.2) by differentiation.

ii) Proof of the scaling property :

f xy Y (ν) x,y (u) = 1 2 y x ν+1 e -a(x,y)u I ν (u) - x y I ν+1 (u)
where a(x, y) :=

x 2 + y 2 2xy = 1 2 x y + y x • Thus, the law of xy Y (ν) x,y
, hence that of Y • We may also prove this scaling property as a direct consequence of the scaling property of the Bessel process. Indeed, we have :

P (Y (ν)
x,y ≤ t xy) = P (ν) x (R txy ≥ y)

= P (ν) x 1 √ xy R txy ≥ y √ xy = P (ν) x √ xy R t ≥ y x
(by scaling of the Bessel process)

= P (ν) √ x y R t ≥ y x •
iii) Proof of (3.6) This formula (3.6) has been obtained during the proof of Theorem 2.1 (point 2.3.2, in the case x = 0). We note that one can deduce (3.3) from (3.6). Indeed, from (3.6), since

P (Y (ν) 0,y ≤ t) = P y 2 2γ ν+1 ≤ t , one has : ∞ 0 e -λt P (Y (ν) 0,y ≤ t)dt = ∞ 0 e -λt P y 2 2γ ν+1 ≤ t dt, i.e. :
E(e -λY (ν) 0,y ) = E(e

-λ 2 y 2 γ ν+1 ) = 1 Γ(ν + 1) ∞ 0 z ν exp - 1 2 2z + λy 2 z dz = 2 1-ν 2 Γ(ν + 1) λ ν+1 2 y ν+1 K ν+1 (y √ 2λ) [Leb], p.119 = 1 2 ν Γ(ν + 1) (y √ 2λ) ν+1 K ν+1 (y √ 2λ)
iv) We now prove (3.8) and (3.7) Since, from the Hirsch-Song formula (2.17), one has :

∂ ∂x p (ν) (t, x, y) = - ∂ ∂y x y p (ν+1) (t, x, y)
then :

∂ ∂x P (Y (ν) x,y ≤ t) = ∂ ∂x P (ν) x (R t ≥ y) = ∞ y ∂ ∂x p (ν) (t, x, z)dz = - ∞ y ∂ ∂z x z p (ν+1) (t, x, z) dz = x y p (ν+1) (t, x, y)
and the relation :

x y p (ν+1) (t, x, y) = x ν + 1 f G (ν+1) x,y (t)
follows from (2.6). Formula (3.7) is obtained by integration of (3.8) (with respect to x) : ν+1) (t, z, y)dz from (3.6) and (3.8). We note that (3.7) may also be written, from (3.8) :

P (Y (ν) x,y ≤ t) = P (Y (ν) 0,y ≤ t) + x 0 ∂ ∂z P (Y (ν) z,y ≤ t)dz = P γ ν+1 ≥ y 2 2t + x 0 z y p (
P (Y (ν) x,y ≤ t) = P γ ν+1 ≥ y 2 2t + x 0 z ν + 1 f G (ν+1)
x,y

(t)dz (3.14)

We also remark that the computation of E(Y

(ν)
x,y ), given by (3.10) when ν > 0, may be obtained from (3.14). Indeed we deduce from (3.14) that :

P (Y (ν) x,y ≥ t) = P γ ν+1 ≤ y 2 2t - x 0 z ν + 1 f G (ν+1) z,y (t)dz
Thus, integrating this relation in t from 0 to ∞, for ν > 0, we obtain :

P (Y (ν) x,y ≥ t) = ∞ 0 dt y 2 2t 0 1 Γ(ν + 1) e -z z ν dz - x 2 2(ν + 1) = 1 Γ(ν + 1) ∞ 0 e -z z ν y 2 2z dz - x 2 2(ν + 1) = y 2 2ν - x 2 2(ν + 1)
v) Formula (3.9) is an immediate consequence of (3.13).

vi) We now prove that Y (ν)

x,y satisfies equation (3.11) Indeed, this follows from a simple application of the Markov property. Let x < y < z. Since the process (R t , t ≥ 0) starting from x needs to pass through y to reach z, we obtain :

P (ν) x (R t ≥ z) = P (ν) x (T y < ∞ ; R t ≥ z) = P (ν) x 1 Ty≤t P (ν) y ( R t-Ty ≥ z)
where in the expression P (ν) y ( R t-Ty ≥ z) the term T y is frozen. Hence, conditioning with respect to T y = u, we obtain :

P (ν) x (R t ≥ z) = t 0 P (ν)
x (T y ∈ du)P (ν) y (R t-u ≥ z)du, i.e. :

P (Y (ν) x,z ≤ t) = t 0 P (ν) x (T y ∈ du)P (Y (ν) y,z ≤ t -u)du hence : Y (ν)
x,z

(law) = T (ν) x,y + Y (ν) y,z
It is also possible to obtain (3.11) by using the Laplace transforms of Y (ν)

x,y and T (ν)

x,y . Indeed, from (3.2) :

E(e -λY (ν) x,z ) = I ν (x √ 2λ) (x √ 2λ) ν (z √ 2λ) ν+1 K ν+1 (z √ 2λ)
whereas, from (3.2) and (2.9) :

E(e -λ(T (ν)

x,y +Y (ν)

y,z ) ) = y x ν I ν (x √ 2λ) I ν (y √ 2λ) • I ν (y √ 2λ) (y √ 2λ) ν (z √ 2λ) ν+1 K ν+1 (z √ 2λ) = I ν (x √ 2λ) (x √ 2λ) ν (z √ 2λ) ν+1 K ν+1 √ 2λ) = E(e -λY (ν) x,z )
vii) The proof of (3.13) hinges on the same arguments as previously and, from (2.7), on :

E(e -λG (ν+1) y ) = 1 2 ν Γ(ν + 1) (y √ 2λ) ν+1 K ν+1 (y √ 2λ)
Indeed, from (2.6), :

E(e -λG (ν+1) y ) = E(e -λ y 2 2γ ν+1 ) = 1 Γ(ν + 1) ∞ 0 e -λy 2 2z e -z z ν dz = 1 2 ν Γ(ν + 1) (y √ 2λ) ν+1 K ν+1 (y √ 2λ)
see [START_REF]Lebedev Special functions and their applications[END_REF]p. 119,(5.10.25) . This ends the proof of Theorem 3.1.

Remark 3.2 i) We note that the r.v.'s (Y

(ν)
x,y , y > x) are not the only ones which satisfy equation (3.11). Indeed, let, for x < y :

A (ν),- x,y := ∞ 0 1 R (ν) x (s)≤y ds = Gy 0 1 R (ν)
x (s)≤y ds where here R (ν)

x (s), s ≥ 0 denotes the Bessel process with index ν starting from x. Then, an application of the Markov property yields, for x < y < z :

A (ν),- x,z (law) = T (ν)
x,y + A (ν),- y,z

(3.15)

We note that, although the r.v.'s (A (ν),-

x,y , y > x) satisfy the same equation as the r.v.'s (Y (ν) x,y , y > x), they do not have the same law. Indeed, for x = 0, from Ciesielski-Taylor [CT], 1962

A (ν),- 0,y (law) = T (ν-1) 0,y (ν > 0) (3.
E(e -λT (ν-1) 0,y ) = 1 2 ν-1 Γ(ν) (y √ 2λ) ν-1 I ν-1 (y √ 2λ)
, from (2.10), whereas, from (3.3) and (3.6),

E(e -λY (ν) 0,y ) = E (ν+1) 0 (e -λGy ) = 1 2 ν Γ(ν + 1) (y √ 2λ) ν+1 K ν+1 (y √ 2λ)
ii) For ν > 0, the equality :

E(Y (ν) x,y ) = y 2 2ν - x 2 2(ν + 1)
may be obtained in a different manner than the previously developed ones, by making use this time of the r.v.'s A (ν),x,y . Indeed, for x < y :

E(Y (ν) x,y ) = ∞ 0 P (Y (ν) x,y ≥ t)dt = E (ν) x ∞ 0 1 (Rt≤y) dt = E(A (ν),- x,y ) (3.17)
Thus, we need to compute E(A (ν),-

x,y ). Itô's formula, for ν > 0, implies, from (2.43) :

(R t ∧ y) 2 = x 2 + 2 t 0 R s 1 (Rs<y) dB s + 2ν + 1 2R s ds -y L y t + t 0 1 (Rs<y) ds (3.18)
where (L y t , t ≥ 0) denotes the local time of (R t , t ≥ 0) at level y. Hence, taking expectation in (3.18) and letting t tend to ∞, we obtain :

y 2 = x 2 + (2ν + 2)E(A (ν),- x,y ) -y E (ν) x (L y ∞ ) (3.19)
However :

E (ν) x (L y ∞ ) = E (ν) y (L y ∞ ) = ∞ 0 p (ν) (t, y, y)dt = ∞ 0 y t exp - y 2 t I ν y 2 t dt from (2.3) = y ∞ 0 I ν (z)e -z dz z after the change of variable y 2 t = z = y ν (3.20)
from the Lipschitz-Hankel formula see (2.40). Finally, from (3.20), (3.19) and (3.17) :

E(Y ν x,y ) = E(A (ν),- x,y ) = 1 2ν + 2 y 2 1 + 1 ν -x 2 = y 2 2ν - x 2 2(ν + 1)
iii) It may be of interest to express the law of the r.v. Y (ν)

x,y in terms of the only process (R t , t ≥ 0 ; P (ν)

x ). Here is our result : Proposition 3.3 Let ν > 0 and h a generic positive Borel function. Then :

E h(Y (ν)
x,y ) =

y 2 2ν E (ν) x h(G y ) G y 1 - x y • exp -ν + 1 2 Gy 0 ds R 2 s (3.21) = y(y -x) 2ν E (ν) x h(G y ) G y + xy 2ν E (ν) x h(G y ) G y 1 -e -(ν+ 1 2 ) R Gy 0 ds R 2 s (3.22)
Proof of Proposition 3.3

The density of the r.v. Y

x,y , f Y (ν)

x,y , equals, from (2.36) :

f Y (ν) x,y (t) = p (ν) (t, x, y) 2t y - I ν+1 I ν xy t = y p (ν) (t, x, y) 2t 1 - x y E (ν) x exp -ν + 1 2 t 0 ds R 2 s R t = y (3.23) from (2.15)
On the other hand, from (2.6), the density f G (ν)

x,y of G y under P (ν) x equals :

f G (ν) x,y (t) = ν p (ν) (t, x, y) y (3.24)
Thus, since for every positive and predictable process H, one has see [FPY] :

E (ν) x (H Gy ) = ∞ 0 P (ν) x (G y ∈ dt)E (ν) x H t |R t = y (3.25)
we derive from (3.22) that :

E h(Y (ν) x,y ) = ∞ 0 h(t) 2t p (ν) (t, x, y) y 1 - x y E (ν) x exp -ν + 1 2 t 0 ds R 2 s R t = y = y 2 2ν ∞ 0 h(t) t P (ν) x (G y ∈ dt) 1 - x y E (ν) x exp -ν + 1 2 t 0 ds R 2 s R t = y = y 2 2ν E (ν) x h(G y ) G y 1 - x y exp -ν + 1 2 Gy 0 ds R 2 s from (3.24)
. This is formula (3.21), which, after some mild rearrangement, yields formula (3.22).

We note that both terms in (3.22) :

y(y -x) 2ν E (ν) x h(G y ) G y and xy 2ν E (ν) x h(G y ) G y 1 -e -(ν+ 1 2 ) R Gy 0 ds R 2 s
are positive measures (viewed via the integration of h). The first one has total mass :

y(y -x) 2ν E (ν) x 1 G y = y(y -x) 2ν ∞ 0 1 ty ν p (ν) (t, x, y)dt from (3.23) = y(y -x) 2ν ν y 2y y 2 -x 2 = y y + x from (2.42)
whereas the second has total mass x y + x

• In particular, for x = 0, we recover :

E h(Y (ν) 0,y = y 2 2ν E (ν) 0 h(G y ) G y
and the second term in (3.22) vanishes.

iv) The r.v.'s Y (-1)

x,y (x > y)

In the same way that point i) of Theorem 2.1 allows to define the r.v.'s Y

x,y , y > x ; ν ≥ -

1 2 ,
point iii) of this Theorem allows to define the positive r.v.'s (Y

(-1)
x,y , x > y) characterised by :

Q (-1) x (M t ≤ y) = P (Y (-1) x,y ≤ t)
where (M t , t ≥ 0) under Q (-1) x is a Bessel square process with index -1 started from x. With arguments similar to those used to prove Theorem 3.1, we obtain :

E(e -λY (-1) x,y ) = K 1 ( √ 2λx) √ 2λx + y √ x √ 2λ I 2 ( 2λy) (3.26) v) Use of

additivity property of squares of Bessel processes

We have shown, in Theorem 2.1, that for ν ≥ -1 2 and for y > 1 :

t → P (ν) 1 (R t ≥ y) is an increasing function of t. It is the distribution function of the r.v. Y (ν)
1,y (the general case, with x < y, may be deduced from this one, by scaling). On the other hand, it is well known see, e.g., [RY], Chap. XI that Q

(d) x * Q (d ′ ) x ′ = Q (d+d ′ ) x+x ′ , where Q (d) x (resp. Q d ′
x ′ ) denotes the law of a squared Bessel process, with dimension d = 2ν + 2, starting from x (resp. with dimension d ′ = 2ν ′ + 2 starting from x ′ ). Hence

P (ν) 1 (R t ≥ √ y) = P (ν) 1 (R 2 t ≥ y) = P R 2 0 (t) + X(t) ≥ y 29
where R 2 0 (t), t ≥ 0 is a squared Bessel process with dimension d = 2ν + 2 starting from 0 and (X t , t ≥ 0) is a squared Bessel process with dimension 0 started from 1. Hence :

P (ν) 1 (R t ≥ √ y) = P γd 2 ≥ 1 2t (y -X t ) + (3.27) since R 2 0 (t) (law)
= t R 2 0 (1) by scaling and R 2 0 (1)

(law) = 2γd
2 . Consequently :

P (ν) 1 (R t ≥ √ y) = 1 Γ( d 2 ) ∞ 0 P 1 2t (y -X t ) + ≤ z z d 2 -1 e -z dz (3.28)
Consequently, from (3.27), Theorem 2.1 hinges only upon the properties of the process (X t , t ≥ 0), which is a squared Bessel process started from 1. For d = 2 (i.e. ν = 0), (3.28) becomes :

P (0) 1 (R t ≥ √ y) = E exp - 1 2t (y -X t ) + (3.29)
since γd 2 = γ 1 is then a standard exponential variable. We note that, Theorem 2.1, applied for ν = 0 > -1 2 ! implies that the RHS of (3.29) is a distribution function (i.e. : that of the r.v. Y

1,y ). For d = 1 i.e. ν = -1 2 , (3.28) becomes :

P (-1 2 ) 1 (R t ≥ y) = 2 π E ∞ q (y-X t ) + t e -u 2 2 du (3.30)
which is also a distribution function with respect to t, that of the r.v. Y

(-1 2 ) 1,y . This remark invites to ask the following questions :

• Which are the probabilities π on R + such that : Proposition 3.4

E π 1 2t (y -X t ) + is a distribution function in t, with π(x) = π ]x,
i) Let x > 0 and V x a positive, 1 2 stable r.v. such that E(e -λVx ) = exp(-x √ 2λ) (λ ≥ 0) If ν > 1 2 , then V x + Y (ν)
x,y (with V x and Y (ν)

x,y independent) is infinitely divisible.

ii) If ν > 1 2 and y x is large enough, then Y (ν)

x,y is infinitely divisible.

Proof of Proposition 3.4

We prove i) : By scaling we may suppose x = 1 and y > 1. From (3.2) we have :

E(e -λY (ν) 1,y ) = y ν+1 √ 2λ I ν ( √ 2λ)K ν+1 (y √ 2λ) (3.31) := exp -h(λ) Hence : h ′ (λ) = - 1 2λ - 1 √ 2λ I ′ ν ( √ 2λ) I ν ( √ 2λ) - y √ 2λ K ′ ν+1 (y √ 2λ) K ν+1 (y √ 2λ) = y √ 2λ K ν K ν+1 (y √ 2λ) - 1 √ 2λ
I ν+1 I ν ( √ 2λ) (3.32) since ∂ ∂x x ν K ν (x) = -x ν K ν-1 (x), ∂ ∂x x -ν I ν (x) = x -ν I ν+1 (x)
; see [Leb], p. 110 . On the other hand see Ismail [I] :

1 √ 2λ K ν ( √ 2λ) K ν+1 ( √ 2λ) = 4 π 2 ∞ 0 1 2λ + t 2 dt t J 2 ν+1 (t) + Y 2 ν+1 (t) (3.33)
and :

I ν+1 ( √ 2λ) I ν ( √ 2λ) = √ 2λ ∞ n=1 2 2λ + j 2 ν,n (3.34)
where J ν+1 , Y ν+1 are the Bessel functions with index ν see [Leb], p. 98 and (j ν,n , n ≥ 1) the increasing sequence of the positive zeroes of J ν . Hence :

1 2 h ′ (λ) = 1 π ∞ 0 y 2 2λy 2 + t 2 2 πt dt J 2 ν+1 (t) + Y 2 ν+1 (t) - n≥1 1 2λ + j 2 ν,n (3.35) Let now Z y := V 1 + Y (ν)
1,y (y > 1). We have, from (3.35) :

E(e -λZy ) = exp -h(λ) - √ 2λ := exp -g(λ)
with :

1 2 g ′ (λ) = 1 2 h ′ (λ) + 1 2 √ 2λ = 1 π ∞ 0 y 2 2λy 2 + t 2 2 πt dt J 2 ν+1 (t) + Y 2 ν+1 (t) + 1 π ∞ 0 dt 2λ + t 2 - ∞ n≥1 1 2λ + j 2 ν,n (3.36) := (1) + (2) we used 1 √ 2λ = 2 π ∞ 0 dt 2λ + t 2 •
One needs, to prove i), to show that g ′ is completely monotone. For (1), it is clear since it is the Stieltjes transform of a positive measure. For (2), one has :

1 π ∞ 0 dt 2λ + t 2 - ∞ n≥1 1 2λ + j 2 ν,n = 1 π ∞ n=1 jν,n j ν,n-1 dt 2λ + t 2 -π ∞ n=1 1 2λ + j 2 ν,n = 1 π ∞ n=1 jν,n j ν,n-1 1 2λ + t 2 - 1 2λ + j 2 ν,n dt + (j ν,n -j ν,n-1 -π) 1 2λ + j 2 ν,n = 1 π ∞ n=1 jν,n j ν,n-1 j 2 ν,n -t 2 (2λ + t 2 )(2λ + j 2 ν,n ) dt + (j ν,n -j ν,n-1 -π) 1 2λ + j 2 ν,n
which is clearly completely monotone since see [I], p. 357 , j ν,nj ν,n-1 > π as ν > 1 2 •

ii) We now prove ii) One needs to show that h ′ (λ), as given by (3.24) is the Laplace transform, for y large enough, of a positive measure. Now, one has ;

1 2 h ′ (λ) = 1 π ∞ 0 1 2λ + v 2 2 πv dv J 2 ν+1 (yv) + Y 2 ν+1 (yv) - n≥1 1 2λ + j 2 ν,n (3.37) after the change of variable t = yv in (3.35) = 1 π ∞ n=1 jν,n j ν,n-1 1 2λ + v 2 - 1 2λ + j 2 ν,n 2 πv dv J 2 ν+1 (yv) + Y 2 ν+1 (yv) +π ∞ n=1 jν,n j ν,n-1 2 πv 1 J 2 ν+1 (yv) + Y 2 ν+1 (yv) -1 dv • 1 2λ + j 2 ν,n (3.38)
But, from Watson [Wat], p. 449 , one has :

J 2 ν+1 (z) + Y 2 ν+1 (z) ∼ y→∞ jν,n j ν,n-1 2 πv πvy 2 dv Hence : jν,n j ν,n-1 2 πv dv J 2 ν+1 (yv) + Y 2 ν+1 (yv) ∼ y→∞ jν,n j ν,n-1 2 πv πvy 2 dv ∼ z→∞ y(j ν,n -j ν,n-1 ) uniformly in n
Thus, for y large enough, and since ν > 1 2 implies j ν,nj ν,n-1 > π, one has :

1 2 h ′ (λ) = 1 π ∞ n=1 jν,n j ν,n-1 j 2 ν,n -v 2 (2λ + v 2 )(2λ + j 2 ν,n ) 2 πv dv J 2 ν+1 (yv) + Y 2 ν+1 (yv) + ∞ n=1 α n (y) 1 2λ + j 2 ν,n
where the α n (y) are all positive. It is then clear that h ′ (λ) is completely monotone. vii) On the negative moments of Y (ν)

x,y (x < y) For every m > 0, we have :

E 1 (Y (ν)
x,y ) m < ∞

(3.39) Indeed, we have, for m > 0 :

E 1 (Y (ν) x,y ) m = 1 Γ(m) ∞ 0 E(e -λY (ν)
x,y )λ m-1 dλ and, from (3.2), E(e -λY (ν)

= | -→ m|, for some -→ m ∈ R d .
Here is a first generalization of our Theorem 2.1 (recovered by letting a → 0). Theorem 4.1 : For ν ≥ -1 2 , and a > 0, the process (R t , t ≥ 0), (P (ν,a) x , x ≥ 0) admits an increasing pseudo-inverse.

Sketch of the proof of Theorem 4.1 (It is very similar to that of Theorem 2.1).

We need to show that, for y > x, ∂ ∂t P (ν,a) x (R t ≥ y) ≥ 0. From Fokker-Planck formula, we get : Hence, we deduce from :

∂ ∂t p(t, x, y) = - 1
I ′ ν (z) I ν (z) = I ν+1 (z) I ν (z) + ν z ,
and from (4.2), that : For d = 2(ν + 1) an integer, this process may be constructed as the square of the Euclidian norm of the d-dimensional Ornstein-Uhlenbeck, with parameter β > 0,that is the solution of :

∂ ∂t P (ν,a) x (R t ≥ y) p(t, x, y) = - 1 2y - a 2 I ν+1 (
X t = x 0 + B t + β t 0 X s ds
where (B t ) denotes a d-dimensional Brownian motion. See, e.g., for results about this family of diffusions.

Theorem 4.2 For ν ≥ -1 2 , and β ≥ 0, the process (R t , t ≥ 0) ; Q (ν,β) x , x ≥ 0) admits an increasing pseudo-inverse.

Sketch of the proof of Theorem 4.2 :

We need to show that ∂ ∂t Q (ν,β) x (R t ≥ y) ≥ 0 for every y ≥ x. From the Fokker-Planck formula, we get : ; x < y, ν ≥ 0, 0 ≤ α ≤ 1)

∂ ∂t Q (ν,β) x (R t ≥ y) = -
In Section 2, we introduced the r.v.'s Y

x,y ), corresponding to α = 1, and we shall describe some of their properties. Let us insist again that variables have nothing to do with those introduced in Remark 2 following Theorem 4.2.

5.1 For ease of the reader, we recall some notation and formulae which we have already used and which shall be useful to us below.

T (ν)

x,y : a r.v. whose law, under P (ν)

x , is that of inf{t ≥ 0 ; R t = y} T x,y : a r.v. whose law, under P (ν)

x , is that of sup{t ≥ 0 ; R t = y} G E(e -λG (ν) y ) = 1 2 ν-1 Γ(ν) (y √ 2λ) ν K ν (y √ 2λ) (0 < y, ν > 0)

(5.4)

E(e -λY (ν)

x,y ) =

I ν (x √ 2λ) (x √ 2λ) ν (y √ 2λ) ν+1 K ν+1 (y √ 2λ) (5.5)
From these formulae, we deduced, in Theorem 3.1 :

T (ν) x + Y (ν)
x,y

= G (ν+1) y (5.6)

Y (ν)
x,z

= T (ν)

x,y + Y (ν)

x,z (x < y < z) (5.7)

On the other hand, it is obvious that :

T (ν) x,y + G (ν) y,z (law)
= G (ν)

x,z (x < y < z) (5.8)

Moreover :

ii)The r.v.'s Y (ν,α)

x,y interpolate between Y (ν)

x,y and G (ν)

x,y , i.e. :

Y (ν,1)

x,y

= Y (ν) x,y , Y (ν,0) x,y

= G (ν)

x,y

(5.19)

iii) Y (ν,α) x,y

= Y (ν) x,y + G (ν+1, 1-α) y

(5.20)

and, for (1ν) + ≤ α ≤ 1 (ν,α) x,y

Y
= Y (ν+α-1)

x,y + T (ν, 1-α)

x

(5.21) iv) T (ν) x + Y (ν,α) x,y

= G (ν+α) y (law)

= T (ν+α-1)

x + Y (ν+α-1)

x,y

(5.22) v) G (ν) x,y = G (ν+α, α) y + Y (ν,α) x,y

(5.23)

In particular, for x = 0 we recover (5.12) : = Y (ν+α-1)

G (ν)
x,y + T (ν,1-α)

x (

= Y (ν,α) x,y ) (5.25)

Observe that, from (5.20) and (5.24), the r.v.'s Y (ν,α)

x,y are stochastically decreasing in α.

Proof of Theorem 5.2 i) We show (5.18). We have :

2 1-α Γ(ν + 1) Γ(ν + α)

I ν (x √ 2λ) (x √ 2λ) ν (y √ 2λ) ν+α K ν+α (y √ 2λ) = I ν (x √ 2λ) (x √ 2λ) ν (y √ 2λ) ν+1 K ν+1 (y √ 2λ) • 2 1-α Γ(ν + 1) Γ(ν + α) (y √ 2λ) α-1 K ν+α (y √ 2λ) K ν+1 (y √ 2λ)
= E(e -λY (ν)

x,y ) • E(e -λG (ν+1, 1-α) y ) from (5.5) and (5.9)

This proves (5.18) and (5.20).

ii) We now show (5.21). We have :

E(e -λY (ν,α)

x,y ) = 2 1-α Γ(ν + 1) Γ(ν + α)

I ν (x √ 2λ) (x √ 2λ) ν (y √ 2λ) ν+α K ν+1 (y √ 2λ) = I ν+α-1 (x √ 2λ) (x √ 2λ) ν+α-1 (y √ 2λ) ν+α K ν+α (y √ 2λ)
• 2 1-α Γ(ν + 1) Γ(ν + α)

1 (x √ 2λ) 1-α I ν (x √ 2λ) I ν+α-1 (x √ 2λ)
= E(e -λY (ν+α-1)

x,y

) • E(e -λT (ν, 1-α)

x ) from (5.5) and (5.10)

This shows (5.21).

iii) The relation T (ν)

x + Y (ν,α)

x,y

= G (ν+α) y follows immediately from (5.2), (5.18) and (5.4).

The relation G = T (ν+α-1)

x + Y (ν+α-1)

x,y follows from (5.6) or from the preceding relation where we replace ν by ν+α-1 and we observe, as is obvious, that Y (ν+α-1)

x,y

= Y (ν+α-1, 1)

x,y

. iv) We now show (5.23). We deduce, from (5.8), (5.24) and (5.22) :

T (ν) x + G (ν)
x,y = G (ν+α, α) y + T (ν) x + Y (ν+α) x,y

Hence :

G (ν)
x,y

= G (ν+α, α) y + Y (ν,α) x,y

We might also have proven this last relation by using (5.3), (5.9) and (5.18).

v) The relation (5.19) obviously holds.

vi) Finally, relation (5.25) follows from (5.21) and (5.20). Indeed, from (5.21) and (5.20) :

Y (ν,α)
x,y

= Y (ν+α-1)

x,y + T (ν, 1-α)

x (law)

= T (ν)

x,y + G (ν+1, 1-α) y

x < y, (1ν) + ≤ α ≤ 1 6 remarks : our tools to show existence of pseudo-inverses 6.1 At the beginning of this paper see point iv) of Remark 1.3 we showed that the geometric Brownian motion exp (B t + νt), t ≥ 0 , with ν > 0, admits an increasing pseudo-inverse. On the other hand, Lamperti's representation of geometric Brownian motion is see [RY], Exerc. 1.28, Chap. XI, p. 452 : This formula (6.1) may induce the reader in thinking that Theorem 2.1 -the existence of an increasing pseudo-inverse for R (ν) (u), u ≥ 0 -may be deduced from (6.1) and point iv) of Remark 1.3. This is not at all clear since the clock (C (ν) (t), t ≥ 0) is random.

exp (B t + νt) = R (ν)
6.2 For the proof of Theorem 2.1, we have used, essentially, among the properties of Bessel processes :

• the Fokker-Planck formula, which allowed us to compute ∂ ∂t P (ν) x (R t ≥ y)

• the Pitman-Yor formula (2.15), which allowed to show that I ν+1 I ν (z) ≤ 1, a key point to prove that P (ν)

x (R t ≥ y) is increasing in t

• the formula due to Hirsch-Song, which allowed to give several expressions of the law the r.v. Y (ν)

x,y

for ν ≥ -1 2

• the formula for the resolvent kernel, thanks to which we were able to compute explicitely the Laplace transform of the r.v.'s Y (ν)

x,y .

  a) log y ∈ dt) + P (G y := sup{t ≥ 0 ; B t + at = log y}, T

  We shall now introduce further positive r.v.'s Y (ν,α) x,y ; x < y, ν ≥ 0, 0 ≤ α ≤ 1 which extends the family (Y

  B s + νs) ds and R (ν) (u), u ≥ 0 a Bessel process with index ν.

•

  It is an immediate consequence of (3.2). The fact that the law of Y

	(ν) x,y xy > 1, x > 0 may also be obtained with the help of (2.38). Indeed, we depends only on deduce from (2.36) and (2.3) that : the ratio y x y x
	P	Y xy (ν) x,y	≤ t = P		xy Y (ν) x,y	≥	1 t	=	0	txy	du 2u	y p (ν) (u, x, y) -	x 2 y	p (ν+1) (u, x, y)
	or, equivalently :
	P	xy Y (ν) x,y	≥ t =	0	xy t	du 2u	y p (ν) (u, x, y) -	x 2 y	p (ν+1) (u, x, y)
	Given (2.3), this formula yields the density f xy Y (ν) x,y	of the r.v.	xy Y (ν) x,y

  ∞[ ? Theorem 2.1 implies that it is the case when π is a

	mixture of gamma laws, with parameter • More generally, which are the properties of the process Z t = d 2 ≥ 1 2 but that is not true if 1 2t (y -X t ) + , t ≥ 0 , with d 2 < 1 2 •
	y > 1 which may explain the above increase property ?
	vi) Around infinite divisibility properties for Y	(ν) x,y

  the density with respect to Lebesgue measure of R t under P

	where p(t, x, y), (ν,a) x	, equals,
	from [Wata], p. 117-118 :
	p(t, x, y) = y	I ν (ay) I ν (ax)	e -a 2 t 2 t	exp -	x 2 + y 2 2t	I ν	xy t	(4.3)
	Hence :
	∂ ∂y	p(t, x, y) = p(t, x, y)	1 y	+ a	I ′ ν (ay) I ν (ay)	-	y t	+	x t	I ′ ν ( xy t ) I ν ( xy t )
			2	∂ ∂y	p(t, x, y) +	2ν + 1 2y	+ a	I ν+1 (ay) I ν (ay)	p(t, x, y)	(4.2)

as β ≥ 0 implies e βt ≥ 1 ≥

I ν (z), from Lemma 2.3

Remarks : 1) When d is an integer, the previous computation may be obtained in a simpler manner, using the fact that R t is the square of the norm of a d-dimensional Gaussian variable with mean x 0 exp(βt), and (common) variance e 2βt -1 2β •

2) We also deduce, from the comparison theorem for SDE's, that, for ν ′ ≥ ν and β ′ ≥ β :

hence, with obvious notation, for x ≤ y, and t ≥ 0 :

x,y ≤ t) ≥ P (Y (ν,β) x,y

≤ t)

which states that the r.v.'s. Y 4.3 A third example : the process :

Here, (B t , t ≥ 0) denotes a 1-dimensional Brownian motion, starting from 0. The process (X (ν) t , t ≥ 0) is easily shown to be Markov, since, thanks to Itô's formula, we obtain :

Theorem 4.3 For ν + 1 ≥ 0, the submartingale (X (ν) t , t ≥ 0) admits an increasing pseudo-inverse. Proof of Theorem 4.3 : Using time reversal (from t), we get : for any fixed t ≥ 0,

Now, the process (A (ν) t , t ≥ 0) is increasing, and, if we denote by (τ (ν) y , y ≥ 0) its inverse, we obtain :

We note that, here, we may define the pseudo-inverse process (Y

y , y > 0) of the process (X (ν) t , t ≥ 0) as a "time" process, precisely the process (τ

5.2 Definition of the r.v.'s (G (ν+θ,ν) y , y > 0, ν, θ ≥ 0) and (T (ν+θ,ν) y , y > 0, ν, θ ≥ 0) From Ismail-Kelker [IK], or Pitman-Yor [START_REF] Pitman | Yor Bessel processes and infinitely divisible laws[END_REF], p. 336, formula (9.b.1), , there exists, for every y > 0, ν, θ ≥ 0, a positive r.v. G (ν+θ,ν) y such that :

Likewise from [IK], or [START_REF] Pitman | Yor Bessel processes and infinitely divisible laws[END_REF], p. 336, formula (9.a.1), there exists, for every y > 0, ν, θ ≥ 0, a positive r.v. T

(ν+θ,ν) y such that :

From the relations (5.2) and (5.10) on one hand, and relations (5.4) and (5.9) on the other hand, we immediately deduce the following Proposition which completes the results of Pitman-Yor [PY1] : Proposition 5.1 i) For every y > 0, ν ≥ -1 2 and θ ≥ 0 :

ii) For every y > 0, ν > 0 and θ ≥ 0 :

Of course, it is desirable to give a "more probabilistic" proof of (5.11) and (5.12). Here is such a proof for relation (5.11).

Another proof of (5.11)

To simplify, we take y = 1. Let (R

t , t ≥ 0) and (R (θ-1) t

, t ≥ 0) two Bessel processes starting from 0, independent, with respective indices ν and θ -1 (i.e. with dimension resp. 2ν + 2 and 2θ). Let R

process with index ν + θ, i.e. with dimension 2ν + 2θ + 2, started at 0. Let T (ν+θ) 1

(5.13)

On the other hand, it follows from the intertwining properties of the Bessel semi-groups see [CPY] or [DRVY] that :

with, on the right-hand side of (5.14) r = R (ν+α) t and β ν+1,θ a beta variable with parameter ν + 1 and θ. We then deduce from (5.14) (which is valid for the R (ν+θ) t stopping time T = T (ν+α) 1

) that :

and that T (ν+θ) 1

and (R (ν)

, u ≥ 0) are independent. It then follows from (5.13) that :

x,1 )dx (5.16)

Taking (5.2) into account, we shall have proven (5.11) once we establish :

(5.17) Now (5.17) is easily established, by using the series expansion :

see [Leb], p. 108 and then integrating term by term on the left hand side of (5.17) and using the formula :

Concerning (5.12) a proof close to the preceding one seems harder since G (ν+θ) 1

is not a stopping time. Nonetheless, it may be possible to use the enlarged filtration (G (ν+θ) t , t ≥ 0), i.e. the smallest filtration which contains the natural filtration (R

, t ≥ 0) stopping time. Explicit computations in this new filtration are made possible from the knowledge of the Azéma supermartingale P (G

) -2(ν+θ) ∧ 1, t ≥ 0 . The enlargement formulae allowing to express the process (R (ν+θ) t , t ≥ 0) in this new filtration lead to complicated formulae which do not seem to provide us with the desired explicit result.

5.3

Here is now the main result of this Section 5, i.e. the existence and the properties of the r.v.'s Y (ν,α) x,y ; x < y ; ν ≥ 0, α ∈ [0, 1] . Theorem 5.2 i) For every ν ≥ 0, x < y and α ∈ [0, 1], there exists a positive r.v. Y (ν,α)

x,y such that :

E(e -λY (ν,α) x,y ) = 2 1-α Γ(ν + 1) Γ(ν + α)

(5.18)


Fokouop Nourdi

M.D. Njutapvoui

M.D. Njutapvoui F Onguéné Nourdi
email: njutapvouinourdi@yahoo.fr

Dr Abessolo O Raphael

Dr Jean Grégoire

Rudant Paul

Pr Sakaros Bogning
email: sakaros.bogning@legos.obs-mip.fr

Dr Thomas Stieglitz

Dr Tomedi E Minette

Dr Lecturer, Dada A Olusegun

Pr Edward Anthony
email: anthony@cerege.fr

Donatus Bapentire Angnuureng
email: donatus.angnuureng@ucc.edu.gh

Rafael Almar
email: rafael.almar@ird.fr

Onguéné Raphael

Abessolo O Grégoire

Jean Paul Rudant
email: rudant@univ-mlv.fr

; Tomedi

E Minette

Fernand Verger

F J ; Ngueutchoua

Sédimentologique Carte

Du

Seasonal to decadal scale wave induced shoreline changes along the Cameroonian coastline

Keywords: Coastline, wave climate, Bay of Bonny, optic, radar Coastline, wave climate, Bay of Bonny, optic, radar

The Bay of Bonny along the Cameroon coastline is home to diverse ecosystems. It is under substantial pressure from large human activity, but remains very little studied, like much of the Gulf of Guinea. In order to understand its long-term shoreline variations and the role the wave regime plays in the evolution of the coastline, a study was conducted on the basis of optical image archives from Landsat 5/7/8 and Sentinel-2A/2B satellite missions acquired between October 1986 and May 2020, coupled with daily ERA-Interim wave re-analysis data covering the period from January 1986 to August 2019. Overall, the results show that the evolution of the coast is highly variable in space and time, as indicated by different levels of erosion (30.55 %), and accretion (27.7 %) on the decadal-scale, with the most significant variations occurring in estuarine areas. Nevertheless, 41.75% of the Cameroon shoreline remains stable during the study period. Three main periods during which the coast underwent significant changes a different location were identified, reaching a retreat rate of up to -10 m/year in the northern section during 1986-1994. The wave regime decreases from the southern part towards the north, with an annual trend of significant wave heights of the five segments (-5.6 to -4.1 mm/year) with wave height maxima of 1.46 ± 0.65 m observed in the summer months (July-August). Monthly shoreline changes inversely correlated with the wave climate in some segments but are in line with the seasonal rhythm of the sub-region-s climate. The relatively strong local influence of erosion and accretion caused by natural dynamics (waves, drift currents) associated with the disproportionate presence of estuaries (Wouri, Rio del Rey, Sanaga) accounts for 76.3% of the observed variability, as indicated by analyses of eigenvalue orthogonal decomposition (mode 1, Empirical Orthogonal Functions).

2

wave regime plays in the evolution of the coastline, a study was conducted on the basis of 24 optical image archives from Landsat 5/7/8 and Sentinel-2A/2B satellite missions acquired 25 between October 1986 and May 2020, coupled with daily ERA-Interim wave re-analysis data 26 covering the period from January 1986 to August 2019. Overall, the results show that the 27 evolution of the coast is highly variable in space and time, as indicated by different levels of A series of digital pre-processing was carried out on all the images, ranging from geo-141 referencing, resetting the images of the various sensors, atmospheric and radiometric 142 correction, contrast enhancement, and the application of a cloud mask to eliminate all images 143 that exceed a certain percentage of cloud cover (include a reference that describes this 144 approach). The quality of the satellite images was improved to obtain optimal detection of 145 the coastline by pan-sharpening (combining the visible band (30 m for Landsat and 10 m for 146 Sentinel 2) with the panchromatic band (15 m for Landsat and 10 m for Sentinel). After this 147 step, all bands were reduced to 15 and 10 m for the pixel size of Landsat 7&8 and Sentinel 2A/B respectively, in order to obtain images with better image quality while decreasing as 149 much as possible the digital precision error. 150 Subsequently, image mosaics were created in order to obtain overall images of the study 151 area. A coloured composition of the bands of each sensor was performed (Tab. 1 column 3), 152 in order to choose the most appropriate channel to extract the indicator used as a coastline 153 marker [START_REF] Robin | Télédétection et modélisation du trait de côte et de sa cinématique[END_REF]. Classically, shoreline detection is carried out either in the near 154 infrared or mid-infrared due to the low reflectance of water in the infrared range of the 155 electromagnetic spectrum. On this basis, in order to best highlight the contact between the 156 dry sediments of the high beach and the vegetation cover a coloured composition of the 157 spectral bands was therefore used to discriminate the high beach/vegetation interface 158 [START_REF] Robin | Télédétection et modélisation du trait de côte et de sa cinématique[END_REF]. Following [START_REF] Moore | Long-term cliff retreat and erosion hotspots along the central 620 shores of the Monterey Bay National Marine Sanctuary[END_REF], the total error related to digitizing (Eq. 1) can be To obtain the errors related to tidal sea level oscillations, it is recommended that along coasts 199 with large tidal ranges and/or flatter beach profiles, all shorelines be adjusted to a standard 200 reference elevation [START_REF] Stafford | Air photo survey of coastal erosion[END_REF]. Since all images are acquired between To ground-truth satellite-based data, satellite-derived results were compared with 290 field data at C2 and C5 (Fig. 3). Data form other sections are not available. The comparison 291 indicates that the satellites are able to well identify the main horizontal coastline changes (r 292 = 0.84 for C5 and 0.79 for C2). [1986][1987][1988][1989][1990][1991][1992][1993][1994], of the order of -3 ± 1.4 319 m. Hs peaked in 1990 and subsequently declined with a trend of -0.0047 m/year. Significant 320 accretion shifts began to occur in the second period (1994 to 2005). In the final period the 321 shoreline oscillated between erosion and accretion. In C2, the first two periods are marked 322 by a regressive trend of the order of -4 ± 0.2 m, followed by a period of accretion during the 323 third period (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020). Segment C3 is marked by a pronounced shoreline retreat of -9.4 324 ± 0.7 m over the 1986-1996 period. As in C1, the second period shows an accretion trend of 325 0.8 ± 0.0021m/year, followed by a period of erosion (1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005). Segment C4 is relatively 326 stable with time, marked by some punctual periods of erosion and accretion. In C5 a 327 maximum interannual progradation of the coastline of 10 ± 1.45 m was reached in 1989 328 followed by a long period of average stability of 0.76 ±1.3 m (1992-2007). A continuous regression of the order of -5 ± 0.23 m was observed during the period 2008 to 2015. Linear 330 regression analysis for shoreline changes from 1986 to 2020 showed a general erosion trend 331 of ~ -6.2 m/year during the first period (1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994) and accretion trends during the second period of 2.4 m/year (Fig. 5). The final period is dominated by low-scale erosion of the order 1986, 1990, 2002 and 2019. 350 Interannual coastal changes were correlated to wave climates and the results are

351
shown on Table 4. The variability of wave climate is more correlated with the coastline two periods with erosive dominance regardless of the time scale of observation (C1 and C3).

366

On the other hand, the Kribi zone experienced a first period (1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000) during which its 367 northern part was slightly eroded, followed by a second period (2013-2019) with a 409 Dada et al., 2015;Kuenzer et al., 2014;[START_REF] Obowu | Spatial and multi-temporal change analysis of the Niger 623 Delta coastline using remote sensing and geographic information system (GIS)[END_REF][START_REF] Adegoke | An assessment of 514 recent changes in the Niger Delta coastline using satellite imagery[END_REF] have 410 suggested that the current coastline in West Africa alternates between erosion and accretion.

423

These regions have been identified as the main hot spots on the Cameroonian coast 424 [START_REF] Abessolo Ondoa | Assessment of the evolution of Cameroon coastline: An overview 509 from 1986 to 2015[END_REF], Fossi et al., 2019). And the coastal currents in these estuaries have 425 been identified as responsible for the variability in sediment inputs that affect not only the 426 dynamics of the coastline but also seasonal sediment transport. For example, the eroded area progressive degradation of the coastal system [START_REF] Rovira | Sediment management options for the lower Ebro River and its 641 Delta[END_REF]Batalla et al., 2004; 455 [START_REF] Sanchez-Arcilla | The Ebro Delta: morphodynamics 647 and vulnerability[END_REF].

456

The interannual wave variability obtained in the present study shows a gradual 457 increase in Hs since 1990. The order of magnitude of the increasing trends in mean Hs Ocean in recent decades (e.g. Bertin et al., 2013;[START_REF] Cozannet | 604 Teleconnection pattern influence on sea-wave climate in the Bay of Biscay[END_REF]Dodet et al., 2010).

  28 erosion (30.55 %), and accretion (27.7 %) on the decadal-scale, with the most significant 29 variations occurring in estuarine areas. Nevertheless, 41.75% of the Cameroon shoreline 30 remains stable during the study period. Three main periods during which the coast 31 underwent significant changes a different location were identified, reaching a retreat rate of 32 up to -10 m/year in the northern section during 1986-1994. The wave regime decreases 33 from the southern part towards the north, with an annual trend of significant wave heights 34 of the five segments (-5.6 to -4.1 mm/year) with wave height maxima of 1.46 ± 0.65 m 35 observed in the summer months (July-August). Monthly shoreline changes inversely 36 correlated with the wave climate in some segments but are in line with the seasonal rhythm 37 of the sub-region-s climate. The relatively strong local influence of erosion and accretion 38 caused by natural dynamics (waves, drift currents) associated with the disproportionate 39 presence of estuaries (Wouri, Rio del Rey, Sanaga) accounts for 76.3% of the observed 40 variability, as indicated by analyses of eigenvalue orthogonal decomposition (mode 1, 41 Empirical Orthogonal Functions). These characteristics explain a complex pattern of 42 coastline changes with an almost continuous retreat shoreline during the period 1986-2012 43 and a reversal in trend towards accretion during the period 2012-2020. EOF mode 2 explains 44 23.7% of variability associated with the gradual overall decline from south to north over the 45 study period. This is attributed to an overall decrease in the supply of sediment via the 46 coastal transport system that prevails in the Gulf of Guinea. 1. Introduction at an integrated coastline management strategy, including the production of historical data 72 on coastal dynamics and risks in coastal areas. Giardino et al. 2018 proposed a large-scale 73 coastal policy for the development of sustainable management practices in coastal areas, but 74 this is difficult to implement. 75 Along the Cameroon coastline, few studies have been carried out on coastal dynamics. 76 Abessolo et al. 2018 have studied coastal erosion/accretion processes through the modelling 77 of barotropic currents in the Wouri estuaries in central Cameroon. The recent work of Fotsi 78 et al. 2019 on the temporal (1948-2012) and dynamic evolution of the Wouri estuary's 79 coastline have demonstrated the influence of the sedimentary regime on the variation of its 80 coastline. These works have shown that the Cameroonian coastline has undergone complex 81 patterns of shoreline modification. However, regional-scale sustainable coastal zone 82 management is hampered by a general lack of continuous and accurate historical data on the 83 position of the coastline. It is difficult to study the impact of hydrodynamic forcing on the 84 Cameroonian coastline due to the lack of morphological data. More importantly, the multi-85 scale nature of wave impact: event-driven (impact of storms and waves, and tides); seasonal 86 (seasonal behaviour of the coastline); interannual to decadal trend (difficult to observe due 87 to lack of adequate data) remains almost unknown on this coast. 88 Given the current lack of data and understanding of coastline dynamics on country 89 scale, the aim of this work is to study coastline variations over the last 34 years (1986 to 90 2020) and the role of wave climate in coastal changes in Bonny Bay (Cameroonian coast).

93

  The study area is located between Rio-del-Rey (Cameroon) to Campo (Equatorial 94 Guinea) and covers entire Cameroonian coastline in the Bay of Biafra (Fig.1). This coast is 95 characterized by different ecosystems. Ten major rivers flow into the ocean, the most 96 important of which are Ntem, Nyong, Sanaga, Wouri and the Rio del rey. Following Giresse 97 et al.(1996) four subsets (S-E) of the coast were defined (Fig.1), as follows: 98 -Subset 1: Composed of the Rio del rey and Cross river basins bordering the east of the Bakassi 99 peninsula in Cameroon, this complex is considered as an extension of the eastern zone of the 100 Niger delta set up during the Paleocene-Eocene. It contains the largest zone of mangroves 101 (125,259 hectares) in the country (Ajonina, 2010; MINEPDED-RCM, 2017) and a group of 102 islets. The Rio del Rey's estuary is open and shallow, with tidal amplitudes reaching up to 3 103 m. 104 -Sub-group 2 (Idenau to Limbé): This section is characterized by small basalt cliffs separated 105 by small bays with pebble and sand beaches on the volcanic base of Mount Cameroon. 106 -Subset 3 (Cameroon estuary): with an area of 103,817 ha (ONEQUIP, 2009), it constitutes a 107 vast wetland, 30 km long and almost as wide; it is an integral part of the Wouri watershed. 108 The relief is formed by a succession of sedimentary plains. This bay is also bordered by the 109 second largest mangrove area (93,549 hectares) (MINEPDED-RCM, 2017), and is the 110 receptacle of the main rivers of the region (Wouri, Mungo and Dibamba). The average tidal 111 range is 2 m with maximum reaching 3 m in fast flowing waters. (Onguéné et al., 2015).112 -Subset 4 (Kribi-Campo): From the mouth of the Nyong to the Ntem River, this coast is 113 characterized by the presence of sandy beaches, which alternate with rocky headlands, 114 lagoons, mangroves and other coastal vegetation. It contains two protected areas (Douala-115 Edea Wildlife Reserve and Campo-Ma'an National Park). The tidal range is about 1.8 m at Campo, 1.5 m at Kribi and 1.2 m at Petit-Batanga. (Giresseet al., 1996).

Figure 1 :

 1 Figure 1: The coast of Cameroon in Bay of Bonny. The study area extends from the mouth of

  123

  Mapper), ETM+ (Enhanced Thematic Mapper), OLI (Operational Land Imager), 138 MSI (Multi-spectral instrument); R: red, G: green, B: blue, NIR: near infrared, SWIR: short wave 139 infrared. B1-B12: spectral bands.

  140

159

  To calculate and map the overall coastline displacement rate, the Digital Shoreline 160 Analysis System (DSAS) plug-in version 4.3 (Thieler et al., 2009), in ArcMap software, 161 shoreline detection and extraction was performed manually by scanning the visible 162 vegetation line and using 1986 as the baseline, the rates of shoreline change for 1999, 2005, 163 2010, 2013, 2015, 2016, 2017, 2018, and 2020 were calculated. All extracted coastlines were 164 projected to the World Geodetic System UTM 32N reference system (WGS 84). After selecting 165 the 1986 reference line (baseline), the transects were projected from this line at 500 m 166 intervals which intercepted each coastline initially vectorised. The Net Shoreline Motion 167 (NSM) and End Point Rate (EPR) statistics incorporated in DSAS were calculated to estimate 168 the shoreline variation for the entire study period. We also specified for each coastline a 169 global uncertainty value (+/-15 m) which accounts for the positional and measurement 170 uncertainties. 171 2.3. Landsat and Sentinel-2 image extraction and processing with the CoastSat Toolbox In order to obtain a longer and less constraining time series, to study the across-shore 173 variability of the coastline, we used a stack of 425 images retrieved from the Google Earth 174 Engine (GEE) platform via the open-source CoastSat toolkit (Python version 3.6), that 175 allowed us to obtain time series of shoreline position with a horizontal accuracy of 10 m any 176 sandy coastline worldwide (Vos et al., 2019). This algorithm allows to retrieve TOA (Top-of-177 Atmosphere) reflectance images from the GEE archives of Landsat 5/7/8 and Sentinel-2 178 satellites; and subsequently perform pre-processing of multi-spectral images (cloud 179 masking, pan-sharpening, classification) as described in (Vos et al., 2019). Unlike DSAS, the 180 Coastsat toolbox performs an automatic extraction from a single reference line (baseline) 181 digitized by the user. It then generates the line for each of the sub-pixel resolution images 182 (van der Walt et al., 2014). However, the results only become usable after applying a tidal 183 correction (Eq. 2), which is not implemented in the toolbox as it is specific to each area (Vos 184 et al, 2019). The user has the option of choosing a profile in which the time series of the cros-185 shore distance along the user-defined normal coastal transects is generated. 186 2.4 Estimation of errors associated with the position of the coastline 187 To constrain coastline change rates, we estimated four sources of error based on the 188 method of Hapke et al. 2010, i.e. uncertainties associated with georeferencing, digitization, 189 satellite position and tidal water level uncertainty. The mean values for each uncertainty 190 term are presented in Table 2. 191 The error in digitizing the coastline depends on several factors including the quality of the 192 reference data; the accuracy of the digitization (operator application); the expertise of the 193 operator in interpreting the reference data.

  195 estimated by the sum of the mean ( ) of the offsets recorded during the repetition of the 196 digitization and 2 standard deviations (2σ) : 197 -198

276 Figure 2

 2762 Figure 2 Map of spatio-temporal changes showing erosion/accretion evolution. (a)

Figure 3 :

 3 Figure 3: difference in across-shore displacement between in situ and satellite

  coastline variability and wave regime 297 The results of the seasonal analysis of shoreline changes and waves are shown in Figure 298 4. It is noted that the seasonal variation of the shoreline presents two different phases, an 299 accretion phase in winter (September to February) and an erosion phase in summer (March 300 to August). During this later phase, the maximum erosion is reached in June except in C5 301 where the maximum is reached in August (-2 ± 0.32 m). The wave regime is characterized 302 by summer waves with a significant monthly mean height between 1.2 m and 1.7 m 303 respectively in segments C1 and C5. Summer Hs are highest near C5 and lowest at C1 to C2 304 (northern part). It is during periods of high turbulence that significant erosional disturbance 305 of the shoreline is observed. Note that a delay is generally observed between high Hs values 306 and shoreline movement. During the winter (September to February), Hs values are almost 307 below the annual average in all segment. (Fig. 4).

Figure 4 :

 4 Figure 4: Monthly variation in coastline and waves by segment (C1 to C5). Mean values and

Figure 5 .

 5 Figure 5. Interannual coastline and wave variability by segment (C1 to C5). The coastline

Figure 6

 6 Figure 6 Maps of contributions to shoreline change: (a) Spatial mode EOF (mode1: 76.3% of

  scale analyses of changes of the Cameroonian coastline between 1986 and 398 2020 reveals on one hand a complex and segmented spatial configuration in five areas (C1 399 to C5) and on the other hand three temporal phases characterized by a significant fluctuation 400 of the coastline with considerable variability between segments. 401 Estuaries are characterised by often complex dynamics as a result of interactions of 402 river flow, waves, tides and coastal drift (Anthony et al., 2019; Abessolo et al., 2018; Dada et 403 al., 2016). The strong erosion trends during the first period (1986-1994) in segments C1 and 404 C3 are chiefly attributable to the cutting of mangroves amplified by an increase in sea level, 405 and a resulting reduction in absorbtion of wave energy reaching the coast (Onguéné et al., 406 2015).Dada et al (2015) noted that shoreline retreat from the Niger Delta coastline, at a rate 407 of 28 m/year, between 1950 and 1987 were mainly related to negative anomalies in rainfall and river flows due to climate change. In addition, recent studies (e.g.[START_REF] Almar | Response of the Bight of Benin (Gulf of Guinea, West Africa) 527 coastline to anthropogenic and natural forcing, Part1: Wave climate variability and 528 impacts on the longshore sediment transport[END_REF]

411

  [START_REF] Almar | Response of the Bight of Benin (Gulf of Guinea, West Africa) 527 coastline to anthropogenic and natural forcing, Part1: Wave climate variability and 528 impacts on the longshore sediment transport[END_REF] demonstrated that this can be attributed to two distinct and simultaneous 412 processes: first order, the decrease in southwest swells generated by storms at mid-latitudes 413 and, second order, the increasing contribution of secondary swells southeast (which will 414 cause coastal transport to the west, thereby reducing net coastal transport to the east) 415 generated by positive trade anomalies. However, it remains difficult to infer that the same 416 processes affect the Cameroonian coast at the bottom of the bay of Bonny. 417 The overall coastline displacement rate shows strong variations at the estuaries 418 entrances (C1 and C3), as might be expected. A source of error could be due to the acquisition 419 time of image-derived estuarine shorelines as the position of the sandbank depends on the 420 seasonal variability of river runoff and oceanic forcings (Onguéné, 2015). However, 421 cameroonian estuaries are characterized by the high flows of the Rio del rey rivers (C1), the 422 Wouri estuary, one of the largest in the sub-region (C3 and C4), and the Sanaga river (Fig 2).

427of

  Cape Cameroon (C3) corresponds to the current flowing out of the estuary. Major port 428 works or protection works, located on a coast affected by a predominantly coastal transit, 429 also disturb the balance of sediment transport, by blocking the materials on the "fed" face of 430 the structures, they cause erosion under the "undernourished" face. 431 Segments C2 and C4 (Limbé and Pointe de Souleyba) are characterized by moderate to 432 significant accretion rates of 2.4 m/year and 18.5 m/year, dominated by natural shoreline 433 modification patterns, such as the Souleyba spit (C4) which is relatively recent feature. observed at Souleyba (C4).. It is interesting to note that even if C4 integrates the shore 437 opposite to C3 in the mouth of the Wouri River, this side of the estuary corresponds to the 438 main source of sand supply from the shore in the estuary: the advance estimated at 69 439 m/year in this sector was observed in the period 2010-2019 (Figure 3) and corresponds to 440 a convergence of sediments eroded on the southwestern part of Manoka Island and those 441 sourced from the ocean. The rate of change of coastline displacement remains relatively 442 stable in section C2, as the lithological nature is characterized by the presence of cohesive 443 sediments and rocky coasts as well as a coast dominated by the presence of mangroves and 444 vegetation that absorb wave energy, which are otherwise very low in density in this segment. 445 In the southern part C5, an abrupt change in dynamics is observed as of late 2012, likely 446 due to the construction of a 1.5 km offshore oriented dike at the port of Kribi (30 km to the 447 south of Kribi) in late 2012. The amount of sediment bypassing the dike may not sufficient 448 to balance the high transport potential in this area (Fig. 2a, b). A number of quarries for 449 extracting beach sand for construction purposes are located in the area, from which a volume 450 of more than 45,000 m3/year of sediment is extracted over the entire Kribi coastal strip 451 (MINEPDED, 2017). The distance from the baseline (reference shoreline) shows fluctuating 452 displacements ± 10 m, especially between 2005 and 2020 (Fig. 2a). Such an alteration (with 453 human intervention) of the balance between sediment input and coastal erosion leads to a 454

( 0 .

 0 0017 to 0.003 m/year; Figure5) found in the present study, particularly between 2005 459 and 2020, are in line with the results of many studies carried out in the Gulf of Guinea sub-460 region. Several ocean wave studies have reported an increase in Hs in the North Atlantic 461

462(

  Bertin et al. (2013) observed an increase in Hs over the entire North Atlantic Ocean based 463 on numerical analysis of retrospective winds (1900 to 2008), reaching 0.01 m/yr and 464 gradually decreasing to less than 0.005 m/yr in the southern North Atlantic. The results 465 obtained by Dada et al, 2016 in the coast of the Niger Delta (upper limit of Bonny Bay), are 466 slightly larger but of the same order of magnitude as that obtained in our study. This could 467 be explained by the fact that the work of Dada et al. 2016 was calculated over a much longer 468 period (last 110 years) with a maximum increase in Hs of 0.0054 m/year. The seasonal 469 rhythm of the coastline indicates a period of erosion corresponding to the passage of the 470 ITCZ over the equatorial zone characterized by heavy rainfall. During this period of heavy 471 agitation (May to September), wave heights increase and the coast undergoes significant 472 disturbance. Almar et al. in 2015 demonstrated that the equatorial fluctuation of the 473 Intertropical Convergence Zone (ITCZ) explains most of the variability in wave-induced 474 transport but also the predominant influence of the Southern Annular Mode (SAM) on wave-475 induced transport. 476 The results of the analysis of EOF variability patterns of spatial-temporal shoreline 477 changes show two modes. The first represents the relatively marked local influence of 478 natural dynamics such as waves, drift currents and river flows at C1 (Rio del rey), C2 and C3 479 Wouri estuary) and along the coasts with the presence of coastal infrastructures, such as at 480 C5 (Fig 6a). The second mode represents larger scale changes. These results mainly show 481 that the changes on this coast are due to several factors, and that the wave regime acts in 482 combination with several other factors at the local and regional level.

  the coastline based on a 34-year time series (1986 to 2020) of satellite 485 image data (Earth Explorer and GEE) and the Era-interim re-analysed wave data improved 486 the understanding of coastline variability and local and regional wave climate in the study 487 area. The shoreline change analysis undertaken here reveals a permanent and localized 488 erosive trend on two segments (C1 and C3) over the 402 km along the Cameroonian coast. 489 Nevertheless, a seasonal signal on coastline variability is evident, corroborated by the strong 490 interannual variability in local and regional wave climate. The EOF analysis has highlighted 491 two modes of variability representing respectively the local spatiotemporal scale influence 492 mode 1, but also a larger scale influence mode 2. While this study suggets that shoreline 493 changes are strongly related to local and regional wave climate variability, it is important to 494 note that shoreline stability is also strongly dependent on a constant supply of sediments, 495 and other processes such as: variability of river flows in drift currents, and anthropogenic 496 effects, may act in combination with changes in wave regime and resulting coastal transport 497 to influence shoreline changes along the Cameroonian coast. French Embassy in Cameroon within the framework of a research mobility grant awarded to 501 Njutapvoui Nourdi, as part of a co-supervision agreement between the University of Douala 502 (UDo_EDSFA) and the University Paris Est Marne la Vallée. We used data provided though 503 the ECMWF ERA Interim dataset (www.ECMWF.Int/research/Era), Landsat images 504 downloaded from earthexplorer.usgs.gov, Google Earth Engine (GEE) and Sentinel images 505 from the ESA Copernicus program (https://scihub.copernicus.eu/dhus/).

  506

Table 1 :

 1 Orbiting on the same sun-synchronous trajectory, but with a phase shift of 180° with respect 132 to each other. The revisit time is 10 days for each satellite (Tab. 1), which corresponds to 5 133 days for the couple S-2A and S2-B. A summary of the different characteristics of the images Summary of the satellite images used in this study.

	124	2.2 Landsat and Sentinel-2 Image Extraction and Processing Under DSAS
	125	We extracted a set of 10 images with less than 10% cloud cover (a common problem in
	126	the tropics), i.e. 7 images from Landsat 5/7/8 missions downloaded from
	127	earthexplorer.usgs.gov and 3 images from the Sentinel 2A/2B satellites downloaded from
	128	the Copernicus platform (https://scihub.copernicus.eu/dhus/) of the ESA (European Space
	129	Agency). All Landsat satellites are orbiting at an altitude of 705 km with a swath of 185 km
	130	/ 185 km (Fernand, 2002). The "Sentinel-2" constellation, on the other hand, consists of two

131 135

Table 3

 3

	281	shows small changes (between 0.43 and 0.65 m/year) in the Kribi (C5) and
		Limbé (C2) segment. Significant and very localised changes are observed with an erosion of
	283

282

-11 to -6 m/year respectively in C1 and C3. Section (see above on section, part, segment….) 285

Table 3 :

 3 Area change and overall observed displacement rate by area of interest

	286
			Kribi	Suleybar pic	Cap Cameroon	Limbe	Rio del rey
			(C5)	(C4)	(C3)	(C2)	(C1)
	Surface (km²)	11.5	8.21	6.4	8.24	25.87
	Deplacement rate	0.65 ± 0.04	5.87 ± 0.78	-6.05 ± 1.43	0.43 ± 0.05	-11.601 ±
	(m/yr)					4.02
	Observation (1986-	Stable	strong accretion	erosion	weak	strong
	2020)				accretion	erosion
	287
	288	3.2 Validation of coastline positions derived from satellite imagery with in-situ
	289	measurements

6b). This observation is even more evident on 10-year time scales (Fig.6b). Fig.6bshows

Acknowledgements

499

This project was funded by the service de coopération et d'Action Culturelle (SCAC) of the

dynamics during the first period (1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994), with a correlation coefficient r = 0.5 in 353 segments C1 to C4. While between 2005 and 2020, the wave regime could be responsible for 354 the erosion observed at C3 (r = -0.65). Overall, these results indicate that waves do not fully 355 explain the high variability along the coastline. However, the involvement of waves makes it 356 possible to understand more than 40% of the coastline dynamics in certain segments. An analysis of localized variability shows a dis-proportionality from one area to another.

357

361

Taking into account the local characteristics (nature of the coastal substrate) allows to 362 highlight the influence of waves on the morphological configuration of the coasts (Fig. 6a).

363

However, the trend is reversed when we move to large time scales from 2012 onwards (Fig. in accretion from south to north over the study period (Fig. 6b), highlighting the influence of 391 regional-scale waves on coastline dynamics.

Bosom, E. and Jiménez, J.A., 2011. Probabilistic coastal vulnerability assessment to storms at

Conflict of Interest


Alain Delaët

Turning Catala into a Proof Platform for the Law

THE CATALA TOOLCHAIN

Legal statutes contain descriptions of algorithmic decision processes. In that regard, certain domains of the law are thus more "algorithmic" than others. For instance, tax law describe the computation of the amount of taxes owed by an individual, depending on her income and family situation. More interestingly, these legally-defined algorithms are usually enforced by computer programs called legal expert systems, that government agencies or large organizations have been using for decades. Merigoux et al. [2021] cast light on the issues related to the maintenance and production of these legal expert systems, and their faithfulness with respect to the corresponding legislative specification. To tackle those issues, a novel domain-specific language, Catala, was described and formalized in the same paper. By design, Catala programs closely follow the logical structure of legal statutes. Thus, Catala enables a pair programming and literate programming methodology that raises the level of assurance of legal expert systems, while providing a usable toolchain amenable for production deployments.

Under the hood, Catala relies on an intermediate representation, called default calculus, which is a 𝜆-calculus augmented with default logic terms. Default logic was introduced by Reiter [1987] and is known to be well-suited for legal reasoning, as shown by [START_REF] Lawsky | Formalizing the Code[END_REF][START_REF] Lawsky | A Logic for Statutes[END_REF]: legal statutes define quantities and computations using a base case/exceptions structure often scattered across multiple articles of law. For instance, in the US Internal Revenue Code, sections 61 to 68 define what is considered as income for an individual in general terms, then sections 71 to 291 refine this definition with exemptions, tax cuts and specific inclusions. Encoding such a structure using only conditionals is uneasy, which is why we introduce a default logic term. Formally, this default term is represented as an expression ⟨𝑒 1 , . . . , 𝑒 𝑛 | 𝑒 𝑗𝑢𝑠𝑡 :-𝑒 𝑐𝑜𝑛𝑠 ⟩ where 𝑒 𝑗𝑢𝑠𝑡 :-𝑒 𝑐𝑜𝑛𝑠 is the base case and 𝑒 1 , . . . , 𝑒 𝑛 are the exceptions. The semantics of the default term is as follows: first, each of the exceptions 𝑒 𝑖 is evaluated; if two or more are valid, a conflict error ⊛ is raised. If exactly one exception 𝑒 𝑖 is valid, the final result is 𝑒 𝑖 . If no exception is valid, we evaluate the base case: if 𝑒 𝑗𝑢𝑠𝑡 evaluates to true then the result is the default consequence 𝑒 𝑐𝑜𝑛𝑠 , otherwise the result is ∅. Using this formalism, we can encode sections 61 to 68 as a base case and sections 71 to 291 as exceptions.

A more in-depth description of Catala, including a full formalization of the language and a part of the section 121 of the US Internal Revenue Code, can be found in the original paper of Merigoux et al. [2021].

PROPERTIES OF CATALA PROGRAMS

In this paper, our goal is to verify properties about Catala programs. These properties mostly are of two different kinds, which we present here.

Well-defined executions. We want to prove that the code generated by the Catala compiler is free of crashes. This includes traditional properties like zero-division errors, but is not limited to it. Indeed, the default-logic term introduces interesting new properties to check. First, we need to ensure the absence of conflict (⊛) errors, i.e., no two exceptions are valid inside the same default Authors' addresses: Alain Delaët, Inria Paris, France and ENS Lyon, France, alain.delaet-tixeuil@ens-lyon.fr; Denis Merigoux, Inria Paris, France, denis.merigoux@inria.fr; Aymeric Fromherz, Inria Paris, France, aymeric.fromherz@inria.fr. term during evaluation. In legal terms, it corresponds to the unambiguous of law: no two legal exceptions can occur at the same time. Conversely, we must ensure that the law covers all cases, i.e., every default term must reduce to a value different from ∅.

Those properties are paramount for the use cases that Catala aims to tackle. For instance, we want the automatic enforcement of tax law to be deterministic and applicable to every possible profiles.

Higher-order properties. Additionally, we want to prove more complex and richer properties that reflect the intent of the legislator. For instance, in 2012, the French Supreme Court ruled against marginal tax rates greater than 70% [French Constitutional Council 2012]. The marginal tax rate can be defined as follows: if the income of a household increases by Δ 𝐼 , then its benefits will decrease by Δ 𝐵 and its income tax will increase by Δ 𝑇 . Using these notations, the marginal tax rate is defined as 𝑇 eff = Δ 𝑇 -Δ 𝐵 Δ 𝐼 . In this case, the property we want to guarantee is that, "for any household, 𝑇 eff ≤ 70%".

In the next section, we propose to extend the Catala toolchain to model and verify such properties.

ARCHITECTURE OF THE VERIFICATION FRAMEWORK

Program simplification by refinement. From the two sorts of properties laid out in the previous section, the higher-level ones are of course more interesting. However, proving those properties on a pristine Catala program might be challenging as they require an analysis of the whole program. The complexity of the proof may arise from the corner cases of the default logic structures that obscure the functional meaning of the code, hence significantly raising program verification costs.

But interestingly, if we first prove that the execution of a Catala program is well-defined, then we can rule out most of these default logic corner cases, and perform semantics-preserving rewritings that simplify the program. For instance, by proving the absence of conflict errors in a default term 𝑒, we can remove the runtime mechanism that counts the number of exceptions triggered and simply use nested conditionals to express 𝑒. Thus, we propose a methodology for proving both the well-definedness of executions and higher-level properties on Catala programs. We summarize this methodology in Figure 1. Starting from a Catala program, we apply a series of program rewritings, which must be proved semantics-preserving. By transitivity, the resulting program will be equivalent to the original one, but without most of the complexity due to the default logic encoding. The objective is to simplify the program to enable the proof of higher-level properties inside proof assistants Syntax-directed automated simplification. In a sense, our methodology is similar to program verification efforts that use proof by refinement to phase out a problematic background semantics issue before hitting the functional core of the proof. Examples of this approach abound in the literature [START_REF] Gu | Certikos: An extensible architecture for building certified concurrent {OS} kernels[END_REF][START_REF] Jacob R Lorch | Armada: low-effort verification of high-performance concurrent programs[END_REF]], but so far each simplification step and intermediate program needs to be manually written down and proved. For Catala, we want to leverage the domain-specific nature of the language to propose a heuristics-based, syntax-directed automated procedure to identify potential simplifications in the program, directly embedded as a Catala compiler pass.

𝑃 0 𝑃 1 𝑃 2 ... 𝑃 𝑛-1 𝑃 𝑛 𝑐 1 𝑐 2 𝑐 3 𝑐 𝑛-1 𝑐 𝑛

Verification condition provers

Most heuristics will be focused on the default logic term and its corner cases, as stated earlier. But additional domain-specific program simplifications could be detected and applied using the syntax-directed matching. For instance, the law often defines a concept or quantity piecewise depending on the current date, e.g., minimum salaries indexed on the yearly inflation. In Catala, the piecewise definitions are stitched together in a single default term where each piece is an exception coming with a triggering precondition comparing the current date to fixed interval bounds. Such a pattern can be globally matched and transformed into a more idiomatic interval check, as long as the sub-intervals form a partition of the input space.

Of course, all of these simplifications are contingent on verification conditions justifying their semantic-preserving nature. Therefore, in addition to the syntax-directed simplification pass, the Catala compiler will incorporate a system for collecting and distributing those verification conditions to various proof backends. We sketch such a system in the next section.

CATALA AS A PROOF PLATFORM

As the syntax-directed simplification procedure described in the last section operates, it will generate verification conditions backing the validity of each proposed simplification. Because the proposed simplifications may depend on each other, these verification conditions will have to be dealt with (proved or rejected) in a specific order, prompting the user for proof input. We aim at providing connections from Catala to multiple proof backends : SMT solvers, symbolic execution engines, and proof assistants. Following the methodology of Why3 [START_REF] Filliâtre | Why3-where programs meet provers[END_REF], we want to let the user control how each verification condition should be discharged. Hence, we might design a tactics language that the programmer could use to give hints about the proof strategy for a specific verification condition. For instance, a user could indicate that a specific variable should be treated symbolically (instead of relying on its definition), or provide a cut proof goal and different strategies for the two resulting verification conditions after the cut split.

By allowing the user to control the proof strategies, we aim to reach a tradeoff between a fully automated proof approach that might fail often and leave the user powerless, and a completely manual approach that would likely require a very high effort. We believe that our approach should lead to dozens of verification conditions for simplifying a medium-sized Catala program, and that most of these could be proven by an automated backend with little to no user guidance. This plan might be altered by the data coming back from real-world case studies; overall, we plan to take a pragmatic, application-oriented approach to this Catala proof platform.

Fig. 1 .

 1 Fig. 1. Overview of the proof by refinement architecture, transforming a program 𝑃 0 into a simplified program 𝑃 𝑛 , under verification conditions 𝑐 1 , . . . , 𝑐 𝑛 .