Keywords: Affine Drift Term, Conclusion

En premier lieu, je remercie Jocelyne Erhel, Antoine Lejay et Géraldine Pichot pour avoir dirigé et encadré cette thèse, mais aussi pour m'avoir supporté au sens propre comme au figuré trois années durant. Ils m'ont tout trois laissé une très grande liberté et m'ont accordé une totale confiance pour mener à bien ce travail.

Je tiens ensuite à remercier Monsieur Gilles Pagès pour avoir accepté de présider le jury de cette thèse. Je sais gré à Messieurs Arturo Kohatsu-Higa et Stéphane Menozzi d'avoir rapporté ce travail et de l'intérêt qu'ils ont manifesté à son propos. Enfin, je souhaite remercier Messieurs Mihai Gradinaru et Miguel Martinez d'avoir fait parti de mon jury de soutenance.

Un thèse n'est rien sans une équipe. Je remercie donc l'ensemble des membres de l'équipe SAGE et les participants de l'ANR H2MNO4. Je remercie aussi toutes les personnes présentes à la soutenance de cette thèse. Je vais m'abstenir de citer des noms car la liste est longue et tout oubli serait inconvenant ! Un thèse est aussi l'aboutissement de plusieurs années d'études. Je remercie mes parents pour leur soutien à la poursuite de mes études. C'est un incroyable cadeau qui provient de l'éducation, de l'attention et des soins qu'ils m'ont prodigué et qu'il me prodigue encore, ainsi que de leur travail acharné qui m'a permis de faire des études dans d'excellentes conditions ! Je remercie aussi mes amis d'enfance Guillaume, Stéphane, Khanh et Florian ainsi qu'Hassen et Brahim rencontrés au cours de mes études. Ils m'ont toujours indéfectiblement encouragé.

Table des matières Introduction Présentation générale de l'étude

Nous considérons les processus aléatoires continus dont le générateur infinitésimal possède comme domaine un sous-ensemble de l'espace C b (R n ) des fonctions continues bornées tendant vers 0 aux infinités et d'un opérateur différentiel de la forme

L = 1 2 div(A(x)∇•) + B(x)∇• (1)
où A est une matrice symétrique définie positive et B est un vecteur, tous deux avec des coefficients étant des fonctions bornées de R dans R ayant des discontinuités sur des hyperplans. Dans la suite, nous référerons souvent à A via le terme diffusion et à B via le terme dérive. Ces processus aléatoires sont généralement appelés diffusions biaisées. En dimension un, leur étude commença dès les années dix-neuf cent soixante avec l'introduction du Mouvement Brownien Biaisé par K. Itô et H. P. McKean [START_REF] Itô | Diffusion processes and their sample paths[END_REF]. De nombreux résultats théoriques ont été publiés [START_REF] Harrison | On skew brownian motion[END_REF][START_REF] Gall | Applications du temps local aux équations différentielles stochastiques unidimensionnelles[END_REF][START_REF] Gall | One-dimensional stochastic differential equations involving the local times of the unknown process[END_REF][START_REF] Walsh | A diffusion with discontinuous local time[END_REF][START_REF] Weinryb | Étude d'une équation différentielle stochastique avec temps local[END_REF] concernant ces processus. Cependant, les premiers schémas numériques sont apparus plus tardivement [START_REF] Étoré | Approximation de processus de diffusion à coefficients discontinus en dimension un et appliations à la simulation[END_REF][START_REF] Étoré | A Donsker theorem to simulate one-dimensional processes with measurable coefficients[END_REF][START_REF] Lejay | Simulating diffusion processes in discontinuous media : A numerical scheme with constant time steps[END_REF][START_REF] Uffink | A random walk method for the simulation of macrodispersion in a stratified aquifer. Relation of Groundwater Quantity and Quality[END_REF]. En outre, ils ne permettent que de simuler des processus sans dérive et à coefficients constant par morceaux Lorsque l'on ajoute une dérive tout en gardant l'ensemble des fonctions constantes par morceaux, des schémas numériques permettant de traiter le cas d'une unique discontinuité ont été développé très récemment [START_REF] Étoré | On the existence of a time inhomogeneous skew Brownian motion and some related laws[END_REF][START_REF] Étoré | Exact simulation of one-dimensional stochastic differential equations involving the local time at zero of the unknown process[END_REF]. La principale raison pour laquelle un tel temps a été nécessaire est l'absence complète d'une expression de la fonction de transition lorsque la dérive est présente. En effet, l'expression la plus avancée obtenue [START_REF] Appuhamillage | Occupation and local times for skew Brownian motion with applications to dispersion across an interface[END_REF][START_REF] Étoré | On the existence of a time inhomogeneous skew Brownian motion and some related laws[END_REF] ne concerne que le cas où la dérive est constante.

En dimension supérieure, on trouve là encore quelques résultats théoriques [START_REF] Portenko | Generalized diffusion processes[END_REF][START_REF] Sznitman | A multidimensional process involving local time[END_REF][START_REF] Takanobu | On the existence of solutions of stochastic differential equations with singular drifts[END_REF][START_REF] Takanobu | On the uniqueness of solutions of stochastic differential equations with singular drifts[END_REF], mais nettement moins nombreux. Cependant, il n'y a dans ce cas aucun schéma de simulation hormis [START_REF] Lejay | Simulation of a stochastic process in a discontinuous layered medium[END_REF]. Ce phénomène est principalement dû à l'absence de véritable analogue du temps local en dimension supérieure à un. En effet, il s'agit d'un problème ouvert en théorie du calcul stochastique.

En somme, la simulation des diffusions biaisées est un problème ardu. Ce n'est cependant pas la seule raison qui pousse à l'étudier. En effet, le panel d'applications des schémas numériques de ces processus est impressionnant puisqu'ils interviennent dans les méthodes de Monte 5

Table des matières

Carlo permettant d'approcher des quantités d'intérêts de l'équation d'advection-dispersion.

Cette équation de la physique intervient dans de nombreux domaines comme la biologie, l'écologie, l'imagerie cérébrale, l'astrophysique, la météorologie ou encore la géophysique. Dans ce dernier domaine, l'équation décrit le transport d'un soluté inerte par un fluide en milieux poreux. La dérive correspond alors à la vitesse du fluide dans le milieu et la diffusion à la diffusion moléculaire et à la dispersion cinématique.

La présence de discontinuités au sein de la dérive, mais aussi de la diffusion, est lié à deux phénomènes. Le premier est physique. Il s'agit de la grande hétérogénéité des milieux poreux qui rend la diffusion et la dérive toutes deux discontinues. Le second est lié aux méthodes numériques nécessaires à l'approximation du champ de vitesse. Plus précisément, les méthodes numériques utilisées pour obtenir le champ de vitesse doivent permettre de garder certaines propriétés physiques comme la conservation de la masse. Ceci oblige à déployer par exemple une méthode d'éléments finis mixtes utilisant une base de fonctions dans l'espace de Raviart-Thomas de plus bas ordre qui donne une approximation discontinue.

Dans cette étude, nous proposons une étude théorique des diffusions biaisées ainsi que quelques schémas numériques pour ces processus lorsqu'ils possèdent une diffusion et une dérive constantes par morceaux. Notre travail s'étale sur quatre parties conçues pour être lues indépendamment les unes des autres.

Nous commençons par aborder le problème en dimension un. Dans ce cas, nous procédons tout d'abord à une étude théorique dans le premier chapitre. Cette étude est assez générale puisqu'elle traite des diffusions biaisées à coefficients différentiables par morceaux. Nous y démontrons un résultat remarquable sur l'expression de la densité de la résolvante de ce type de processus. Ce résultat permet notamment d'obtenir une méthode de calcul de ladite densité au moins sous la forme d'une série entière.

Dans le cas des coefficients constants par morceaux ayant une seule discontinuité en 0, la méthode de calcul nous permet de récupérer une expression fermée de la densité de la résolvante. Pour certains jeux de diffusions et de dérives, nous parvenons à l'aide d'une inversion de Laplace à obtenir une expression de la fonction de transition. Cependant, nous n'obtenons pas d'expression de cette fonction dans le cas général.

Dans un second chapitre, nous donnons deux schémas de simulation des diffusions biaisées à coefficients constants par morceaux utilisant non pas la fonction de transition mais la densité de la résolvante. Ils exploitent donc directement les travaux du chapitre 1. Enfin si ces schémas sont construits pour traiter le cas de coefficients possédant une unique discontinuité en 0, leur extension au cas complet est implicite.

Le passage à la dimension supérieure lorsque l'on désire étudier les diffusions biaisées est particulièrement difficile. Nous parvenons dans le cas d'un matrice A diagonale et d'un vecteur B dont aucun coefficient n'est nul à obtenir des fonctionnelles des processus marginaux. Ce travail est fait dans le troisième chapitre. Il constitue un premier travail dans l'optique de prolonger les recherches effectuées dans [START_REF] Lejay | Simulation of a stochastic process in a discontinuous layered medium[END_REF].

Finalement, dans le dernier chapitre, nous considérons une stratégie de parallélisation des méthodes particulaires. En effet, les schémas que nous avons construits dans les deuxième et troisième chapitres pour simuler certaines diffusions biaisées doivent être intégrés au sein de méthodes particulaires pour être utilisés en géophysique ou dans tout autre domaine. La lenteur de ces schémas est alors problématique puisqu'ils sont utilisés un grand nombre de fois pour chaque particule à simuler. En outre, le nombre de particules nécessaires pour une bonne convergence de la méthode particulaire peut être très grand. Comme il est difficile d'accélérer les algorithmes composant les schémas de simulation des diffusions biaisées, ceci explique notre idée de mettre au point un moyen d'accélérer les méthodes particulaires via la parallélisation.

Résumé du chapitre 1 / Expression de la resolvante des diffusions biaisées en dimension un

Nous considérons les processus aléatoires continus dont le générateur infinitésimal possède comme domaine un sous-ensemble Dom(L) de l'espace C b (R) des fonctions continues bornées tendant vers 0 aux infinités et d'un opérateur différentiel de la forme

L = ρ(x) 2 div(a(x)∇•) + b(x)∇• (2)
où a et ρ sont des fonctions mesurables positives et bornées et b est une fonction mesurable et bornée. Ces processus sont souvent appelés diffusions biaisées unidimensionnelles. Nous proposons de faire une étude théorique de ces processus en partant des travaux de [START_REF] Feller | Diffusion processes in one dimension[END_REF][START_REF] Stroock | Diffusion semigroups corresponding to uniformly elliptic divergence form operators[END_REF]. Nous fixons comme objectif de traiter le cas où a, ρ et b sont différentiables par morceaux avec une unique discontinuité en 0. Les idées développées et les résultats obtenus sont les suivants.

En réinterprétant et modifiant la théorie des problèmes de Sturm-Liouville développée dans [START_REF] Kodaira | The eigenvalue problem for ordinary differential equations of the second order and heisenberg's theory of s-matrices[END_REF][START_REF] Titchmarsh | Eigenfunction Expansions, part 1[END_REF][START_REF] Weyl | Über gewöhnliche Differentialgleichungen mit Singularitäten und die zugehörigen Entwicklungen willkürlicher Funktionen[END_REF], nous obtenons deux propositions majeures portant sur la résolvante r(λ, x, y) du générateur infinitésimal (L, Dom(L)). La première est une agrégation de résultats sur les problèmes de Sturm-Liouville obtenus principalement dans [START_REF] Feller | Generalized second order differential operators and their lateral conditions[END_REF][START_REF] Feller | On the intrinsic form for second order differential operator[END_REF][START_REF] Kodaira | The eigenvalue problem for ordinary differential equations of the second order and heisenberg's theory of s-matrices[END_REF].

Proposition. Soit y = 0 fixé. Soit I l'un des intervalles (-∞, 0∧y), (0∧y, y ∨0) et (y ∨0, +∞) et L I l'opérateur dont les coefficients sont ceux de L restreints à I. Alors pour chaque I, il existe u I et u I solutions de (λ -L I )u = 0 continues en x et holomorphes sur le demi-plan complexe H + des nombres de partie réelle strictement positive. De plus, 1. u I pour I = (y ∨ 0, +∞) et u I pour I = (-∞, 0 ∧ y) sont intégrables quel que soit λ ∈ H + , 2. u I et u I sont respectivement croissante et décroissante sur R quel que soit λ réel positif.

La deuxième proposition donne une expression de la densité de la résolvante r(λ, x, y). Elle fournit par la même occasion une nouvelle méthode de calcul de r(λ, x, y). En outre, elle est facilement adaptable au cas des diffusions biaisées présentant de multiples discontinuités ou encore a des diffusions à bord accessibles [START_REF] Feller | Diffusion processes in one dimension[END_REF].

Table des matières

Proposition. Pour y > 0, posons q 1 (λ, x) = u (-∞,0) (λ, x), q 2 (λ, x) = u (0,y) (λ, x), q 3 (λ, x) = u (0,y) (λ, x) et q 4 (λ, x) = u (y,∞) (λ, x), et pour y < 0, q 1 (λ, x) = u (-∞,y) (λ, x), q 2 (λ, x) = u (y,0) (λ, x), q 3 (λ, x) = u (y,0) (λ, x) and q 4 (λ, x) = u (0,∞) (λ, x).

Alors, pour chaque y ∈ R \ {0}, r(λ, x, y) = k 1 (λ, y) q 1 (λ, x)1(x < 0) + k 2 (λ, y) q 2 (λ, x)1(x ∈ (0, y))

+ k 3 (λ, y) q 3 (λ, x)1(x ∈ (0, y)) + k 4 (λ, y) q 4 (λ, x)1(x ≥ y),

pour tout λ ∈ H + et tout x ∈ R avec la fonction k = (k 1 , k 2 , k 3 , k 4 ) solution de
M (λ, y)k(λ, y) T = 0 0 0 1 T , où pour y > 0,

M (λ, y) =    
q 1 (λ, 0) -q 2 (λ, 0) -q 3 (λ, 0) 0 a(0-)q 1 (λ, 0) -a(0+)q 2 (λ, 0) -a(0+)q 3 (λ, 0) 0 0 q 2 (λ, y) q 3 (λ, y) -q 4 (λ, y) 0 ρ(y)a(y) 2 q 2 (λ, y) ρ(y)a(y) 2 q 3 (λ, y) -ρ(y)a(y) 2 q 4 (λ, y)

    ,
alors que pour y < 0,

M (λ, y) =     0 q 2 (
λ, 0) q 3 (λ, 0) -q 4 (λ, 0) 0 a(0-)q 2 (λ, 0) a(0-)q 3 (λ, 0) a(0+)q 4 (λ, 0) q 1 (λ, y) -q 2 (λ, y) -q 3 (λ, y) 0 ρ(y)a(y) 2 q 1 (λ, y) -ρ(y)a(y) 2 q 2 (λ, y) -ρ(y)a(y)

2 q 3 (λ, y) 0     .
En utilisant la méthode de calcul développée, nous obtenons l'expression de la densité de la résolvante pour une diffusion biaisée dont les coefficients (a, ρ, b) sont de la forme

a(x) = β 1(x ≥ 0) + (1 -β) 1(x < 0), ρ(x) = β -1 1(x ≥ 0) + (1 -β) -1 1(x < 0) (3) et b(x) = b + 1(x ≥ 0) + b -1(x < 0) où β ∈ [0, 1].
Cette expression, donnée dans la proposition qui suit, est celle utilisée dans le chapitre suivant conjointement à un changement de variable pour simuler les diffusion biaisées en dimension un.

Proposition. Soit (a, ρ, b) les coefficients de la forme (3) et g(γ, λ, x, y) = 1

γ 2 + 2λ
e (-γ+ √ γ 2 +2λ)(x-y) si x < y, e (-γ- √ γ 2 +2λ)(x-y) si x ≥ y.

Pour y ≥ 0 et λ ∈ C Re λ ≥ 0 ou Im λ = 0, r(λ, x, y) = g(b + , λ, x, y)1(x ≥ 0) + A -+ (λ, y)g(b -, λ, x, y)1(x < 0) + A ++ (λ, y)g(b + , λ, x, -y)1(x ≥ 0) alors que pour y < 0 et λ ∈ C Re λ ≥ 0 ou Im λ = 0, r(λ, x, y) = g(b -, λ, x, y)1(x < 0) + A --(λ, y)g(b -, λ, x, -y)1(x < 0) + A +-(λ, y)g(b + , λ, x, y)1(x ≥ 0) avec

A -+ (λ, y) = Θ(λ, b + , b -) -1 2 β b 2 -+ 2λ e (b + -b -+ √ b 2 -+2λ- √ b 2 + +2λ)y , A ++ (λ, y) = Θ(λ, b + , b -) -1 β( b 2 + + 2λ -b + ) -(1 -β)( b 2 -+ 2λ -b -) e 2b + y , A --(λ, y) = Θ(λ, b + , b -) -1 -β ( b 2 + + 2λ + b + ) + (1 -β) ( b 2 -+ 2λ + b -) e 2b -y , A +-(λ, y) = Θ(λ, b + , b -) -1 2( β -1) b 2 + + 2λe (b --b + + √ b 2 -+2λ- √ b + + 2λ)y ,
où le dénominateur commun est

Θ(λ, b + , b -) = β b 2 + + 2λ + b + + (1 -β) b 2 -+ 2λ -b -.
Résumé du chapitre 2 / Deux schémas de simulation des diffusions biaisées basés sur la résolvante

Nous considérons les processus aléatoires continus dont le générateur infinitésimal est de la forme que celui du chapitre 1 où a et ρ sont des fonctions constantes par morceaux positives et b est une fonction constantes par morceaux. Ces processus sont souvent appelés diffusions biaisées constantes unidimensionnelles.

Nous proposons de construire deux schémas de simulation numérique utilisant la densité de la résolvante calculée dans le premier chapitre et d'établir sa consistance. Les deux schémas sont construits pour des processus n'admettant qu'une seule discontinuité en 0. Ils sont cependant directement extensibles au cas de multiples discontinuités à la différence de [START_REF] Étoré | On the existence of a time inhomogeneous skew Brownian motion and some related laws[END_REF][START_REF] Étoré | Exact simulation of one-dimensional stochastic differential equations involving the local time at zero of the unknown process[END_REF]. Les idées à la base de ces travaux sont les suivantes.

Data: La position X 0 au temps 0. Result: Une trajectoire approchée de X jusqu'au temps T .

for i = 1, • • • , n do Générer X iT n en utilisant la densité n T r( n T , X (i-1)T n , y). end return (X 0 , X iT n , • • • , X T ).
Algorithm 1: Un schéma utilisant la résolvante Soit X = (X t ) t≥0 un processus de Feller conservatif en dimension un dont la résolvante r(λ, x, y) est connue. Alors, ce processus peut être simulé de manière quasi-exacte en utilisant l'Algorithme 1.

La consistance de ce schéma est établie dans la proposition ci-dessous qui établit la convergence en loi de la marche aléatoire construite dans le schéma vers le processus X.

Proposition. Soit (X n t ) 0≤t≤T la suite de marche aléatoires définies pour chaque n ≥ 0 comme l'interpolation linéaire du processus discret en n-étapes donné dans l'Algorithme 1. Alors, (X n t ) n≥0 converge en loi vers X.

A l'aide de ce schéma et de la proposition, nous pouvons donc simuler les diffusions biaisées constantes unidimensionnelles. Nous proposons deux schémas possibles. Le premier est une application directe de l'Algorithme 1. Le second est plus subtil. Il utilise une partie de l'Algorithme 1 et ceci uniquement lorsque le point de discontinuité est franchi.

Plus précisément, le second schéma considère un pas de temps et l'utilise pour simuler une nouvelle position de manière exacte via un mouvement Brownien avec dérive. Si pour atteindre cette nouvelle position le processus a traversé le point de discontinuité, alors le temps qu'il a mis pour atteindre la discontinuité est décompté et il est mis sur la discontinuité. Une fois sur la discontinuité, le processus est simulé sur le restant du pas de temps à l'aide de l'Algorithme 1. Cette idée vient directement du schéma de simulation du Mouvement Brownien Biaisé développé dans [START_REF] Lejay | Simulating diffusion processes in discontinuous media : A numerical scheme with constant time steps[END_REF].

Résumé du chapitre 3 / Quelques fonctionnelles pour les diffusions biaisées multidimensionelles

Nous considérons les diffusions biaisées en dimension n dans le cas où A est une matrice diagonale dont les coefficients admettent une discontinuité sur l'hyperplan formé par les vecteurs dont la dernière coordonnée est nulle et B est un vecteur dont les coefficients sont du même type que ceux de A.

Nous proposons une étude de ces processus par une approche analytique. En effet, le temps local étant difficile à définir en dimension supérieur à un, nous pensons que cette approche est plus à même de permettre la mise au point de schémas numériques ainsi que l'obtention de formule du type Feynman-Kac pour ces processus.

Notre procédons en trois temps. En premier lieu, nous commençons par calculer la transformée de Fourier partielle dans les coordonnées parallèle à l'hyperplan de discontinuité de la résolvante. Puis à l'aide de cette dernière, nous obtenons la résolvante du processus marginal orthogonal à l'hyperplan de discontinuité et identifions ce processus. Enfin, nous calculons la transformée de Laplace de la fonction caractéristique du processus marginal parallèle à l'hyperplan de discontinuité.

Résumé du chapitre 4 / Stratégie de parallélisation des schémas numériques

Nous considérons le problème de la parallélisation des méthodes de Monte Carlo particulaires au sein desquels les processus aléatoires en jeu n'interagissent pas entre eux. Pour y répondre, nous proposons une stratégie de parallélisation qui a été implementée en C++ 2011 [START_REF]The c++ programming language[END_REF] sous la forme d'une librairie [START_REF] Lenôtre | Palmtree library[END_REF] basée sur le générateur de nombres aléatoires RngStream [START_REF] L'ecuyer | An object-oriented random-number package with many long streams and substreams[END_REF]. Les idées principales de cette stratégie sont les suivantes.

Dans un premier temps, nous analysons la manière d'implémenter les méthodes particulaires dans un langage orienté objet. Nous déterminons ensuite un certain nombre de règles. Nous observons alors que le C++ et certaines de ses particularités permettent d'observer ces règles. Ce langage de programmation offre aussi d'autres possibilités dont nous tirons partie. Dans un second temps, nous analysons les différentes manières de paralléliser un générateur de nombres aléatoires. Nous déterminons une bonne manière de faire pour le générateur de nombres aléatoires RngStream [START_REF] L'ecuyer | An object-oriented random-number package with many long streams and substreams[END_REF]. Après cela, nous regardons les différentes méthodes qui permettent de paralléliser les méthodes particulaires. Nous en adoptons une et argumentons notre choix. Notre critère principal est la possibilité de refaire rapidement une partie précise de la simulation. Enfin, nous testons notre implémentation afin de mettre en avant le bon fonctionnement et les performances de notre stratégie.

Introduction

Despite being introduced in the '60, the topic of skew diffusions is now more and more popular for its practical applications in many domains where interfaces and discontinuous coefficients are encountered. Among others, let us cite geophysics, ecology, brain imaging, astrophysics, meteorology or molecular dynamics (See [START_REF] Lejay | Simulating diffusion processes in discontinuous media : A numerical scheme with constant time steps[END_REF][START_REF] Ramirez | Advection-dispersion across interfaces[END_REF] for a list of references).

There are several ways to construct a so-called one-dimensional skew diffusion. One could define it through a scale function and a speed measure. This construction was given for the Skew Brownian motion (SBM) by K. Itô and H.P. McKean [START_REF] Itô | Diffusion processes and their sample paths[END_REF]. Another formulation consider a skew diffusion as the unique strong solutions of stochastic differential equations (SDE) involving local time. This was used by J.M. Harrison and L.A. Shepp [START_REF] Harrison | On skew brownian motion[END_REF] for their construction of the SBM. Among others possibilities (see [START_REF] Lejay | On the constructions of the skew brownian motion[END_REF] for a full panel), N. Portenko choose to see the SBM as a process with a singular drift and give a purely analytic method for its construction (See e.g. [START_REF] Portenko | Generalized diffusion processes[END_REF]).

In the '80, J.-F. Le Gall generalized the notion of Stochastic Differential Equations (SDE) involving local time and provides us, among others, with strong existence and uniqueness results [START_REF] Gall | One-dimensional stochastic differential equations involving the local times of the unknown process[END_REF]. In the same time, J. Groh linked this type of SDEs through scale functions and speed measures [START_REF] Groh | Feller's one-dimensional diffusions as weak solutions to stochastic differential equations[END_REF].

These last studies were of great help for the understanding of diffusions processes whose infinitesimal generators are divergence form operators with piecewise regular coefficients. Indeed, in a series of work [START_REF] Étoré | Approximation de processus de diffusion à coefficients discontinus en dimension un et appliations à la simulation[END_REF][START_REF] Lejay | On the constructions of the skew brownian motion[END_REF][START_REF] Lejay | A scheme for simulating one-dimensional diffusion processes with discontinuous coefficients[END_REF], ..., the authors show that these processes are solutions of SDEs involving local time. This endows their importance to understand the diffusion phenomenon in media with interfaces and therefore justifies the search for numerical simulation algorithms (See [START_REF] Lejay | Simulating diffusion processes in discontinuous media : A numerical scheme with constant time steps[END_REF] and references within).

This article has three goals. Firstly, it aims at completing the study of a large class of SDE with local time by characterizing their infinitesimal generators not only in term of divergence form operators but also of classical differential operators following the approach of W. Feller [START_REF] Feller | Generalized second order differential operators and their lateral conditions[END_REF][START_REF] Feller | On the intrinsic form for second order differential operator[END_REF]. Going back to Sturm-Liouville problem, which was a starting point for W. Feller, we characterize the domain of the operator. This provides us with an alternative formulation to the so-called diffraction problem (see e.g. [START_REF] Ladyženskaja | Equations aux dérivées partielles de type elliptique[END_REF]) with minimal regularity assumption on the coefficients. We also show a stability property of this class of functions under a family of piecewise regular one-to-one space transforms. Such a result was also obtained using the Itô-Tanaka formula in [START_REF] Étoré | On random walk simulation of one-dimensional diffusion processes with discontinuous coefficients[END_REF].

Secondly, using the Weyl-Titchmarsh-Kodaira theory ( [START_REF] Kodaira | The eigenvalue problem for ordinary differential equations of the second order and heisenberg's theory of s-matrices[END_REF][START_REF] Titchmarsh | Eigenfunction Expansions, part 1[END_REF][START_REF] Weyl | Über gewöhnliche Differentialgleichungen mit Singularitäten und die zugehörigen Entwicklungen willkürlicher Funktionen[END_REF] and an analytic machinery, we provide a closed-form for the resolvent kernel of the operator in free space. Since this approach relies on a appropriate combination of known resolvents in free space, it can be easily generalized to deal with boundary conditions and solve related problems such as the elastic skew Brownian motion and the Feynman-Kac formula. This also leads to expressions for the Laplace transform of the local time at 0 or the occupation time of some segment or half-line. Although our framework presents similarities with the one developed by B. Gaveau et al. [START_REF] Gaveau | Second order differential operators and dirichlet integrals with singular coefficients, i. functional calculus of one-dimensional operators[END_REF], we differ by our application of the Weyl-Titchmarsh-Kodaira theory. Besides, the drift term, which was our primarily purpose, can only be taken into account in [START_REF] Gaveau | Second order differential operators and dirichlet integrals with singular coefficients, i. functional calculus of one-dimensional operators[END_REF] at the price of using limiting arguments on rather complex formulae.

We point out that obtaining explicit expressions for the resolvent or the density is of higher importance to set-up numerical schemes as in e.g. [START_REF] Lejay | Simulating diffusion processes in discontinuous media : A numerical scheme with constant time steps[END_REF], since the dynamic of a process solution to a SDE with local time is difficult to apprehend.

Finally, we consider the skew Brownian motion with a possibly piecewise constant drift which takes two values b + on R + and b -on R -. For several practical situations, we give an explicit closed-form of the density. Although it seems hard to perform this inversion for general values of b + and b -, we provide an explicit expression for a new case, called constant Péclet, which possesses a natural interpretation in fluid mechanics.

While the density of the SBM is easily found [START_REF] Walsh | A diffusion with discontinuous local time[END_REF], a heavy probabilistic machinery was set up for the same purpose when some particular values of the drift are encountered [START_REF] Appuhamillage | Occupation and local times for skew Brownian motion with applications to dispersion across an interface[END_REF][START_REF] Étoré | Exact simulation of one-dimensional stochastic differential equations involving the local time at zero of the unknown process[END_REF][START_REF] Fernholz | Two brownian particles with rank-based characteristics and skew-elastic collisions[END_REF][START_REF] Gairat | Density of skew Brownian motion and its functionals with application in finance[END_REF]. We consider our work as a first incursion in alternative methods to overcome the difficulties raised by approaches based on stochastic calculus.

Outline

In Section 1.1, we show the link between processes generated by divergence form operators, those defined in the sense of Feller through their scale functions and their speed measures, and SDE with local time. We also exhibit the behavior of the resolvent kernel. In Section 1.2, we provide a method for the explicit computation of the resolvent. Finally, in Section 1.3, densities are computed through the Laplace inversion of resolvents in some particular cases.

Operators, semi-groups and diffusion processes

This section is mostly theoretical and pursues multiple goals: 1. To identify the stochastic process X generated by a divergence form operator (in the sense defined by D. Stroock in [START_REF] Stroock | Diffusion semigroups corresponding to uniformly elliptic divergence form operators[END_REF]) with the one defined through its scale function and speed measure. This complements and extends the results obtained in [START_REF] Lejay | A scheme for simulating one-dimensional diffusion processes with discontinuous coefficients[END_REF]. 2. For coefficients that are piecewise smooth with discontinuities of the first kind, to identify the domain of the generator of X = (X t ) t≥0 and state a result on the regularity of the function given by the probabilistic representation

(t, x) → E x [f (X t )]
. This extends the result provided in [START_REF] Martinez | One-dimensional parabolic diffraction equations : pointwise estimates and discretization of related stochastic differential equations with weighted local times[END_REF]. 3. To obtain some analytic properties on the resolvent kernel associated with the process X. 4. To relate the process X when its coefficients are as in 2. above with an SDE involving local time in the sense defined by J.-F. Le Gall in [START_REF] Gall | One-dimensional stochastic differential equations involving the local times of the unknown process[END_REF]. When I is a union of disjoint open abutting intervals, we mean by f in C k (I, R) a function of class C k in each interval composing I such that f and its derivatives up to order k have limits at the inner separation of the intervals.

We set R = [-∞, ∞] and for K = R or C,

C b (R, K) = {f : R → K | f is continuous and bounded}, C 0 (R, K) = {f ∈ C b (R, K) | f (±∞) = 0}.
Any continuous function f : R → K vanishing at infinity is naturally identified with a function in C 0 (R, K).

Let 0 < λ < Λ be two fixed real numbers. We consider the sets

C = {a = (a, ρ, b) measurable from R 3 → R | λ ≤ a, ρ ≤ Λ and |b| ≤ Λ}, S = {a = (a, ρ, b) ∈ C | a, ρ, b ∈ C ∞ (R, R)} and P = (a, ρ, b) ∈ C For a set J without cluster points,a ∈ C 1 (R \ J, R), b, ρ ∈ C(R \ J, R) .

Diffusion process associated with a divergence-form operator

Here we construct a Feller process X associated with a divergence form operator ρ 2 ∇(a∇•) + b∇• when a is discontinuous. Unlike the solutions of SDEs, such a process is not always a semimartingale. In Subsection 1.1.6, we however characterize the situations where X solves a SDE involving local time.

This section relies on [START_REF] Stroock | Diffusion semigroups corresponding to uniformly elliptic divergence form operators[END_REF] (See also [START_REF] Lejay | A scheme for simulating one-dimensional diffusion processes with discontinuous coefficients[END_REF]). We use the framework of the variational theory of PDE where the solutions are primary sought in subspaces of L 2 (R), the space of square integrable functions on R. We base our construction on a Gaussian estimate on the fundamental solution.

To each a = (a, ρ, b) ∈ C is associated a family of bilinear form

E α (u, v) = 1 2 a(x)∇u(x)∇v(x) dx -b(x)∇u(x)v(x) dx ρ(x) + α u(x)v(x) dx ρ(x) , ∀(u, v) ∈ H 1 (R) × H 1 (R),
where H 1 (R) is the Sobolev space of square of functions f with weak derivatives ∇f such that f, ∇f ∈ L 2 (R).

Clearly, E α is a Gårding form for α ≥ α 0 large enough (depending only on λ and Λ). By an application of the Lax-Milgram theorem, one may associate with (E α ) α>0 a resolvent (G α ) α≥α 0 on L 2 (R).

One-dimensional skew diffusions: explicit expressions of densities and resolvent

When the coefficients a are smooth, any solution u to the variational problem

E α (u, v) = f, v 1/ρ def = R f (x)v(x) dx ρ(x)
is smooth when f is. As a result, u becomes a classical solution to

(α -L)u(x) = f (x) with Lu(x) = 1 2 ρ(x)∇(a(x)∇u(x)) + b(x)∇u(x) (1.1.1)
and the domain Dom(L) = G α (L 2 (R)) for any α > 0 where G α = α -L) -1 is the resolvent. But it works even for non-smooth coefficients. Indeed, such a (L, Dom(L)) exists and satisfies

(α -L)u, v 1/ρ = E α (u, v), ∀u ∈ Dom(L), ∀v ∈ H 1 (R).
From the Hille-Yosida theorem, the operator L is the infinitesimal generator of a contraction semi-group (P t ) t≥0 . A result from J. Nash and J. Moser extended by D. Aronson [START_REF] Aronson | Non-negative solutions of linear parabolic equations[END_REF] and D. Stroock [START_REF] Stroock | Diffusion semigroups corresponding to uniformly elliptic divergence form operators[END_REF] shows that there exists a density q with respect to the measure dx/ρ(x) such that for any

f ∈ L 2 (R) ∩ L ∞ (R), u(t, x) = R q(t, x, y)f (y) dy ρ(y) , t > 0, x ∈ R is a version of P t f (x) ∈ L 2 (R)
. Moreover, q is Hölder continuous and satisfies a Gaussian upper and lower bounds. The main feature is that the constants appearing in the bound and in the Hölder modulus of continuity depend only on λ and Λ. This proves the existence of a diffusion as well as some stability properties (See [START_REF] Étoré | A Donsker theorem to simulate one-dimensional processes with measurable coefficients[END_REF][START_REF] Lejay | A scheme for simulating one-dimensional diffusion processes with discontinuous coefficients[END_REF][START_REF] Rozkosz | Weak convergence of diffusions corresponding to divergence form operators[END_REF] Of course, for this diffusion process X a ,

P X a t f (x) def = E x [f (X a t )] = R p(t, x, y)f (y) dy a.e. = P t f (x), ∀t ≥ 0, (1.1.3) and G X a α f (x) def = +∞ 0 e -αt E x [f (X a t )] dt a.e. = G α f (x), ∀α > 0, (1.1.4)
for any x ∈ R and f ∈ L 2 (R). Even if P t and G α are only operators from L 2 (R) to L 2 (R), we now identify P t and G α to P X a t and G X a α . These operators P X a t and G X a α also act on C b (R, R). Since X a is a Feller process, P X a t and G X a α also map C 0 (R, R) to C 0 (R, R).

The infinitesimal generator of X a

In the previous Subsection, a process X a has been associated with a ∈ C. The domain of its infinitesimal generator may be described as the range of the resolvent G α (L 2 (R)). Here we give an alternative description of the infinitesimal generator of X a using the approach from W. Feller on scale functions and speed measures where the underlying Banach space is the one of continuous bounded functions [START_REF] Feller | Generalized second order differential operators and their lateral conditions[END_REF][START_REF] Feller | On the intrinsic form for second order differential operator[END_REF].

We point out that every one-dimensional regular diffusion process Y is characterized by its scale function S and its speed measure M (for classical references on this subject, see e.g. [START_REF] Breiman | of Classics in Applied Mathematics[END_REF][START_REF] Rogers | Diffusions, Markov processes, and martingales[END_REF][START_REF] Mandl | Analytical treatment of one-dimensional Markov processes, volume 151 of Die Grundlehren der mathematischen Wissenschaften[END_REF]). We also recall that the speed measure M is a right-continuous, increasing function, which can be identified with a measure on R and the scale function S, which is unique up to an affine transformation, is a convex, continuous and increasing function.

For a continuous, increasing function g : R → R, let us define

D + g f (x) = lim y x f (y+) -f (x-) g(y) -g(x) and D - g f (x) = lim y x f (x+) -f (y-) g(y) -g(x)
.

(1.1.5)

From now, we only consider the case of a diffusion Y having R as a state of space and possessing a continuous speed measure M . The infinitesimal generator of such a process Y associated with (S, M ) is (1.1.8)

Af = 1 2 D M D S f where D S f (x) = D + S f (x) = D - S f (x) and D M = D + M = D - M whose domain is [23, Theorem 8.2] for R = [-∞, ∞], Dom(A) = f ∈ C b (R, R) D M D S f ∈ C b (R, R) . ( 1
Thus D M f (x) = m(x) -1 D + f (x) and D S f (x) = a(x) exp(h(x))D + f (x)
, where D + f is the right derivative operator.

Proposition 1.1.2. For a ∈ C, the infinitesimal generator of X a is A = 1 2 D M D S with (S, M ) defined by (1.1.7)-(1.1.8) and Dom(A) given by (1.1.6).

Remark 1.1.1. The infinitesimal generator 1 2 D M D S is the same for (a, ρ, b) ∈ C as well as for ( a, ρ, 0) with a(x) = a(x) exp(h(x)) and ρ(x) = exp(-h(x))ρ(x). However ( a, ρ, 0) does not necessarily belong to C as the coefficients could decrease to 0 or increase to ∞ for large values of |x|. But, since h is continuous whatever b, this explains why the domain of A does not depend on the regularity of b.

Proof. Let Y a be the diffusion process with infinitesimal generator A. When a ∈ S, it is easily seen that Y a is equal in distribution to X a associated with a using Proposition 1.1.1. From (1.1.7)-(1.1.8), we deduce that the pointwise convergence of a n to a implies the one of the associated scale function and speed measure. Thus, (Y an ) n≥0 converges in distribution to Y a through the results of [START_REF] Freidlin | Necessary and sufficient conditions for weak convergence of one-dimensional Markov processes[END_REF] on the convergence of the scale function and speed measure. Now Proposition 1.1.1 implies that (X an ) n≥0 converges in distribution to X a . Hence X a is equal in distribution to Y a by using for a n a regularization of the coefficients.

The resolvent kernel

Let us assume that a ∈ C. Throughout this section, we consider the scale function S and the speed measure M given by (1.1.7) and (1.1.8), and A the associated operator A = 1 2 D M D S with domain given by (1.1.6).

For an interval I of R, let AC(I, C) be the space of absolutely continuous functions from I to C and S a (I, C) the space of functions u in AC(I, C) with D S u ∈ AC(I, C). The next proposition concerns the existence and the regularity of solutions to (λ -A)u = 0. This proposition is the cornerstone of the Sturm-Liouville theory and its extension by H. Weyl, K. Kodaira and E.C. Titchmarsh using complex analysis. Proof. For any λ ∈ C, any solution Z : R → C 2 to the first order differential equation

Z (x) = 0 s(x) 2λm(x) 0 Z(x), Z(0) = α β ∈ C 2 for any x ∈ R (1.1.9)
is such that (λ -A)u(x) = 0 when Z(x) = (u(x), D S u(x)). Regarding the spatial regularity of the solution, let u and f be continuous functions from R to C such that (λ -A)u = f . Then integrating between x and y against the measure M leads to

λ y x u(z)m(z) dz -D S u(x) + D S u(y) = 2 y x f (z)m(z) dz. (1.1.10)
In particular, D S u ∈ AC(R, C) and u ∈ S a (R, C).

The Wronskian of two functions u and v in S a (R, C) is

Wr[u, v](x) def = u(x)D S v(x) -v(x)D S u(x).
If u and v solve (λ -A)u = 0, then Wr[u, v](x) is constant over x (See (1.1.14) below). The solutions u and v are independent if and only Wr[u, v](x) = 0. In this case, any solution to (λ -A)w = 0 on some interval I is a linear combination of u and v on I.

We introduce

H + = {λ ∈ C | Re(λ) > 0}, h + = {λ ∈ C | Im(λ) = 0, Re(λ) > 0}, and h -= {λ ∈ C | Im(λ) = 0, Re(λ) ≤ 0}. Lemma 1.1.1. Let τ x = inf{t > 0 | X a t =
x} be the first hitting time of x ∈ R of the process X a . For any y < x, the function

φ(λ, x) = E x [exp(-λτ y )] is holomorphic on H + . Besides, for any λ ∈ H + , φ(λ, x) is solution to (λ -A)φ(λ, x) = 0 for x ≥ y and φ(λ, y) = 1. Proof. Whatever λ ∈ H + , 0 ≤ |φ(λ, x)| ≤ φ(Re(λ), x) ≤ 1 so that φ is well defined on H + . In addition, from the Feynman-Kac formula [37, § 2.6], (λ -A)φ(λ, x) = 0, ∀x ≥ y.
We decide to show that φ(λ, x) can be developed as a power series on a neighborhood of any λ ∈ H + since it is equivalent of being holomorphic for complex function.

Notice that φ(λ, x) = 1 -λE x τy 0 exp(-λs) ds , Now using the series representation of the exponential, for σ ∈ C,

φ(λ + σ, x) = 1 -λE x τy 0 exp(-λs) exp(-σs) ds -σE x τy 0 exp(-λs) exp(-σs) ds = φ(λ, x) - k≥1 (-1) k σ k c k (λ, x) with c k (λ, x) = E x τy 0 λ s k k! - s k-1 (k -1)! exp(-λs) ds . (1.1.11) Since +∞ 0 s k exp(-λs) ds = k!/λ k+1 and τ y ≥ 0, |c k (λ, x)| ≤ 1+1/|λ| for each k.
Hence, the series in the right-hand side of (1.1.11) converges absolutely for any σ ∈ H + , |σ| < min 1, Re λ and the result is proved.

The function λ → φ(λ, x) is not necessarily analytic around λ = 0. Otherwise, it would mean that τ y has finite moments of all orders. This is not true for example for the Brownian motion.

Proposition 1.1.4. For any λ ∈ H + , the functions

u (λ, x) = E x [exp(-λτ 0 )], u (λ, x) = 1 E 0 [exp(-λτ x )] if x < 0, (1.1.12) u (λ, x) = 1 E 0 [exp(-λτ x )] , u (λ, x) = E x [exp(-λτ 0 )] if x > 0 (1.1.13) are solutions to (λ -A)u = 0. For λ on the half-line h + , x → u (λ, x) is increasing from 0 to +∞, x → u (λ, x) is decreasing from +∞ to 0. Moreover, u (λ, 0) = u (λ, 0) = 1. For λ ∈ H + , |u(λ, x)| converges to 0 as |x| → ∞. Besides, for any x ∈ R, λ → u(λ, x) and λ → D S u(λ, x) for u = u , u are holomorphic on H + .
Proof. We follow the classical probabilistic construction of u and u [START_REF] Breiman | of Classics in Applied Mathematics[END_REF][START_REF] Rogers | Diffusions, Markov processes, and martingales[END_REF] for λ ∈ R, which extends to λ ∈ H + . In the order, 1. From the Feynman-Kac formula, u and u are solutions to (λ -A)u = 0 with u(0) = 1.

Using the strong Markov property

, E x [e -λτy ] = E x [e -λτz ]E z [e -λτy ] for any x ≤ y ≤ z. From this, one easily checks that E x [e -λτy ] = u(x)/u(z) for u = u (λ, •) or u = u (λ, •)
. Since these functions are positive, Lemma 1.1.1 implies that they are holomorphic on H + .

3. Using for a basis of solutions for (λ -A)u = 0 the functions u 1 and

u 2 with u 1 (λ, 0) = D S u 2 (λ, 0) = 1 and D S u 1 (λ, 0) = u 2 (λ, 0) = 0, u may be written u (λ, x) = u 1 (λ, x) + α(λ)u 2 (λ, x). Since u 1 and u 2 are holomorphic from Proposition 1.1.3, α(λ) is meromorphic. Since u (•, x) is holomorphic on H + from Lemma 1.1.1, α is necessarily holomorphic on H + . With Proposition 1.1.3, D S u (•, x) is holomorphic on H + for any x ∈ R since α, D S u 1 (•, x)
and D S u 2 (•, x) are holomorphic on H + (actually, α(λ) = D S u(λ, 0) from our choice of u 1 and u 2 ). The same argument holds for u . We have then proved that D S u and D S u are holomorphic on H + .

4. For λ ∈ h + , it is easily checked that the functions u and u defined by (1.1.12) and (1.1.13) are positive, monotone, and range respectively from 0 to ∞ and from ∞ to 0.

Corollary 1.1.1. For λ ∈ H + , let u in S a (R + , C) (resp. S a (R -, C)) be a solution to (λ-A)u = 0 on R + (resp. on R -) which vanishes at +∞ (resp. -∞) and satisfies u(0) = 0. Then u = 0.
Proof. For a continuous, bounded function f on R + , set

H λ f (x) def = E x τ 0 0 e -λs f (X a s ) ds .
Using the strong Markov property of X a ,

H λ f (x) = G X a λ f (x) -u (λ, x)G X a λ f (0). With the Feynman-Kac formula, (λ -A)H λ f (x) = f (x) on R + . In addition, H λ f (0) = 0.
Finally, H λ maps bounded functions to bounded functions. Besides, if f vanishes at infinity, then G X a λ f and u (λ, •) vanish at infinity, so does H λ f . This proves that H λ is a bounded inverse of λ -A for λ ∈ H + . Therefore, there exists a unique solution to (λ -A)u = f when f is bounded, continuous on R + , and this solution vanishes at infinity when f does. Hence, for f = 0, u = 0.

Remark 1.1.2. Proposition 1.1.4 and Corollary 1.1.1 could have been treated by a purely analytic method with the theory of limit circle and limit points that was initiated by H. Weyl [START_REF] Weyl | Über gewöhnliche Differentialgleichungen mit Singularitäten und die zugehörigen Entwicklungen willkürlicher Funktionen[END_REF] (See also [START_REF] Kodaira | The eigenvalue problem for ordinary differential equations of the second order and heisenberg's theory of s-matrices[END_REF][START_REF] Titchmarsh | Eigenfunction Expansions, part 1[END_REF]). It could be proved that u and u are actually meromorphic on C \ h -.

The functions u and u serve to build the resolvent. For our purpose, it is important that they are holomorphic on H + so that we could identify the resolvent kernel with another analytic expression to recover the densities through Laplace inversion (See Proposition 1.2.1 below).

The holomorphic property is the key to prove results such as Proposition 1.1.5 below on the regularity of the density (See also [START_REF] Hille | The abstract Cauchy problem and Cauchy's problem for parabolic differential equations[END_REF]) as well as the existence of a diffusion process. This is not our purpose here. Therefore, we took benefit of the known probabilistic representation of the resolvent to refine our knowledge on the resolvent kernel.

[φ(x)Aψ(x) -ψ(x)Aφ(x)]m(x) dx = 1 2 Wr[φ, ψ](b) - 1 2 Wr[φ, ψ](a). (1.1.14) Since u (λ, 0) is solution to (λ-A)u = 0, and when u (λ, •) = u (λ, •), (1.1.14) with φ = u and ψ = u yields 2i Im b a λu (λ, x)u (λ, x)m(x) dx = 0 or Im(λ) = 0.
Therefore, Wr[u , u ](λ, x) = 0 implies that Im(λ) = 0. On the other hand, since u (λ, •) increases whereas u (λ, •) decreases for λ ∈ h + while both functions are positive implies that

Wr[u , u ](λ, x) = 0 only on h -.
We give two examples that are of great importance later.

Example

1.1.1. For s(x) = m(x) = 1 so that A = 1 2 D 2 x , u (λ, x) = exp( √ 2λx) and u (λ, x) = exp(- √ 2λx).

They are holomorphic on

C \ h -. Besides, Wr[u , u ](λ, x) = 0 if and only if λ = 0. Example 1.1.2. For s(x) = exp(-γx) and m(x) = exp(γx) so that A = 1 2 D 2 x + γD x , u (γ, λ, x) = exp((-γ + γ 2 + 2λ)x)
and u (γ, λ, x) = exp((-γγ 2 + 2λ)x).

(1.1.15)

They are holomorphic on C \ {λ ∈ C | Im(λ) = 0, Re(λ) < -|γ|/2}. Besides, Wr[u , u ](λ, x) = γ 2 + 2λ and Wr[u , u ](λ, x) = 0 if and only if λ = -γ 2 /2.
Let us set for λ ∈ H + and the density m of the speed measure M ,

r(λ, x, y) = m(y) Wr[u , u ](λ, y) u (λ, x)u (λ, y) if x < y, u (λ, x)u (λ, y) if x > y. The resolvent R λ = G X a λ |C 0 (R,R
) can be extended for λ ∈ H + as a bounded (by 1/ Re(λ)) operator. The resolvent is actually (λ -A) -1 for λ ∈ H + . From standard results (See e.g. [START_REF] Breiman | of Classics in Applied Mathematics[END_REF]Theorem 16.75] or [START_REF] Feller | Generalized second order differential operators and their lateral conditions[END_REF]), for any λ ∈ H + ,

R λ f (x) = R r(λ, x, y)f (y) dy, f ∈ C 0 (R, R).
Exploiting (1.1.10), we easily see that r(λ, x, y) satisfies D S r(λ, y-, y) -D S r(λ, y+, y) = 2m(y), for any y ∈ R, λ ∈ H + .

(1.1.16)

The density p(t, x, y) of the process X a generated by (A, Dom(A)) or equivalently by (L, Dom(L)) is related to r(λ, x, y) by the Laplace transform L. The next proposition is due to H.P. McKean. It is actually valid for any process with generators of type 1 2 D M D S . Proposition's proof exploits the fact that (A, Dom(A)) is a self-adjoint operator to provide the existence of a density p which is analytic 1 .

1. approach, L is self-adjoint and the semi-group (P t ) t>0 is analytic and maps L 2 (R) to Dom(L). The nontrivial part is the existence of a density. In [START_REF] Mckean | Elementary solutions for certain parabolic partial differential equations[END_REF], H.P. McKean uses a spectral representation on the kernel of the resolvent to show the existence of a density transition function and analyse its properties, in particular to show that p(t, •, y) belongs to Dom(A n ) for any n ∈ N. These two approaches are parallel, one at the level of the operators, the second on their kernels. Many of the properties of p exhibited in [START_REF] Mckean | Elementary solutions for certain parabolic partial differential equations[END_REF] are not quoted in the subsequent literature such as [START_REF] Breiman | of Classics in Applied Mathematics[END_REF][START_REF] Itô | Diffusion processes and their sample paths[END_REF].

Proposition 1.1.5 (H.P. McKean [START_REF] Mckean | Elementary solutions for certain parabolic partial differential equations[END_REF]). The density p(t, x, y) of the stochastic process X a for a ∈ C satisfies p(t, x, y) = L -1 (r(λ, x, y)). In addition, for any n ≥ 1, ∂ n t p exists and is continuous on p defined by (1.1.2). In addition, it satisfies ∂ n t p = A n p on p and ∂ n t p(t, •, y) vanishes at infinity for any t > 0, y ∈ R.

Remark 1.1.3. The above theorem is done with S(x) = x. Nevertheless we could always reduce to this case (see also the results in Section 1.1.5).

Diffusion with piecewise regular coefficients

We now focus on the process associated with a ∈ P, that is a Feller process with piecewise regular coefficients. Notation 1.1.1. For K = R or C, given a ∈ P and I an interval of R, we denote by A a (I, K)

the set of functions u ∈ C(I, K) ∩ C 2 (I \ Dis(a), K) such that a(x-)u (x-) = a(x+)u (x+), ∀x ∈ I ∩ Dis(a).
(1.1.17)

We show how close the domains of A and L are.

Proposition 1.1.6. For a ∈ P,

Dom(A) ⊂ A a (R, R) and A a (R, R) ∩ {f ∈ L 2 (R) | Af ∈ L 2 (R)} ⊂ Dom(L).
Besides, the latter inclusion is dense in Dom(L) with respect to the L 2 -norm.

Proof. The inclusion Dom(A) ⊂ A a (R, R) is straightforward using the fact that Dom(A) ⊂ S a (R, R) and the continuity of m and s out of Dis(a). Let u ∈ A a (R, R). The set R \ Dis(a) may be expressed as the disjoint union of intervals of type

(x i , x i+1 ) with Dis(a) = {x i } i∈J with x i < x i+1 , i ∈ J ⊂ Z. Thus, for φ ∈ C 1 c (R, R), 2 R Au(x)φ(x) dx ρ(x) = i∈Z x i+1 x i e -h(x) 1 s(x) u (x) φ(x) dx = i∈Z 1 s(x i+1 -) u (x i+1 -)φ(x i+1 )e -h(x i+1 ) - 1 s(x i ) u (x i )φ(x i )e -h(x i ) - R a(x)u (x)φ (x) dx + 2 R b(x)u (x)φ(x) dx ρ(x) .
Using the interface condition (1.1.17

) at x i , R (α -A)u(x)φ(x) dx ρ(x) = E α (u, φ), ∀φ ∈ C 1 c (R, R). (1.1.18) If u ∈ L 2 (R) is such that Au ∈ L 2 (R), it follows from (1.1.18) that u = (α -L) -1 (α -A)u ∈ Dom(L). The space C c (R, R) is dense in L 2 (R). For f ∈ C c (R, R), u A = (α -A) -1 F ∈ Dom(A) ⊂ A a (R, R) and u L = (α -L) -1 f ∈ Dom(L). With (1.1.18), we deduce that u A is a version of u L . This proves that (α -A) -1 (C c (R, R)) ⊂ Dom(A) ⊂ A a (R, R) is dense in Dom(L) since C c (R, R) is.
We characterize the resolvent kernel.

Proposition 1.1.7. The resolvent kernel r is such that 1. For any λ ∈ H + , y ∈ R, x → r(λ, x, y) belongs to A a (R \ {y}, C) and vanishes at infinity.

The map

r(λ, •, •) is continuous from R 2 to C for any λ ∈ H + and r(•, x, y) is holomorphic on H + for any (x, y) ∈ R 2 .
3. For any y ∈ Dis(a),

D x r(λ, y-, y) -D x r(λ, y+, y) = 2 a(y)ρ(y) . (1.1.19)
Besides, if for some λ ∈ H + and y ∈ R, g ∈ A a (R \ {y}, C) vanishes at infinity and satisfies (λ -A)g(x) = 0 for x = y and (1.1.19) at y, then g(x) = r(λ, x, y) for any x ∈ R.

Proof. The direct part is an immediate consequence of Proposition 1.1.6 and (1.1.16).

Concerning the uniqueness, for any y ∈ Dis(a), the difference u(x) = g(x)r(λ, x, y) is solution to (λ -A)u(x) = 0 for any x ∈ Dis(a) ∪ {y} and vanishes at infinity. Because both g and r(λ, •, y) satisfy (1.1.19), u belongs to A a (R, C). It also satisfies (λ -A)u(x) = 0 for x ∈ Dis(a). Since Dis(a) is discrete and due to (1.1.10), (λ -A)u(x) = 0 for every x ∈ R. As R λ is a bounded operator for λ ∈ H + , u vanishes everywhere because it vanishes at infinity. Using Proposition 1.1.5, we weaken the conditions of applications of Theorem 3.1 in [START_REF] Martinez | One-dimensional parabolic diffraction equations : pointwise estimates and discretization of related stochastic differential equations with weighted local times[END_REF].

Corollary 1.1.3. If f ∈ L ∞ (R), let us set u(t, x) = R p(t, x, y)f (y) dy = E x [f (X a t )].
For each t > 0, u(t, x) ∈ A a (R). Hence, it is a classical solution to ∂ t u(t, x) = Lu(t, x) for t > 0, x ∈ R \ Dis(a) satisfying the transmission condition (1.1.17) at any point of Dis(a).

If in addition, f ∈ L 2 (R), then u(t, x) is a variational solution in L 2 ([0, T ], H 1 (R))∩C([0, T ], L 2 (R)) to ∂ t u(t, x) = Lu(t, x).

Stability under space transforms

Let us consider a continuous, increasing one-to-one function φ and define for every

f : R → R the transform φ * f = f • φ -1 .
For a continuous increasing function g, it is easily checked that Clearly,

D φ * g (φ * f ) = φ * (D g f ) when D g f or D φ * g (φ * f ) exist
f ∈ Dom(A) if and only if φ * f ∈ Dom(A φ ) since A φ (φ * f ) = 1 2 D φ * M D φ * S (φ * f ) = 1 2 φ * (D M D S f ) = φ * (Af ). (1.1.20) Besides, (f, F ) satisfies (λ -A)f = F if and only if (φ * f, φ * F ) satisfies (λ -A φ )φ * f = φ * F .
By an identification of the domains, we easily found the infinitesimal generator of φ(X).

Proposition 1.1.8. Let X be the process with infinitesimal generator A. Then X φ = φ(X) has infinitesimal generator A φ .

As a result, a large class of transforms keep the family {X a } a∈C invariant.

Proposition 1.1.9. Assume that φ is continuous, locally absolutely continuous on R and that for some η > 0, η < φ (x) < η -1 for almost every

x ∈ R. Let us set φ g(x) = g•φ -1 •φ •φ -1 .
Then for a = (a, ρ, b) ∈ C, φ a ∈ C (the constants λ and Λ may be changed) and φ(X a ) = X φ a φ where φ a = (φ a, φ ρ, φ b).

Remark 1.1.4. This choice is such that g(x)∇(φ * f (x)) = φ * (g∇f )(x) for almost every x ∈ R.

Proof. From the conditions on φ, lim x→±∞ φ(x) = ±∞ so that φ is one-to-one. Moreover, we have seen in Proposition 1.1.8 that the scale function and the speed measure of X φ are φ * S and φ * M . From a change of variable, φ * S and φ * M are indeed identified with the scale function and speed measure of φ a. Thus, the result follows from Proposition 1.1.2.

Corollary 1.1.4. If φ is continuous and piecewise of class C 2 , where the points of discontinuity of φ have no cluster points (hence absolutely continuous), and 0 < η < φ (x) < η -1 for almost every x. If a ∈ P, then φ a ∈ P.

Remark 1.1.5. We could have also prove the Proposition 1.1.9 using the Itô-Tanaka formula as in [START_REF] Étoré | Approximation de processus de diffusion à coefficients discontinus en dimension un et appliations à la simulation[END_REF][START_REF] Étoré | On random walk simulation of one-dimensional diffusion processes with discontinuous coefficients[END_REF].

Remark 1.1.6. Using a well chosen transform is actually the key for reducing the problem to a simpler one. It was heavily used with different transforms in [START_REF] Lejay | A scheme for simulating one-dimensional diffusion processes with discontinuous coefficients[END_REF][START_REF] Lejay | Simulating diffusion processes in discontinuous media : A numerical scheme with constant time steps[END_REF] and in [START_REF] Martinez | Interprétations probabilistes d'opérateurs sous forme divergence et analyse de méthodes numériques probabilistes associées[END_REF][START_REF] Martinez | One-dimensional parabolic diffraction equations : pointwise estimates and discretization of related stochastic differential equations with weighted local times[END_REF].

The case we consider in the sequel is the following: if a ∈ P is piecewise continuous, and φ is piecewise linear, then a φ is also piecewise continuous.

Stochastic differential equations

Here we prove that a diffusion X a associated with a ∈ P is the unique strong solution to a SDE involving local time. For this purpose, we introduce

G = (σ, ν) σ is of bounded variation over R and λ ≤ σ(x) ≤ Λ, ∀x ∈ R ν is a finite measure over R with |ν({x})| < 1, ∀x ∈ R .
The measure ν is associated with a function s of bounded variation through

ν( dx) = s ( dx) s(x) + s(x-)
where s ( dx) is the measure associated with s. Now let us consider the SDE given by

X t = x + t 0 σ(X s ) dB s + R ν( dy)L y t (X) (1.1.21)
where B is a Brownian motion and L x (X) the symmetric local time of X at x. From [46, Theorem 2.3], this SDE has a unique strong solution. Now using [46, Proposition 2.2], we easily deduce that X is a diffusion process with scale function S = s and speed measure with density m( dx) = s(x) -1 σ -2 (x) dx.

Considering a = (a, ρ, b) ∈ P, the expressions (1.1.7)-(1.1.8) for the scale function and the speed measure leads us to associate to a the function σ(x) = a(x)ρ(x) as well as the measure

ν( dx) = b(x) a(x)ρ(x) dx + a (x) 2a(x) dx + x∈Dis(a) a(x+) -a(x-) a(x+) + a(x-) , (1.1.22)
where a (x) is the derivative of a on the intervals on which a is C 1 and a = 0 on the set of discontinuities of a of zero Lebesgue measure. Because ν is not necessarily a finite measure over R, σ and ν do not necessarily belong to G. However, using a localization argument, there is no problem to set (a, ρ, b) = (1, 1, 0) outside a given compact K so that ν becomes a finite measure over R and then to let K grow. Since a one-dimensional diffusion process is uniquely characterized by its scale function and speed measure, X a is a semi-martingale solution to (1.1.21). [START_REF] Feller | Diffusion processes in one dimension[END_REF], it is easily seen that X a and X aκ are solutions to the same SDE and equal in distribution.

For a ∈ C, the diffusion X a has been characterized as a diffusion with infinitesimal generator (L, Dom(L)) obtained through the quadratic form (E α , H 1 (R)) as well as a diffusion (A, Dom(A)) with A = 1 2 D M D S with (S, M ) given by (1.1.7)-(1.1.8).

Proposition 1.1.11. Let σ be piecewise C 1 and ν be a finite measure of type

ν( dx) = d(x) dx + i∈Z α x i δ x i , |α x i | ≤ 1 -, i ∈ J ⊂ Z,
for an increasing family {x i } i∈J with no cluster points, some > 0, and a bounded function d.

Then there exists a = (a, ρ, b) ∈ P such that X a is solution to (1.1.21).

Proof. Select an interval [x i , x i+1 ]. Set a i (x) = β i = 1, ρ i (x) = σ(x) and b(x) = d(x)a i (x)ρ i (x) for x ∈ (x i , x i+1
). Now, at the endpoints x i+1 , find β i+1 such that

β i+1 -β i β i+1 + β i = α i . Set a i+1 (x) = β i , ρ i (x) = σ(x)/β i and b(x) = d(x)a i (x)ρ i (x) for x ∈ (x i+1 , x i+2
). Hence, we construct (a, ρ, b) iteratively on the intervals (x j , x j+1 ) for j > i. The same construction could be performed for j < i by going downward instead of upward. It is easily checked that a = (a, ρ, b) ∈ P with Dis(a) = {x i } i∈J .

Resolvent in the presence of one discontinuity 1.2.1 A computation method

We now assume that the coefficients a = (a, ρ, b) ∈ P are discontinuous only at one point, say 0. We fix y = 0 and λ ∈ C.

We denote by I the intervals (-∞, 0 ∧ y), (0 ∧ y, y ∨ 0) and (y ∨ 0, +∞) and introduce the operator A I whose coefficients a I are such that a I = a on I.

Let u I and u I be the two functions on C × R given by Proposition 1.1.4 for each A I . When y > 0, we set

q 1 (λ, x) = u (-∞,0) (λ, x), q 2 (λ, x) = u (0,y) (λ, x), q 3 (λ, x) = u (0,y) (λ, x) and q 4 (λ, x) = u (y,∞) (λ, x),
whereas, for y < 0,

q 1 (λ, x) = u (-∞,y) (λ, x), q 2 (λ, x) = u (y,0) (λ, x), q 3 (λ, x) = u (y,0) (λ, x) and q 4 (λ, x) = u (0,∞) (λ, x).
One-dimensional skew diffusions: explicit expressions of densities and resolvent Proposition 1.2.1. When a ∈ P is only discontinuous at 0, for each y ∈ R \ {0},

r(λ, x, y) = k 1 (λ, y) q 1 (λ, x)1(x < 0) + k 2 (λ, y) q 2 (λ, x)1(x ∈ (0, y)) + k 3 (λ, y) q 3 (λ, x)1(x ∈ (0, y)) + k 4 (λ, y) q 4 (λ, x)1(x ≥ y), (1.2.1)
for any λ ∈ H + and any x ∈ R, with the function k

= (k 1 , k 2 , k 3 , k 4 ) solution to M (λ, y)k(λ, y) T = 0 0 0 1 T (1.2.2)
where for y > 0,

M (λ, y) =     q 1 (λ, 0) -q 2 (λ, 0) -q 3 (λ, 0) 0 a(0-)q 1 (λ, 0) -a(0+)q 2 (λ, 0) -a(0+)q 3 (λ, 0) 0 0 q 2 (λ, y) q 3 (λ, y) -q 4 (λ, y) 0 ρ(y)a(y) 2 q 2 (λ, y) ρ(y)a(y) 2 q 3 (λ, y) -ρ(y)a(y) 2 q 4 (λ, y)     while for y < 0, M (λ, y) =     0 q 2 (λ, 0) q 3 (λ, 0)
-q 4 (λ, 0) 0 a(0-)q 2 (λ, 0) a(0-)q 3 (λ, 0) a(0+)q 4 (λ, 0) q 1 (λ, y) -q 2 (λ, y) -q 3 (λ, y) 0 ρ(y)a(y) 2 q 1 (λ, y) -ρ(y)a(y) 2 q 2 (λ, y) -ρ(y)a(y)

2 q 3 (λ, y) 0     .
Proof. When λ ∈ C, any solution to (λ -A)u = 0 on some interval I should be sought as a linear combination of u (λ, •) and u (λ, •) provided that Wr[u , u ](λ) = 0. This is the case for λ ∈ H + (See Corollary 1.1.2). Our choice of q i (λ, •) in the decomposition 1.2.1 ensures the correct vanishing behaviour at infinity for λ ∈ H + . The system (1.2.2) is an algebraic transcription of the compatibility condition at 0 given by (1.1.17) and y given by (1.1.19) in Proposition 1.1.7. Hence, with any solution k(λ, y) to (1.2.2), the right-hand-side of (1.2.1) is indeed equal to r(λ, x, y) for any λ ∈ H + .

When M (λ, y) is invertible for y ∈ R and λ in a domain H + then λ → M (λ, y) -1 is holomorphic on H + since it is composed of sums, products and ratios of the terms of M (•, y) which are themselves holomorphic on H + .

Let us prove now that M (λ, y) is invertible for λ ∈ H + . Since q 2 (λ, y) = q 4 (λ, y) and q 2 (λ, y) = q 4 (λ, y) for any λ ∈ C, the determinant of

M (λ, y) is for y > 0 det M (λ, y) = ρ(y) exp(-h(y)) 2 Wr[q 2 , q 3 ](λ, y)
× {a(0-)q 1 (λ, 0)q 2 (λ, 0)a(0+)q 1 (λ, 0)q 2 (λ, 0)}

where Wr[q 2 , q 3 ](λ, •) denotes the Wronskian of q 2 (λ, •) and q 3 (λ, •) in x. The above equality follows from the relation a(x)D x = exp(-h(x))D S . From Corollary 1.1.2, Wr[q 2 , q 3 ](λ, y) = 0 for any λ ∈ H + . Let us consider the solutions u and u associated with the coefficients a given by Proposition 1.1.4. From our choice of q 1 and q 2 , q 1 (λ, 0) = q 2 (λ, 0) = u (λ, 0) = u (λ, 0) = 1 and the functions q 1 and u vanish at -∞ whereas q 2 and u vanish at +∞. With Corollary 1.1.1, then u (λ, x) = q 1 (λ, x) for x ≤ 0 and u (λ, x) = q 2 (λ, x) for x ≥ 0. Thus,

a(0-)q 1 (λ, 0)q 2 (λ, 0) -a(0+)q 1 (λ, 0)q 2 (λ, 0) = exp(-h(0)) Wr[u , u ](λ, 0).
This quantity cannot vanish on H + according to Corollary 1.1.2. Thus, det M (λ, y) = 0 for y > 0, λ ∈ H + and M (λ, y) is holomorphic on H + . The case y < 0 is similar. Remark 1.2.1. Of course, this approach is also feasible when the state space is bounded, one should consider up to 6 functions with the appropriate boundary conditions. This is useful for example to provide formulae related to the first exit time or first passage time as in [START_REF] Lefebvre | First passage problems for asymmetric Wiener processes[END_REF][START_REF] Abundo | First-passage problems for asymmetric diffusions and skew-diffusion processes[END_REF]. One could also use an expansion relying on eigenvalues and eigenfunctions. For practical purposes, especially for numerical simulation which was our original goal, a large number of eigenfunctions should be considered (See [START_REF] Guerrero | Analytical solution for the advection-dispersion transport equation in layered media[END_REF]). Actually, the eigenvalue formulation is the best for estimating the density for large times.

Skew Brownian motion with a piecewise constant drift

We now apply the computation method developed in the previous section to skew diffusions with a piecewise constant drift. Therefore, let us consider the operator

L = ρ(x) 2 d dx a(x) d dx + b(x) d dx (1.2.3) with for some β ∈ (0, 1), a = (a, ρ, b) = a(β, b + , b-) defined in Table 1.1. for x ≥ 0 for x < 0 a(x) β 1 -β ρ(x) β -1 (1 -β) -1 b(x) b + b - Table 1.1 -Coefficients a(β, b + , b -) ∈ P.
Remark 1.2.2. Out of 0, ρ(x)a(x) = 1. Thus, from (1.1.21) and (1.1.22), the associated diffusion process X a is solution to

X a t = x + B t + (2β -1)L 0 t (X a ) + t 0 b(X a s ) ds for t ≥ 0. (1.2.4)
For some constant γ ∈ R, we consider the two functions u and u given by (1.1.15) in Example 1.1.2 and remark that the Green function of L γ = 1 2 D 2 x + γD x has a density with respect to the Lebesgue measure containing these functions and which is defined by g(γ, λ, x, y) = 1

γ 2 + 2λ e (-γ+ √ γ 2 +2λ)(x-y) if x < y, e (-γ- √ γ 2 +2λ)(x-y) if x ≥ y. (1.2.5) Proposition 1.2.2. The resolvent kernel of L with coefficients a(β, b + , b -) is for y ≥ 0 and λ ∈ C \ (-∞, 0), r(λ, x, y) = g(b + , λ, x, y)1(x ≥ 0) + A -+ (λ, y)g(b -, λ, x, y)1(x < 0) + A ++ (λ, y)g(b + , λ, x, -y)1(x ≥ 0) (1.2.6)
and for y < 0,

r(λ, x, y) = g(b -, λ, x, y)1(x < 0) + A --(λ, y)g(b -, λ, x, -y)1(x < 0) + A +-(λ, y)g(b + , λ, x, y)1(x ≥ 0) (1.2.7) with A -+ (λ, y) = Θ(λ, b + , b -) -1 2 β b 2 -+ 2λ e (b + -b -+ √ b 2 -+2λ- √ b 2 + +2λ)y , A ++ (λ, y) = Θ(λ, b + , b -) -1 β( b 2 + + 2λ -b + ) -(1 -β)( b 2 -+ 2λ -b -) e 2b + y , A --(λ, y) = Θ(λ, b + , b -) -1 -β ( b 2 + + 2λ + b + ) + (1 -β) ( b 2 -+ 2λ + b -) e 2b -y , A +-(λ, y) = Θ(λ, b + , b -) -1 2( β -1) b 2 + + 2λe (b --b + + √ b 2 -+2λ- √ b + + 2λ)y ,
where the common denominator is

Θ(λ, b + , b -) = β b 2 + + 2λ + b + + (1 -β) b 2 -+ 2λ -b -. (1.2.8)
Remark 1.2.3. Using the expressions of g(b ± , λ, x, y) and A ±,± (λ, y), the inverse transform L -1 can be applied to each terms of the decompositions (1.2.6) and (1. Proof. It is an application of Proposition 1.2.1. Since u (γ, λ, 0) = u (γ, λ, 0) = 1, we set

q 1 (λ, x) = u (b -, λ, x), q 3 (λ, x) = u (b(y), λ, x), q 2 (λ, x) = u (b(y), λ, x) and q 4 (λ, x) = u (b + , λ, x).
The solution of M (λ, y)k(λ, y) = 0 0 0 1

T is for y ≥ 0,

                         k 1 (λ, y) = -2 β (β-1)(-b -+ √ b 2 -+2λ)+β (-b + - √ b 2 + +2λ) e -(-b + + √ b 2 + +2λ)y , k 2 (λ, y) = (1-β)(-b -+ √ b 2 -+2λ)-β(-b + + √ b 2 + +2λ) √ b 2 + +2λ[(β-1)(-b -+ √ b 2 -+2λ)+β (-b + - √ b 2 + +2λ)] e -(-b + + √ b 2 + +2λ)y k 3 (λ, y) = 1 √ b 2 + +2λ e -(-b + + √ b 2 + +2λ)y , k 4 (λ, y) = (1-β)(-b -+ √ b 2 -+2λ)-β(-b + + √ b 2 + +2λ) √ b 2 + +2λ[(β-1)(-b -+ √ b 2 -+2λ)+β (-b + - √ b 2 + +2λ)] e -(-b + + √ b 2 + +2λ)y , + 1 √ b 2 + +2λ e -(-b + - √ b 2 + +2λ)y , while for y < 0,                          k 1 (λ, y) = (1-β) (-b -- √ b 2 -+2λ)-β (-b + - √ b 2 + +2λ) √ b 2 -+2λ[(β-1)(-b -+ √ b 2 -+2λ)+β (-b + - √ b 2 + +2λ)] e -(-b -- √ b 2 -+2λ)y , + 1 √ b 2 -+2λ e -(-b -+ √ b 2 -+2λ)y , k 2 (λ, y) = 1 √ b 2 -+2λ e -(-b -- √ b 2 -+2λ)y , k 3 (λ, y) = (1-β) (-b -- √ b 2 -+2λ)-β (-b + - √ b 2 + +2λ) √ b 2 -+2λ[(β-1)(-b -+ √ b 2 -+2λ)+β (-b + - √ b 2 + +2λ)] e -(-b -- √ b 2 -+2λ)y k 4 (λ, y) = 2( β-1) (β-1)(-b -+ √ b 2 -+2λ)+β (-b + - √ b 2 + +2λ) e -(-b -- √ b 2 -+2λ)y .
Finally, it remains to plug the coefficients at their rightful places to obtain the solutions in each case. After some simplifications, we get (1.2.6) and (1.2.7).

Skew diffusion with piecewise constant coefficients

Let us now consider the diffusion whose infinitesimal generator is (L, Dom(L)) given by (1.2.3) and coefficients a = (a, ρ, b) are given in Table 1.2.

In order to set up the diffusion coefficient √ aρ to 1, the one-to-one transform φ(x) =

x/ a(x)ρ(x) is convenient. If X a is the diffusion process associated with a = (a, ρ, b), then X φ a is the diffusion process associated to

φ a = √ a √ ρ , √ ρ √ a , b √ aρ . for x ≥ 0 for x < 0 a(x) a + a - ρ(x) ρ + ρ - b(x) b + b - Table 1.2 -Coefficients a = (a, ρ, b) constants on R + and R -.
With Remark 1.1.7, we can consider κ = ( a + /ρ + + a -/ρ -) -1 and see that X φ a has the same distribution as the SBM with parameter β = κ a + /ρ + ∈ (0, 1) and piecewise

constant drift b(x) = b(x)/ a(x)ρ(x), that is φ a = a(β, b + / √ a + ρ + , b -/ √ a -ρ -).
Hence, the computations of Section 1.2.1 could be applied to get the resolvent of X φ a . Now using a change of variable,

r a (λ, x, y) = r φ a (λ, φ(x), φ(y)) a(y)ρ(y) ,
where r a (resp. r φ a ) is the resolvent kernel of X a (resp. X φ a ).

Some explicit expressions for the density

We introduce the functions

H(y) = 1(y ≥ 0), H -(y) = 1(y < 0) = 1(y ∈ R) -H(y), sgn(y) = 1(y ≥ 0) -1(y < 0), erf(x) = 2 √ π x 0 e -v 2 dv and erfc(x) = 1 -erf(x) = 2 √ π +∞ x e -v 2 dv.

The Skew Brownian Motion

The SBM of parameter β is given by the choice of the coefficients a(β, 0, 0). Thus,

A -+ (λ, y) = 2β, A ++ (λ, y) = 2β -1, A --(λ, y) = 1 -2β, A --(λ, y) = 1 -2β.
To recover the density of the SBM, we need to invert the Green function of the Laplace operator since none of the A ±± (λ, y) depend on λ. It is well known that it gives the Gaussian kernel.

More precisely, from [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables[END_REF], and

L -1 1 √ λ e -k √ λ = 1 √ πt e -k 2 4t , k ≥ 0, (1.3
L -1 (f (cλ + d)) = 1 c e -d c t L -1 (f ) t c , c > 0. (1.3.2)
Thus, after an application of (1.3.1)-(1.3.2) to (1.2.6) and (1.2.7) as well as some rewriting

p(t, x, y) = 1 √ 2πt e -(y-x) 2 2t + sgn(y)(2β -1) 1 √ 2πt e -(|y|+|x|) 2 2t
.

This density was obtained in [START_REF] Walsh | A diffusion with discontinuous local time[END_REF] through a probabilistic argument. Some plots are given in Figure 1.1.

The cumulative distribution function P : R + × R × R → (0, 1) of p is

P (t, x, y) = 1 2 1 + erf y -x √ 2t - 2β -1 2 1 -sgn(y) erf y + sgn(y)|x| √ 2t . and A +-(λ, y) = 2(β -1) √ b 2 + 2λ √ b 2 + 2λ + b(2β -1) .
Let us compute the Laplace inverse of A ±± (λ, y)g(b, λ, x, y). From [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables[END_REF], and

L -1 e -k √ λ d + √ λ = 1 √ πt e -k 2 4t -de dk e d 2 t erfc d √ t + k 2 √ t , k ≥ 0, (1.3.3)
L -1 e -k √ λ √ λ(d + √ λ) = e dk e d 2 t erfc d √ t + k 2 √ t , k ≥ 0. ( 1
p(t, x, y) = 1 √ 2πt e -(bt-(y-x)) 2 2t + sgn(y)(2β -1) 1 √ 2πt e -(|y|+|x|) 2 2t e -b 2 t 2 e b(y-x) + sgn(y)(H -(y) -β)b(2β -1)e b(2β-1)(|y|+|x|) e b(y-x) e 2b 2 β(β-1)t × erfc b(2β -1) t 2 + |y| + |x| √ 2t .
We recover the density obtained in [START_REF] Appuhamillage | Occupation and local times for skew Brownian motion with applications to dispersion across an interface[END_REF][START_REF] Étoré | Exact simulation of one-dimensional stochastic differential equations involving the local time at zero of the unknown process[END_REF] up to a conversion of erfc to its probabilistic counterpart.

The Bang-Bang Skew Brownian Motion

The Bang-Bang SBM is the diffusion corresponding to the choice of the coefficients a(β, b, -b) for some b = 0. This process was introduced in [

p(t, x, y) = 1 √ 2πt e -(bt-sgn(y)(y-x)) 2 2t + sgn(y)(2β -1) 1 √ 2πt e -(|y|+|x|) 2 2t e -b 2 t 2 e b(|y|-|x|) + sgn(y)(H -(y) -β)be 2b|y| erfc b t 2 + |y| + |x| √ 2t .
This is the expression given in [START_REF] Fernholz | Two brownian particles with rank-based characteristics and skew-elastic collisions[END_REF] up to a conversion of erfc into its probabilistic counterpart.

The constant Péclet case

Computing the transition density in all situations seems to be a very difficult problem. There are however other cases were simplifications occur. Among these situations, we consider the case where b + and b -are linked by

b -= µb + with µ = β 1 -β or equivalently (1 -β)b -= βb + , β = 1 2 . (1.3.5) Remark 1.3.1.
This assumption on b + and b -is natural. If we consider the diffusion X a with piecewise constant coefficients a as in Section 1.2.3, for φ(x) = x/ a(x)ρ(x), then φ(X) is a SBM with a piecewise constant drift. Is it easily checked that (1.3.5) is satisfied when

b + ρ + = b - ρ - .
When a = 1, the ratio b/ρ is called the Péclet number. It is a dimensionless quantity which plays a very important role in fluid mechanics by characterizing the effect of the convection against the diffusion and vice versa.

When (1.3.5) holds,

Θ(λ, b + , µb -) = β b 2 + + 2λ + (1 -β) b 2 -+ 2λ.
Hence,

1 Θ(λ, b + , µb -) = β b 2 + + 2λ -(1 -β) b 2 -+ 2λ β 2 (b 2 + + 2λ) -(1 -β) 2 (b 2 -+ 2λ) = β b 2 + + 2λ -(1 -β) b 2 -+ 2λ 2λ(2β -1)
.

Therefore, for x = 0,

r(λ, y) = β e b + y 2β -1 β b 2 + + 2λ -(1 -β) b 2 -+ 2λ λ e - √ b 2 + +2λy 1(y ≥ 0) - (1 -β) e b -y 2β -1 β b 2 + + 2λ -(1 -β) b 2 -+ 2λ λ e √ b 2
-+2λy 1(y < 0).

Lemma 1.3.1. For a ∈ R, L -1 √ a + se -√ s y (t, y) = 1 π a 0 sin(y √ r) √ a -r e -rt dr - 1 π +∞ 0 cos(y √ r + a) √ r e -(r+a)t dr.
Proof. Inspired by [START_REF] Puri | Some inverse laplace transforms of exponential form[END_REF], we the use the Bromwich formula with the contour Γ illustrated in Figure 1.4.

Since the integrals on the outer and inner arcs as well as half-circles are null,

L -1 √ a + s e -√ s y (r, y) = 1 2iπ
p(t, 0, y) = 2 (H -(y) -β) 2 (2β -1) √ 2πt e -(b(x) 2 t-y) 2 2t + b(x) (H -(y) -β) 2 2(2β -1) erfc y √ 2t -b(x) t 2 -e 2b(x) y erfc y √ 2t + b(x) t 2 - β (1 -β) e b(x) y 2(2β -1)π t 0 e -b(x) 2 τ 2 b 2 (-x)-b 2 (x) 0 sin(y √ r) b 2 (-x) -b 2 (x) -re -r τ 2 dr dτ β (1 -β) e b(x) y 2(2β -1)π t 0 e -(b(-x) 2 ) τ 2 +∞ 0 cos(y r + b 2 (-x) -b 2 (x)) √ re -r τ 2
dr dτ.

Proof. First,

L -1 (r)(t, y) = β 2 e b + y 2β -1 L -1 b 2 + + 2λ λ e - √ b + +2λy (t, y)1(y ≥ 0) - β (1 -β) e b + y 2β -1 L -1 b 2 -+ 2λ λ e - √ b + +2λy (t, y)1(y ≥ 0) - β (1 -β) e b -y 2β -1 L -1 b 2 + + 2λ λ e √ b -+2λy (t, y)1(y < 0) + (1 -β) 2 e b -y 2β -1 L -1 b 2 -+ 2λ λ e √ b -+2λy
(t, y)1(y < 0).

By remarking that

L -1 b 2 + + 2λ λ e - √ b + +2λy (t, y) = L -1 b 2 + λ b 2 + + 2λ e - √ b + +2λy
(t, y)

+ 2 L -1 1 b 2 + + 2λ e - √ b + +2λy
(t, y)1(y ≥ 0), and using formulas (1.3.1), (1.3.2) and

L -1 e -ak √ s (s -d) √ s = e dt 2d √ de -ak √ d erfc( ak 2 √ t - √ dt) - e dt 2d √ de ak √ d erfc( ak 2 √ t + √ dt)
which can be found in [START_REF] Lindstrom | A note on a laplace transform pair associated with mass transport in porous media and heat transport problems[END_REF],

L -1 b 2 + + 2λ λ e - √ b + +2λy (t, y) = 1 2 b + e -b + y erfc y √ 2t -b + t 2 - 1 2 b + e b + y erfc y √ 2t + b + t 2 + √ 2 √ πt e -b 2 + t 2 e -y 2 2t .
Then, using Lemma 1.3.1,

L -1 b 2 -+ 2λ λ e - √ b + +2λy (t, y) = 1 2π t 0 e -b 2 + τ 2 b 2 --b 2 + 0 sin(y √ r) b 2 --b 2 + -r e -r τ 2 dr dτ - 1 2π t 0 e -b 2 + τ 2 +∞ 0 cos(y r + b 2 --b 2 + ) √ re -(r+b 2 --b 2 + )τ 2 dr dτ.
Finally, when y ≥ 0,

L -1 (r)(t, y) = b + β 2 2(2β -1) erfc y √ 2t -b + t 2 1(y ≥ 0) - b + β 2 e 2b + y 2(2β -1) erfc y √ 2t + b + t 2 1(y ≥ 0) + 2 β 2 (2β -1) √ 2πt e -(b + t-y) 2 2t - β (1 -β) e b + y 2(2β -1)π t 0 e -b 2 + τ 2 b 2 --b 2 + 0 sin(y √ r) b 2 --b 2 + -re -r τ 2 dr dτ 1(y ≥ 0) + β (1 -β) e b + y 2(2β -1)π t 0 e -b 2 -τ 2 +∞ 0 cos(y r + b 2 --b 2 + ) √ re -r τ 2
dr dτ 1(y ≥ 0).

The result proved using a symmetry argument when y < 0.

Conclusion

We have proposed a new method to compute the resolvent kernel of skew diffusions. In some cases, it leads to explicit expressions of their density.

Our approach is different from the one of B. Gaveau et al. [START_REF] Gaveau | Second order differential operators and dirichlet integrals with singular coefficients, i. functional calculus of one-dimensional operators[END_REF]. Drift term can be taken in consideration, even when the drift is discontinuous. Not only known expressions are easily recovered but also new ones are computed. Hence, our approach overcomes the difficulties of a purely stochastic one.

At first glance, this approach seems to be rooted in the Feller theory of one-dimensional diffusions. Some generalizations are possible in particular cases of multi-dimensional diffusions. In addition, it leads to the development of new simulation techniques for the diffusions relying on exponential time steps instead of constant time step. This way, the expressions of the resolvent kernel are used as they are simpler than the ones of the density. These two directions of research are currently under investigation.

Introduction

We consider the skew diffusions in dimension one, that is the continuous stochastic processes whose transition functions are the solutions to

∂ ∂t p(t, x, y) = ρ(x)
2 div (a(x)∇p(t, x, y))

+ b(x)∇p(t, x, y), ∀x ∈ R, 0 < t ≤ T, p(0, x, y) = δ y (x), y ∈ R (2.0.1)
when (a, ρ, b) are piecewise constant. We decide to call them constant skew diffusions.

Recently [START_REF] Dereudre | An explicit representation of the trnaistion densities of the skew brownian motion with drift and two semipermeable barriers[END_REF][START_REF] Étoré | On the existence of a time inhomogeneous skew Brownian motion and some related laws[END_REF][START_REF] Étoré | Exact simulation of one-dimensional stochastic differential equations involving the local time at zero of the unknown process[END_REF] new algorithms were developed to sample the paths of these processes. But, none of them actually give an answer for any combination of piecewise constant coefficients (a, ρ, b). Indeed, [START_REF] Étoré | On the existence of a time inhomogeneous skew Brownian motion and some related laws[END_REF][START_REF] Étoré | Exact simulation of one-dimensional stochastic differential equations involving the local time at zero of the unknown process[END_REF] provide very interesting exact schemes but only suited for a single discontinuity while [START_REF] Dereudre | An explicit representation of the trnaistion densities of the skew brownian motion with drift and two semipermeable barriers[END_REF] is designed to handle multiple discontinuities but does not work for a discontinuous b.

In this note, we propose to give two new simulation schemes for one-dimensional constant skew diffusions with one discontinuity at 0 which can be extended in case of multiple discontinuities as they are of the Euler type [START_REF] Lejay | Simulating diffusion processes in discontinuous media : A numerical scheme with constant time steps[END_REF]. Besides, both of our schemes can be mixed with a classical scheme as in [START_REF] Lejay | Simulating diffusion processes in discontinuous media : A numerical scheme with constant time steps[END_REF] to simplify the simulation of the skew diffusion when it is far from the discontinuities.

The outline of this note is as follows. In §2.1, using the idea of [START_REF] Jansons | Efficient numerical solution of stochastic differential equations using exponential timestepping[END_REF][START_REF] Jansons | Exponential timestepping with boundary test for stochastic differential equations[END_REF], we build a simulation scheme for Feller processes that uses their resolvent density . Then, §2.2 we prove a Donsker invariance principle to argue for the consistency of this scheme. In §2.3 we quickly show that a constant skew diffusion is a Feller process and possesses a resolvent density. We also provide a change of variable allowing to simplify the coefficients (a, ρ, b) and to reduce the problem. After that, we introduce always in §2.3 a closed-form of the resolvent density for the constant skew diffusions with simplified coefficients that were computed in [START_REF] Lejay | One-dimensional skew diffusions : explicit expressions of densities and resolvent kernel[END_REF]. Finally in §2.4 and §2.5 we detail respectively the first and second schemes for the reduced problem.

A numerical scheme for Feller processes using their resolvent density

Let X = (X t ) 0≤t≤T be a one-dimensional Feller process starting at x, p(t, x, y) be its transition function, r(λ, x, y) be its resolvent density and τ be an exponentially distributed random variable with parameter λ which is independent of X. Proposition 2.1.1. The random variable X τ possesses the density λ r(λ, x, y).

for each n ≥ 0 the n-steps random walk (Y n t ) 0≤t≤1 and (Z n t ) 0≤t≤1 from the linear interpolation of respectively Y k/n = X n τ k for each k = 1, . . . , n and (X τ 0 , . . . , X τn ). From their definition, (Y n t ) 0≤t≤1 and (X n t ) 0≤t≤1 have the same law for each n ≥ 0. Thus, if (Y n t ) 0≤t≤1 converges in distribution to X, then (X n t ) 0≤t≤T converges also to X. By [8, Theorem 3.1], (Y n t ) 0≤t≤1 converges in distribution to X if we show (i) d(Y n , Z n ) converges almost surely 0 which reduced here from the interpolation properties to prove sup 0≤k≤n |τ kk/n| converges a.s to 0.

(ii) (Z n ) n≥0 converges in distribution to X.

By the construction above, (τ kk/n) 0≤k≤n is a martingale. Therefore, using the Doob's martingale inequality and the fact that the sums of exponential random variables are Gamma distributed,

P sup 0≤k≤n (τ k -k/n) ≥ K ≤ 1 K p n p 2 .
(2.2.1)

and we can choose p such that, for any K > 0,

+∞ n=1 P sup 0≤k≤n (τ k -k/n) ≥ K < +∞. (2.2.2)
Now from the properties of Gamma variables, one derives that P sup 0≤k≤n (τ kk/n) ≥ K converges in distribution to 0. Thus, using (2.2.2) and the well-known Borel-Cantelli lemma, sup 0≤k≤n |τ kk/n| converges a.s to 0 and (i) is proved.

Since Z n is the linear interpolation of (X 0 , X τ 1 , • • • , X τn ),

|Z n t -X t | ≤ |Z n t -X τ nt | + |X t -X τ nt | ≤ |X t -X τ nt +1 | + |X t -X τ nt | + |X t -X τ nt |.
Hence,

P sup 0≤t≤1 |Z n t -X t | ≥ K ≤ P sup 0≤t≤1 |X t -X τ nt +1 | ≥ K + 2 P sup 0≤t≤1 |X t -X τ nt | ≥ K , (2.2.3)
and (Z n ) n≥0 converges in probability to X if we show that both terms in the right hand side of (2.2.3) go to 0. The proof is similar for both terms of (2.2.3). Therefore, we only deal with the second one. As the paths of X are continuous [START_REF] Feller | Generalized second order differential operators and their lateral conditions[END_REF][START_REF] Feller | On the intrinsic form for second order differential operator[END_REF], for any K > 0 and any ε > 0, it exists δ > 0 such that

P sup |t-s|<δ |X t -X t | ≥ K ≤ ε.
(2.2.4)

Two numerical schemes for the simulation of skew diffusions using their resolvent Subsequently, since

P sup 0≤t≤1 |X t -X τ nt | ≥ K ≤ P sup |t-s|<δ |X t -X t | ≥ K + P sup 0≤k≤n (τ k -k/n) ≥ δ , (2.2.5)
(2.2.1) and (2.2.4) immediately provide that (2.2.5) goes to 0. As a result, (ii) is proved.

Existence of the resolvent density, problem reduction and formulae

We have the following fundamental result on skew diffusions.

Proposition 2.3.1. Let the coefficients in (2.0.1) be such that

(a, ρ, b) are measurable from R 3 → R with λ ≤ a, ρ ≤ Λ and |b| ≤ Λ.
Then, (2.0.1) possesses a unique solution

p(t, x, y) ∈ L 2 ([0, T ], H 1 (R)) ∩ C([0, T ], L 2 (R)) (2.3.1)
which is the transition function of a conservative Feller process (X t , F t , P x ) t≥0,x∈R with continuous paths. Besides, X possesses the infinitesimal generator

A = ρ(x) 2 div(a(x)∇•) + b(x)∇ • with Dom(A) = f ∈ C b (R, R) | Af ∈ C b (R, R) ,
where C b (R, R) = f : R → R | f is continuous and bounded and the resolvent of (A, Dom(A)) defined for any λ ∈ C such that Re[λ] > 0 admits a density with respect to Lebesgue measure

r(λ, x, y) = +∞ 0 p(t, x, y)e -λt dt, ∀(x, y) ∈ R 2 ,
that is the Laplace transform of p(t, x, y).

Proof. First, the existence of the solution p(t, x, y) and its regularity is proved in [START_REF] Aronson | Non-negative solutions of linear parabolic equations[END_REF]. Secondly, the existence of the continuous stochastic process X is deducible through an adaptation of [START_REF] Feller | Generalized second order differential operators and their lateral conditions[END_REF][START_REF] Feller | On the intrinsic form for second order differential operator[END_REF] as well as the infinitesimal generator with [START_REF] Lejay | One-dimensional skew diffusions : explicit expressions of densities and resolvent kernel[END_REF]. At last, the result on the density of the resolvent associated with the infinitesimal generator can be found in [START_REF] Hille | Functional analysis and semi-groups[END_REF][START_REF] Pazy | Semigroups of linear operators and applications to partial differential equations[END_REF].

Remark 2.3.1. The above result is classical when the coefficients belongs to C ∞ (R) and the Feller process obtained is a semi-martingale. This is in general not true when the coefficients are just measurable. With this result, Algorithm 2 can be used to simulate skew diffusions assuming a closedform of the resolvent density is known. Now let X be a constant skew diffusion having only one discontinuity at 0, that is with

a(x) = a + 1(x ≥ 0) + a -1(x < 0), ρ(x) = ρ + 1(x ≥ 0) + ρ -1(x < 0) and b(x) = b + 1(x ≥ 0) + b -1(x < 0), (2.3.2) with (a + , a -, ρ + , ρ -) ∈ R * + × R * + × R * + × R * + and (b + , b -) ∈ R 2 .
In [START_REF] Étoré | Approximation de processus de diffusion à coefficients discontinus en dimension un et appliations à la simulation[END_REF][START_REF] Lejay | One-dimensional skew diffusions : explicit expressions of densities and resolvent kernel[END_REF], it is shown that X = φ -1 (X ) where φ(x) = x/ a(x) ρ(x) and X is the constant skew diffusion whose coefficients are

â(x) = β 1(x ≥ 0) + (1 -β) 1(x < 0), ρ(x) = β -1 1(x ≥ 0) + (1 -β) -1 1(x < 0) with β = a + /ρ + a + /ρ + + a -/ρ - (2.3.3) and b(x) = b + √ a + ρ + 1(x ≥ 0) + b - √ a -ρ - 1(x < 0) = b+ 1(x ≥ 0) + b-1(x < 0).
The constant skew diffusions with coefficients of type (2.3.3) are well-known. They are the unique strong solution to the stochastic differential equation involving local time

Y t = Y 0 + B t + t 0 b(Y s ) ds + (2β -1)L 0 t (Y ), X 0 = x. (2.3.4)
where B = (B t ) t≥0 is a Brownian motion and Y = (Y t ) t≥0 is the unknown process [START_REF] Gall | One-dimensional stochastic differential equations involving the local times of the unknown process[END_REF]. Besides, a closed-form of their resolvent density was recently computed in [START_REF] Lejay | One-dimensional skew diffusions : explicit expressions of densities and resolvent kernel[END_REF].

Proposition 2.3.2 (A. Lejay et al., [START_REF] Lejay | One-dimensional skew diffusions : explicit expressions of densities and resolvent kernel[END_REF]). Let g(γ, λ, x, y) = 1

γ 2 + 2λ e (-γ+ √ γ 2 +2λ)(x-y) if x < y, e (-γ- √ γ 2 +2λ)(x-y) if x ≥ y.
(2.3.5)

The resolvent density of a skew diffusion X with coefficients satisfying

(2.3.3) is r(λ, x, y) = g( b+ , λ, x, y)1(x ≥ 0) + A -+ (λ, y)g( b-, λ, x, y)1(x < 0) + A ++ (λ, y)g( b+ , λ, x, -y)1(x ≥ 0). (2.3.6)
and for y < 0,

r(λ, x, y) = g( b-, λ, x, y)1(x < 0) + A --(λ, y)g( b-, λ, x, -y)1(x < 0) + A +-(λ, y)g( b+ , λ, x, y)1(x ≥ 0) (2.3.7) with A -+ (λ, y) = Θ(λ, b+ , b-) -1 2 β b2 -+ 2λ e ( b+ -b-+ √ b2 -+2λ- √ b2 + +2λ)y , A ++ (λ, y) = Θ(λ, b+ , b-) -1 β( b2 + + 2λ -b+ ) -(1 -β)( b2 -+ 2λ -b-) e 2 b+ y , A --(λ, y) = Θ(λ, b+ , b-) -1 -β ( b2 + + 2λ + b+ ) + (1 -β) ( b2 -+ 2λ + b-) e 2 b-y , A +-(λ, y) = Θ(λ, b+ , b-) -1 2( β -1) b2 + + 2λe ( b--b+ + √ b2 -+2λ- √ b+ +2λ)y ,
where the common denominator is

Θ(λ, b+ , b-) = β b2 + + 2λ + b+ + (1 -β) b2 -+ 2λ -b-. (2.3.8)
With this last result and φ, we are finally able to simulate any constant skew diffusions with coefficients of type (2.3.2) using Algorithm 2. Furthermore, any scheme that sample X with (â, ρ, b) given by (2.3.3) simulate Y with (a, ρ, b) given by (2.3.2) after the application of φ -1 .

The first scheme

We treat the case where X has coefficients given by (2.3.3) through a series of algorithms which once nested one within the others give the scheme when the starting point of X is non negative. By the symmetry of (2.3.5), these algorithms can be easily adapted when the starting point is negative, thus providing a complete scheme.

We start with a few notations. E(α) denotes an exponential distribution with parameters α and B(p) a Bernoulli distribution with parameters p. The notations Q(α, x) and T (α, x) stand for the distribution whose densities are respectively α e -α(y-x) 1(y ≥ x ≥ 0) and α 1e αx e αy 1(0 ≤ y ≤ x)

From the expression in the Proposition 2.3.2,

P[X τ X 0 < 0|X 0 ≥ 0] = 2λ (β -1) (β -1)(-b-+ b2 -+ 2λ) + β (-b+ -b2 + + 2λ) × e (-b+ - √ b2 + +2λ)x R e -(-b-- √ b2 -+2λ)y dy.
Now using the densities of the exponential laws with parameters

α -= -b--b2 -+ 2λ and α + = -b+ + b2 + + 2λ,
the probabilities simplifies to

P[X τ X 0 < 0|X 0 ≥ 0] =      2λ (β-1) (-b-- √ b2 -+2λ) -1 (β-1)(-b-+ √ b2 -+2λ)+β (-b+ - √ b2 + +2λ) e (-b+ - √ b2 + +2λ)x , if x ≥ 0, -2λ β (-b+ + √ b2 + +2λ) -1 (β-1)(-b-+ √ b2 -+2λ)+β (-b+ - √ b2 + +2λ) e (-b-+ √ b2 -+2λ)x , if x < 0. Hence, P[X τ ∈ dy|X τ < 0, X 0 ≥ 0] P[X τ X 0 < 0|X 0 ≥ 0] = α + e -α + y dy.
and the Algorithm 3 for sampling X τ given X 0 ≥ 0 follows, thus giving the basis of the scheme.

Data: An initial position

X 0 = x ≥ 0 at time 0. Result: The position X τ at time τ . Sample U ∼ B(P[X τ X 0 < 0|X 0 ≥ 0]). if U = 1 then X τ = E + ∼ E(α + ). else Sample X τ using (1 -P[X τ X 0 < 0|X 0 ≥ 0]) -1 λ r(λ, x, y) 1(y ≥ 0) 1(x ≥ 0). end return X τ .
Algorithm 3: Sampling X τ .

Now using the densities of

Q(α + , X 0 ), (1 -P[X τ X 0 < 0|X 0 ≥ 0]) b2 -+ 2λ (-b+ + b2 -+ 2λ) λ e -( b+ - √ b2 + +2λ)X 0 P[X τ ≥ X 0 |X 0 ≥ 0] = (1 -β)(-b-+ b2 -+ 2λ) -β (-b+ -b2 + + 2λ) (β -1)(-b-+ b2 -+ 2λ) + β (-b+ -b2 + + 2λ) e (-b+ - √ b2 + +2λ)X 0 + e (-b+ + √ b2 + +2λ)X 0 .
Henceforth, with a classical splitting of (1 -P[X τ X 0 < 0|X 0 ≥ 0]) -1 λ r(λ, x, y) 1(y ≥ 0) 1(x ≥ 0), Algorithm 4 for sampling of X τ given X τ ≥ 0 comes naturally and complete Algorithm 3. From (2.3.6),

(1 -P[X τ X 0 < 0|X 0 ≥ 0]) -1 (1 -P[X τ ≥ X 0 |X 0 ≥ 0]) 1 λ r(λ, x, y) 1(x > y ≥ 0) 1(x ≥ 0)
is bounded density with its whole mass in the (0, x). Thus, we can use a classical rejection sampling [START_REF] Devroye | Non-Uniform Random Variate Generation[END_REF] to complete Algorithm 3 and obtain the scheme.

Data: An initial position X 0 = x ≥ 0 at time 0.

Result: The position X τ at time τ given X τ ≥ 0.

Sample U ∼ B(P[X τ ≥ X 0 |X 0 ≥ 0]). if U = 1 then X τ = F ∼ Q( b+ -b2 + + 2λ, X 0 ). else Sample X τ using (1 -P[X τ X 0 < 0|X 0 ≥ 0]) -1 (1 -P[X τ ≥ X 0 |X 0 ≥ 0]) 1 λ r(λ, x, y) 1(x > y ≥ 0) 1(x ≥ 0). end return X τ . Algorithm 4: Sampling of X τ given X τ ≥ 0. Remark 2.4.1. When S = (1 -β)(-b-+ b2 -+ 2λ) -β(-b+ + b2 + + 2λ) > 0, since b2 + + 2λ P[X τ ∈ dy|X 0 ≥ 0, 0 ≤ X τ < X 0 ] (1 -P[X τ X 0 < 0|X 0 ≥ 0]) -1 (1 -P[X τ ≥ X 0 |X 0 ≥ 0]) -1 = (1 -β)(-b-+ b2 -+ 2λ) -β(-b+ + b2 + + 2λ) (β -1)(-b-+ b2 -+ 2λ) + β (-b+ -b2 + + 2λ) e -b+ (x-y) e - √ b2 + +2λ(x+y) + e (-b+ - √ b2 + +2λ)(x-y) ,
a dichotomy provides Algorithm 5 which gives an alternative to the rejection sampling.

This scheme has shown relevant performances in a few numerical experiments. Basically, it competes with the exact scheme for the Skew Brownian Motion given in [START_REF] Lejay | Simulating diffusion processes in discontinuous media : A numerical scheme with constant time steps[END_REF] (Figure 2.2 and Table 2.1) and gives a reliable approximation of the Skew Brownian Motion with a constant drift (Figure 2.1) which is the most complicated constant skew diffusion with transition function that can be drawn.

The second scheme

We consider a scheme for X with (a, ρ, b) satisfying (2.3.3) which is exact when the X stays on the positive or the negative half line and approximated when a crossing occurs. More precisely, we use an exact Euler scheme when the process starts from a point that is not 0 and Data: The position X 0 = x ≥ 0 at time 0. Result: The position X τ at time τ knowing that 0 ≤ X τ < X 0 and X 0 ≥ 0. Sample U ∼ B(P ). if U = 1 then

X τ = G 1 with E + ∼ T (-b+ -b2 + + 2λ, X 0 ). else X τ = G 2 with E + ∼ T (-b+ + b2 + + 2λ, X 0 ). end return X τ .
Algorithm 5: Sampling of X τ given 0 ≤ X τ < X 0 and S > 0. . check if the process crosses 0. If 0 is crossed, then we estimate the time s the process takes to reach 0 and sample a position using the resolvent for the time length δts. This idea actually comes from the two steps scheme developed in [START_REF] Lejay | Simulating diffusion processes in discontinuous media : A numerical scheme with constant time steps[END_REF] for the skew Brownian motion.

We start by proving the consistency of this approach. Since the local time in 0 of a semimartingale grows only when the process crosses x [73, Chapter VI, Proposition 1.3], (2.3.4) provides that X behaves exactly like a Brownian motion with drift when it stays on one side of 0.

As in §2.4, we gives a series of algorithms which must also be nested one within the others, cover only the case where X starts from non negative position and can be easily adapt in the other case.

The general structure of the scheme is given in Algorithm 6 where B b,x is Brownian Motion with drift and τ x is its first hitting time of 0. At first glance, it appears that we must sample a stopped and pinned Brownian motion with drift. This is actually very easy to do since the exact scheme given in [START_REF] Lejay | Simulating diffusion processes in discontinuous media : A numerical scheme with constant time steps[END_REF] for the case of a Brownian motion is completely suited. In effect, a bridge of a Brownian motion with drift is a simple Brownian bridge. As evidence of this statement, using a Doob's h-transform, the density of a bridge of a Brownian motion with drift is given by

q T (t, x, y) = p b (t, x, y)p b (T -t, y, z) p b (T, x, z) = √ T 2πt(T -t) exp (T (x -y) + t(z -x)) 2 2t(t -T )T = p 0 (t, x, y)p 0 (T -t, y, z) p 0 (T, x, z) ,
where p b (t, x, y) is the transition function of B b,x . For the sake of completeness, we recalled the scheme in Algorithm 8 where IG(ν, γ) denote the inverse gaussian distribution.

Data: An initial position x and a time δt > 0.

Result: A couple (τ x , 0) or (δt, B b,x δt ). Sample a realization y of a B b,x with a Euler scheme.

if xy < 0 then Generate ξ ∼ IG( |x| |y| , x 2 δt ). Set τ x = δt ξ 1+ξ . return (τ, 0) else Generate U ∼ U([0, 1]). if U < e -2xy δt then Generate ξ ∼ IG( |x| |y| , x 2 δt ). Set τ x = δt ξ
1+ξ . return (τ, 0) else return (δt, y). end end Algorithm 8: Simulation of (δt ∧ τ x , B b,x δt∧τx ). . addition, it can be interesting to compare which of these schemes is the best particularly when they are extended to handle a large number of discontinuities. This is certainly a good topics for a future work together with possible schemes which use the resolvent and are exact.

Chapter 3

Some functionals of n-dimensional skew diffusions Abstract

We propose some functionals of n-dimensional skew diffusions with piecewise constant coefficients that possess a discontinuity on the hyperplane whose vectors have their last coordinate equal to 0. Among these functionals, we give the Laplace transform of the transition function of the marginal process orthogonal to the hyperplane of discontinuity and the Laplace transform of the characteristic function of the other marginal processes. For the sake of simplicity, we make the computation in dimension two.

Introduction

We consider the continuous stochastic processes generated by the divergence form operator

L = 1 2 div(A(x)∇•) + B(x)∇• (3.0.1)
where

A =    a 1 (x 2 ) • • • 0 . . . . . . . . . 0 • • • a n (x n )    and B =    b 1 (x 2 ) . . . b n (x n )    , for x =    x 1 . . . x n    with the coefficients (a 1 , . . . , a n , b 1 , . . . , b n ) satisfying a i (x n ) = a + i 1(x n ≥ 0) + a - i 1(x n < 0), a + i , a - i > 0, i = 1, . . . , n b i (x 2 ) = b + i 1(x 2 ≥ 0) + b - i 1(x 2 < 0), (b + i , b - i ) ∈ R 2 , i = 1, . . . , n.
In one dimension, these processes are studied since the '60 through the study of the Skew Brownian Motion [START_REF] Itô | Diffusion processes and their sample paths[END_REF]. Along the years, a few numerical schemes appears. The first generation was only able to treat the second order term [START_REF] Étoré | Approximation de processus de diffusion à coefficients discontinus en dimension un et appliations à la simulation[END_REF][START_REF] Lejay | Simulating diffusion processes in discontinuous media : A numerical scheme with constant time steps[END_REF][START_REF] Martinez | Interprétations probabilistes d'opérateurs sous forme divergence et analyse de méthodes numériques probabilistes associées[END_REF][START_REF] Uffink | A random walk method for the simulation of macrodispersion in a stratified aquifer. Relation of Groundwater Quantity and Quality[END_REF], while the next one managed to deal with the first order term but only with one discontinuity [START_REF] Étoré | On the existence of a time inhomogeneous skew Brownian motion and some related laws[END_REF][START_REF] Étoré | Exact simulation of one-dimensional stochastic differential equations involving the local time at zero of the unknown process[END_REF]. Only very recently, schemes were built to handle multiple discontinuities [START_REF] Dereudre | An explicit representation of the trnaistion densities of the skew brownian motion with drift and two semipermeable barriers[END_REF][START_REF] Lenôtre | A strategy for parallel implementations of stochastic lagrangian simulation[END_REF]. Each of them share a common point: they relies on an analytical approach of stochastic processes following the seminal work of W. Feller [START_REF] Feller | Diffusion processes in one dimension[END_REF] in spite of stochastic calculus.

For the multidimensional processes, a bunch of theoretical results appears [START_REF] Portenko | Generalized diffusion processes[END_REF][START_REF] Sznitman | A multidimensional process involving local time[END_REF][START_REF] Takanobu | On the existence of solutions of stochastic differential equations with singular drifts[END_REF][START_REF] Takanobu | On the uniqueness of solutions of stochastic differential equations with singular drifts[END_REF]], but almost no numerical schemes exist. At our knowledge, we only find [START_REF] Lejay | Simulation of a stochastic process in a discontinuous layered medium[END_REF] which concern discontinuous coefficients along an hyperplane and cannot be used when b n comes into play.

In this note, we propose to investigate the case where the coefficients are discontinuous along an hyperplane using a purely analytical approach. We try to compute the marginal distribution of the process in order to ease the construction of numerical scheme for these processes and also extend the algorithm in [START_REF] Lejay | Simulation of a stochastic process in a discontinuous layered medium[END_REF]. We actually only succeed to obtain the distribution of the marginal process orthogonal to the hyperplane where the discontinuity is. Nevertheless, we compute the Laplace transform of the characteristic function of this marginal process in the general case. Finally, some of our computations in this note can clearly lead to some kind of Feynman-Kac formulae.

The reason for an analytical approach is pretty simple. In the light of [START_REF] Dereudre | An explicit representation of the trnaistion densities of the skew brownian motion with drift and two semipermeable barriers[END_REF][START_REF] Lenôtre | A strategy for parallel implementations of stochastic lagrangian simulation[END_REF], we think this kind of approaches are more suited to obtain the necessary material to build a numerical scheme, particularly in dimension more than one where the theory of local time rise problem. Moreover, we are almost convince that this work can lead after further investigation to a numerical scheme for the general case. But we also think that the approach given here can be developed to solve the case where the A is symmetric with coefficients discontinuous along an hyperplane.

This note starts by recalling in §3.1 that the stochastic processes generated by (3.0.1) are Feller processes and possess a resolvent density which is the Laplace transform of its transition function. We also show that the resolvent density is the solution of an elliptic partial differential equation. Then, for the sake of clarity, we start working only in dimension two, since the extension to the dimension n is pretty straightforward.

In §3.2 using the elliptic equation, we compute the partial Fourier transform in the coordinate parallel to the hyperplane where the discontinuity is of the resolvent density for a twodimensional process. The computational method is an adaptation of the one developed in [START_REF] Ramm | Fundamental solutions to some elliptic equations with discontinuous senior coefficients and an inequality for these solutions[END_REF] and heavily detailed in [START_REF] Elschner | Uniqueness in inverse transmission scattering problems for multilayered obstacles[END_REF] for more classical elliptic partial differential equations.

The dynamic of a process in dimension two is determined by the dynamic of the marginal processes. In §3.3 the dynamic of the process in the direction orthogonal to the hyperplane. Then, we study the dynamic of the process parallel to the hyperplane when its starts from the hyperplane in §3.4. This is fairly enough since the process behaves like a standard diffusion when it stays on one side and numerical scheme can clearly be build to obtain the mean hitting time of 0 which allows to use strategies similar to those in [START_REF] Lejay | Simulating diffusion processes in discontinuous media : A numerical scheme with constant time steps[END_REF][START_REF] Lenôtre | Two numerical schemes for the simulation of skew diffusions using their resolvent kernel[END_REF].

Divergence form operator, stochastic processes and resolvents

We consider the infinitesimal generator (L, Dom(L)) where L is defined by (3.0.1) and Dom(L) is the closure in the space C b (R n ) of continuous bounded functions vanishing when x → +∞ endowed with the supremum norm of the set

{f ∈ C 1 b (R n \p) | a + 2 (x 2 )∂ x 2 f (x 1 , 0+) = a - 2 (x 2 )∂ x 2 f (x 1 , 0), ∀x 1 ∈ R} (3.1.1)
where

p = {x ∈ R n | x = (x 1 , . . . , x n-1 , 0)} and C 1 b (R n \p) is the subspace of C b (R n )
consisting of the functions differentiable on R n \p. The existence of a Feller process associated with (L, Dom(L)) easily follows from Gaussian estimates on the semi-group [START_REF] Stroock | Diffusion semigroups corresponding to uniformly elliptic divergence form operators[END_REF] as well as the existence of the transition function p(t, x, y) solution to

∂ t p(t, x, y) = L p(t, x, y), ∀x ∈ R n , 0 ≤ t ≤ T, p(0, x, y) = δ y (x), y ∈ R n , (3.1.2)
where δ y (x) is the Dirac mass at the position y. Now the Gaussian estimates and regularity properties provided in [START_REF] Stroock | Diffusion semigroups corresponding to uniformly elliptic divergence form operators[END_REF] also allows to define the Laplace transform [START_REF] Widder | The Laplace transform[END_REF] of p(t, x, y) for any λ ∈ C such that Re λ > 0, that is r(λ, x, y) = +∞ 0 e -λt p(t, x, y) dt.

(3.1.3)

In one dimension, we easily see that r(λ, x, y) is the density of the resolvent associated with (L, Dom(L)) [START_REF] Feller | Diffusion processes in one dimension[END_REF] and the unique solution in the space H 1 (R n ) of distributions which are square integrable and whose generalized derivatives are also square integrable to λr(λ, x, y) = Lr(λ, x, y), ∀x ∈ R n , r(λ, x, y) = δ y (x), y ∈ R n .

(3.1.4)

for any λ > 0. However, in higher dimension, it is clearly less evident. Consequently, we give a sketch of the proof.

From the variational formulation of (3.1.2), the definition (3.1.3) and the estimates on p(t, x, y), we deduce that r(λ, x, y) is the unique solution for any λ > 0 to (3.1.4) belonging to H 1 (R n ).

From the regularity properties [START_REF] Evans | of Graduate studies in mathematics[END_REF][START_REF] Ladyženskaja | Equations aux dérivées partielles de type elliptique[END_REF] of the solution to (3.1.4), we can use r(λ, x, y) to define for any λ > 0 the integral operator

R λ f (x) = R n r(λ, x, y) f (y) dy, ∀f ∈ L 2 (R n ). (3.1.5)
which has an analytical continuation to the half plane in C formed by the numbers with a positive real part [START_REF] Titchmarsh | Eigenfunction Expansions, part 2[END_REF]. Now, using the Hille-Yosida theorem, we derive that the analytical prolongation of (R λ ) λ>0 is the resolvent of a semi-groups on L 2 (R n ) associated with (3.0.1). Finally, from the density of

C ∞ 0 in both L 2 (R n ) and C b (R n ),
it is easy to identify the infinitesimal generators of the process and (3.1.5) and conclude that r(λ, x, y) is the resolvent density.

Computation of the partial resolvent density

For the sake of simplicity, we give the details in two dimensions. The extension is straightforward. Moreover, we only compute the solution when a + 1 = a + 2 = a + and a - 1 = a - 2 = a -. In effect, computing the solution r(λ, x, y) to (3.1.4) for

â+ 1 = â+ 2 = â+ = a + 1 , â- 1 = â- 2 = â-= a - 1 , b1 (x) = b 1 (x) and b2 (x) = a 2 (x)b 2 (x) a 1 (x) , (3.2.1)
then using

Φ(x) = φ 1 (x 1 , x 2 ) = x 1 , φ 2 (x 1 , x 2 ) = x 2 a 1 (x) a 2 (x) , (3.2.2)
one recover the solution r(λ, x, y)

to (3.1.4) when (a 1 , a 2 , b 1 , b 2 ) by setting r(λ, x, y) = 1 det |Φ(x)| r(λ, Φ(x), Φ(y)).
From [44, §16], the unique solution in H 1 (R 2 ) to (3.1.4) is also the unique solution in

H 1 (R 2 ) to                λr(λ, x 1 , x 2 , y 1 , y 2 ) -L + r(λ, x 1 , x 2 , y 1 , y 2 ) = 0, x 1 ∈ R, x 2 ≥ 0, λr(λ, x 1 , x 2 , y 1 , y 2 ) -L -r(λ, x 1 , x 2 , y 1 , y 2 ) = 0, x 1 ∈ R, x 2 < 0, r(λ, x 1 , x 2 , y 1 , y 2 ) = δ (y 1 ,y 2 ) (x 1 , x 2 ), (y 1 , y 2 ) ∈ R 2 , r(λ, x 1 , 0+, y 1 , y 2 ) = r(λ, x 1 , 0-, y 1 , y 2 ), ∀x 1 ∈ R, a + d dxn r(λ, x 1 , 0+, y 1 , y 2 ) = a - d dxn r(λ, x 1 , 0-, y 1 , y 2 ), ∀x 1 ∈ R, (3.2.3)
where

L + = a + 2 (∂ 2 x 1 + ∂ 2 x 2 ) + b + 1 ∂ x 1 + b + 2 ∂ x 2 and L -= a - 2 (∂ 2 x 1 + ∂ 2 x 2 ) + b - 1 ∂ x 1 + b - 2 ∂ x 2 .
Consequently, depending on the sign of y 2 , r(λ, x 1 , x 2 , y 1 , y 2 ) is either the unique solution in

H 1 (R 2 ) to          λr(λ, x 1 , x 2 , y 1 , y 2 ) -L + r(λ, x 1 , x 2 , y 1 , y 2 ) = δ (y 1 ,y 2 ) (x 1 , x 2 ), x 1 ∈ R, x 2 ≥ 0, λr(λ, x 1 , x 2 , y 1 , y 2 ) -L -r(λ, x 1 , x 2 , y 1 , y 2 ) = 0, x 1 ∈ R, x 2 < 0, r(λ, x 1 , 0+, y 1 , y 2 ) = r(λ, x 1 , 0-, y 1 , y 2 ), ∀x 1 ∈ R, a + d dxn r(λ, x 1 , 0+, y 1 , y 2 ) = a - d dxn r(λ, x 1 , 0-, y 1 , y 2 ), ∀x 1 ∈ R, (3.2.4) or 
        λr(λ, x 1 , x 2 , y 1 , y 2 ) -L + r(λ, x 1 , x 2 , y 1 , y 2 ) = 0, x 1 ∈ R, x 2 ≥ 0, λr(λ, x 1 , x 2 , y 1 , y 2 ) -L -r(λ, x 1 , x 2 , y 1 , y 2 ) = δ (y 1 ,y 2 ) (x 1 , x 2 ), x 1 ∈ R, x 2 < 0, r(λ, x 1 , 0+, y 1 , y 2 ) = r(λ, x 1 , 0-, y 1 , y 2 ), ∀x 1 ∈ R, a + d dxn r(λ, x 1 , 0+, y 1 , y 2 ) = a - d dxn r(λ, x 1 , 0-, y 1 , y 2 ), ∀x 1 ∈ R, (3.2.5)
Since the solutions to (3.2.4) and (3.2.5) can be computed similarly, we treat only (3.2.4) which concern y 2 ≥ 0. In this case, with the regularity of the coefficients on each side of p in (3.2.4), we can perform a partial Fourier transform and it follows the new system

         λ + (ξ)r(λ, ξ 1 , x 2 , y 1 , y 2 ) -L + 2 r(λ, ξ 1 , x 2 , y 1 , y 2 ) = e iξ 1 y 1 δ y 2 (x 2 ), x 1 ∈ R, x 2 ≥ 0, λ -(ξ)r(λ, x 1 , x 2 , y 1 , y 2 ) -L - 2 r(λ, x 1 , x 2 , y 1 , y 2 ) = 0, x 1 ∈ R, x 2 < 0, r(λ, x 1 , 0+, y 1 , y 2 ) = r(λ, x 1 , 0-, y 1 , y 2 ), ∀x 1 ∈ R, a + d dxn r(λ, x 1 , 0+, y 1 , y 2 ) = a - d dxn r(λ, x 1 , 0-, y 1 , y 2 ), ∀x 1 ∈ R, (3.2.6)
where

λ + (ξ) = (λ + a + ξ 2 1 2 -ib + 1 ξ 1 ), λ -(ξ) = (λ + a -ξ 2 1 2 -ib - 1 ξ 1 ), L + 2 = a + 2 ∂ 2 x 2 + b + 2 ∂ x 2 and L - 2 = a - 2 ∂ 2 x 2 + b - 2 ∂ x 2 .
We assume that the unique solution in H 1 (R) to (3.2.6) can be written as

r(λ, ξ 1 , x 2 , y 1 , y 2 ) = g δ (λ, ξ 1 , x 2 , y 1 , y 2 ) + α(λ, ξ 1 )g 0 (λ, ξ 1 , x 2 , y 1 , y 2 ) (3.2.7)
where

g δ (λ, x , x n , y n ) is a solution in H 1 (R) to      λ + (ξ)g δ (λ, ξ 1 , x 2 , y 1 , y 2 ) -L + 2 g δ (λ, ξ 1 , x 2 , y 1 , y 2 ) = δ y 2 (x 2 ), x 1 ∈ R, x 2 ≥ 0, λ -(ξ)g δ (λ, x 1 , x 2 , y 1 , y 2 ) -L - 2 g δ (λ, x 1 , x 2 , y 1 , y 2 ) = 0, x 1 ∈ R, x 2 < 0, g δ (λ, x 1 , 0+, y 1 , y 2 ) = g δ (λ, x 1 , 0-, y 1 , y 2 ), ∀x 1 ∈ R, (3.2.8) g 0 (λ, ξ, x n ) is a solution in H 1 (R) to      λ + (ξ)g 0 (λ, ξ 1 , x 2 , y 1 , y 2 ) -L + 2 g 0 (λ, ξ 1 , x 2 , y 1 , y 2 ) = 0, x 1 ∈ R, x 2 ≥ 0, λ -(ξ)g 0 (λ, x 1 , x 2 , y 1 , y 2 ) -L - 2 g 0 (λ, x 1 , x 2 , y 1 , y 2 ) = 0, x 1 ∈ R, x 2 < 0, g 0 (λ, x 1 , 0+, y 1 , y 2 ) = g 0 (λ, x 1 , 0-, y 1 , y 2 ), ∀x 1 ∈ R, (3.2.9)
and α(λ, ξ 1 ) which allows r(λ, ξ 1 , x 2 , y 1 , y 2 ) to be such that

a + d dx n r(λ, x 1 , 0+, y 1 , y 2 ) = a - d dx n r(λ, x 1 , 0-, y 1 , y 2 ), ∀x 1 ∈ R.
The continuity of r(λ, ξ 1 , x 2 , y 1 , y 2 ) naturally holds since both g δ (λ, x , x n , y n ) and g 0 (λ, ξ, x n ) are continuous by Sobolev's injection of

H 1 (R) in C(R).
For g δ (λ, x , x n , y n ), we take the Green function of the differential equation defined by the extension to R of the differential operator in (3.2.8) when x 2 ≥ 0 and restrain it to the half line x 2 ≥ 0. Using the computation method in [START_REF] Lejay | One-dimensional skew diffusions : explicit expressions of densities and resolvent kernel[END_REF],

g δ (λ, ξ 1 , x 2 , y 1 , y 2 )1(x 2 ≥ 0) = exp b + 2 -(b + 2 ) 2 -a + (-2λ + 2b + 1 iξ 1 -a + ξ 2 1 ) a + y 2 × exp -b + 2 + √ (b + 2 ) 2 -a + (-2λ+2b + 1 iξ 1 -a + ξ 2 1 ) a + x 2 (b + 2 ) 2 -a + (-2λ + 2b + 1 iξ 1 -a + ξ 2 1 ) e iξ 1 y 1 1(0 ≤ x 2 < y 2 ) + exp b + 2 + (b + 2 ) 2 -a + (-2λ + 2b + 1 iξ 1 -a + ξ 2 1 ) a + y 2 × exp - b + 2 + √ (b + 2 ) 2 -a + (-2λ+2b + 1 iξ 1 -a + ξ 2 1 ) a + x 2 (b + 2 ) 2 -a + (-2λ + 2b + 1 iξ 1 -a + ξ 2 1 ) e iξ 1 y 1 1(x 2 ≥ y 2 )
Then, we complete with the unique solution in H 1 (R) [START_REF] Evans | of Graduate studies in mathematics[END_REF][START_REF] Ladyženskaja | Equations aux dérivées partielles de type elliptique[END_REF] to

     λ + (ξ)g(λ, ξ 1 , x 2 , y 1 , y 2 ) -L + 2 g(λ, ξ 1 , x 2 , y 1 , y 2 ) = 0, ∀ξ 1 ∈ R, x 2 < 0 g(λ, ξ 1 , 0-, y 1 , y 2 ) = g δ (λ, ξ 1 , 0+, y 1 , y 2 ), ∀ξ 1 ∈ R, g(λ, ξ 1 , x 2 , y 1 , y 2 ) = 0, ∀x 2 ≥ 0, where g + (λ, ξ 1 , x 2 , y 1 , y 2 ) = exp b + 2 -(b + 2 ) 2 -a + (-2λ + 2b + 1 iξ 1 -a + ξ 2 1 ) a + y 2 × exp -b + 2 + √ (b + 2 ) 2 -a + (-2λ+2b + 1 iξ 1 -a + ξ 2 1 ) a + x 2 (b + 2 ) 2 -a + (-2λ + 2b + 1 iξ 1 -a + ξ 2 1 ) e iξ 1 y 1 1(0 ≤ x 2 < y 2 ) + exp b + 2 + (b + 2 ) 2 -a + (-2λ + 2b + 1 iξ 1 -a + ξ 2 1 ) a + y 2 × exp - b + 2 + √ (b + 2 ) 2 -a + (-2λ+2b + 1 iξ 1 -a + ξ 2 1 ) a + x 2 (b + 2 ) 2 -a + (-2λ + 2b + 1 iξ 1 -a + ξ 2 1 ) e iξ 1 y 1 1(x 2 ≥ y 2 ) + exp b + 2 -(b + 2 ) 2 -a + (-2λ + 2b + 1 iξ 1 -a + ξ 2 1 ) a + y 2 × exp -b - 2 + √ (b - 2 ) 2 -a -(-2λ+2b - 1 iξ 1 -a -ξ 2 1 ) a - x 2 (b + 2 ) 2 -a + (-2λ + 2b + 1 iξ 1 -a + ξ 2 1 )
e iξ 1 y 1 1(x 2 < 0)

g -(λ, ξ 1 , x 2 , y 1 , y 2 ) = exp b - 2 -(b - 2 ) 2 -a -(-2λ + 2b - 1 iξ 1 -a -ξ 2 1 ) a - y 2 × exp -b - 2 + √ (b - 2 ) 2 -a -(-2λ+2b - 1 iξ 1 -a -ξ 2 1 ) a - x 2 (b - 2 ) 2 -a -(-2λ + 2b - 1 iξ 1 -a -ξ 2 1 )
e iξ 1 y 1 1(x 2 < y 2 )

+ exp b - 2 + (b - 2 ) 2 -a -(-2λ + 2b - 1 iξ 1 -a -ξ 2 1 ) a - y 2 × exp - b - 2 + √ (b - 2 ) 2 -a -(-2λ+2b - 1 iξ 1 -a -ξ 2 1 ) a - x 2 (b - 2 ) 2 -a -(-2λ + 2b - 1 iξ 1 -a -ξ 2 1 ) e iξ 1 y 1 1(y 2 ≤ x 2 < 0) + exp b - 2 + (b - 2 ) 2 -a -(-2λ + 2b - 1 iξ 1 -a -ξ 2 1 ) a - y 2 × exp - b + 2 + √ (b + 2 ) 2 -a + (-2λ+2b + 1 iξ 1 -a + ξ 2 1 ) a + x 2 (b - 2 ) 2 -a -(-2λ + 2b - 1 iξ 1 -a -ξ 2 1 ) e iξ 1 y 1 1(x 2 ≥ 0) h(λ, ξ 1 , x 2 , y 1 , y 2 ) = exp - b + 2 + (b + 2 ) 2 -a + (-2λ + 2b + 1 iξ 1 -a + ξ 2 1 ) a + x 2 1(x 2 ≥ 0) + exp -b - 2 + (b - 2 ) 2 -a -(-2λ + 2b - 1 iξ 1 -a -ξ 2 1 ) a - x 2 1(x 2 < 0)
where

g 2 + (λ, x 2 , y 2 ) = exp b + 2 -(b + 2 ) 2 + 2a + λ a + y 2 exp -b + 2 + √ (b + 2 ) 2 +2a + λ a + x 2 (b + 2 ) 2 + 2a + λ 1(0 ≤ x 2 < y 2 ) + exp b + 2 + (b + 2 ) 2 + 2a + λ a + y 2 exp - b + 2 + √ (b + 2 ) 2 +2a + λ a + x 2 (b + 2 ) 2 + 2a + λ 1(x 2 ≥ y 2 ) + exp b + 2 -(b + 2 ) 2 + 2a + λ a + y 2 exp -b - 2 + √ (b - 2 ) 2 +2a -λ a - x 2 (b + 2 ) 2 + 2a + λ 1(x 2 < 0) g 2 -(λ, x 2 , y 2 ) = exp b - 2 -(b - 2 ) 2 + 2a -λ a - y 2 exp -b - 2 + √ (b - 2 ) 2 +2a -λ a - x 2 (b - 2 ) 2 + 2a -λ 1(x 2 < y 2 ) + exp b - 2 + (b - 2 ) 2 + 2a -λ a - y 2 exp - b - 2 + √ (b - 2 ) 2 +2a -λ a - x 2 (b - 2 ) 2 + 2a -λ 1(y 2 ≤ x 2 < 0) + exp b - 2 + (b - 2 ) 2 + 2a -λ a - y 2 exp - b + 2 + √ (b + 2 ) 2 +2a + λ a + x 2 (b - 2 ) 2 + 2a -λ e iξ 1 y 1 1(x 2 ≥ 0) h 2 (λ, x 2 , y 2 ) = exp - b + 2 + (b + 2 ) 2 + 2a + λ a + x 2 1(x 2 ≥ 0) + exp -b - 2 + (b - 2 ) 2 + 2a -λ a - x 2 1(x 2 < 0) with Θ + (λ) = (-b + 2 + b - 2 -(b - 2 ) 2 + 2a -λ + (b + 2 ) 2 + 2a + λ) exp b + 2 - √ (b + 2 ) 2 +2a + λ a + y 2 (b + 2 ) 2 + 2a + λ (b + 2 -b - 2 + (b - 2 ) 2 + 2a -λ + (b + 2 ) 2 + 2a + λ) , Θ -(λ) = (-b + 2 + b - 2 + (b - 2 ) 2 + 2a -λ -(b + 2 ) 2 + 2a + λ) exp b - 2 + √ (b - 2 ) 2 +2a -λ a - y 2 (b - 2 ) 2 + 2a -λ(b + 2 -b - 2 + (b - 2 ) 2 + 2a -λ + (b + 2 ) 2 + 2a + λ .
Proof. Playing with the adjoint, it is actually easy to check that if we invert the sign of the coefficient in B, then r(λ, x, y) is the resolvent of another skew diffusions Then, using the characterization of the resolvent given in [55, §1], we obtain that λ r(λ, x, y) is a probability density function for any λ > 0 and x ∈ R. As a result, we can introduce a random vector (R λ,y 1 , R λ,y 2 ) whose distribution has the density λ r(λ, x, y) dx. Now using the definition of (R λ,y 1 , R λ,y 2 ),

E[e iξ 1 R λ,y 1 f (R λ,y 2 )] = R R e iξ 1 x 1 λ r(λ, x 1 , x 2 , y) dx 1 dx 2 = R λ r(λ, ξ 1 , x 2 , y) dx 2
for any f measurable and bounded. Taking r(λ, ξ 1 = 0, x 2 , y), the formula in the proposition with -b 1 and -b 2 follows. Since the resolvent is uniquely characterized by its value for λ > 0, the formula is the one of the resolvent starting from (y 1 , y 2 ). By switching the sign one more time x 2 and y 2 , we obtain the result.

We clearly see that the process X 2 does not depend on the process X 1 . Now using the formula and the change of variable provided by [START_REF] Lejay | One-dimensional skew diffusions : explicit expressions of densities and resolvent kernel[END_REF], we deduce that the process X 2 is a one dimensional skew diffusion with constant coefficients. As a result, X 2 can be simulated independently using the numerical scheme provided in [START_REF] Dereudre | An explicit representation of the trnaistion densities of the skew brownian motion with drift and two semipermeable barriers[END_REF][START_REF] Étoré | On the existence of a time inhomogeneous skew Brownian motion and some related laws[END_REF][START_REF] Étoré | Exact simulation of one-dimensional stochastic differential equations involving the local time at zero of the unknown process[END_REF][START_REF] Lejay | Simulating diffusion processes in discontinuous media : A numerical scheme with constant time steps[END_REF][START_REF] Lenôtre | Two numerical schemes for the simulation of skew diffusions using their resolvent kernel[END_REF] depending on the value of the coefficients. This is a part of the result obtain in [START_REF] Lejay | Simulation of a stochastic process in a discontinuous layered medium[END_REF].

The partial dynamic of the parallel process

We suppose here that the marginal process parallel to the hyperplane starts from the hyperplane. To compute the Laplace transform of the characteristic function r(λ, ξ 1 , y 1 ), we just have to set x 2 = 0 in the formula given in Proposition 3.2.1 and to integrate in y 2 .

After some simplifications, one derives

r(λ, ξ 1 , 0, y 1 , y 2 ) = g 1 + (λ, ξ 1 , y 1 , y 2 )1(y 2 ≥ 0) + (g 1 -(λ, ξ 1 , y 1 , y 2 )1(y 2 < 0).
where

g 1 + (λ, ξ 1 , y 1 , y 2 ) = 2 Θ -1 (λ, ξ 1 ) exp b + 2 -(b + 2 ) 2 -a + (-2λ + 2b + 1 iξ 1 -a + ξ 2 1 ) a + y 2 e iξ 1 y 1 g 1 -(λ, ξ 1 , y 1 , y 2 ) = 2 Θ -1 (λ, ξ 1 ) exp b - 2 + (b - 2 ) 2 -a -(-2λ + 2b - 1 iξ 1 -a -ξ 2 1 ) a - y 2 e iξ 1 y 1 with Θ 1 (λ, ξ) = b + 2 -b - 2 + s -+ s + ,
where

s + = (b + 2 ) 2 -a + (-2λ + 2b + 1 iξ 1 -a + ξ 2 1 ) and s -= (b - 2 ) 2 -a -(-2λ + 2b - 1 iξ 1 -a -ξ 2 1
).

Then integrating r(λ, ξ 1 , 0, y 1 , y 2 ) in y 2 , one derives

r(λ, ξ 1 , y 1 ) = 2 a + Θ -1 (λ, ξ 1 ) -b + 2 + s + e iξ 1 y 1 + 2a -Θ -1 (λ, ξ 1 ) b - 2 + s - e iξ 1 y 1 where Θ(λ, ξ 1 ) = b + 2 -b - 2 + s -+ s + .

Conclusion

We do not known of to inverse both the Laplace and Fourier transforms in the Laplace transform of the characteristic function. Nevertheless, we obtain interesting functionals. We also think that further work can lead to interesting Feynman-Kac functionals.

Chapter 4

A parallelization strategy for stochastic Lagrangian methods

Abstract

Here we present some investigations on the parallelization of stochastic Lagrangian simulations. The challenge is the proper management of the random numbers. We review two different object-oriented strategies: to draw the random numbers on the fly within each MPI's process or to use a different random number generator for each simulated path. We show the benefits of the second technique which is implemented in the PALMTREE software developed by the Project-team Sage of Inria Rennes. The efficiency of PALMTREE is demonstrated on two classical examples.

Introduction

Monte Carlo simulation is a very convenient method to solve problems arising in physics like the advection-diffusion equation with a Dirichlet boundary condition

       ∂ ∂t c(x, t) = div(σ(x) • ∇c(x, t)) -v(x)∇c(x, t)), ∀(x, t) ∈ D × [0, T ], c(x, 0) = c 0 (x), ∀x ∈ D, c(x, t) = 0, ∀t ∈ [0, T ] and , ∀x ∈ ∂D, (4.0.1)
where, for each x ∈ D, σ(x) is a d-dimensional square matrix which is definite, positive, symmetric, v(x) is a d-dimensional vector such that div(v(x)) = 0, D ⊂ R d is a regular open bounded subset and T is a positive real number. In order to have a well-posed problem [START_REF] Evans | of Graduate studies in mathematics[END_REF][START_REF] Friedman | Partial Differential Equations of Parabolic Type[END_REF] and to be able to use later the theory of stochastic differential equations, we required that σ satisfies an ellipticity condition and has its coefficients at least in C 2 (D), and that v is bounded and in C 1 (D).

Interesting computations involving the solution c(t, x) are the moments

M k (T ) = D x k c(T, x) dx, ∀k ≥ 1 such that M k (T ) < +∞.
One possibility for their computation is to perform a numerical integration of an approximated solution of (4.0.1). Eulerian methods (like Finite Difference Method, Finite Volume Method or Finite Element Method) are classical to obtain such an approximated solution. However, for advection-diffusion problems, they can induce numerical artifacts such as oscillations or artificial diffusion. This mainly occurs when advection dominates [START_REF] Hundsdorfer | Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations[END_REF].

An alternative is to use Monte Carlo simulation [START_REF] Gardiner | A Handbook for the Natural and Social Sciences[END_REF][START_REF] Zheng | Applied Contaminant Transport Modelling[END_REF] which is really simple. Indeed, the theory of stochastic processes implies that there exists X = (X t ) t≥0 whose law is linked to (4.0.1) and is such that

M k (T ) = E[X k T ]. (4.0.2)
The above expectation is nothing more than an average of the positions at time T of particles that move according to a scheme associated with the process X. This requires a large number of these particles to be computed. For linear equations, the particles do not interact with each other and move according to a Markovian process.

The great advantage of the Monte-Carlo method is that its rate of convergence is not affected by the curse of dimensionality. Nevertheless, the slowness of the rate caused by the Central-Limit theorem can be considered as a drawback. Precisely, the computation of the moments requires a large amount of particles to achieve a reliable approximation. Thus, the use of supercomputers and parallel architectures becomes a key ingredient to obtain reasonable computational times. However, the main difficulty when one deals with parallel architectures is to manage the random numbers such that the particles are not correlated, otherwise a bias in the approximation of the moments is obtained.

In this paper, we investigate the parallelization of the Monte Carlo method for the computation of (4.0.2). We will consider two implementation's strategies where the total number of particles is divided into batches distributed over the Floating Point Units (FPUs):

1. SAF: the Strategy of Attachment to the (FPUs) where each FPU received a Virtual Random Number Generator (VRNG) which is either different independent Random Number Generators (RNGs) or copies of the same RNG in different states [START_REF] L'ecuyer | Random numbers for parallel computers : Requirements and methods, with emphasis on gpus[END_REF]. In this strategy, the random numbers are generated on demand and do not bear any attachment to the particles.

2. SAO: the Strategy of Attachment to the Object where the particles carries their own Virtual Random Number Generator.

Both schemes clearly carry the non correlation of the particles assuming that all the drawn random numbers have enough independence which is a matter of RNGs. Sometimes particles with a singular behavior are encountered and the examination of the full paths of such particles is necessary. With the SAF, a particle replay requires either to re-run the simulation with a condition to record only the positions of this particle or to keep track of the random numbers used for this particle. In both cases, it would drastically increase the computational time and add unnecessary complications to the code. On the contrary, a particle replay is straightforward with the SAO.

The present paper is organized in two sections. The first one describes SAF and SAO. It also treat of the work done in PALMTREE, a library we developed with the generator RNGStreams [START_REF] L'ecuyer | An object-oriented random-number package with many long streams and substreams[END_REF] and which contains an implementation of the SAO. The second section presents two numerical experiments which illustrate the performance of PALMTREE [START_REF] Lenôtre | Palmtree library[END_REF] and the SAO. Characteristic curves like speedup and efficiency are provided for both experiment.

Parallel and Object-Oriented Implementations in Monte Carlo

All along this section, we assume that we are able to simulate the transition law of particles undergoing a Markovian dynamics such that there is no interactions between them. As a result, the presentation below can be applied to various Monte Carlo schemes involving particle tracking where the goal is to compute moments. Moreover, this shows the flexibility of the implementation we choose.

An Object-Oriented Design for Monte Carlo

C++ offers very interesting features which are of great help for a fast execution or to treat multidimensional processes. In addition, a consistent implementation of MPI is available in this language. As a result, it becomes a natural choice for PALMTREE.

In what follows, we describe and motivate the choices we made in the implementation of PALMTREE. We refer to a FPU as a MPI's process.

We choose to design an object called the Launcher which conducts the Monte Carlo simulation. Roughly speaking, it collects all the generic parameters for the simulation (the number of particles or the repository for the writing of outputs). It also determines the architecture of the computer (cartography of the nodes, number of MPI's process, etc.) and is responsible for the parallelization of the simulation (managing the VRNGs and collecting the result on each MPI's process to allow the final computations).

Some classical designs introduce an object consisting of a Particles Factory which contains all the requirements for the particle simulations (the motion scheme or the diffusion and advection coefficients). The Launcher's role is then to distribute to each MPI's process a factory with the number of particles that must be simulated and the necessary VRNGs. The main job of the factory is to create objects which are considered as the particles and to store them. Each one of these objects contains all the necessary information for path simulation including the current time-dependent position and also the motion simulation algorithm.

This design is very interesting for interacting particles as it requires the storage of the path of each particle. For the case we decide to deal with, this implementation suffers two major flaws: a slowdown since many objects are created and a massive memory consumption as a large number of objects stay instantiated.

As a result, we decide to avoid the above approach and to use a design based on recycling. In fact, we choose to code a unique object that is similar to the factory, but does not create redundant particle objects. Let us call this object the Particle.

In few words, the recycling concept is the following. When the final position at time T is reached for each path, the Particle resets to the initial position and performs another simulation. This solution avoids high memory consumption and allows complete management of the memory. In addition, we do not use a garbage collector which can cause memory leaks.

Another thing, we adopt in our design, is the latest standards in the C++11 library [START_REF]The c++ programming language[END_REF] which offers the possibility to program an object with a template whose parameter is the spatial dimension of the process we want to simulate. Thus, one can include this template parameter into the implementation of the function governing the motion of the particle. If it is, the object is declared with the correct dimension and automatically changes the function template. Otherwise, it checks the compatibility of the declared dimension with the function.

Such a feature allows the ability to preallocate the exact size required by the chosen dimension for the position in a static array. Subsequently, we avoid writing multiple objects or using a pointer and dynamic memory allocation, which provoke slowdown. Moreover, templates allow for a better optimization during the compilation. Now a natural parallel scheme for a Monte Carlo simulation consists in the distribution of a particle on the different MPI's processes. Then, a small number of paths are sequentially simulated on each MP. When each MPI's process has finished, the data is regrouped on the master MPI process using MPI communications between the MPI's processes. Thus, the quantities of This scheme is typically embarrassingly parallel and can be used with both shared or distributed memory paradigm. Here we choose the distributed memory paradigm as it offers the possibility to use supercomputers based on SGI Altix or IBM Blue Gene technologies. Furthermore, if the path of the particles needs be recorded, the shared memory paradigm can not be used due to a very high memory consumption.

Random Number Generators

The main difficulty with the parallelization of the Monte Carlo method is to ensure the independence of all the random numbers split on the different MPI's processes. To be precise, if the same random numbers are used on two different processes, the simulation will end up with non-independent paths and the targeted quantities will be erroneous.

Various recognized RNGs such as RNGStreams [START_REF] L'ecuyer | An object-oriented random-number package with many long streams and substreams[END_REF], SPRNG [START_REF] Mascagni | Algorithm 806 : Sprng : A scalable library for pseudorandom number generation[END_REF] or MT19937 [START_REF] Matsumoto | Mersenne twister : A 623-dimensionally equidistributed uniform pseudo-random number generator[END_REF] offer the possibility to use VNRGs an can be used on parallel architectures. Recently, algorithms have been proposed to produce advanced and customized VRNGs with MRG32k3a and MT19937 [START_REF] Bradley | Parallelization techniques for random number generations[END_REF].

In PALMTREE, we choose RNGStreams which possesses the following two imbricated subdivisions of the backbone generator MRG32k3a:

1. Stream: 2 127 consecutive random numbers 2. Substream: 2 76 consecutive random numbers and the VRNGs are just the same MRG32k3a in different states (See Figure 4.1). Moreover, this RNG has already implemented VRNGs [START_REF] L'ecuyer | An object-oriented random-number package with many long streams and substreams[END_REF] and passes several statistical tests which can be found in TestU01 that ensure the independence of random numbers [START_REF] Ecuyer | Testu01[END_REF]. Now a possible strategy with RNGStreams is to use a stream for each new simulation of a moment as we must have a new set of independent paths and to use the 2 51 substreams contain in each stream to allocate VRNGs on the FPU or to the objects for each moment simulation. This decision clearly avoids the need to store the state of the generator after the computations. An implementation of SAF with RNGStreams and the C++ design proposed in Subsection 4.1.1 is very easy to perform as the only task is to attach a VRNG to each MPI's process in the Launcher. Then the particles distributed on each MPI's process are simulated, drawing the random number from the attached VRNG.

Sometimes a selective replay may be necessary to capture some singular paths in order to enable a physical understanding or for debugging purposes. However, recording the path of every particle is a memory intensive task as keeping the track of the random numbers used by each particle. This constitutes a major drawback for this strategy. SAO is preferred in that case.

Strategy of Object-Attachment (SAO) and PALMTREE

Here a substream is attached to each particle which can be considered as an object and all that is needed to implement this scheme is a subroutine to quickly jump from the first substream to the nth one. We show why in the following example: suppose that we need 1,000,000 paths to compute the moment and have 5 MPI's processes, then we distribute 200,000 paths to each MPI's process, which therefore requires 200,000 VRNGs to perform the simulations (See Figure 4.2).

The easiest way to solve this problem is to have the mth FPU that starts at the (m -1) × 200, 000 + 1st substream and then to jump to the next substream until it reaches the m × 200, 000th substream.

RNGStreams possesses a function that allows to go from one substream to the next one (See Figure 4.3). Thus the only problem is to go quickly from the first substream to the (m -1) × 200.000 + 1st substream so that we can compete with the speed of the SAF.

A naive algorithm using a loop containing the default function that passes through each substream one at a time is clearly too slow. As a result, we decide to modify the algorithm for MRG32k3a proposed in [START_REF] Bradley | Parallelization techniques for random number generations[END_REF].

The current state of the generator RNGStreams is a sequence of six numbers, suppose that {s 1 , s 2 , s 3 , s 4 , s 5 , s 6 } is the start of a substream. With the vectors Y 1 = {s 1 , s 2 , s 3 } and Y 2 =

X 1 = A 1 × Y 1 mod m 1 and X 2 = A 2 × Y 2 mod m 2
with X 1 and X 2 the states providing the first number of the next substream.

As we said above, it is too slow to run these computations n times to make a jump from the 1st-substream to the nth-substream. Subsequently, we propose to use the algorithm developed in [START_REF] Bradley | Parallelization techniques for random number generations[END_REF] based on the storage in memory of already computed matrix and the decomposition s = k j=0 g j 8 j , for any s ∈ N.

Since a stream contains 2 51 = 8 17 substreams, we decide to only store the already computed matrices

A i A 2 i • • • A 7 i A 8 i A 2
A s i Y i = k j=0 A g j 8 j i Y i mod m i
This solution provides a process that can be completed with a complexity less than O(log 2 p) which is much faster [START_REF] Bradley | Parallelization techniques for random number generations[END_REF] than the naive solution. The Figure 4.4 illustrates this idea. In effect, we clearly see that the second FPU receive a stream and then performs a jump from the initial position of this stream to the first random number of the n + 1 substream of this exact same stream.

Experiments with the Advection-Diffusion Equation

The Advection-Diffusion Equation

In physics, the solution c(x, t) of (4.0.1) is interpreted as the evolution at the position x of the initial concentration c 0 (x) during the time interval [0, T ]. The first moment of c is often called the center of mass.

Let us first recall that it exists a unique regular solution of (4.0.1). Proofs can be found [START_REF] Friedman | Partial Differential Equations of Parabolic Type[END_REF][START_REF] Nash | Continuity of solutions of parabolic and elliptic equations[END_REF]. This clearly means, as we said in the introduction, that we deal with a well-posed problem.

The notion of fundamental solution [?, [START_REF] Evans | of Graduate studies in mathematics[END_REF][START_REF] Friedman | Partial Differential Equations of Parabolic Type[END_REF][START_REF] Nash | Continuity of solutions of parabolic and elliptic equations[END_REF] which is motivated by the fact that c(x, t) depends on the initial condition plays an important role in the treatment of the advection- which is easy to compute since div(v(x)) = 0 for every x ∈ R.

Assuming that σ and v satisfy the hypotheses settled in (4.0.1), then using the Feynman-Kac formula [START_REF]Stochastic Differential Equations[END_REF] and (4.2.2), we can define the process X as the unique strong solution of the Stochastic Differential Equation dX t = v(X t ) dt + σ(X t ) dB t , (4.2.3) starting at the position y and killed on the boundary D. Here, (B t ) t≥0 is a d-dimensional Brownian motion with respect to the filtration (F t ) t≥0 satisfying the usual conditions [START_REF] Revuz | Continuous martingales and Brownian motion[END_REF].

The path of such a process can be simulated step-by-step with a classical Euler scheme. Therefore a Monte Carlo algorithm for the simulation of the center of mass simply consists in the computation until time T of a large number of paths and the average of all the final positions of every simulated particle still inside the domain.

As we are mainly interested in computational time and efficiency, the numerical experiments that will follow are performed in free space. Working on a bounded domain would only require to set the accurate stopping condition, which is a direct consequence of the Feynman-Kac formula that is to terminate the simulation of the particle when it leaves the domain.

Brownian Motion Simulation

Let us take an example in dimension one. We suppose that the drift term v is zero and that σ(x) is constant. We then obtain the renormalized Heat Equation whose solution is the standard Brownian Motion. and 4.3, the nth particle is simulated by a certain FPU using the nth substream. As a result, it is easy to replay the nth particle since we just have to use the random numbers of the nth substream. The point is that the parameters must stay exactly the same particularly the time step. Otherwise, the replay of the simulation will use the same random numbers but not for the exact same call of the generator during the simulation.

Conclusion

The parallelization of Stochastic Lagrangian solvers relies on a careful and efficient management of the random numbers. In this paper, we proposed a strategy based on the attachment of the Virtual Random Number Generators to the Object.

The main advantage of our strategy is the possibility to easily replay some particle paths. This strategy is implemented in the PALMTREE software. PALMTREE use RNGStreams to benefit from the native split of the random numbers in streams and substreams.

We have shown the efficiency of PALMTREE on two examples in dimension one: the simulation of the Brownian motion in the whole space and the simulation of an advection-diffusion problem with an affine drift term. Independence of the paths was also checked.

Our current work is to perform more tests with various parameters and to link PALMTREE to the platform H2OLAB [START_REF]Project team Sage. H2olab software[END_REF], dedicated to simulations in hydrogeology. In H2OLAB, the drift term is computed in parallel so that the drift data are split over MPI's processes. The challenge is that the computation of the paths will move from one MPI's process to another which raises issues about communications, good work load balance and an advanced managemement of the VRNGs in PALMTREE.

Conclusion

In this study, we tried to give a few contributions to the problem of simulating the sample paths of skew diffusions. For the records, the skew diffusions are the stochastic processes generated by the divergence form operator

L = 1 2 div(A(x)∇•) + B(x)∇• (4.2.5)
where A is a symmetric matrix and B a vector which are both composed of functions from R to R which possess discontinuities on hyperplanes. Our first contribution concerned the dimension one and is twofold. The first part [START_REF] Lejay | One-dimensional skew diffusions : explicit expressions of densities and resolvent kernel[END_REF] furnishes a result on the resolvent density of any Feller processes [START_REF] Feller | Diffusion processes in one dimension[END_REF] and provides a computation method. Precisely, this method allows to obtain the resolvent density as at least a power series and sometimes as a closed form. For skew diffusions with piecewise constant coefficients and only one discontinuity, the method provide a closed form [START_REF] Lejay | One-dimensional skew diffusions : explicit expressions of densities and resolvent kernel[END_REF]. The second part [START_REF] Lenôtre | Two numerical schemes for the simulation of skew diffusions using their resolvent kernel[END_REF] gives two schemes based on the resolvent density. In effect, we were not able to perform a Laplace inversion of the resolvent, thus recovering the transition function of the skew diffusions with piecewise constant coefficients and one discontinuity. Now the computation method can be used to compute the resolvent density of any type of Feller processes. However, it must be done on a case-by-case basis. Consequently, a first perspective is clearly to find a family of functions such as Hermite polynomials or else which allow to universally compute the finite basis of functions required in the method. When this is done, we have a numerical scheme for simulating any Feller processes. In effect, one of the scheme we built is not restricted to skew diffusions with piecewise constant coefficients. It can be used for any Feller processes assuming the resolvent density is known as a closed form or a well truncated power series. On balance, this means that it can be built a complete automaton for the simulation of any Feller processes.

Before we move forward, we must point out that the generic algorithm [START_REF] Lenôtre | Two numerical schemes for the simulation of skew diffusions using their resolvent kernel[END_REF] is just an approximation scheme not an exact one. Subsequently, a prospect may be to construct an exact scheme using the resolvent density and some informations that can be derive from it by proving a kind of Feynman-Kac formula [START_REF] Kac | On distributions of certain wiener functionals[END_REF].

Another perspective is to extend the result on resolvent density to higher dimensions. Such a thing would be of major importance not only for the stochastic analysis but also for the study 91 of differential operators. We also think in the light of [START_REF] Boutet De Monvel | Comportement d'un opérateur pseudo-différentiel sur une variété à bord, I La propriété de transmission[END_REF][START_REF] Boutet De Monvel | Comportement d'un opérateur pseudo-différentiel sur une variété à bord, II Pseudo-noyaux de poisson[END_REF] that it will provide help for the understanding of some pseudo-differential operators.

In dimension more than one, we obtained a few results [START_REF] Lenôtre | Some functionals of n-dimensional skew diffusions[END_REF] when A is diagonal and the coefficients of both A and B have discontinuities on hyperplanes parametrized by the last coordinate. We primarily get the partial Fourier transform in the coordinates parallel to the hyperplanes where discontinuities are of the resolvent density. It is then used to compute some functionals of the marginal processes.

We decided to go for an analytical approach for a pretty simple reason. Indeed, we think that this approach can be refined to obtain a simulation algorithm even when the matrix A is full. This is an other prospect we can brought to light.

Whatever the dimension, the numerical scheme for the path simulation we developed are particularly heavy. This is clearly a problem for practical applications where a huge number of sample paths must be computed to approximate some quantities of interest such as the moments of the processes at a given time T . An answer we found to this problem consists in the parallelization of the path simulation [START_REF] Lenôtre | A strategy for parallel implementations of stochastic lagrangian simulation[END_REF]. This was a hard task mainly because of the necessity to distribute subtly the random numbers on the floating-point units to conserve their independence.

While the parallelization strategy we developed can handle the basic Monte Carlo simulation, it actually cannot take care of what we call hybrid simulation. Roughly speaking, this is when Monte Carlo and deterministic methods are brought together to solve a complex problem such as the transport of an inert solute in geophysics. As a consequence, the strategy in [START_REF] Lenôtre | A strategy for parallel implementations of stochastic lagrangian simulation[END_REF] is only a beginning. There is clearly much progress to accomplish in this wide open area.

For a piecewise smooth

  function f : R → R d , we denote by Dis(f ) the set of discontinuities of the first kind. When we extend this notation to a = (a, ρ, b), we obtain Dis(a) = Dis(a) ∪ Dis(ρ) ∪ Dis(b).

.1. 6 )

 6 To a = (a, ρ, b) ∈ C, we associate h(x) = 2 x 0 b(y) a(y)ρ(y) dy, S(x) = x 0 s(y) dy where s(x) = exp(-h(x)) a(x) (1.1.7) and M (x) = x 0 m(y) dy where m(x) = exp(h(x)) ρ(x) .

Proposition 1 . 1 . 3 .

 113 For any (λ, α, β) ∈ C 3 , there exists a unique function u ∈ S a (R, C) that solves (λ -A)u = 0 with u(0) = α and D S u(0) = β. Besides, (λ, α, β) → (u(x), D S u(x)) is holomorphic for any x ∈ R on C 3 .

Corollary 1 . 1 . 2 .

 112 The map λ → Wr[u , u ](λ, x) is holomorphic on H + and vanishes only on the closed half-line h -. Proof. First, Wr[u , u ](λ, x) does not depend on x. Second, from the very definition of u and u , Wr[u , u ](λ, x) = 0 for x ∈ h + and is holomorphic in λ. Third, Wr[u , u ](λ, x) = 0 for some λ ∈ C if and only of u = αu for α ∈ C\{0}. Necessarily, α = 1 since u (λ, 0) = u (λ, 0). Using an integration by parts on an interval [a, b], for any φ, ψ ∈ Dom(A), b a

  and are continuous. Now let us introduce the operator A φ = 1 2 D φ * M D φ * S associated to the scale function φ * S and the speed measure φ * M .

Proposition 1 . 1 . 10 .

 1110 When a ∈ P, X a is the unique strong solution to (1.1.21) with ν given by (1.1.22) and σ = √ aρ. Remark 1.1.7. For a = (a, ρ, b) and κ > 0, let us set a κ = (κa, κ -1 ρ, b). With (1.1.

2 . 7 )

 27 , as they are holomorphic whenever Θ(λ, b + , b -) = 0 and 2λ = -b 2 ± (only for λ ∈ R -) and decreasing to 0 as |λ| → ∞. Remark 1.2.4. When x = 0, (1.2.6) and (1.2.7) simplify to r(λ, 0, y) = 2βΘ(λ, b + , b -) -1 g(b + , λ, 0, y) for y > 0, 2(β -1)Θ(λ, b + , b -) -1 g(b -, λ, 0, y) for y < 0, (1.2.9) where Θ(λ, b + , b -) is given by (1.2.8).

. 1 )Figure 1 . 1 -

 111 Figure 1.1 -The density p(t, x, y) of the SBM for β = 3 4 at time t = 1 for the initial positions: x = -1, x = 0 and x = 1 (from left to right).

Figure 1 . 2 -

 12 Figure 1.2 -The density p(t, x, y) of the SBM for β = 3 4 with a constant drift b = 2 at time t = 1 for the initial positions: x = -1, x = 0 and x = 1 (from left to right).

Figure 1 . 3 -

 13 Figure 1.3 -The density p(t, x, y) of the SBM for β = 3 4 with a Bang bang drift b = 2 at time t = 1 for the initial positions x = -1, x = 0 and x = 1 (from left to right).

  .[START_REF] Abundo | First-passage problems for asymmetric diffusions and skew-diffusion processes[END_REF]

Figure 1 . 4 -Proposition 1 . 3 . 1 .

 14131 Figure 1.4 -The contour Γ.

Figure 2 . 1 -

 21 Figure 2.1 -Histograms using the first algorithm and 100,000 final positions for a skew Brownian motion and a skew Brownian motion with constant drift b = 2 with the parameter β = 0.7, 0 as starting point and T = 1 for both of them.

Figure 2 . 2 -

 22 Figure 2.2 -A Q-Q plot between the first algorithm and the SBM's exact algorithm for a skew Brownian motion with β = 0.7, 0 as starting point and T = 1..

  (a) Skew Brownian motion (b) Skew Brownian motion with constant drift

Figure 2 . 3 -

 23 Figure 2.3 -The histograms using the second algorithm and 100,000 final positions for a skew Brownian motion and a skew Brownian motion with constant drift b = 2 with the parameter β = 0.7, 0 as starting point and T = 1 for both of them.

Figure 2 .

 2 Figure 2.4 -A Q-Q plot between the second algorithm and the SBM's exact algorithm for a skew Brownian motion with β = 0.7, 0 as starting point and T = 1..

Figure 4 . 1 -

 41 Figure 4.1 -The structure of RNGStreams

Figure 4 . 2 -

 42 Figure 4.2 -Distribution of 200, 000 particles to each FPU

Figure 4 . 3 -

 43 Figure 4.3 -Distribution of VRNGs or substreams to each FPU

Figure 4 . 4 -

 44 Figure 4.4 -Illustration of the stream repartition on FPUs

  diffusion equation. It is the unique solution Γ(x, t, y) of , t, y) = div x (σ(x) • ∇ x Γ(x, t, y))v(x)∇ x Γ(x, t, y), ∀(x, t, y) ∈ D × [0, T ] × D, Γ(x, 0, y) = δ y (x), ∀(x, y) ∈ D × D, Γ(x, t, y) = 0, ∀t ∈ [0, T ], ∀y ∈ D, ∀x ∈ ∂D.

( 4 . 2 . 1 )

 421 This parabolic partial differential equation derived from (4.0.1) is often called the Kolmogorov Forward equation or the Fokker-Planck equation. The probability theory provides us with the existence of a unique Feller process X = (X t ) t≥0 whose transition function density is the solution of the adjoint of (4.2.1), that is , t, y) = div y (σ(y)• ∇ x Γ(x, t, y)) + v(y)∇ y Γ(x, t, y), ∀(x, t, y) ∈ D × [0, T ] × D, Γ(x, 0, y) = δ x (y), ∀(x, y) ∈ D × D, Γ(x, t, y) = 0, ∀t ∈ [0, T ], ∀x ∈ D, ∀y ∈ ∂D,(4.2.2)

Figure 4 . 5 -

 45 Figure 4.5 -Brownian motion: (a) The dash line represents the linear acceleration and the black curve shows the speedup. (b) The dash line represents the 100% efficiency and the black curve shows the PALMTREE's efficiency.

Figure 4 . 6 -

 46 Figure 4.6 -χ 2 test: (a) between substreams i and i + 1 for i = 0 . . . 40000 of the first stream. (b) between substreams i of the first and second streams for i = 0 . . . 10000.

  Moreover, if a sequence {a n } n∈N of elements of C converges pointwise to a ∈ C, then X an converges in distribution to X a under P x for any starting point x.

	where similar arguments are
	deployed).
	Let us set
	(1.1.2) Proposition 1.1.1. To each element a = (a, ρ, b) ∈ C is associated a conservative Feller pro-p = {(t, x, y) ∈ R 3 | t > 0, x ∈ R, y ∈ R}. cess (X a , (F t ), P x ) t≥0,x∈R with continuous paths. Its transition function admits a density with
	respect to the Lebesgue measure denoted by p(t, x, y) which satisfies a Gaussian upper bound
	and is jointly continuous on p defined by (1.1.2).

Table 2 .

 2 1 -Statistical quantities collected from simulations of the first algorithm and the Exact SBM's algorithm of 100,000 final position of a skew Brownian motion with β = 0.7, 0 as starting point and T = 1.

		approximated exact
	max	4.71137	4.80811
	75 % quantile 0.927567	0.922331
	median	0.374352	0.365781
	25 % quantile -0.176939	-0.218652
	min	-4.59835	-4.27381

Table 2 .

 2 2 -Statistical quantities collected from simulations of the first algorithm and the Exact SBM's algorithm of 100,000 final position of a skew Brownian motion with β = 0.7, 0 as starting point and T = 1.

		approximated 1 approximated 2 exact
	max	4.71137	4.07464	4.80811
	75 % quantile 0.927567	0.923473	0.922331
	median	0.374352	0.364082	0.365781
	25 % quantile -0.176939	-0.217547	-0.218652
	min	-4.59835	-4.24103	-4.27381

Remerciements

Two numerical schemes for the simulation of skew diffusions using their resolvent Proof. The result follows directly from the computation below. Remark 2.1.1. Looking at this proof, the extension to other dimensions is pretty straightforward.

From the definition of an exponential random variable, the larger λ is, the closest to its mean m = 1/λ its mass is concentrated. Hence, using Proposition 2.1.1, we can consider for λ large enough that Algorithm 2 gives a discrete approximation of a path until the time T of X.

Data: An initial position x at time 0. Result: An approximated path of X in n steps until the time T . Set X 0 = x. for i = 1, • • • , n do Sample X iT n using the density n r(n, X (i-1)T n , y).

end return (X 0 , X iT n , • • • , X T ).

Algorithm 2: A scheme using the resolvent.

An invariance principle for the numerical scheme

We denote here by a n-steps random walk on [0, T ] a sequence of n random variables indexed by n elements of [0, T ] and linearly interpolated on [0, T ]. In order to show that Algorithm 2 produces an approximation of X, we give an invariance principle for the n-steps random walk it produces. Proposition 2.2.1. Let (X n t ) 0≤t≤T be for each n ≥ 0 the n-steps random walk defined by the random variables generated by Algorithm 2. Then, (X n t ) n≥0 converges in distribution to X.

Proof. For the sake of clarity, the proof is given assuming T = 1.

For each n ≥ 0, we put τ 0 = 0 and τ k = k i=1 θ i for k = 1, . . . , n with θ 1 , ..., θ n a finite sequence of independent exponential variables of parameter n. Using this sequence, we define Now to complete the skeleton given in Algorithm 6, we set x = 0 in (2.3.5). It comes

and Algorithm 7 follows.

As for the first scheme, we perform some numerical experiments given in Data: An initial position X 0 = x > 0 and a time δt > 0.

Result: An approximated position X δt at time δt. Generate a realization (S, Y ) of (δ t ∧ τ x , B b,x δt∧τx ). if S = δt then X δt = Y . else Sample X δt using the resolvent for an exponential time step of mean 1/(δ ts) end return X δt .

Algorithm 6: The skeleton of the second scheme Data: An initial position X 0 = 0 and a time δt > 0.

Result: An approximated position X δt at time δt.

Algorithm 7: Sampling from 0.

Conclusion

While both schemes clearly answer to the problem of simulating constant skew diffusions, they are not exact scheme. Besides, the error of both schemes remains actually unknown. In which is

e iξ 1 y 1 1(x 2 < 0).

It is easy to check that our solution has square integrable derivatives and therefore belongs to

For g 0 (λ, x , x n , y n ), the computation above leads directly to

follows immediately.

where

Proof. It is enough to check that r(λ, ξ, x 2 , y 1 , y 2 ) satisfies (3.2.6) and is holomorphic on the half complex plane formed by the numbers with positive real part.

The dynamic of the orthogonal process

We denote by X = (X t ) t≥0 the stochastic process whose transition function is (3.1.2) and by X 2 = (X t ) t≥0 and X 2 = (X t ) t≥0 its associated marginal process.

We easily see in (3.1.2) that the marginal process which really matters is X 2 . We give here the Laplace transform of its transition function. Then, we compare it to a formula provided in [START_REF] Lejay | One-dimensional skew diffusions : explicit expressions of densities and resolvent kernel[END_REF] to conclude that this process is a one-dimension skew diffusion with piecewise constant coefficients.

with ∆B n = B t n+1 -B tn . In this case, the Euler scheme presents the advantage of being exact.

Since the Brownian motion is easy to simulate, we choose to sample 10,000,000 paths starting from the position 0 until time T = 1 with 0.001 as time step. We compute the speedup S and the efficiency E which are defined as

where T 1 is the sequential computational time with one MPI's process and T p is the time in parallel using p MPI's process. The speedup and efficiency curves together with the values used to plotted them are respectively given in Figure 4.5 and Table 4.1. The computations were realized with the supercomputer Lambda from the Igrida Grid of INRIA Research Center Rennes Bretagne Atlantique. This supercomputer is composed of 11 nodes with 2 × 6 Intel Xeon(R) E5647 CPUs at 2.40 Ghz on Westmere-EP architecture. Each node possesses 48 GB of Random Access Memory and is connected to the others through infiniband. We choose GCC 4.7.2 as C++ compiler and use the MPI library OpenMPI 1.6 as we prefer to use opensource and portable software. These tests include the time used to write the output file for the speedup computation so that we also show the power of the HDF5 library.

The Table 4.1 clearly illustrates PALMTREE's performance. It appears that the SAO does not suffer a significant loss of efficiency despite it requires a complex preprocessing. Moreover, the data show that the optimum efficiency (89.26%) is obtained with 24 MPI's processes.

As we mentioned in Subsection 4.1.2, the independence between the particles is guaranteed by the non correlation of random numbers generated by the RNG. Moreover, Figure 4.6 shows that the sum of the squares of the positions of the particles at T = 1 follow a χ 2 distribution in two different cases: (a) between substreams i and i + 1 for i = 0, . . . , 40000 of the first stream. (b) between substreams i of the first and second streams for i = 0, . . . , 10000.

Advection-Diffusion Equation with an Affine Drift Term

We now consider the advection-diffusion equation whose drift term v is an affine function, that is for each x ∈ R, v(x) = a x + b and σ is a constant. We simulate the associated stochastic process X through the exact scheme

where N (0, 1) is a standard Gaussian law [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF].

For this scheme with an initial position at 0 and the parameters σ = 1, a = 1, b = 2 and T = 1, we give the speedup and efficiency curves represented in Figure 4.7 based on the simulation of hundred millions of particles. The Table 4.2 provides the data resulting from the simulation and used for the plots.

Whatever the number of MPI's processes involved, we obtain the same empirical expectation E = 3.19 and empirical variance V = 13.39 with a standard error S.E. = 0.0011 and a confidence interval C.I. = 0.0034. Moreover, a good efficiency (89.29%) is obtained with 60 MPI's processes.

In this case, the drift term naturally pushes the particles out of 0 relatively quickly. If this behavior is not clearly observed in a simulation, then the code has a bug and a replay of a selection of few paths can be useful to track it in spite of reviewing all the code. This can clearly save time.

With the SAO, this replay can be easily performed since we know which substreams is used by each particle as it is shown in Figure 4.4. Precisely, in the case presented in Figure 4.2

Résumé

Nous considérons les processus de diffusion biaisés et leur simulation. Notre étude se divise en quatre parties et se concentre majoritairement sur les processus à coefficients constants par morceaux dont les discontinuités se trouvent le long d'un hyperplan simple. Nous commençons par une étude théorique dans le cas de la dimension un pour une classe de coefficients plus large. Nous donnons en particulier un résultat sur la structure des densités des résolvantes associées à ces processus et obtenons ainsi une méthode de calcul. Lorsque cela est possible, nous effectuons une inversion de Laplace de ces densités et donnons quelques fonctions de transition. Nous nous concentrons ensuite sur la simulation des processus de diffusions baisées. Nous construisons un schéma numérique utilisant la densité de la résolvante pour tout processus de Feller. Avec ce schéma et les densités calculées dans la première partie, nous obtenons une méthode de simulation des processus de diffusions biaisées en dimension un. Après cela, nous regardons le cas de la dimension supérieure. Nous effectuons une étude théorique et calculons des fonctionnelles des processus de diffusions biaisées. Ceci nous permet d'obtenir entre autre la fonction de transition du processus marginal orthogonal à l'hyperplan de discontinuité. Enfin, nous abordons la parallélisation des méthodes particulaires et donnons une stratégie permettant de simuler de grand lots de trajectoires de processus de diffusions biaisées sur des architectures massivement parallèle. Une propriété de cette stratégie est de permettre de resimuler quelques trajectoires de précédentes simulations.

Abstract

We consider the skew diffusion processes and their simulation. This study are divided into four parts and concentrate on the processes whose coefficients are piecewise constant with discontinuities along a simple hyperplane. We start by a theoretical study of the one-dimensional case when the coefficients belong to a broader class. We particularly give a result on the structure of the resolvent densities of these processes and obtain a computational method. When it is possible, we perform a Laplace inversion of these densities and provide some transition functions. Then we concentrate on the simulation of skew diffusions process. We build a numerical scheme using the resolvent density for any Feller processes. With this scheme and the resolvent densities computed in the previous part, we obtain a simulation method for the skew diffusion processes in dimension one. After that, we consider the multidimensional case. We provide a theoretical study and compute some functionals of the skew diffusions processes. This allows to obtain among others the transition function of the marginal process orthogonal to the hyperplane of discontinuity. Finally, we consider the parallelization of Monte Carlo methods. We provide a strategy which allows to simulate a large batch of skew diffusions processes sample paths on massively parallel architecture. An interesting feature is the possibility to replay some the sample paths of previous simulations.