Discrétisation de processus stochastiques, estimées de

densités et applications
Stephane Menozzi

To cite this version:

Stephane Menozzi. Discrétisation de processus stochastiques, estimées de densités et applications.
Mathématiques [math]. Université Paris-Diderot - Paris VII, 2010. �tel-00533333�

HAL Id: tel-00533333

[Link]
Submitted on 5 Nov 2010

HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est

archive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents
entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non,
lished or not. The documents may come from émanant des établissements d’enseignement et de
teaching and research institutions in France or recherche français ou étrangers, des laboratoires
abroad, or from public or private research centers. publics ou privés.
 Habilitation à Diriger des Recherches en Sciences

Spécialité : Mathématiques Appliquées

présentée par

Stéphane Menozzi

à l’Université Paris 7-Denis Diderot

Discrétisation de processus stochastiques,

estimées de densité et applications.

Soutenue le 10 Novembre 2010 devant le jury composé de :

Vlad BALLY Professeur à l’Université de Marne la Vallée-Paris Est

Patrick CATTIAUX Professeur à l’Université Toulouse 3-Paul Sabatier (Rapporteur)
Josselin GARNIER Professeur à l’Université Paris 7-Denis Diderot (Rapporteur)
Emmanuel GOBET Professeur à l’Ecole Polytechnique
Sylvie MÉLÉARD Professeur à l’Ecole Polytechnique
Gilles PAGÈS Professeur à l’Université Paris 6-Pierre et Marie Curie
Alexander VERETENNIKOV Professeur à l’Université de Leeds (Rapporteur)
Présentation du manuscrit

Ce mémoire présente et met en perspective l’ensemble des travaux de recherche effectués pendant ma thèse,
de 2001 à 2004, à l’Université Pierre et Marie Curie (Paris VI) et à l’Ecole Polytechnique, sous la direction de
Vlad Bally et Emmanuel Gobet, et par la suite depuis mon recrutement en 2005 à l’Université Denis Diderot
(Paris VII).
Ces travaux se sont développés autour de deux axes principaux d’intersection non vide : le premier relatif
à la discrétisation de processus stochastiques, le second à l’estimation de densité de processus dégénérés et à
certaines applications associées.

Les quatre premiers chapitres relèvent majoritairement du premier axe. Le Chapitre 1 est introductif. Nous
y présentons les principaux types de problèmes associés à la discrétisation de processus, les résultats existants
ainsi que les outils utilisés pour les obtenir. Nous détaillons plus particulièrement les techniques employées pour
analyser l’approximation faible.
Le Chapitre 2 regroupe les résultats des articles [1,3,5,7]. Nous y étudions principalement l’approximation
faible d’un processus de diffusion tué/stoppé à la sortie d’un domaine régulier [1,7]. Sont également abordés le
cas du cône pour le mouvement Brownien [3], ainsi qu’une extension de l’analyse de l’erreur à des processus
non Markoviens [5].
Au Chapitre 3, nous rappelons les principes de la méthode Parametrix qui permet d’obtenir une représentation
de la densité de transition d’une équation différentielle stochastique (EDS) comme somme d’une série où cha-
cun des termes est la convolution de la densité d’une EDS à coefficients constants et des convolutions itérées
d’un noyau régularisant. En appliquant cette approche à temps discret nous en déduisons un théorème limite
local pour l’approximation par chaı̂ne de Markov d’un processus de diffusion dégénéré [8] et un développement
d’erreur sur les densités pour une EDS dirigée par un processus de Lévy stable symétrique et approchée par son
schéma d’Euler [9].
Après une brève présentation des Equations Différentielles Progressives Rétrogrades (EDSR) et de leurs
connexions avec certaines équations aux dérivées partielles non linéaires, nous donnerons au Chapitre 4 un
aperçu des méthodes numériques développées pour leur approximation. Nous détaillerons plus particulièrement
les articles [2,4] qui traitent du cas quasi-linéaire (EDSR fortement couplées) qui constitue le cadre le plus
général pour les applications en physique ou en finance. Enfin, nous exposerons certaines propriétés de schémas
de semi-discrétisation pour des EDSR tuées à la sortie d’un domaine, article [6].
Le Chapitre 5 relève lui majoritairement du deuxième axe, i.e. l’estimation de densité. Il débute par une
partie introductive destinée à présenter une synthèse des techniques intervenant dans l’estimation quantitative
de densités de processus stochastiques ou en termes analytique de solutions fondamentales d’équations aux
dérivées partielles. Nous y présentons ensuite les résultats de l’article [10]. Nous obtenons des estimées de type
Aronson pour la densité de transition d’une équation différentielle stochastique dont la dynamique peut se voir
comme celle d’une équation différentielle ordinaire (EDO) dont la première composante serait perturbée par
un bruit Brownien. Bien sûr, pour que la densité existe, les coefficients de l’EDO doivent vérifier une condition
de non dégénérescence de type Hörmander faible dans la mesure où seule la première composante est non
dégénérée. Ces résultats font intervenir contrôle stochastique et techniques analytiques, méthode parametrix
présentée au Chapitre 3. Nous concluons le chapitre par une application des estimées de densité : les contrôles
non asymptotiques d’intervalles de confiance dans les méthodes de Monte Carlo. Nous montrons pour des
schémas d’Euler, éventuellement dégénérés, des estimées d’Aronson et nous utilisons la domination Gaussienne
ainsi que des techniques de concentration de la mesure : inégalités de Sobolev logarithmique, argument de Herbst
pour arriver au résultat, cf. [11].
La note [12] établit une condition suffisante pour avoir l’unicité en loi pour les systèmes dégénérés considérés
au Chapitre 5 et abordés dans les articles [10],[11]. Les idées de ce travail sont présentées à l’Appendice A. En
particulier, ces résultats apportent des éléments de réponse à une question soulevée par le rapport d’Alexander
Veretennikov.

1
 Travaux présentés pour l’habilitation

[1] E. Gobet et S. Menozzi.“Exact Approximation rate of killed hypoelliptic diffusions using the discrete
Euler scheme”. Stochastic Processes and their Applications, 112-2 (2004), 201-223.
[2] F. Delarue et S. Menozzi. “A forward backward algorithm for quasi-linear PDEs”.
Annals of Applied Probability, 16-1 (2006), 140-184.
[3] S. Menozzi. “Improved simulation for the killed Brownian motion in a cone”.
SIAM Journal on Numerical Analysis, 44-6 (2006), 2610-2632.
[4] E. Gobet et S. Menozzi. “Discrete Sampling of functionals of Itô processes”.
Séminaire de Probabilités, XL (2007), 355-375.
[5] F. Delarue et S. Menozzi. “An interpolated Stochastic Algorithm for Quasi-Linear PDEs”
Mathematics of Computation, 261-77 (2008), 125-158.
[6] B. Bouchard and S. Menozzi. “Strong Approximations of BSDEs in a domain”.
Bernoulli, 15-4 (2009), 1117-1147.
[7] E. Gobet et S. Menozzi. “Stopped diffusion processes : overshoots and boundary correction”.
Stochastic Processes and their Applications, 120-2 (2010), 130-162.
[8] V. Konakov, S. Menozzi et S. Molchanov. “Explicit parametrix and local limit theorems for some de-
-generate diffusion processes”. Annales de l’Institut Henri Poincaré, Série B, Vol. 46-4 (2010), 908-923.
[9] V. Konakov et S. Menozzi.“Weak Error for stable driven SDEs : expansion of the densities”.
A paraı̂tre dans Journal of Theoretical Probability, 25 pages.
[10] F. Delarue et S. Menozzi. “Density Estimates for a Random Noise Propagating through a Chain of
Differential Equations”. Journal of Functional Analysis, 259-6 (2010), 1577-1630.
[11] V. Lemaire et S. Menozzi. “On some Non Asymptotic Bounds for the Euler Scheme”.
Electronic Journal of Probability, Volume 15 (2010), 1645-1681.
Travail non présenté
[12] S. Menozzi. “Parametrix techniques and martingale problem for some degenerate Kolmogorov’s equations”.

2
Chapitre 1

Discrétisation de processus : un bref

tour d’horizon

Dans ce chapitre nous présentons brièvement différents problèmes associés à la discrétisation d’un processus
stochastique de “diffusion-sauts” en distiguant tout d’abord approximation forte et faible. Ces approximations
sont de nature très différentes, tant en terme des résultats de convergence auxquels elles permettent d’aboutir,
que dans les techniques mathématiques utilisées pour les établir.
Nous abordons premièrement le cas de l’approximation forte qui est celle associée trajectoriellement à la
différence entre le processus initial et une approximation donnée. Les contrôles que l’on souhaite obtenir pour
cette erreur sont soit presque sûrs soit en norme Lp . Pour les démontrer, on a recours aux outils de l’analyse
stochastique : développements d’Itô-Taylor, inégalités de martingales de type Burkholder Davies Gundy (BDG),
associés à des arguments de stabilité semblables à ceux utilisés pour l’approximation d’équations différentielles
ordinaires.
Nous poursuivons ensuite par l’approximation faible. Dans ce cadre, on ne souhaite plus avoir un contrôle
trajectoriel, mais proposer un schéma d’approximation pour l’espérance d’une fonctionnelle d’un processus.
Cette fonctionnelle peut être par exemple une fonction du processus observé à une date donnée fixée ou à un
temps d’arrêt, une fonction du maximum de l’une des composantes ou de l’intégrale du processus sur la période.
Le travail fondateur de Talay et Tubaro [TT90] a mis en évidence l’importance fondamentale pour l’analyse de
ce type d’erreur dans le cas Markovien, du lien entre la quantité à approcher et l’équation aux dérivées partielles
dont cette quantité n’est autre que la représentation probabiliste de la solution, formule de Feynman-Kac.
Nous avons essayé de rassembler quelques résultats majeurs de ces vingt dernières années qui seront utiles à
la compréhension des sujets abordés dans les chapitres suivants. Nous insistons particulièrement sur l’approxi-
mation faible et les techniques d’étude de l’erreur associée en distinguant les deux types d’hypothèses requises :
régularité de la fonctionnelle et des coefficients ou non dégénerescence. En effet, cette thématique sera l’objet
d’étude des Chapitres 2 (cas des processus stoppés) et 3 (processus dégénérés et de certains processus à sauts).
Nous conclurons par une section dédiée à quelques autres problèmes de discrétisation : l’approximation en
loi, on entend par là la caractérisation en terme de dynamique d’une renormalisation convenable de la différence
entre le processus initial et le schéma de discrétisation, l’approximation de mesures invariantes et de dynamiques
non linéaires où les contrôles sur l’erreur forte interviennent de nouveau.

1.1 Modèle initial et discrétisation

Soit une équation différentielle stochastique (EDS) à valeurs dans Rd de dynamique
Z t Z t
Xt = x + b(s, Xs− )ds + σ(s, Xs− )dZs , t ≥ 0, (1.1)
0+ 0+

où les coefficients b : R+ × Rd → Rd , σ : R+ × Rd → Rd ⊗ Rr sont mesurables bornés en temps, Lipschitziens en

espace et où (Zt )t≥0 est un processus de Lévy r-dimensionnel de décomposition de Lévy-Khintchine
 Z 
|u|2

E[exp(ihu, Zt i)] = exp(tψ(u)), ψ(u) = exp − + exp(ihu, zi) − 1 − ihu, ziI|z|≤1 ν(dz) . (1.2)
2 Rr

Nous supposons sans perte de généralité que la partie de dérive est nulle dans (1.2) dans la mesure où nous
l’avons introduite explicitement dans la dynamique de (1.1). De même, la partie diffusive non normalisée et/ou

3
anisotrope de Z peut être prise en compte
R dans le coefficient σ de (1.1). Enfin, la mesure de Lévy ν de Z vérifie
la condition d’intégrabilité usuelle Rr (1 ∧ |z|2 )ν(dz) < +∞.
Ainsi l’équation (1.1) est dirigée par un processus Z qui s’exprime comme la somme d’un mouvement Brow-
nien standard W et d’un processus de saut pur N de mesure de Lévy ν.

L’EDS (1.1) intervient dans de nombreux domaines d’applications. En particulier, des liens profonds existent
entre équations différentielles stochastiques et équations aux dérivées partielles (EDP) via les formules de
représentation de Feynman Kac, voir Freidlin [Fre85] ou Friedman [Fri75]. Ainsi, dans de nombreux cas détaillés
ci-après, cf. Section 1.3, il est possible d’exprimer la solution d’une EDP linéaire, parabolique ou elliptique,
comme espérance d’une fonctionnelle d’un processus de dynamique (1.1).
Dans le cas où le processus directeur Z = W , le processus X solution de (1.1) interviendra dans les
représentations probabilistes des équations de diffusion de la chaleur en milieu inhomogène ou des équations
de la neutronique. Le cas où Z est un processus de saut pur N intervient quant à lui de façon naturelle dans
la modélisation du trafic dans des réseaux informatiques ou routiers. Enfin, en mathématiques financières,
l’équation (1.1) décrit la dynamique d’un actif dont les variations peuvent faire intervenir composantes brow-
niennes et de sauts, ces dernières intervenant de façon significative dans la modélisation en cette période de
crise. Dans de nombreux cas, le prix d’une option sur l’actif X peut là aussi s’exprimer comme l’espérance d’une
fonctionnelle de X.

Il est donc naturel d’introduire une procédure d’approximation de (1.1) qui de façon générale ne peut être
simulée de façon exacte. Mentionnons à cet égard dans le cas monodimensionnel continu non dégénére, i.e.
d = r = 1, Z = W et σ ≥ σ0 > 0, le travail de Beskos et al. [BPR06] qui permet une simulation exacte.
En toute généralité l’approximation la plus naturelle de (1.1) est celle du schéma d’Euler, qui par ailleurs
fournit dans certains cas un procédé constructif de résolution de (1.1), cf. Revuz et Yor [RY99] pour le cas
continu. Pour un pas de temps fixé h > 0 on pose pour i ∈ N, ti := ih, et l’on définit
Z t Z t
Xth = x + h
b(φ(s), Xφ(s) )ds + h
σ(φ(s), Xφ(s) )dZs , t ≥ 0, (1.3)
0 0

où φ(s) = sup{ti : ti ≤ s}. En d’autres termes, pour une grille temporelle donnée (ti )i∈N∗ , le schéma d’Euler
est obtenu en gelant sur chaque pas de temps la dynamique spatiale et temporelle aux points de discrétisation
en temps. Précisons que dans le cas où Z = W + N possède une composante de sauts, le schéma d’Euler n’est
pas toujours exactement simulable, cf. Section 1.3.2. Il le reste en revanche pour de nombreux cas importants
pour les applications : si N est un processus stable, Γ, ...
Pour la suite nous allons distinguer deux grands types de résultats d’approximation qui font intervenir
(Xt )t≥0 solution de (1.1) et son approximation (Xth )t≥0 solution de (1.3). Historiquement, l’approximation
de processus remonte à Stroock et Varadhan qui à la suite de l’étude des problèmes de martingales s’étaient
intéréssés à des approximations en loi des processus de diffusions avec des contrôles “à la Donsker”, i.e. conver-
gence sans vitesse de l’approximation. Pour un aperçu de cette approche, voir le Chapitre 11 de [SV79] pour les
approximations dans tout l’espace ou l’article [SV71] pour les problèmes de Dirichlet.
Il a fallu attendre les travaux de Milstein [Mil74, Mil76] pour que des approximations de type Itô-Taylor
permettent de quantifier plus précisément l’erreur forte trajectorielle. C’est ensuite dans les années 80 que
Milstein [Mil85] et Talay [Tal86] ont montré que l’erreur faible était d’ordre h sous de bonnes hypothèses.

1.2 Erreur forte

On désigne par erreur forte associée au schéma d’Euler la quantité supt∈[0,T ] |Xt − Xth | pour T > 0 fixé.
Dans le cas d’un processus de diffusion, si l’on suppose que les coefficients b, σ dans (1.1) sont uniformément
Lipshitz continus en espace et 1/2-Hölder en temps, des développements d’Itô-Taylor permettent de déduire
que pour tout p ≥ 1, ∃C := C(p, T, b, σ),
√
E[ sup |Xt − Xth |p ]1/p ≤ C h, (1.4)
t∈[0,T ]

où la dépendance en b, σ dans la constante C ci-dessus intervient via les modules de continuité en temps et
espace des coefficients. Par rapport à l’analyse de la discrétisation d’une EDO, ce sont les mêmes arguments de
stabilité qui sont utilisés avec une étape supplémentaire pour contrôler les intégrales stochastiques à l’aide des
inégalités BDG. C’est lors de cette étape qu’apparaı̂t l’ordre de l’erreur en h1/2 qui correspond à la norme Lp de
l’accroissement Brownien sur un pas de temps. L’utilisation du Lemme de Gronwall pour conclure nous indique

4
ensuite qu’en toute généralité ces estimations ne seront pertinentes que pour un horizon de temps T > 0 fini,
du fait de la dépendance exponentielle de C en ce paramètre. Un corollaire immédiat de (1.4) et du lemme de
Borel Cantelli est la convergence presque sûre de h−α supt∈[0,T ] |Xt − Xth | vers 0 avec h pour tout α < 1/2.
C’est pour compenser cette vitesse que Milstein a introduit son schéma [Mil76] qui permet de retrouver la
vitesse usuelle en h. En revanche, ce schéma est difficilement simulable en dimension d > 1, sans condition
supplémentaire de commutation sur les colonnes de la matrice de diffusion identifiées à des champs de vecteurs.
Il existe de nombreux schémas d’ordre élévés pour l’erreur forte. Nous renvoyons à l’ouvrage de Kloeden et
Platten [KP92] pour un panorama de ces schémas. Mentionnons toutefois que leur simulation reste beaucoup
plus complexe que celle du schéma d’Euler et que ces schémas ne sont pas forcément les plus appropriés dans
de nombreux cas, voir Section 1.3.
Des arguments similaires peuvent être utilisés pour caractériser l’erreur forte, avec la même vitesse, lorsque
Z est un processus de Lévy de carré intégrable, voir Kohatsu-Higa et Protter [KHP94]. Mentionnons enfin
les travaux de Bruti Liberati et Platten [BLP07] relatifs à la construction de schémas d’ordre plus élévé pour
l’erreur forte dans le cas mixte Z = W + N ou de saut pur Z = N lorsque l’intensité du processus de saut
est finie, i.e. ν(Rr ) < +∞ dans (1.2). Ces auteurs utilisent en particulier l’idée d’une grille de discrétisation
temporelle adaptée aux instants de saut.

1.3 Erreur faible

Dans de nombreuses applications : approximation probabiliste de la solution d’une EDP linéaire, valorisation
ou couverture d’option en finance, il n’est toutefois pas nécéssaire d’avoir un contrôle trajectoriel sur l’approxi-
mation du processus. En effet, dans ces cas là les quantités à estimer sont de la forme E[Ψ ((Xt )t≥0 )], où
Ψ : C(R+ , Rd ) → R est une fonctionnelle de la trajectoire du processus. La représentation en terme d’espérance
de la quantité à estimer tend naturellement à faire envisager une méthode de Monte Carlo pour l’approxima-
tion sous réserve d’intégrabilité. Par ailleurs, outre le problème déjà évoqué de l’impossibilité de simuler de
h
façon exacte la diffusion, s’ajoute
 ici celui
 de l’approximation de la fonctionnelle Ψ. Si l’on désigne par Ψ une
approximation de Ψ et par (Xth,i )t≥0 une suite i.i.d. de trajectoires du schéma d’Euler, l’estimateur de
i∈N∗  
1
PM h,i
Monte Carlo naturel s’écrit pour tout M ≥ 1, ΘhM := M i=1 Ψ h
(X t )t≥0 . Ainsi, l’erreur globale entre la
quantité explicitement simulable sur ordinateur et la quantité exacte se décompose comme :
h





EM = ΘhM − E[Ψh (Xth )t≥0 ] + E[Ψh (Xth )t≥0 ] − E[Ψ ((Xt )t≥0 )]
S
:= EM,h + EhD . (1.5)
S


Le terme EM,h est un terme d’erreur statistique. Si suph>0 Ψh (Xth )t≥0 ∈ L2 (P), ce terme sera contrôlé par le
√

théorème central limite à l’ordre σ h / M avec forte probabilité, où l’écart-type σ h := Var(Ψh (Xth )t≥0 )1/2 est
en général assez stable, ou tout du moins uniformément dominé, en h. Nous reviendrons au Chapitre 5.3 sur les
contrôles non asymptotiques d’intervalles de confiance.
L’étude de l’erreur faible consiste à analyser le comportement du terme EhD . Ce terme correpond à l’erreur
de discrétisation : du processus, et de la fonctionnelle. Pour évaluer la pertinence de l’approximation choisie, ou
mettre en évidence ses limites comme pour le problème de Dirichlet, cf. Section 1.3.3 et Chapitre 2, la première
étape consiste à quantifier en fonction de h la convergence de EhD vers 0. Le découpage (1.5) permet ensuite de
déterminer la précision globale de l’approximation numérique. Il est alors possible d’équilibrer de façon optimale
les paramètres M de Monte Carlo et h de discrétisation de sorte à atteindre un seuil fixé.
Si l’on arrive par ailleurs à obtenir un développement asympotique de EhD en fonction de h, il est alors
possible, à complexité équivalente, d’utiliser des techniques usuelles et purement algébriques d’extrapolation
pour accélérer de façon significative la convergence de EhD . Cette démarche, connue sous le nom d’extrapolation
de Richardson-Romberg est présentée en détails dans l’article [TT90]. Nous la rappelons en Section 1.3.1.
Dans la suite de ce chapitre nous allons rassembler des résultats relatifs à trois grands types de fonctionnelles :
1. ψ ((Xt )t≥0 ) = g(XT ) pour T > 0 fixé et une fonction g donnée.
2. ψ ((Xt )t≥0 ) = g(τ ∧ T, Xτ ∧T ), τ := inf{t ≥ 0 : Xt 6∈ D}, pour D domaine de Rd .
RT
3. ψ ((Xt )t≥0 ) = g(T −1 0 Xt dt, XT ).
La technique eprouvée dans le cas Markovien pour l’analyse de l’erreur EhD , consiste à identifier la quantité
à estimer E[Ψ ((Xt )t≥0 )] à la solution d’une EDP linéaire avec conditions aux limites appropriées en fonction
du problème considéré. Dans le cas 1. ci-dessus nous serons dans le cadre du problème de Cauchy. Dans le cas
2. du problème de Cauchy-Dirichlet et dans le cas 3. toujours dans le cadre du problème de Cauchy mais pour

5
un opérateur dégénéré. L’intégrale est vue comme une variable supplémentaireR t à laquelle ne sera pas associé de
terme d’ordre 2 dans le générateur infinitésimal du couple (Xt , Yt )t≥0 , Yt := 0 Xs ds.

1.3.1 Problème de Cauchy : cas continu

Nous considérons ici d’abord le cas où X dans (1.1) est un processus de diffusion à coefficients réguliers, i.e.
b, σ sont supposés C ∞ et à dérivées bornées et ψ ((Xt )t≥0 ) = g(XT ) où g est une fonction C ∞ à croissance poly-
nomiale. L’idée originale des articles [Tal86] et [TT90], consiste tout d’abord à observer que sous les hypothèses
précédentes, la fonction u(t, x) := E[g(XTt,x )] := E[g(XT )|Xt = x] définie pour (t, x) ∈ [0, T ] × Rd satisfait l’EDP
parabolique 
(∂t + Lt )u(t, x) = 0, (t, x) ∈ [0, T ) × Rd ,
(1.6)
u(T, x) = g(x),
où Lt ϕ(x) désigne le générateur infinitésimal de (1.1). Pour ϕ ∈ C0∞ (Rd ), Lt ϕ(x) = hb(t, x), ∇x ϕ(x)i +
1 2 ∗
2 tr(a(t, x)Dx ϕ(x)), a(t, x) := σσ (t, x). Ces hypothèses permettent ensuite de déduire des contrôles sur les
dérivées de u à l’aide d’estimées Lp sur les dérivées du flot de (1.1). Précisément, pour tout multi-indice de
dérivation α, il existe kα (T ), Cα (T ) > 0 tels que pour tout (t, x) ∈ [0, T ] × Rd , |∂α u(t, x)| ≤ Cα(T ) (1 + |x|kα (T ) ).
C’est à dire que les dérivées gardent une croissance polynomiale en espace, ce qui permet en particulier de déduire
des contrôles usuels sur le schéma d’Euler que pour tout p ≥ 1, i ∈ [[0, N ]], T := N h, |∂α u(ti , Xthi )| ∈ Lp (P).
L’étape suivante consiste à utiliser (1.6) pour réécrire l’erreur de discrétisation comme :
N
X −1
EhD = Ex [g(XTh )] − Ex [g(XT )] = Ex [u(T, XTh )] − u(0, x) = E[u(ti+1 , Xthi+1 ) − u(ti , Xthi )].
i=0

A partir de cette dernière formulation, l’idée est ensuite d’appliquer un développement de Taylor en temps-
espace autour de (ti , Xthi ) pour chaque i ∈ [[0, N − 1]] grâce à la régularité a priori de u. En exploitant que u
satisfait (1.6), on obtient que les termes d’ordre h s’annulent dans E[u(ti+1 , Xthi+1 ) − u(ti , Xthi )] := EhD,i . Pour
établir ensuite que EhD,i = O(h2 ) il est nécessaire de poursuivre le développement à l’ordre deux en temps et
quatre en espace, en exploitant que l’on peut éliminer les moments d’ordre 3 par centrage. Ainsi, si l’on a a
priori que u ∈ C 2,4 ([0, T ] × Rd ), on déduit que ED
h
= O(h) par simple sommation.
Notons que dans la procédure décrite ci-avant, il n’est pas nécessaire d’avoir un aléa Gaussien dans le
schéma d’Euler (1.3). En effet, le développement ne fait intervenir que les propriétés de centrage des moments
d’ordre 1 et 3 et que la covariance de l’aléa est la matrice identité. Cette dernière observation permet de gagner
en complexité en utilisant par exemple des lois de Bernoulli convenablement normalisées en lieu et place des
accroissements Browniens.
Si l’on poursuit maintenant les développements à l’ordre 3 en temps et 6 en espace, en tuant toujours par
centrage les moments impairs des accroissements, on aboutit au développement d’erreur :

EhD = Ch + O(h2 ), (1.7)

où la constante C := C(T, a, b, u, x), cf. Théorème 1 de [TT90] pour une expression explicite. Du développement
D
précédent il est aisé de voir que si l’on considère le schéma d’Euler de pas h/2, il vient : Eh/2 = C h2 + O(h2 ).
Une simple manipulation algébrique donne ensuite
n o
D h/2
2Eh/2 − EhD = 2Ex [g(XT )] − Ex [g(XTh )] − Ex [g(XT )] = O(h2 ). (1.8)
nP  o
M h/2,i
Cette identité suggère d’utiliser l’estimateur de Monte Carlo Θ̂hM = 1
M i=1 2g(XT ) − g(XTh,i ) , où les
h/2,i h/2,i
(XT , XTh,i )i∈[[1,M]] sont des réalisations i.i.d. et pour tout i ∈ [[1, M ]] les XT , XTh,i sont simulés à partir des
mêmes accroissements Browniens. Ce procédé est connu sous le nom d’extrapolation de Richardson-Romberg.
Notons que sous les hypothèses précédentes,
P la régularité de u est suffisante pour obtenir un développement
asymptotique à l’ordre n, i.e. Ehd = ni=1 Ci hi + O(hn+1 ) à partir duquel il serait possible d’itérer le procédé
d’extrapolation et d’obtenir une vitesse d’ordre hn+1 en multipliant la complexité par n + 1. Mentionnons
toutefois que cette procédure peut accroı̂tre la variance de l’estimateur Θ̂hM . On considère en effet n + 1 fois
plus de variables aléatoires.
L’analyse de l’erreur faible s’effectuera toujours suivant le schéma présenté ci-dessus, avec les spécificités
associées à chaque problème, voir Sections 1.3.2 et 1.3.3. Le point clé est d’avoir de la régularité sur la fonc-
tion u et des contrôles sur ses dérivées, afin d’effectuer des développements d’Itô-Taylor. Sous les hypothèses

6
 précédentes, ces contrôles s’obtiennent par des techniques usuelles de calcul stochastique. On a cependant be-
soin pour cela de beaucoup de régularité, sur b, σ et g. Or, dans de nombreux cas cette régularité n’est pas
nécessaire. C’est en particulier vrai si l’on peut bénéficier de l’effet régularisant des noyaux de la chaleur. En
d’autres termes, il est possible de relaxer la régularité de la fonction g si par exemple XRT admet une densité
de transition régulière dont on contrôle bien les dérivées. En effet dans ce cas u(t, x) = Rd g(y)p(t, T, x, y)dy
où p(t, T, x, .) désigne la densité de XTt,x . Dans le cas de coefficients homogènes en temps nous noterons pour
simplifier l’écriture p(t, T, x, .) = p(T − t, x, .).
Néanmoins, prouver que XTt,x possède une telle densité nécéssite une hypothèse de non-dégénérescence de
l’équation (1.1). On considère deux jeux d’hypothèses qui garantissent que l’on obtient les propriétés souhaitées :
l’hypoellipticité à la Hörmander (hypothèse (H)) et l’uniforme ellipticité de la matrice de diffusion a (hypothèse
(UE)).
Hypothèse P (H) : Supposons que les coefficients b, σ sont C ∞ , bornés, homogènes en temps et réécrivons dXt =
r
A0 (Xt )dt + i=1 Ai (Xt ) ◦ dWti (i.e. les champs de vecteurs (Ai )i∈[[1,r]] correspondent aux colonnes de la matrice
de diffusion, A0 = b − 21 σ.∇σ et ◦dW désigne l’intégrale de Stratonovitch). Si l’espace vectoriel engendré par

A1 , · · · , Ar , [Ai , Aj ]0≤i,j≤r , [Ai , [Aj , Ak ]]0≤i,j,k≤r , · · ·

au point initial x est Rd (où [., .] désigne le crochet de Lie entre champs de vecteurs), alors pour tout t > 0,
Xt0,x admet une densité, voir par exemple le Chapitre 2.3 de Nualart [Nua95].
Si l’on suppose de plus qu’il existe L ∈ N∗ tel que que pour tout x ∈ Rd l’espace Rd est engendré en
ne considérant que les crochets de longueur inférieure ou égale à L, alors les travaux de Kusuoka et Stroock
[KS85], donnent les contrôles suivants sur la densité. Soient m, k des entiers et des multi-indices α, β tels que
2m + |α| ≤ k. ll existe des constantes C(T ), c, q, Q dépendant de L, m, k, α, a, b, telles que
 
C(T )(1 + |x|Q ) |x − y|2
|∂tm ∂xα p(t, x, y)| ≤ q exp −c . (1.9)
t (1 + |y − x|)k t(1 + |x|)2

Par ailleurs, l’inverse de la matrice de Malliavin que nous noterons Γt (x) est contrôlée dans tous les Lp . Pour
tout p ≥ 1 il existe C(T ), µ t.q.
1 + |x|µ
E[|Γt (x)|p ]1/p ≤ C(T ) dL . (1.10)
t
Le propos du travail de Bally et Talay [BT96a] est d’exploiter les contrôles (1.9), (1.10) pour étendre le
développement d’erreur (1.7) au cas d’une diffusion hypoelliptique sans régularité a priori sur g, que l’on peut
ici supposer mesurable et à croissance exponentielle.
Rappelons que lors des développements d’Itô-Taylor précédents apparaissaient pour tout i ∈ [[0, N − 1]], des
quantités du type Qhi := E[∂xα u(ti , Xthi )G(ti , Xthi )], où G fait intervenir les coefficients b, σ et leurs dérivées. Le
contrôle sur la densité (1.9) permet de contrôler les Qhi pour i ≤ N/2. Pour les i > N/2, l’idée naturelle est
d’utiliser une intégration par parties de Malliavin. Notons Qi = E[∂xα u(ti , Xti )G(ti , Xti )], i.e. on remplace Xthi
par Xti . L’intégration par parties pour Qi s’écrit Qi = E[u(ti , Xti )Hα (Xti , G(ti , Xti ))]. Les inégalités de Meyer
permettent alors d’obtenir qu’il existe p ≥ 1 tel que |Qi | ≤ CE[|Γti (x)|p ]1/p , où C dépend de normes des dérivées
de Malliavin de Xti .
On pourra faire l’intégration par parties dans Qhi sur les événements de forte probabilité pour lesquels
l’inverse de la matrice du schéma d’Euler Xthi est “proche” de l’inverse de la matrice de Malliavin Γti (x) de la
diffusion (1.1) et contrôler la probabilité des événements complémentaires, cf. [BT96a].
L’idée sous jacente globale est que l’hypothèse d’hypoellipticité permet de contrôler les dérivées enR utilisant en
temps petit les estimées ponctuelles sur la densité qui sont pertinentes dans la mesure où u(t, x) = R g(y)p(T −
t, x, y)dy, puis d’exploiter qu’en temps plus long c’est la matrice inverse (Γt )t≥T /2 qui est non-dégénérée. Quant
à la croissance de g, c’est le contrôle (1.9) qui rend cette condition suffisante.
Mentionnons que si la non-dégénérescence de l’inverse de la matrice de diffusion est bien contrôlée dans le
cas de la diffusion X, il n’en est pas de même pour son schéma d’Euler X h . C’est ce qui a motivé l’introduction
des Qi ci-dessus pour présenter l’intégration par parties. Il n’y a à ma connaissance pas de preuve d’existence
de la densité pour le schéma d’Euler hypoelliptique. C’est d’ailleurs ce dernier point qui a amené Bally et Talay
à considérer un schéma d’Euler légèrement perturbé pour obtenir un développement d’erreur sur les densités,
cf. [BT96b].
Hypothèse (UE) : On suppose la matrice de diffusion a uniformément elliptique et bornée, i.e. ∃Λ ≥ 1, ∀(t, x) ∈
[0, T ] × Rd , ∀ξ ∈ Rd , Λ−1 |ξ|2 ≤ ha(t, x)ξ, ξi ≤ Λ|ξ|2 . Pour b et σ Lipschitziens en espace et b borné, il est bien
connu que l’unique solution forte de (1.1) admet une densité vérifiant une estimée supérieure Gaussienne, i.e.

7
 2
C
∀T > 0, ∃(C, c) := (C, c)(T, |b|∞ , Λ), ∀(t, x, y) ∈ (0, T ] × Rd × Rd , p(t, x, y) ≤ td/2 exp(−c |x−y|
t ). Ce contrôle
s’obtient aisément à partir de la représentation en série parametrix de la densité, voir Friedman [Fri64] ou
McKean et Singer [MS67]. L’idée de ce type de méthode est d’écrire :
+∞
X
p(t, x, y) = p̃(t, x, y) + p̃ ⊗ H (i) (t, x, y). (1.11)
i=1

où p̃(t, x, y) est une densité Gaussienne, dans laquelle on a gelé les cofficients de (1.1) au point d’arrivée y et
où les termes p̃ ⊗ H (i) sont des convolutions temps-espace de densités Gaussiennes gelées et de noyaux H (i)
toujours plus régularisants en temps avec i. En d’autres termes, en temps petit, la contribution des p̃ ⊗ H (i) est
de plus en plus négligeable avec i. Nous reviendrons en détail sur la méthode parametrix au Chapitre 3.
L’approche développée dans [MS67] s’est en fait révelée très adaptée aux discrétisations. Considérons le
schéma d’Euler généralisé de (1.1) où l’on remplace les accroissements Browniens (Wti −Wti−1 )i∈[[1,N ]] dans (1.3)
par les réalisations i.i.d. (ξih )i∈[[1,N ]] d’une loi à densité plus générale ayant mêmes moments que l’accroissement
Brownien jusqu’à un certain ordre. Si ce dernier admet une densité pN , il est possible de la représenter aux
instants de discrétisation sous la forme
+∞
X
pN (t, x, y) = p̃N (t, x, y) + p̃N ⊗N H N,(i) (t, x, y). (1.12)
i=1

Dans l’équation ci-dessus p̃N designe la densité de la somme convenablement normalisée des ξih , ⊗N désigne une
convolution temps-espace discrète en temps et H N,(i) une version discrète des noyaux régularisants de (1.11).
Cette écriture et une analyse de la différence entre les représentations (1.11) et (1.12) a permis à Konakov et
Mammen [KM00] d’obtenir un théorème limite local à l’ordre h1/2 pour la différence (p − pN )(T, x, y). Cette
vitesse correspond à la vitesse d’approximation classique du théorème limite local en régime Gaussien, cf.
Bhattacharya et Rao [BR76]. C’est celle qui contrôle la convergence de p̃N vers p̃. Le résultat entre p et pN
découle ensuite d’une longue analyse de type stabilité. Nous présenterons au Chapitre 3 une extension de ce
résultat pour certains processus dégénérés, article [8].
Dans le cas particulier du schéma d’Euler (1.3), i.e. lorsque les (ξih )i∈[[1,N ]] sont déjà Gaussiens, on retrouve la
vitesse d’ordre h comme dans le cas de l’erreur faible mais sur les densités [KM02] et leur dérivées, cf. [Guy06] :

(p − pN )(1, x, y) = π(1, x, y)h + r(1, x, y)h2 ,

∃(c, C) := (c, C)(1, |b|∞ , Λ), ∀(x, y), (|π| + |r|)(1, x, y) ≤ C exp −c|x − y|2 . (1.13)

Sous des hypothèses plus fortes que celles courantes de régularité sur les coefficients, Gobet et Labart [GL08]
ont précisément quantifié l’ordre du premier terme pour T éventuellement petit. Ils ont établi |p − pN |(T, x, y) ≤
ChT −(d+1)/2 exp(−c|x−y|2 /T ), où les constantes C, c sont uniformes sur un compact en temps. C’est l’utilisation
du calcul de Malliavin qui impose une certaine régularité sur les coefficients (b, σ supposés Cb1,3 et ∂t σ ∈ Cb0,1 ).
L’ordre de la singularité en temps du terme de reste est néanmoins intuitivement clair : les développements
parametrix (1.11) et (1.12) font intervenir dans les noyaux H et H N des termes d’ordre deux en espace. Lorsque
l’on fait la différence de ses deux développements on a une sensibilité qui fait donc intervenir des termes d’ordre
trois. Les contrôles usuels sur les dérivées des noyaux de la chaleurs elliptiques, cf. Friedman [Fri64] donnent
que ces termes auront une singularité d’ordre 3/2 en temps. Celle-ci doit néanmoins être intégrée en temps
dans la convolution ce qui, après intégrations par parties, donne bien une singularité globale additionnelle en
T −1/2 correspondant au résultat de Gobet et Labart. Le même résultat peut se retrouver à partir des travaux
de Konakov et Mammen [KM09] qui précisent les résultats de [KM00] pour T éventuellement petit toujours
pour une chaı̂ne de Markov générale approchant (1.1).
L’avantage des développements de type (1.13) pour l’étude de l’erreur faible, par rapport aux approches
précédentes, est de pouvoir précisément relâcher les hypothèses sur la fonction g. Guyon montre en particulier
que le développement (1.7) reste valable si g est une distribution tempérée, incluant en particulier les masses
de Dirac et les dérivées de masse de Dirac. Ce dernier point est en particulier important pour les applications
au calcul de sensibilité, ou d’estimation des dérivées de l’EDP (1.6), qui sont cruciales en finance.
Enfin, l’approche parametrix de [KM00], [KM09] pour l’approximation permettrait d’étendre les résultats
mentionnés à des équations à coefficient de diffusion a uniformément Hölder en espace, borné, mesurable en
temps et dérive b bornée dans le cadre du problème de martingale classique. On retrouverait de la sorte les
hypothèses naturelles en analyse pour obtenir les contrôles sur les noyaux de la chaleur associés à des opérateurs
sous forme non divergence elliptiques, cf. Friedman [Fri64] ou Sheu [She91] pour une approche probabiliste de
ces résultats.

8
 En conclusion de cette présentation, notons que l’approche originelle de Talay et Tubaro ne requiert pas de
condition de non-dégénérescence. Elle est robuste dès lors que l’on a de la régularité. Pour se passer de régularité
sur g, une hypothèse de non dégénérescence (hypoellipticité ou uniforme ellipticité) est nécessaire. Cela permet
de faire porter les dérivations spatiales intervenant lors de l’analyse de l’erreur, sur la densité du processus que
l’on sait alors contrôler.
Mentionnons également la recherche active ces dernières années sur les schémas d’erreur faible élevée. Ces
derniers sont issus des travaux de Kusuoka [Kus01] dont le schéma a ensuite été développé par Lyons et Victoir
[LV04] et Ninomiya et al. [NV08], [NN09]. L’idée est toujours de procéder à des développements d’Itô-Taylor
dans les intégrales stochastiques et de remplacer ensuite les intégrales stochastiques itérées apparaissant par des
variables aléatoires ayant les mêmes moments jusqu’à un certain ordre en utilisant des formules de quadrature
pour ces intégrales Browniennes multiples. Le gros du travail est ensuite de déterminer des bonnes fonctions de
quadrature, ce qui requiert un travail algébrique conséquent.
Ces schémas imposent néanmoins beaucoup de régularité sur les coefficients b, σ ce qui n’est pas a priori
nécessaire pour avoir un développement (1.7) (cf. [KM02]) à partir duquel l’extrapolation de Richardson Rom-
berg est toujours possible.

1.3.2 Problème de Cauchy : cas discontinu

L’approche initiale de [TT90] est suffisamment robuste pour s’étendre naturellement au cas où le processus
Z dans (1.1) comprend une partie à saut. L’idée sous jacente est la même et consiste à contrôler les dérivées de
u qui peut dans ce cadre se voir comme la solution de l’EDP intégro différentielle :


 ∂ u(t, x) + hb(t, x), Dx u(t, x)i + 21 tr(a(t, x)Dx2 u(t, x))
 tZ

+ u(t, x + σ(t, x)z) − u(t, x) − h∇u(t, x), σ(t, x)ziI|z|≤1 ν(dz) = 0, (t, x) ∈ [0, T ) × Rd , (1.14)

 Rr
 u(T, x) = g(x), x ∈ R . d

Sous des hypothèses similaires à celles du cas continu, i.e. régularité des coefficients (supposés homogènes) et
de la condition terminale g, en imposant par ailleurs des conditions d’intégrabilité assez fortes sur la mesure
de Lévy ν, Protter et Talay [PT97] ont réussi à obtenir un développement de l’erreur faible similaire à (1.7)
lorsque (1.1) est approchée par son schéma d’Euler (1.3), cf. leur Théorème 2.2. L’idée reste celle de dérivation
du flot. Les hypothèses d’intégrabilité sur ν sont nécessaires pour contrôler en norme Lp , où p est relié à l’indice
de dérivation, les dérivées du flot.
Dans le cas d’EDS dirigées par des processus de Lévy généraux Z se pose également la question de simulation
du schéma d’Euler. En effet, bien que les accroissements de processus de Lévy soient explictement simulables
dans certains cas : processus stables, Gamma, Poisson composé, ..., ils ne le sont pas pas en toute généralité.
L’analyse de l’erreur faible associée à un schéma “simulable” dérivant du schéma d’Euler est effectuée dans
Jacod et al [JKMP05]. Les auteurs introduisent le schéma
h
X̄thi+1 = X̄thi + b(ti , X̄thi )h + σ(ti , X̄thi )ζi+1 , i ∈ [[0, N − 1]], X̄0h = x,

où les (ζi )i≥1 sont des variables i.i.d. approchant en loi les accroissements (Zti+1 − Zti )i≥1 . Toujours sous
des hypothèses fortes concernant l’intégrabilité de la mesure de Lévy ν, ils obtiennent une majoration d’er-
reur de la forme |E[g(X̄Th ) − g(XTh )]| ≤ C(γh ∨ h). Dans l’inégalité précédente γh est tel que pour H ∈ C04 (Rd ),
|E[H(ζh )−H(Zh )]| ≤ Chγh . Sous des hypothèses toujours plus fortes de régularité et d’intégrabilité sur b, σ, g, ν,
un développement d’erreur est également obtenu.

L’hypothèse d’intégrabilité de la mesure ν requise dans les travaux précédents, et cruciale du point de
vue des techniques utilisées, est en revanche génante pour de nombreuses applications qui font intervenir des
processus à queues lourdes. Le cas des processus stables d’indice α ∈ (0, 2), i.e. on exclut le cas Brownien,
illustre bien cet aspect. Ils interviennent dans de nombreux champs d’applications : physique mathématique
[IP06], télécommunications [SK74], finance [JMW96] et ne vérifient pas du tout les conditions d’intégrabilité
précédentes. Par ailleurs, le schéma d’Euler est dans ce cas parfaitement simulable, cf. [PT97] ou [ST94]. Dans
ce cadre, Hausenblas [Hau02] a obtenu dans le cas scalaire une majoration de l’erreur faible à l’ordre 1 en h
pour g ∈ L∞ en utilisant le calcul de Malliavin.
Au chapitre 3, nous développerons comment à partir d’une représentation en série parametrix pour la densité
de l’EDS (1.1) dirigée par un processus stable, ou en termes analytiques de la solution fondamentale de l’EDP
(1.14), on peut obtenir un développement d’erreur entre la densité de la diffusion et de son schéma d’Euler,
article [9]. En utilisant que les queues de la densité ont un comportement asymptotique en espace de l’ordre

9
de (1 + |x|d+α )−1 , qui est optimal au sens où une minoration du même ordre existe, nous pourrons également
relaxer les hypothèses de bornitude de g pour permettre des croissances en |g(x)| ≤ C(1 + |x|β ), β < α.
Mentionnons dans le cas général pour Z le travail de Kohatsu Higa et Tankov [KHT10] qui étudie l’erreur
faible pour des schémas de pas adaptés à l’amplitude des sauts avec une amélioration de l’erreur de convergence
usuelle en h du schéma d’Euler. L’approche reste néanmoins celle de dérivation du flot et requiert donc de la
régularité.

1.3.3 Autres fonctionnelles

Le problème de Cauchy Dirichlet : cas continu
Supposons que l’équation (1.1) soit dirigée par un mouvement Brownien. Considérons pour un domaine

D ⊂ Rd (i.e. un ouvert connexe) τ := inf{t ≥ 0 : Xt 6∈ D}. Si l’on considère la fonctionnelle Ψ (Xt )t∈[0,T ] =


g(T ∧τ, XT ∧τ ), la question de l’approximation de Ex [Ψ (Xt )t∈[0,T ] ] est liée à des problématiques de valorisation
d’options à barrières en finance et de façon plus générale à la résolution numérique de problèmes de Cauchy-
Dirichlet pour l’équation de la chaleur.
Dans la pratique, une question supplémentaire se pose par rapport à l’étude du problème sans condition
de bord, celle de
l’approximation de Ψh . Un choix naturel et peu coûteux numériquement est de prendre
Ψ (Xt )t∈[0,T ] = g(T ∧ τ , ΠD̄ (XTh∧τ h )), où τ h := inf{ti := ih, Xthi 6∈ D} et ΠD̄ désigne la projection sur
h h h

l’adhérence du domaine. En effet, ce choix permet de ne vérifier qu’aux instants de discrétisation si le schéma
d’Euler se trouve ou non dans le domaine. L’analyse de l’erreur de discrétisation

EhD = Ex [g(T ∧ τ h , ΠD̄ (XTh∧τ h ))] − Ex [g(T ∧ τ, XT ∧τ )] = E[u(T ∧ τ h , ΠD̄ (XTh∧τ h ))] − u(0, x),

s’effectue ensuite suivant le même schéma que dans le cas sans bord en utilisant les propriétés de régularité de
u solution de 
(∂t + Lt )u(t, x) = 0, (t, x) ∈ [0, T ) × D,
(1.15)
u(t, x) = g(t, x), (t, x) ∈ [0, T ) × ∂D ∪ {T } × D̄.
Des hypothèses de non dégénérescence similaires à celles présentées à la section précédente (hypoellipticité ou
uniforme ellipticité) permettent d’avoir des estimées si le domaine D est suffisamment régulier. Dans ce cadre,
Gobet [Gob00] a montré que EhD = O(h1/2 ) avec certaines conditions de support et/ou régularité sur g. Il est
à noter que ce résultat n’est pas dû à l’approximation du processus de diffusion mais est inhérent au choix de
l’approximation discrète du temps d’arrêt. En effet, dans le cas monodimensionnel, pour Xt = Wt et qu’il n’y a
donc pas d’erreur de discrétisation du processus, on peut déduire des travaux de Siegmund et Yuh [SY82] que
pour l’approximation discrète de la fonction de répartition du processus stoppé, i.e. g(t, x) = It=T Ix>A , A fixé,
EhD = Ch1/2 + o(h1/2 ). Ce résultat a motivé les travaux [1] et [7] qui ont permis d’établir respectivement une
h
borne inférieure et un développement d’erreur à l’ordre 1/2 pour ED .
Par ailleurs, cette mauvaise vitesse en comparaison du cas sans bord a donné lieu à une recherche active
pour l’amélioration de la convergence des méthodes numériques : approximation plus fine du temps de sortie
par ponts Browniens [Gob00, Gob01], techniques de grandes déviations pour les asymptotiques de probabilités
de sortie [Bal95, BCI99], correction de domaine [BGK97]. Détaillons ici cette dernière. Dans le cas du modèle
de Black et Scholes monodimensionnel et pour un demi-espace, Broadie et al. [BGK97], ont introduit l’idée de
baisser convenablement en fonction de h la barrière de sorte à compenser la surestimation des trajectoires prises
en compte du fait de l’estimation discrète. Cette méthode a été généralisée à une large classe de domaines non
cylindriques réguliers dans [7] pour des processus de diffusions.
Nous reviendrons en détail sur ces aspects au Chapitre 2 dédié aux approximations de processus stoppés.

Les modèles cinétiques et options asiatiques : cas continu

RT
Nous considérons ici le cas où ψ (Xt )t∈[0,T ] = g(T −1 0 Xt dt, XT ). Typiquement, si l’on considère le
RT RT
cas d’une option asiatique en finance g(T −1 0 Xt dt, XT ) = (T −1 0 Xt dt − XT )+ pour un call. Le facteur
de normalisation en T −1 devant l’intégrale permet aux deux quantités d’être à la même échelle globale et
la comparaison entre valeur d’un actif à l’instant T et sa moyenne sur la période. Pour l’analyse de l’erreur
Z T Z T
faible EhD := Ex [g(T −1 h
Xφ(t) dt, XTh )] − Ex [g(T −1 Xt dt, XT )], le découpage précédent peut s’effectuer en
0 0
considérant l’EDP : 
(∂t + Lt )u(t, x) = 0, (t, x) ∈ [0, T ) × R2d ,
(1.16)
u(T, x) = g(x), x ∈ R2d ,

10


où l’on note x = (x, y) et Lt ϕ(x) = hb(t, x), ∇x ϕ(x)i + 12 tr a(t, x)Dx2 ϕ(x) + hx, ∂y ϕ(x)i. Cette équation
Ry
est dégénérée et le générateur Lt est celui du couple (Xt , Yt )t≥0 , Yt = y + 0 Xs ds. Lorsque les coefficients b, σ
dépendent également de y alors l’équation (1.16) peut être associée à un système Hamiltonien stochastique. Pour
2 2
un Hamiltonien donné de la forme H(x, y) = V (y) + |x|2 où V est un potentiel et |x|2 l’énergie cinétique d’une
particule de masse unité, le système Hamiltonien stochastique associé correspondrait à b(Xs , Ys ) = −(∂y V (Ys )+
F (Xs , Ys )Xs ), où F ≥ 0 est un terme de friction.
Dans le cas d’un mouvement Brownien géométrique à coefficient de diffusion aléatoire, l’analyse de EhD
a été effectuée par Temam [Tem01]. L’approximation par chaı̂ne de Markov de la solution fondamentale de
(1.16), avec b, σ dépendant de x = (x, y), est l’objet du travail [8] et sera présentée au Chapitre 3. A l’aide de
développements parametrix nous obtiendrons un théorème limite local en régime Gaussien pour la différence
entre les densités du processus et de la chaı̂ne de Markov. L’étude de systèmes dégénérés de type (1.16) sera
développée au Chapitre 5.2. Nous y présenterons les résultats de [10] qui donnent en particulier des estimées
“à la Aronson” pour la solution fondamentale de (1.16).

1.4 Quelques autres problèmes usuels d’approximation de processus

Erreur en loi normalisée
Lorsque l’on s’intéresse à l’approximation de fonctionnelles de la trajectoire, une autre question naturelle est
celle de la limite en loi (au sens des processus) de l’erreur normalisée, i.e. la limite de U h := ρh (X − X h ) pour
une renormalisation ρh → +∞ à préciser. La dynamique de cette limite a été étudiée initialement par Kurtz
h→0
et Protter [KP91] pour des processus de diffusion puis dans un cadre plus général par Jacod et Protter [JP98].
Ce travail inclut en particulier le cas déterministe et de certains processus de Lévy.
Pour les processus de diffusion, la bonne renormalisation est ρh = h−1/2 , ce qui pouvait se prévoir à partir
des contrôles sur l’erreur forte. Si b = 0, Z = W dans (1.1), la dynamique limite de l’erreur s’écrit
Z t Z t
1
Ut = σ ′ (Xs )Us dZs − √ (σσ ′ )(Xs )dW̃s , (1.17)
0 2 0

où W̃ est un mouvement Brownien indépendant de W . Cette dynamique peut intuitivement se déduire du
développement au premier ordre de la dynamique de la différence :
Z t Z t
′
h
Ut = h
σ (Xs )Us dWs + σσ ′ (Xφ(s)
h
)h−1/2 (Ws − Wφ(s) )dWs + Rth ,
0 0

où Rth est un terme de reste.

Le mouvement Brownien indépendant du deuxième terme dans (1.17) apparaı̂t grâce à l’indépendance √ asymp-
totique des deux intégrales stochastiques dans la dynamique de U h . Le terme de normalisation en 1/ 2 provient
quant à lui de la caractérisation à l’aide de la variation quadratique.
Le cas des processus de Lévy est plus délicat. Mentionnons le travail de Jacod [Jac04] qui relie le paramètre
de normalisation ρh à la concentration de la mesure de Lévy ν de Z autour de 0 (cf. (1.2)). La dynamique limite
de U est similaire à celle de (1.17) mais l’intégration pour le deuxième terme ne fait plus intervenir un Brownien
mais un autre processus de Lévy V , lié au processus directeur Z mais toujours indépendant de ce dernier. En
particulier, le cas stable est entièrement caractérisé.
Le problème de caractérisation de la dynamique limite de l’erreur a été abordé par Rubenthaler [Rub03] dans
le cas d’une version simulable du schéma de discrétisation de l’EDS dirigée par un procéssus de Lévy général.

Mesures invariantes
Lorsque le processus X vérifie certaines propriétés de stabilité de type Lyapunov, la question de l’approxi-
mation de mesures invariantes se pose de façon naturelle.
Des schémas à pas décroissants ont été introduits dans le cas continu par Lamberton et Pagès [LP02, LP03].
Le raffinement du pas de temps permet de converger vers la mesure invariante du processus initial. Cette
approche présente par ailleurs l’avantage de pouvoir s’appliquer au cas de processus dégénérés. Elle a été
généralisée par Panloup pour des d’EDS dirigées par des processus de Lévy [Pan08] et dans des cas continus à
dérive singulière par Lemaire [Lem07].
Dans le cadre des systèmes Hamiltoniens, i.e. lorsque le terme de friction F > 0 dans l’exemple ci-dessus,
Talay [Tal02], Mattingly et al. [MSH02] ont étudié différents schémas d’approximation de la mesure invariante

11
par chaı̂ne de Markov à pas constants. Ces derniers présupposent une meilleure connaissance a priori du semi-
groupe de la diffusion. En effet dans ces travaux c’est la mesure invariante de la chaı̂ne qui est approchée et il
est ensuite nécessaire d’avoir des contrôles précis sur l’écart avec la mesure invariante du processus initial.

Transferts de propriétés du modèle continu au modèle discret

Nous mentionnerons ici les travaux de Veretennikov et ses coauteurs qui se sont intéressés au transfert de
propriétés de mixing [KV06] et de principes de grandes déviations du modèle continu de diffusion dirigée par
un Brownien à un schéma numérique discret [Ver03].

Vers des modèles non linéaires

L’essor de la théorie des équations différentielles stochastiques rétrogrades (EDSR) initiée par Pardoux et
Peng [PP90], a permis d’étendre les représentations de Feynman Kac à une certaine classe de problèmes non
linéaires : équations semi/quasi-linéaires, problèmes de frontière libre.
Ce lien a donné lieu à une série de travaux sur l’approximation de telles équations. Citons Zhang [Zha04]
pour l’analyse de l’erreur d’un schéma de semi-discrétisation d’une EDSR, puis Gobet et Labart [GL07] pour
un développement de l’erreur. Il ressort de ces travaux que l’erreur globale en norme L2 commise sur le triplet
(X, Y, Z) est guidée par l’erreur forte associée à la composante progressive. Si la composante rétrograde est
un peu régulière, on propage simplement en temps rétrograde l’erreur usuelle forte associée à l’approximation
de la composante progressive. On retrouve ce même phénomène pour les équations couplées, cf. [2] et [5], les
approximations d’EDSR dans un domaine [6], ou les problèmes de frontière libre, voir Bally, Pagès et al. [BP03].
Nous reviendrons en détail sur ces aspects au Chapitre 4.

12
Chapitre 2

Approximation de processus
tués/stoppés

Ce chapitre vise à faire une synthèse des travaux [1], [3], [4], [7] écrits en partie avec Emmanuel Gobet de
l’ENSIMAG et relatifs à différents problèmes d’approximation pour des processus tués et/ou stoppés. L’idée
est tout d’abord d’établir une vitesse de convergence pour l’erreur faible lorsque l’on fait le choix de stopper le
schéma d’approximation à temps discret. Par opposition au cas de l’espace tout entier, cette vitesse est d’ordre
1/2 en h si h désigne le pas de temps uniforme de discrétisation. Cet ordre correpond à celui de l’overshoot du
schéma d’Euler, i.e. la distance au bord du domaine lorsque le schéma sort à temps discret. Cet ordre est encore
celui de l’accroissement Brownien.
Le travail initial [1] a permis d’isoler l’overshoot comme terme dominant dans l’erreur et d’obtenir un
encadrement de l’erreur à l’ordre 1/2 pour des diffusions hypoelliptiques et des domaines cylindriques réguliers
vérifiant une condition de frontière non caractéristique.
Des théorèmes limites sur l’overshoot existent depuis longtemps dans le cadre du Brownien monodimension-
nel tué à temps discret. En particulier, Siegmund et Yuh [SY82] ont établi l’indépendance asymptotique du
temps de sortie discret et de l’overshoot renormalisé. C’est précisément ce type de résultat qui est à la base
du développement de l’erreur faible. Une grande partie de [7] a consisté à caractériser la loi limite du triplet
temps de sortie discret, position de sortie et overshoot renormalisé pour le schéma d’Euler d’une diffusion uni-
formément elliptique stoppée à la sortie d’un domaine régulier non cylindrique. Il s’agit là de l’un des résultats
principaux de ce chapitre. Dans le cas particulier de cônes formés par l’intersection de demi-espaces un résultat
similaire avait été obtenu pour le Brownien dans [3]. Ce théorème limite a dans les deux cas mentionnés permis
d’aboutir au développement de l’erreur faible (avec une extension au problème stationnaire pour un domaine
régulier).
Comme nous l’avons indiqué au Chapitre 1.3.1, un développement d’erreur permet d’utiliser des méthodes
d’accélération de la convergence pour l’approximation numérique. En particulier, la méthode d’extrapolation de
Romberg est toujours possible, puisque purement algébrique, sans connaissance a priori de la constante dans le
développement limité. Dans le cas de l’approximation faible de processus tués, une caractérisation “explicite” de
cette constante permet d’utiliser une technique alternative d’accélération de convergence. Celle-ci consiste à tuer
le processus à temps discret à la sortie d’un domaine plus petit que celui initial que l’on peut choisir de sorte à
tuer le terme dominant de l’erreur. La construction de ce domaine dépend du pas de temps h, et explicitement de
l’expression de la constante dans le développement limité de l’erreur. L’idée est d’utiliser la notion de sensibilité
du problème de Dirichlet par rapport au domaine, qui avait d’abord été introduite en optimisation de formes,
cf. Sokoloski et Zolesio [SZ92], puis développée dans un cadre probabiliste par Costantini et al. [CKG06].
Nous conclurons ce chapitre en présentant les résulats de [4] où l’on montre comment étendre l’analyse de
l’erreur de discrétisation au cas non Markovien des processus d’Itô pour lesquels il n’est plus possible d’utiliser
une représentation de l’erreur via une EDP.

2.1 Modèle général et lien avec le problème de Dirichlet

Soit (Xt )t≥0 un processus de diffusion d-dimensionnel de dynamique
Z t Z t
Xt = x + b(s, Xs )ds + σ(s, Xs )dWs (2.1)
0 0

13
où W est un mouvement Brownien de Rd , et les coefficients b, σ sont Lipschitz en espace et localement bornés
en temps.
On se donne une famille (Dt )t≥0 , dépendant du temps de domaines réguliers bornés de Rd . On suppose
également cette famille régulière en temps. Pour un instant T > 0 déterministe fixé, on associe à la famille
précédente le domaine temps-espace (voir Figure 2.1)
[
D= {t} × Dt = {(t, x) : 0 < t < T, x ∈ Dt } ⊂]0, T [×Rd.
0<t<T

Les domaines cylindriques sont des cas particuliers des domaines temps-espace précédents. On a dans ce cas

Dt
DT
D0

temps
0 t T
Rd

Figure 2.1 – Domaine temps espace et ses sections en temps.

D =]0, T [×D où D est un domaine de Rd (i.e. Dt = D pour tout t).

Définissons τ = inf{t > 0 : Xt 6∈ Dt }. Ainsi, τ ∧ T est le premier temps de sortie de (t, Xt )t≥0 du do-
maine temps-espace D. Pour des fonctions continues g, f, k : D̄ → R, nous présenterons des résultats relatifs à
l’approximation numérique de quantités de la forme
Z τ ∧T Z s
Q(T, g, f, k, x) := Ex [g(τ ∧ T, Xτ ∧T )Zτ ∧T + Zs f (s, Xs )ds], Zs = exp(− k(r, Xr )dr). (2.2)
0 0

Sous de bonnes hypothèses de régularité des fonctions précédentes et de non-dégénerescence du coefficient de

diffusion dans un voisinage du bord de D, l’expression ci-dessus n’est autre que la représentation Feynman-Kac
de la solution de l’EDP (
(∂t u + Lt u − ku + f )(t, x) = 0, (t, x) ∈ D,
(2.3)
u(t, x) = g(t, x), (t, x) ∈ PD,
où Lt ϕ(x) = hb(t, x), ∇ϕ(x)i + 12 tr(a(t, x)Dx2 ϕ(x)) est le générateur infinitésimal de Xt , PD := ∂D\[{0} × D0 ]
désigne la frontière parabolique du domaine D et les fonctions k et f correspondent respectivement à des termes
potentiels et source dans l’équation de la chaleur ci-dessus.
L’approximation numérique la plus simple de la quantité (2.2) consiste à utiliser le schéma d’Euler stoppé
à temps discret. Pour un pas de temps h = T /m > 0, m ∈ N∗ fixé, on définit les instants de discrétisation
ti = ih = iT /m)i≥0 et pour t ≥ 0, φ(t) := ti si ti ≤ t < ti+1 . Le schéma d’Euler s’écrit alors :
Z t Z t
h h h
Xt = x + b(φ(s), Xφ(s) )ds + σ(φ(s), Xφ(s) )dWs . (2.4)
0 0

Nous associons au schéma d’Euler son temps de sortie discret τ h := inf{ti > 0 : Xthi ∈
/ Dti }, ce qui conduit à
approcher (2.2) par :
Z τ h ∧T Rt h
Qh (T, g, f, k, x) := Ex [g(τ h ∧ T, Xτhh ∧T )Zτhh ∧T + h
Zφ(s) h
f (φ(s), Xφ(s) )ds], Zth = e− 0
k(φ(r),Xφ(r) )dr
. (2.5)
0

L’erreur faible associée à cette procédure de discrétisation est la différence des équations (2.2) et (2.5), EhD :=
(Qh − Q)(T, g, f, k, x). Précisons que dans l’approximation (2.5) ci-dessus la fonction g est presque sûrement
évaluée hors de l’adhérence du domaine. Cette approximation peut à première vue sembler grossière. Elle
n’affecte pas néanmoins la vitesse de convergence et évite numériquement l’étape potentiellement coûteuse de
projection sur D̄.
Une première hypothèse de travail est celle de frontière non caractéristique du domaine. Précisément, si
F (t, x) désigne la distance signée de x ∈ Rd au bord du domaine Dt 1 , on suppose que
1. Dans le cas d’un domaine cylindrique F (t, x) = F (x), distance signée au bord du domaine D.

14
 S
(FNC) Pour un r0 > 0 il existe a0 > 0 tel que ∇F (t, x)[σσ ∗ ](t, x)∇F (t, x)∗ ≥ a0 pour tout (t, x) ∈ 0≤t≤T {t} ×
V∂Dt (r0 ), V∂Dt (r0 ) := {x ∈ Dt : d(x, ∂Dt ) ≤ r0 }.
Cette condition est suffisante pour garantir la convergence de EhD vers 0, cf. Proposition 3 de [1]. Elle apparaı̂t
comme minimale dans la mesure où des contre-exemples existent lorsqu’elle n’est pas vérifiée, cf. [Gob00]. Elle
exprime que dans un voisinage du bord le processus de distance au bord est non dégénéré, ou dit autrement
que la matrice de diffusion est non dégénérée le long des normales au bord des sections Dt du domaine D.
Comme le processus de distance au bord est mono-dimensionnel, il est en quelque sorte similaire à un
mouvement Brownien changé de temps. C’est la formalisation de cette observation heuristique qui a permis de
généraliser au schéma d’Euler les résultats asymptotiques sur l’overshoot à partir des résultats Browniens, voir
Section 2.2.1.
Comme nous l’avons indiqué en Section 1.3.1, un outil essentiel pour l’analyse de l’erreur faible est le contrôle
des dérivées de la solution u de (2.3). Nous avons considéré deux types d’hypothèses :
(A) Le domaine est cylindrique, i.e. D :=]0, T [×D. Les coefficients b, σ sont homogènes en temps et de classe Cb∞
(bornés à dérivées bornées). Le coefficient de diffusion a vérifie (FNC) et une hypothèse de type Hörmander
forte, i.e. seuls les crochets associés à la partie diffusive engendrent l’espace.
Dans ce cadre le noyau de la chaleur tué possède de bonnes propriétés, cf. Cattiaux [Cat91]. Ainsi, lorsque
g|PD\[{T }×DT ] = 0 et que le support de g(T, .) est à une distance positive du bord, il est possible de contrôler
les dérivées de u en adaptant au cas tué les intégrations par parties de Malliavin utilisées par Bally et Talay
[BT96a] (Article [1]).
(B) Le domaine D, les coefficients b, σ, f, g sont suffisamment réguliers et a est uniformément elliptique.
Dans ce cas il est possible d’utiliser des estimées de type Schauder pour les équations dans des domaines non
cylindriques, cf. Lieberman [Lie96] (Article [7]).

Précisons que sous les hypothèses (A) nous avons étudié l’erreur faible sans terme potentiel ni source, i.e.
f = k = 0. En revanche, sous (B) nous avons considéré le cas général de l’équation (2.3). Dans les deux cas
mentionnés, les formules d’Itô-Taylor nécessaires pour contrôler EhD font néanmoins intervenir un développement
en semi-martingale de ΠD̄ (Xth ). En particulier, le contrôle des restes, voir équation (2.6), requiert des estimées
assez fines sur la mesure d’occupation d’un voisinage du bord pour le schéma et sur son temps local au bord, i.e.
temps local de F (t, Xth )t≥0 en 0. Nous renvoyons à l’annexe A de [7] et à la Section 3 de [1] pour des détails.
Sous les hypothèses (A), pour lesquelles g(s, x) = Is=T g0 (x), (s, x) ∈ PD, nous avons établi dans [1] que :
"Z #
T
EhD := Ex [g0 (XTh )Iτ h >T ] − E[g0 (XT )Iτ >T ] = Ex ∂n u(t, Xth )d[F − (Xt∧τ
h
h )] + O(h). (2.6)
0

Ce contrôle illustre l’importance de l’overshoot pour le développement d’erreur. Pour g0 positive, nous avons
déduit de (2.6) dans [1] un encadrement de l’erreur pour h suffisamment petit, i.e. ∃C ≥ 1, C −1 h1/2 ≤ EhD ≤
Ch1/2 , en utilisant la stricte positivité de la dérivée normale. C’est la caractérisation de la loi limite de l’overshoot
qui manquait pour obtenir le développement. C’est précisément cette loi que nous avons pu établir dans [7]
pour des domaines temps-espace vérifiant la condition de frontière non caractéristique (FNC).

2.2 Contrôles sur l’erreur faible

2.2.1 Théorème limite sur l’overshoot
Le résultat principal de cette section est le théorème limite sur la loi jointe du triplet formé par le temps de
sortie discret, la position de sortie et l’overshoot renormalisé du schéma de discrétisation :

(τ h , Xτhh , h−1/2 F − (τ h , Xτhh )).

Les lois limites d’overshoot ont été étudiées dans de nombreux cadres. Mentionnons les travaux de Siegmund
[Sie79] et Siegmund et Hu [SY82] dans le cas de la marche simple scalaire, ceux d’Alsmeyer [Als94] ou Fuh
et Lai [FL01] pour les chaı̂nes de Markov ergodiques, et enfin Doney et Kyprianou [DK06] pour les processus
de Lévy. Le point commun à toutes ces approches est l’utilisation d’arguments de renouvellement. Une théorie
du renouvellement non linéaire a même été développée afin d’étudier les overshoot de marches aléatoires par
rapport à des frontières courbes, cf. [Woo82] et [Zha88]. Nous étendons ce type de résultats au schéma d’Euler.

15
Théorème 2.2.1 (Lois jointes limites associées à l’overshoot) Supposons que (FNC) est vérifiée, que le
domaine D est de classe C 2 et les coefficients b, σ ∈ C (1+θ)/2,1+θ (D̄) où θ ∈]0, 1[ avec les notations usuelles pour
les espaces de Hölder 2 . Soit ϕ une fonction continue à support compact. Pour tout t ∈ [0, T ], x ∈ D0 , y ≥ 0,


Ex [Iτ h ≤t Zτhh ϕ(Xτhh )IF − (τ h ,X h )≥y√h ] −→ Ex Iτ ≤t Zτ ϕ(Xτ ) 1 − H(y/|∇F σ(τ, Xτ )|) ,
τh h→0
Ry Pn
où H(y) := (E0 [sτ + ])−1 0 P0 [sτ + > z]dz et s0 := 0, ∀n ≥ 1, sn := i=1 Gi , les Gi étant des lois normales
centrées standards i.i.d., τ + := inf{n ≥ 0 : sn > 0}.
En d’autres termes, le triplet (τ h , Xτhh , h−1/2 F − (τ h , Xτhh )) converge donc faiblement vers (τ, Xτ , |∇F σ(τ, Xτ )|Y )
où Y est une variable aléatoire indépendante de (τ, Xτ ) et de fonction de répartition H. En fait, Y est la loi
asymptotique de l’overshoot Brownien, voir [Sie79]. La moyenne de l’overshoot est une quantité importante
E0 [s2τ + ]
qui interviendra dans la procédure de correction de domaine. Notons que E(Y ) = 2E0 [sτ + ] := c0 , et d’après
Siegmund [Sie79]
ζ(1/2)
c0 = − √ = 0.5826... (2.7)
2π
où ζ désigne la fonction de Riemann.
L’idée clé pour prouver le Théorème 2.2.1 est de pratiquer une approximation de la frontière du domaine
par un demi-espace et du schéma d’Euler par un processus à coefficients gelés dans un voisinage du bord. Une
double localisation est toutefois nécessaire, voir Figure 2.2. En effet, l’approximation à domaine et coefficients
constants sera pertinente lorsque le schéma d’Euler se trouve dans un voisinage de taille hα du bord, α > 0. Mais
pour pouvoir appliquer les arguments de renouvellement, il faut également que le schéma d’Euler se trouve hors
d’un voisinage de taille h1/2−ε < hα . Par rapport à l’échelle Brownienne typique de h1/2 , le schéma se trouve
suffisamment loin du bord. Le théorème résulte ensuite de contrôles assez précis sur l’approximation évoquée,
Rd
∂D 
1
}∆ 2 −ε
∆α

temps
0 t T
1
Figure 2.2 – Les deux voisinages de localisation avec α < 2 − ε.

avec en particulier des estimées sur les échelles typiques pour les temps de sortie en fonction de l’échelle du
voisinage, voir Section 3 de [7].

2.2.2 Application au développement d’erreur

Le Théorème 2.2.1 était l’argument manquant pour obtenir le développement d’erreur dans [1]. Nous avons
utilisé ce résultat dans [7] pour établir le :
Théorème 2.2.2 (Développement au premier ordre) Supposons que les hypothèses (B) soient vérifiées.
Alors, pour h suffisamment petit
√ √
EhD = c0 hEx [Iτ ≤T Zτ (∇u − ∇g)(τ, Xτ ) · ∇F (τ, Xτ )|∇F σ(τ, Xτ )|] + o( h), (2.8)

où c0 est défini en (2.7).

Nous avons indiqué dans (B) qu’une régularité suffisante était requise. La régularité minimale considérée dans
[7] présuppose que le domaine est C 2 , les fonctions g, k, f sont C (1+θ)/2,1+θ (D̄) de même que les coefficients
b, σ.
Dans le cas particulier du jeu d’hypothèses (A), abordé dans [1], i.e. cas hypoelliptique tué sans second
membre ni potentiel, on déduirait sans trop de difficultés le développement à partir du Théorème 2.2.1 et de
l’équation (2.6), après une intégration par parties pour faire explicitement apparaı̂tre l’overshoot F − (XTh∧τ h )
et non plus ses variations.
2. On renvoie pour cela au livre de Friedman [Fri64] ou Lieberman [Lie96]

16
2.2.3 Correction de domaine
Le Théorème 2.2.2 caractérise précisément la “mauvaise” convergence de l’approximation discrète Qh . De
nombreuses techniques avaient été utilisées numériquement pour y remédier. Mentionnons les approches par
grandes déviations de Baldi et al., [Bal95, BCI99], ou celles par ponts Browniens, pour affiner l’estimation du
temps de sortie entre deux instants de discrétisation, cf. [Gob00, Gob01]. Ces dernières sont néanmoins coûteuses
pour des géométries complexes dans la mesure où il faut procéder à des approximations par plan tangent.
D’autres méthodes numériques ont été développées. Indiquons celles par marches aléatoires sur des sphères,
cf. Milstein [Mil98], ou sur des rectangles, cf. Milstein et Tretyakov [MT99] et Deaconu et Lejay [DL06], qui
sont d’autant plus efficaces que la géométrie du domaine est proche de leur structure naturelle. Nous pensons
néanmoins que notre approche reste une des plus directes.
Par ailleurs, le développement d’erreur du Théorème 2.2.2 ouvre de nouvelles perspectives. Il permet déjà
d’effectuer l’extrapolation de Romberg. Néanmoins, la connaissance explicite de la constante c0 ainsi que la forme
du développement, qui fait intervenir les dérivées normales de la solution de l’EDP (2.3) et le gradient sur la
frontière parabolique de la condition de bord, nous a amenés à proposer une technique alternative d’accélération
de convergence. Celle-ci a par ailleurs un moindre impact sur la variance que l’extrapolation de Romberg.
Il s’agit simplement de stopper à temps discret le processus lorsque ce dernier sort d’un domaine plus petit
que celui initial que l’on déforme suivant une amplitude et une direction qui permettent de supprimer le premier
terme du développement dans (2.8). Précisément, introduisons Dh ⊂ D dont les sections en temps sont définies
par Dth = {x ∈ Dt : d(x, ∂Dt ) > c0 h1/2 |∇F σ(t, x)|}, voir Figure 2.3. Notons τ̂ h = inf{ti > 0 : Xthi ∈ Dthi } le

∂Dt
1111
0000 Dth
0000
1111
0000
1111
0000
1111
0000
1111
0000
1111

Figure 2.3 – Frontière ∂Dt et domaine restreint Dth .

temps de sortie discret de Dh . Cette procédure de correction de frontière avait initialement été introduite dans
le cas mono-dimensionnel Brownien par Broadie et al. [BGK99]. Elle amène à simuler la quantité
Z τ̂ h
g(Xτ̂hh )Zτ̂hh + h
Zφ(s) h
f (Xφ(s) )ds. (2.9)
0

R τ̂ h
Notons Q̂h (T, g, f, k, x) := E[g(Xτ̂hh )Zτ̂hh + 0
h
Zφ(s) h
f (Xφ(s) )ds]. Nous avons alors le résultat suivant :

Théorème 2.2.3 (Correction de frontière) Sous les hypothèses du Théorème 2.2.2 et si ∇F (.)|∇F σ(.)| est
de classe C 1,2 , alors pour h suffisamment petit :
√
ÊhD := (Q̂h − Q)(T, g, f, k, x) = o( h).

L’hypothèse supplémentaire sur ∇F (.)|∇F σ(.)| permet d’assurer que le domaine modifié reste de classe C 2 .
L’idée de la preuve repose sur le développement :
Z h Z !
τ̂ τDh ∧T
ÊhD := E[g(Xτ̂hh )Zτ̂hh + h
Zφ(s) h
f (Xφ(s) )ds] − Ex [g(τDh ∧ T, XτDh ∧T )ZτDh ∧T + Zs f (s, Xs )ds]
0 0
Z Z !
τDh ∧T τ ∧T
+ Ex [g(τDh ∧ T, XτDh ∧T )ZτDh ∧T + Zs f (s, Xs )ds] − Ex [g(τ ∧ T, Xτ ∧T )Zτ ∧T + Zs f (s, Xs )ds]
0 0

:= ÊhD,1 + ÊhD,2 ,

où τDh := inf{t > 0 : Xt 6∈ Dth }. Le terme ÊhD,1 est de même nature que EhD . En particulier, les estimées de
type Schauder sur l’EDP associée au domaine Dh sont uniformes pour h petit et la solution de l’EDP dans le

17
domaine perturbé et son gradient vont converger respectivement vers u et ∇u. Le développement du Théorème
2.2.2 va donc rester valable pour ÊhD,1 .
Le terme ÊhD,2 peut être relié à la sensibilité du problème de Dirichlet, cf. Theorème 2.2 dans [CKG06]. Assez
intuitivement, cette sensibilité fait intervenir les gradients de u et g sur la frontière parabolique, l’intensité de
la pertubation, −c0 h1/2 |∇F σ| dans la mesure où l’on rétrécit le domaine, et la direction normale ∇F qui fait
rentrer “le plus vite” dans le domaine. On obtient :
ÊhD,2 = −c0 h1/2 Ex [Iτ ≤T Zτ (∇u − ∇g)(Xτ ) · ∇F (Xτ )|∇F σ(Xτ )|] + o(h1/2 ),
ce qui va exactement compenser le terme du développement du Théorème 2.2.2 et donner l’estimée du Théorème
2.2.3.
Insistons sur le fait qu’il était crucial pour pouvoir effectuer cette correction de connaı̂tre explicitement la
constante c0 dans le développement (2.8). C’est la forme de ce développement qui permet le parallèle avec
la sensibilité du problème de Dirichlet. La connaissance précise de la constante est nécessaire pour prendre la
bonne intensité de perturbation dans le schéma corrigé. Cette technique reste donc très liée au type de problème
considéré.
Nous conclurons en indiquant que le développement d’erreur et la correction associée ont pu être démontrés
sous des hypothèses de type (B), i.e. domaine et coefficients réguliers, uniforme ellipticité, pour un problème
stationnaire, voir Section 4 de [7] et les résultats numériques ci-dessous (Section 2.3). Un autre cadre est celui du
mouvement Brownien tué à la sortie d’un cône, cf. [3]. Ce dernier cas est important en finance où les domaines de
validité des options à barrières sont souvent définis en terme d’espaces produits. Un théorème limite similaire au
Théorème 2.2.1 a été obtenu pour une intersection de demi-espaces. Le développement et la correction n’ont en
revanche été établis en toute généralité que pour des cônes bidimensionnels lorsque g est nulle sur un voisinage
de la singularité du cône. La difficulté est ici de contrôler le gradient de la solution de l’EDP jusqu’au bord,
y compris au voisinage de la singularité. Pour ce dernier point, nous obtenons la régularité gâce à l’hypothèse
sur le support de g et à partir d’un contrôle radial sur les dérivées de la densité du processus tué. Celle-ci se
décompose en une série dont chaque terme fait intervenir une partie sphérique et radiale : la partie sphérique
provient de la décomposition spectrale du Laplacien pour le problème de Dirichlet associé à la trace du cône sur
la sphère, la partie radiale fait intervenir des fonctions de Bessel d’indices associés à la dimension et aux valeurs
propres du Laplacien. Dans le cas bidimensionnel, l’expression explicite de cette densité sous forme de série est
bien connue et s’obtient comme extension de la méthode des images, cf. Carslaw et Jaeger [CJ59]. Dans le cas
d ≥ 2, la formulation générale a été obtenue par Bañuelos et Smits [BS97].

2.3 Résultats numériques

Pour des exemples numériques dans le cas du Brownien à la sortie d’un cône nous renvoyons à la Section 3
de [3]. Nous choisissons de présenter ici un exemple pour un processus tridimensionnel stoppé à la sortie d’une
boule, problème elliptique.
Soit (Xt )t≥0 de dynamique :
∗
dXt = b(Xt )dt + σ(Xt )dWt , ∀x ∈ R3 , b(x) = (x2 x3 x1 ) ,
 
(1 + |x3 |)1/2 0 0
 1

3 1/2 
σ(x) =  2 (1 + |x1 |)1/2 4 (1 + |x1 |)1/2 0 ,
1 1/2


3 1/2 1/2
0 2 (1 + |x2 |) 4 (1 + |x2 |)
(2.10)
et X0 appartenant à D = B(0, 2). On se donne une fonction u(x) = x1 x2 x3 définie sur D̄. On obtient l’EDP
de type (2.3), en version stationnaire, associée à (2.10) satisfaite par u en prenant g = u|∂D , f = −Lu, où L
désigne le générateur infinitésimal de X dans (2.10) et k = 0. Ce procédé permet donc d’avoir une expression
explicite de la solution à approcher.
Pour x0 tel que (xi0 )1≤i≤3 ∈ {−.7, −.3, .3, .7}, nous avons pris NMC = 106 trajectoires pour la simulation
de Monte Carlo et et nous avons fait varier h dans {.01, .05, .1}. Pour tous les calculs effectués, la taille de
l’intervalle de confiance à 95% était comprise dans [1.5 × 10−3 , 2 × 10−3 ]. Nous présentons au tableau 2.1 les
résultats obtenus pour le supremum de la valeur absolue des erreurs absolues et relatives pour les 3 × 43 = 192
points de la grille spatiale. Ces résultats semblent indiquer que pour la correction le terme de reste en o(h1/2 )
dans le Théorème 2.2.3 est en fait un O(h).
Nous concluons cette section en illustrant numériquement le résultat du Théorème 2.2.2. Pour x0 = (−.7, .3, .7)
et les précédentes valeurs de h, nous représentons en Figure 2.4 la quantité − log(ErrMC ) en fonction de − log(h),

18
 h Sans correction Correction de domaine
.1 0.169 (199%) 0.0220 (24.4%)
.05 0.114 (133%) 0.0115 (13.1%)
.01 0.0471 (54.7%) 0.0026 (2.98%)

Table 2.1 – Supremum de la valeur absolue du schéma d’Euler (erreur relative en % entre parenthèses)

( Z !)
MC
X τ h,i
où ErrMC := 1
MC g(Xτh,i
h,i ) + h,i
f (Xφ(s) )ds − u(x0 ). La courbe ainsi obtenue est assez proche d’une
i=1 0
droite de pente 1/2 telle que l’indique le théorème.

-Log(Err_MC)) in function of -Log(Delta)

4.5
-Log(Err_MC)
.5*x +1

3.5

2.5

2
2 2.5 3 3.5 4 4.5 5
-Log(Delta)

Figure 2.4 – Erreur de la méthode de Monte Carlo (sans correction) en fonction de h (échelle logarithmique).
Evaluation en x0 = (−.7, .3, .7).

2.4 Echantillonage discret de processus d’Itô

En ce qui concerne l’étude de l’erreur faible pour des processus stoppés/tués, le dernier axe abordé est celui
de l’observation discrète de processus d’Itô.
Le point principal dans ce cadre non Markovien est que l’on ne dispose plus d’une EDP de référence (2.3)
pour analyser l’erreur. L’idée développée dans [4] pour un domaine cylindrique régulier, D :=]0, T [×D, consiste à
exploiter que le processus Ut := E[g(T ∧τt , ΠD̄ (XT ∧τt ))|Ft ], τt := inf{s ≥ t : Xs 6∈ D}, t ≥ 0 est une martingale.
Ce processus est un objet naturel, en effet dans le cas Markovien, pour f = k = 0, on a bien Ut := u(t, Xt ) où
u est solution de (2.3). La projection sur D̄ est introduite pour la cohérence avec les observations discrètes. En
écrivant les accroissements de cette martingale entre les instants de discrétisation, on arrive à se ramener à une
discussion sur les temps de sortie proche de celles de [1],[7].
Sous des hypothèses de type (A), de frontière non caractéristique, bornitude des coefficients b, σ, de continuité
en probabilité de l’intégrale stochastique à l’ordre usuel 1/2 en temps, pour un domaine D de classe C 2 est une
fonction g ∈ C 1,2 ([0, T ] × Rd ), nous retrouvons le contrôle
√
|Ex [g(T ∧ τ h , ΠD̄ (XT ∧τ h ))] − Ex [g(T ∧ τ, ΠD̄ (XT ∧τ ))]| ≤ C h,

où τ h := inf{ti := ih : Xti 6∈ D}. Dans le cas tué, i.e. lorsque g(s, x) = Is=T g0 (x), (s, x) ∈ [0, T ] × D̄, et
que g0 bornée vérifie une condition de support, d(supp(g0 ), ∂D) ≥ ε > 0, le résultat précédent s’étend à une
intersection de domaines réguliers.

Extensions et perspectives
Le bilan des travaux [1],[3],[4],[7] est que nous avons pu obtenir dans le cas des processus de diffusion un
développement de l’erreur et une procédure d’accélération de convergence associée en contrôlant précisément
l’overshoot du schéma d’Euler.
Une question naturelle, tant pour la modélisation financière que pour l’approximation numérique d’équations
intégro-différentielles, est celle du comportement de l’erreur faible lorsque le processus dirigeant l’EDS (2.1)

19
comporte une partie de saut. Intuitivement il semble naturel de penser que dans le cas où l’on ajoute dans la
dynamique différentielle les variations d’un processus de Lévy d’intensité finie, il sera possible de distinguer dans
l’analyse de l’erreur les sorties dues aux grands sauts de celles dues à l ’observation discrète du schéma. Par
ailleurs, les propriétés de la solution de l’EDP associée ont été étudiées par Garroni et Menaldi [GM02] dans le
cas elliptique ou Bensoussan et Lions dans le cas parabolique [BL82], cf. Chapitre 3 Section 3. La difficulté reste
de bien contrôler le gradient ponctuellement jusqu’au bord du domaine. Dans ce cadre, l’erreur faible devrait
néanmoins rester du même ordre.
Les résultats numériques de la Section 2.3 suggèrent que sous de bonnes hypothèses de régularité le terme
de reste du Théorème 2.2.3 est un O(h). La question de la poursuite du développement de l’erreur et de la
correction associée aux ordres supérieurs reste ouverte.
Un dernier développement possible de cette thématique est associé à la caractérisation de la régularité des
noyaux de la chaleur dans des domaines non-réguliers. Des techniques variationnelles ont été développées, voir
par exemple Grisvard [Gri92] ou Kozlov et al. [KMR97]. L’approche par décomposition spectrale proposée dans
[BS97] et les contrôles sur le problème elliptique en dimension 2, que l’on peut penser être de même nature que
ceux associés à la trace d’un cône tridimensionnel sur la sphère, devraient permettre de passer à la dimension
3 et de procéder par suite à une sorte de récurrence sur la dimension. Indiquons également que pour des cônes
convexes des techniques de fonctions barrières permettent d’obtenir assez directement la Lipschitz continuité
jusqu’au bord, cf. Chapitre 4, Section 4.4.

20
 Chapitre 3

La méthode parametrix : théorèmes

limites et développement d’erreur

Nous présentons dans ce chapitre les résultats des articles [8] (en collaboration avec V. Konakov de l’Univer-
sité de Moscou et S. Molchanov de l’Université de Charlotte) et [9] (en collaboration avec V. Konakov). Le point
commun à ces deux travaux est l’utilisation de la représentation de la densité de certains processus en terme
de série parametrix. L’idée de cette représentation est ancienne et remonte pour les opérateurs paraboliques
d’ordre deux non dégénérés sous forme non divergence à E. Levi, cf. Friedman [Fri64] pour une présentation.
Elle se base sur une idée de continuité ou de méthode perturbative. En temps petit, la densité d’une EDS à
coefficients variables est a priori proche de la densité d’une EDS à coefficients constants, pour laquelle on a de
bons contrôles de densité. L’idée de la méthode est ensuite d’utiliser les équations de Kolmogorov satisfaites par
ces deux densités pour précisément estimer la différence.
Outre l’approche de Levi, d’autres versions de la méthode parametrix ont été développées. En particulier,
c’est l’approche proposée par McKean et Singer [MS67] pour obtenir des développements asymptotiques du
spectre du Laplacien sur une variété en fonction de sa courbure, qui dans la mesure où elle admet un pendant
discrétisé direct, est la plus propice pour analyser les erreurs associées aux schémas d’approximation. C’est cette
approche qui a permis à Konakov et Mammen d’obtenir dans le cadre des processus de diffusion non dégénérés
inhomogènes un théorème limite local sur les densités pour l’approximation par chaı̂ne de Markov [KM00] ainsi
qu’un développement d’erreur pour le schéma d’Euler [KM02].
Nous présenterons d’abord dans le cadre abstrait des processus de “diffusion-saut” la méthode parametrix de
Mc Kean et Singer, cf. Section 3.1. Nous illustrerons ensuite la robustesse de cette approche pour les problèmes
de discrétisation à travers deux exemples :
- l’approximation par chaı̂ne de Markov de processus dégénérés, qui interviennent entre autres dans des modèles
cinétiques et en finance (options asiatiques), cf. Section 3.2 et article [8].
- le développement d’erreur pour des EDS dirigées par des processus stables. cf. Section 3.3 et article [9].

3.1 Parametrix à la Mc Kean et Singer pour les processus de “diffu-

sion-saut”
3.1.1 Présentation de la méthode
Soit (Xt )t≥0 un processus de Markov à valeurs dans Rd de générateur de la forme :
Z  
1 2

hy, ∇ϕ(x)i
Lt ϕ(x) = hb(t, x), ∇ϕ(x)i + tr a(t, x)Dx ϕ(x) + ϕ(x + y) − ϕ(x) − ν(t, x, dy),(3.1)
2 Rd 1 + |y|2

pour t ≥ 0, a(t, .) matrice de diffusion semi-définie positive, b(t, .) coefficient de dérive borné, ν mesure de Lévy
et ϕ ∈ C0∞ (Rd ) fonction test. Les opérateurs de la forme 3.1 vérifient le principe du maximum faible et sont
quasi locaux, i.e. si supp(ϕ) ∈ B(x, ε), |Lt ϕ(x)| ≤ Cε |ϕ|∞ , cf. Hille and Philips [HP57].
Intuitivement, pour h > 0 petit, l’incrément Xt+h − Xt est “proche” de l’incrément du processus à accrois-
sements indépendants et stationnaires (PAIS) de caractéristiques gelées en b(t, Xt ), a(t, Xt ), ν(t, Xt , .). C’est là
l’idée qui est à la base de l’approche de Levi. Supposons désormais que le processus X admet pour tout t > 0

21
une densité p par rapport à la mesure de Lebesgue et que l’on s’intéresse à l’estimation de p(t, t + h, Xt , x′ ),
i.e. densité en t + h au point x′ . Il est alors également naturel de dire que l’accroissement Xt+h − Xt est
proche de celui du PAIS de caractéristiques gelées en espace au point final où l’on considère la densité, i.e.
b(t, x′ ), a(t, x′ ), ν(t, x′ , .). Cette autre approximation est à la base du travail de Mc Kean et Singer.
Pour obtenir une représentation de la densité p à partir de ces approximations heuristiques, il faut ensuite
quantifier précisément l’erreur d’approximation. Comme l’on a supposé que X admettait une densité pour tout
′
t > 0 il en est de même pour le processus à caractéristiques gelées. Notons p̃x cette densité, où l’indice x′
′
indique explicitement le point de congélation, et L̃x son générateur :
Z  
′ 1
hy, ∇ϕ(x)i
L̃xt ϕ(x) = hb(t, x′ ), ∇ϕ(x)i + tr a(t, x′ )Dx2 ϕ(x) + ϕ(x + y) − ϕ(x) − 2
ν(t, x′ , dy).
2 R d 1 + |y|
(3.2)

Sous de bonnes hypothèses, la densité p vérifie les équations de Kolmogorov progressives et rétrogrades, i.e

∂t p(s, t, x, y) = L∗t p(s, t, x, y), t > s, (x, y) ∈ (Rd )2 , p(s, t, x, .) → δx (.),

t↓s
d 2
∂s p(s, t, x, y) = −Ls p(s, t, x, y), t > s, (x, y) ∈ (R ) , p(s, t, ., y) → δy (.), (3.3)
s↑t

où L∗t désigne l’adjoint formel de Lt (nous n’avons en effet pas précisé d’hypothèses de régularité sur les
′
coefficients dans (3.1)). La densité p̃x vérifie également l’équation rétrograde :
′ ′ ′ ′
∂s p̃x (s, t, x, y) = −L̃xs p̃x (s, t, x, y), t > s, (x, y) ∈ (Rd )2 , p̃x (s, t, ., y) → δy (.). (3.4)
s↑t

Ainsi,
Z t Z
x′ ′ ′
(p − p̃ )(s, t, x, x ) = du∂u dzp(s, u, x, z)p̃x (u, t, z, x′ )
s Rd
Z t Z  
′ ′
= du dz ∂u p(s, u, x, z)p̃x (u, t, z, x′ ) + p(s, u, x, z)∂u p̃x (u, t, z, x′ )
s Rd
Z t Z  
′ ′ ′
= du dz L∗u p(s, u, x, z)p̃x (u, t, z, x′ ) − p(s, u, x, z)L̃xu p̃x (u, t, z, x′ )
s Rd
Z t Z
′ ′
= du dzp(s, u, x, z)(Lu − L̃xu )p̃x (u, t, z, x′ ), (3.5)
s Rd

où l’on utilise la convergence vers la masse de Dirac pour la première égalité, les équations de Kolmogorov (3.3)
et (3.4) pour la troisième et le passage à l’adjoint pour la dernière. Nous avons aussi formellement dérivé sous
l’intégrale en espace pour la deuxième égalité. Z Z t
Introduisons la notation f ⊗ g(s, t, x, x′ ) = du dzf (s, u, x, z)g(u, t, z, x′) pour la convolution temps-
s Rd
′ ′
espace et posons p̃(s, t, x, x′ ) := p̃x (s, t, x, x′ ), H(s, t, x, x ) := (Ls − L̃xs )p̃(s, t, x, x′ ). L’équation (3.5) se réécrit
′

alors :
(p − p̃)(s, t, x, x′ ) = p ⊗ H(s, t, x, x′ ).
A partir de cette expression, l’idée est ensuite d’itérer cette procédure pour p(s, u, x, z) dans (3.5) en utilisant
cette fois ci la densité de transition d’un processus de caractéristiques gelées en z variable d’intégration. On fait
de la sorte apparaı̂tre les convolutions itérées du noyau H pour obtenir le développement formel :
∞
X
p(s, t, x, x′ ) = p̃ ⊗ H (r) (s, t, x, x′ ), (3.6)
r=0

où H (0) = I, H (r) = H ⊗ H (r−1) , r ≥ 1. Déduire des estimées sur p à partir de l’expression formelle (3.6)
suppose d’avoir de bons contrôles sur le membre de droite. De façon générale, on utilisera les contrôles a priori
connus sur la densité gelée et des propriétés de régularisation du noyau de convolution H. Nous renvoyons à la
Section 3.1.2 pour une illustration dans le cas des diffusions.
Nous avons jusqu’à présent supposé l’existence de la densité p. Lorsque les coefficients dans (3.1) sont réguliers
(en particulier Lipschitz continus en espace), il est facile de voir Xt comme la solution forte d’une EDS. Sous
des hypothèses de non dégénérescence des coefficients, le calcul de Malliavin permet d’obtenir l’existence et la

22
régularité de la densité, voir Norris [Nor86] pour le cadre de diffusions hypoelliptiques homogènes, Cattiaux et
Mesnager [CM02] pour le cas inhomogène, et Bichteller et al. [BGJ87] pour le cas avec sauts.
Réciproquement, il est toujours possible de partir du terme de droite de l’équation (3.6) résultant du
précédent développement formel. Si la série ainsi que ses dérivées par rapport aux variables de temps et d’espace
sont bien définies, il vient bien par construction que p est une solution des équations de Kolmogorov. Le principe
du maximum donne l’unicité. Si cette procédure aboutit, on construit bien ainsi l’unique solution fondamentale
de l’EDP parabolique intégro-différentielle associée à (3.1). L’identification à la densité du processus de Markov
de générateur (3.1) se fait ensuite par un argument à la Dynkin, cf. Théorème 2.3 de [Dyn63] .
C’est ce que nous allons faire dans la section suivante dans le cadre de processus de diffusion inhomogènes
à coefficient de diffusion uniformément Hölder continu en espace, mesurable en temps et dérive bornée.

3.1.2 L’exemple des processus de diffusion non dégénérés

Considérons le cas ν = 0 dans (3.1) et supposons également que la matrice de diffusion est uniformément
elliptique bornée de constante Λ ≥ 1 (i.e. Λ−1 Id ≤ a(t, x) ≤ ΛId , ∀(t, x) ∈ R+ × Rd ) et uniformément η-Hölder
continue en espace, η ∈ (0, 1]. Le coefficient b est supposé mesurable borné.
Le problème de martingale est alors bien posé, cf. Théorème 7.2.1 de Stroock et Varadhan [SV79]. On
′
peut définir p comme en (3.6). La densité gelée p̃(s, t, x, x′ ) n’est autre que celle du processus Gaussien X̃tx =
Rt R t
x+ s b(u, x′ )du+ s σ(u, x′ )dWu , où σσ ∗ (u, x′ ) = a(u, x′ ), au point de congélation x′ . Les hypothèses précédentes
sur les coefficients a, b donnent en particulier que pour T > 0 fixé, il existe C := C(T, Λ, |b|∞ ), c := c(Λ) telles
que pour tout multi-indice de dérivation α, |α| ≤ 4, et pour tout 0 < s < t < T , (x, x′ ) ∈ (Rd )2 ,
 d/2  
C c |x − x′ |2
|∂xα p̃(s, t, x, x′ )| ≤ pc (t − s, x, x′ ), pc (t − s, x, x′ ) := exp −c . (3.7)
(t − s)|α|/2 2π(t − s) t−s

Le noyau de convolution s’écrit :

′
H(s, t, x, x′ ) = (Ls − L̃xs )p̃(s, t, x, x′ )
1

= hb(s, x) − b(s, x′ ), ∇x p̃(s, t, x, x′ )i + tr (a(s, x) − a(s, x′ ))Dx2 p̃(s, t, x, x′ ) . (3.8)
2
En particulier, le contrôle (3.7) donne :
 
′ 2|b|∞ |a(s, x) − a(s, x′ )|
|H(s, t, x, x )| ≤ C + pc (t − s, x, x′ )
(t − s)1/2 (t − s)
η
≤ C((t−1/2 + (t − s)−1+ 2 )pc (t − s, x, x′ ),

où l’on utilise la η-Hölder continuité de a pour compenser la singularité en temps non intégrable des termes
d’ordre 2 en modifiant le c pour la dernière inégalité. Le noyau H est donc régularisant, au sens où lorsqu’on le
convole on obtient :
∃C := C(T, Λ, |b|∞ , η), c := c(Λ, η), ∀r ≥ 1,
r+1
Y  
′ (i − 1)η η
|p̃ ⊗ H (r)
(s, t, x, x )| ≤ C r+1 rη/2
(t − s) B 1+ , pc (t − s, x, x′ ), (3.9)
i=1
2 2
R1
où B(m, n) := 0 sm−1 (1 − s)n−1 ds désigne la fonction β. Pour t − s petit, ce terme est donc de plus en plus
négligeable lorsque r croı̂t.
De (3.9) et (3.6) on déduit directement des asymptotiques de la fonction β une borne supérieure de type Aron-
son pour p ainsi définie, i.e. ∃C := C(T, Λ, |b|∞ , η), c := c(Λ, η), ∀0 ≤ s < t ≤ T, (x, x′ ) ∈ (Rd )2 , p(s, t, x, x′ ) ≤
Cpc (t − s, x, x′ ).
Par ailleurs, l’équation (3.6) permet de justifier que p vérifie l’équation de Kolmogorov rétrograde (∂s +
Ls )p(s, t, x, x′ ) = 0, 0 ≤ s < t ≤ T . La propriété régularisante en temps de H, cf. (3.9), permet ensuite
de déduire p(s, t, ., x′ ) → δx′ (.). On a donc bien construit une solution fondamentale de l’EDP parabolique.
s↑t
L’identification avec la densité de transition du processus X solution du problème de martingale se fait de façon
standard, cf. Chapitre 3 de [SV79] ou Dynkin [Dyn63].

23
 3.1.3 Représentation parametrix pour un schéma de discrétisation
Considérons un processus de diffusion-saut (Xt )t≥0 à valeurs dans Rd de dynamique :
Z t Z t
Xt = x + b(s, Xs− )ds + σ(s, Xs− )dZs , (3.10)
0 0

pour (Zt )t≥0 processus de Lévy dont le générateur peut se mettre sous la forme (3.1) et admet une densité.
Supposons également que Xt possède pour tout t ∈ (0, T ], T > 0 fixé, une densité que l’on peut exprimer sous
la forme (3.6).
Soit h > 0 un pas de discrétisation tel que T /h = N ∈ N∗ et {(ti = ih)i∈[[0,N ]]} la grille temporelle associée.
Une approximation naturelle de (3.10) par chaı̂ne de Markov est donnée par le schéma
h
Xthi+1 = Xthi + b(ti , Xthi )h + σ(ti , Xthi )ξi+1 , i ∈ [[0, N − 1]], X0h = x, (3.11)

(loi)
où les (ξih )i∈N∗ sont des variables à densité i.i.d. qui approchent l’incrément Zti+1 − Zti = Zh (i.e. Z0 = 0). Si
le coefficient σ vérifie une hypothèse de non-dégénérescence, alors pour tout i ∈ [[1, N ]], Xthi admet une densité.
Notons ph cette densité.
La méthode parametrix de Mc Kean et Singer présentée en Section 3.1 présente le grand avantage de
pouvoir assez directement se transposer au cas discret et d’établir une représentation de ph similaire à (3.6). Si
l’on s’intéresse à la représentation de la densité ph (tj , tj ′ , x, x′ ), 0 ≤ j < j ′ ≤ N , l’idée analogue au cas continu
′
est d’introduire une chaı̂ne de Markov X̃ h,x à coefficients gelés au point d’arrivée x′ où l’on considère la densité.
Sa dynamique s’écrit,
′ ′ ′
X̃th,x
i+1
= X̃th,x
i
+ b(ti , x′ )h + σ(ti , x′ )ξi+1
h
, i ∈ [[j, j ′ − 1]], X̃th,x
j
= x. (3.12)
′
On notera p̃h,x sa densité. Introduisons désormais les analogues discrets des générateurs définis en (3.1) et (3.2).
Soit pour ϕ ∈ C0∞ (Rd ) et j ∈ [[0, j ′ ),
′ ′
E[ϕ(Xthj +h )|Xthj = x] − ϕ(x) ′
e h,x )|X
E[ϕ(X eth,x = x] − ϕ(x)
e h,x tj +h
Lhtj ϕ(x) =
j
, et L tj ϕ(x) = .
h h
′
En utilisant la notation peh (tj , tj ′ , x, x′ ) = peh,x (tj , tj ′ , x, x′ ), on définit de façon similaire à (3.8) le noyau discret
H h par
 ′

H h (tj , tj ′ , x, x′ ) = Lhtj − L e h,x
tj peh (tj + h, tj ′ , x, x′ )
Z h i
′ ′
= h−1 ph − peh,x (tj , tj+1 , x, z)e ph,x (tj+1 , tj ′ , z, x′ )dz, 0 ≤ j < j ′ ≤ N.
Rd

Une simple manipulation algébrique, cf. Lemme 3.6 de [KM00], permet ensuite d’établir que pour 0 ≤ tj <
tj ′ ≤ T ,
′
jX −j  
′
h
p (tj , tj ′ , x, x ) = peh ⊗h H h,(r) (tj , tj ′ , x, x′ ), (3.13)
r=0
P ′ −j−1 R
où f ⊗h g(tj , tj ′ , x, x′ ) = jk=0 h Rd g(tj , tj+k , x, z)f (tj+k , tj ′ , z, x′ )dz est l’opérateur de convolution discret
h,(0) h,(r)
et H = I et H = H ⊗h H h,(r−1) , r ≥ 1.
h

D’après (3.6) et (3.13), les densités respectives du processus X et de son schéma d’Euler X h peuvent se
représenter sous des formes similaires. Les étapes pour contrôler l’écart (p − ph )(0, T, x, x′ ) sont les suivantes :
1. Contrôler les queues de la série dans (3.6), i.e. les termes r ≥ N .
2. Quantifier l’impact de la discrétisation de l’opérateur de convolution dans (3.6), i.e. lorsque l’on remplace ⊗
par ⊗h dans le développement tronqué en r = N .
3. Contrôler la convergence de p̃h vers p̃.
4. Appliquer des arguments fins de stabilité pour quantifier les différences entre les noyaux de convolution continus
H et discrets H h .

24
 Dans [KM00], Konakov et Mammen ont pu appliquer ce programme pour des processus de diffusion inhomogènes
à coefficients bornés Lipschitziens et non dégénérés, i.e. la matrice de diffusion a est uniformément elliptique.
Dans ce cadre, c’est l’étape 3. qui donne la vitesse de l’erreur d’approximation qui est de h1/2 et qui est celle
du théorème limite local sur les densités en régime Gaussien, cf. Bhattacharya et Rao [BR76]. Les auteurs
′
−x S −1
établissent que |(p − ph )(0, T, x, x′ )| ≤ h1/2 π(T, x, x′ ), π(T, x, x′ ) ≤ CT −d/2 (1 + | xT 1/2 | ) où l’exposant S est
fonction des propriétés d’intégrabilité des variables (ξih )i≥1 et C est une constante dépendant des coefficients,
uniforme sur les compacts en temps. Ces mêmes techniques leur ont permis d’établir des asymptotiques précises
de l’erreur en temps T petit [KM09].
h
Lorsque les ξi+1 ont même loi que l’accroissement Zti+1 − Zti , i.e. X h est le schéma d’Euler, l’étape 3.
ci-dessus n’est pas nécessaire. On retrouve alors que l’erreur dominante est celle associée à la discrétisation de
l’opérateur de convolution qui est intuitivement de l’ordre de h. C’est ce qu’ont montré Konakov et Mammen
[KM02] pour des processus de diffusion homogènes non dégénérés à coefficients bornés et Lipschitziens, i.e.
pour T = 1, (p − ph )(1, x, x′ ) = hπ(1, x, x′ ) + h2 r(1, x, x′ ), (|π| + |r|)(1, x, x′ ) ≤ C exp(−c|x − x′ |2 ), C :=
C(1, Λ, |b|∞ ), c := c(Λ). L’extension au cas inhomogène ne présente pas de difficultés et se ferait suivant la
même procédure.
Un développement au même ordre a été obtenu dans le cadre hypoelliptique par Bally et Talay [BT96b] pour
un schéma d’Euler perturbé. L’une des difficultés dans ce cadre est en effet de prouver l’existence d’une densité
pour le schéma d’approximation ce que permet la perturbation. Nous avons retrouvé cette difficulté dans [8]
que nous présentons dans la section suivante.

3.2 Théorème limite local en régime Gaussien dans un cas dégénéré

Nous présentons dans cette section les résultats de [8]. Nous avons pu établir un théorème limite local pour
l’approximation par chaı̂ne de Markov d’un processus de diffusion dégénéré (Xt1 , Xt2 )t≥0 à valeurs dans Rd × Rd
de dynamique :
Z t Z t
Xt1 = x1 + b(Xs1 , Xs2 )ds + σ(Xs1 , Xs2 )dWs ,
0 0
Z t
Xt2 = x2 + Xs1 ds, (3.14)
0

où W est un mouvement Brownien de Rd et où l’on suppose

(S) b, σ sont bornés, Lipschitz continus et de classe C 1 .
(UE) La matrice de diffusion a : R2d → Rd ⊗ Rd est uniformémement elliptique, ce qui associé à l’hypothèse
(S) donne : ∃Λ ≥ 1, ∀x = (x1 , x2 ) ∈ R2d , ∀ξ ∈ Rd , Λ−1 |ξ|2 ≤ ha(x)ξ, ξi ≤ Λ|ξ|2 .
Comme nous l’avons mentionné au Chapitre 1, ce type d’équations intervient dans de nombreuses applications,
en finance pour les problèmes de valorisation/couverture d’options Asiatiques, en physique pour les modèles
cinétiques ou les systèmes Hamiltoniens.

3.2.1 Contrôles sur l’objet continu

Sous (S), (UE), le théorème de Hörmander donne directement l’existence de la densité pour Xt , t > 0.
Pour comprendre plus précisément la structure de cette densité, considérons d’abord le cas de coefficients b, σ
constants pour lequel Xt est un processus Gaussien. Il suffit alors pour exprimer sa densité de calculer sa
moyenne et sa matrice de covariance. L’expression de cette densité remonte au travail de Kolmogorov [Kol34].
Dans ce cas dégénéré très simple on remarque déjà l’échelle intrinsèque non diffusive due à la dégénérescence, i.e.
la première composante vit à l’échelle de temps caractéristique du Brownien en t1/2 , la deuxième composante
qui est l’intégrale de la première vit elle à l’échelle t3/2 . Un autre aspect remarquable est dû à la non bornitude
du coefficient de la deuxième composante, qui induit un transport de la condition initiale x1 dans Xt2 , i.e.
Rt Rs
E[Xt2 ] = x2 + tx1 + 0 E[ 0 b(Xu1 , Xu2 )du]ds.
La solution fondamentale de l’EDP parabolique associée au générateur de (3.14) a, sous des hypothèses
supplémentaires de régularité des coefficients, été étudiée par Weber [Web51]. L’approche est celle du parametrix
à la Levi. L’auteur obtient une borne supérieure sur la densité p de (3.14) similaire à celle de l’exemple de

25
Kolmogorov, soit pour T > 0 fixé,
∃C := C(T, Λ, |b|∞ ), c := c(Λ), ∀(t, x, x′ ) ∈ (0, T ] × Rd × Rd , p(t, x, x′ ) ≤ Cpc (t, x, x′ ),
√ !d/2   ′ x +x′ 
′ 3c |x1 − x1 |2 |x′2 − x2 − 1 2 1 t|2
pc (t, x, x ) := exp −c + 3 . (3.15)
2π(t − s)2 4t t3

A la Section 3.1.2 nous avons gelé les coefficients au point d’arrivée où l’on souhaitait estimer la den-
sité. Au vu des contrôles (3.7) et de la définition de l’opérateur de convolution (3.8), cela permettait en
effet bien d’équilibrer la singularité associée au termes d’ordre deux. Le point clé est qu’il y avait dans ce
cas “compatibilité” entre le terme apparaissant dans l’exponentielle exp(−c|x′ − x|2 /t) et celui résultant de
′
−x|η −1+η/2
|a(t, x) − a(t, x′ )|t−1 ≤ C |xtη/2 t , i.e. cette singularité était compensée par l’exponentielle. Dans le cas
présent ce choix de processus gelé n’est pas satisfaisant dans la mesure où il ne permet pas de prendre en compte
le transport de la condition initiale pour la deuxième composante qui apparaı̂t dans l’expression explicite (3.15).
Pour adapter la méthode de Mc Kean et Singer au cas de processus à coefficient de dérive non borné, l’idée est
de construire un processus Gaussien à coefficients gelés qui permette d’avoir lors de la différence des générateurs
un terme homogène à la composante dans l’exponentielle. Le candidat naturel est :
Z s
e 1,t,x′
Xs = x1 + σ(x′1 , x′2 − x′1 (t − u))dWu + b(x′1 , x′2 )s,
0
Z s
Xes2,t,x′ = x2 + eut,x′ du, s ∈ [0, t].
X (3.16)
0

En effet, les parties dégénérées des deux générateurs sont les mêmes et se simplifient. Par ailleurs, la différence
entre les termes d’ordre 2 des générateurs de (3.14) et (3.16) est dominée par |a(x′1 , x′2 −x′1 t)−a(x)|t−1 p̂c (t, x, x′ ) ≤
x′ +x
C(|x′1 − x1 |(1 + 2t ) + |x′2 − x2 − 1 2 1 t|)t−1 p̂c (t, x, x′ ) ≤ Ct−1/2 pc (t, x, x′ ) grâce à l’hypothèse de Lipschitz conti-
nuité de a et avec une modification de c pour la dernière inégalité. Ce qui permet bien de rattrapper la singularité.
Nous pouvons ainsi retrouver sous nos hypothèses courantes, plus faibles, le contrôle de type Aronson obtenu
dans [Web51].
Proposition 3.2.1 Soit T > 0 fixé. Sous (UE), (S), il existe des constantes C := C(T, Λ, |b|∞ ), c := c(Λ)
telles que pour tout (t, x, x′ ) ∈ (0, T ] × (Rd )2 ,
p(t, x, x′ ) ≤ Cpc (t, x, x′ ),
où pc est défini en (3.15).
En fait, de façon similaire à la Section 3.1.2, il suffit pour obtenir cette borne de considérer le coefficient a,
η-Hölder continu en espace et b borné. De même une dépendance en temps des coefficients n’induirait aucune
difficulté additionnelle.
Pour conclure cette partie, mentionnons que le choix des paramètres de congélation qui permet de contrôler
le développement (3.6) est fortement lié à la métrique dérivant du problème de contrôllabilité déterministe
associé à l’EDS de départ. Le terme dans l’exponentielle apparaı̂t alors comme l’énergie du contrôle optimal qui
permet de joindre x à x′ en un temps t. Nous reviendrons en détails sur ce point au Chapitre 5.2.

3.2.2 Approximation par chaı̂ne de Markov

On introduit tout d’abord un schéma de discrétisation du processus gelé (3.16). Soit h0 > 0  un pas de temps
′ 2d 2d h0 ,t,x′ e h0 ,1,t,x′ h0 ,2,t,x′
“microscopique”, on définit pour (x, x ) ∈ R ×R , t > 0 fixés, X0 = X0 , X0 = x = (x1 , x2 )
et pour i ≥ 0,
′ ′ p
Xeth0 ,1,t,x = X eth0 ,1,t,x + b(x′ )h0 + σ(x′1 , x′2 − tx′1 ) h0 ξei+1 ,
i+1 i
′ ′ ′
eth0 ,2,t,x
X = X eth0 ,1,t,x h0
eth0 ,2,t,x + X
i+1 i i+1
′ ′
= X e h0 ,1,t,x + h2 b(x′ ) + h3/2 σ(x′ , x′ − tx′ )ξei+1 ,
e h0 ,2,t,x + h0 X (3.17)
ti ti 0 0 1 2 1

où les (ξei )i∈N∗ sont des variables centrées i.i.d. de densité q, de covariance l’identité. Si l’on itère n fois le schéma
il vient
′ p
Xeth0 ,1,t,x = x1 + (nh0 )b(x′ ) + σ(x′1 , x′2 − x′1 t) nh0 ξen(1) ,
n

eth0 ,2,t,x = x2 + (nh0 )x1 + γn (nh0 )2 b(x′ ) + (nh0 )3/2 σ(x′1 , x′2 − x′1 t)ξen(2) ,
′
X (3.18)
n
2

26
 (1) Pn e ξen(2) = Pn n−(i−1) e
où γn = (1 + n1 ) et ξen = √1
n i=1 ξi ,
√1
n i=1 ( n )ξi .
Sous des hypothèses d’intégrabilité et régularité de la densité q des (ξei )i∈N∗ , on peut montrer en calcu-
(1) (2)
lant explicitement la fonction caractéristique que pour n suffisamment grand (ξen , ξen ) admet une densité qn
possédant des propriétés 
de régularité et intégrabilité 
similaires à celles de q, cf. Proposition 4.2 de [8].Par
(1) (2) 1 γn /2 1 1/2
ailleurs Cov(ξen , ξen ) = qui est “proche” pour n grand de =
γn /2 1/3γn (1 + 1/(2n)) 1/2 1/3
R1
Cov(W1 , 0 Ws ds).
Si l’on pose maintenant h := nh0 , pas de temps “macroscopique”, et que l’on s’intéresse à l’approximation
de la densité en tj := jh, j ∈ (0, N ]], N h = T , on introduit le schéma gelé agrégé, dérivant de (3.18) et résultant
du regroupement de n itérations de (3.17) où l’on prend t = tj − ti au ième instant de discrétisation de l’échelle
′
macroscopique. Précisément, pour (x, x′ ) ∈ R2d × R2d , j ∈ (0, N ]] on définit X e h,tj ,x = x et pour i ∈ [[0, j − 1]],
0
′ ′ √ 1
eth,1,tj ,x
X eth,1,tj ,x
= X + b(x′ )h + σ(x′1 , x′2 − x′1 (tj − ti )) he
ηi+1 ,
i+1 i
n √ 2 o
eth,1,tj ,x + γn b(x′ )h + σ(x′1 , x′2 − x′1 (tj − ti )) he
′ ′ ′
eth,2,tj ,x
X eth,2,tj ,x
= X + X ηi+1 h, (3.19)
i+1 i i
2
où les variables i.i.d. (η̃i1 , η̃i2 )i∈(0,j]] ont pour densité qn (.).
Le schéma précédent guide le choix de la dynamique de la chaı̂ne de Markov qui va approcher X solution
de (3.14). Soit pour (x, x′ ) ∈ R2d × R2d , X0h = x et pour i ∈ [[0, N − 1]],
√ 1
Xth,1i+1
= Xth,1
i
+ b(Xthi )h + σ(Xthi ) hηi+1 ,
γ n
√ 2
Xth,2 = Xth,2 + (Xth,1 + b(Xthi )h + σ(Xthi ) hηi+1 )h, (3.20)
i+1 i i
2
où les (ηi1 , ηi2 )i∈(0,N ]] ont également pour densité qn .
En résumé, l’itération du schéma “naturel” gelé (3.17) permet d’obtenir une densité après suffisamment de
pas de temps “microscopiques”. Ceci donne la forme du schéma gelé à l’échelle “macroscopique”, cf. (3.19).
La structure de ce schéma donne ensuite celle de la chaı̂ne de Markov associée à l’approximation du processus
initial. A partir des équations (3.20), (3.19), on peut obtenir pour la densité ph de la chaı̂ne de Markov X h une
représentation de type (3.13) où p̃h est la densité associée à (3.19). Cette représentation sera comparable au
développement (3.6) obtenu pour le processus continu à la Section 3.2.1 à partir de (3.16). En suivant ensuite
le programme détaillé en fin de Section 3.1.3 on obtient :
Théorème 3.2.1 (Théorème limite local pour les densités) . R
Supposons (UE) et (B). Si pour tout multi-indice de dérivation ν, |ν| ≤ 4, Rd |u|S |Duν qn (u)|du < +∞, S >
4(d + 1) + 2d2 , alors ∃c > 0,
"   ′ #−1
′ ′ ′ ′ x1 + x1
sup (1 + |x1 | + |x1 |) sup pc (T (1 + δ), x, x ) + χT 1/2 x1 − x1 , x2 − x2 − T
x,x′ ∈R2d δ∈[0,1] 2
×|ph (T, x, x′ ) − p(T, x, x′ )| = O(h1/2 ),
où pc est défini en (3.15) et ∀u = (u1 , u2 ) ∈ R2d ,
 (S−4)/4d−1 !−1
−2d u1 2 u2 2
χT 1/2 (u) = T 1+ + 3/2 .
T 1/2 T

La vitesse de convergence reste celle du théorème limite local Gaussien dans la mesure où la densité p̃ du
processus gelé continu (3.16) reste Gaussienne. L’ordre de h1/2 est celui de la convergence de p̃h vers p̃. Le
terme en pc provient d’une correction technique nécessaire pour traiter le coefficient non borné dans l’analyse
de stabilité.
Pour illustrer ce résultat on peut considérer l’exemple de l’approximation numérique de Px [(T −1 XT2 −XT1 )+ >
K] pour K > 0 donné. Cette quantité correspond en finance au prix d’une option Call digital de maturité T
lorsque la dynamique de l’actif est donnée par la première composante de (3.14). Si l’on associe maitenant les
composantes du système (3.14) à la vitesse et à la position d’une particule, Px [(T −1 XT2 − XT1 )+ > K] est un
quantile de la partie positive de l’écart entre vitesse moyenne et vitesse instantanée en T . Le théorème limite
local 3.2.1 permet de directement relier les densités des objets discrets et continus de dynamiques respectives
données en (3.14) et (3.20). Cette approche est préférable à une discrétisation de la composante non-dégénérée
associée une estimation numérique de l’intégrale dans la mesure où dans ce cas là le schéma d’approximation
peut ne pas avoir de densité.

27
3.3 Approximation d’EDS dirigées par un stable : développement
d’erreur des densités
On présente ici les résultats de l’article [9]. Nous considérons le cas où X à valeurs dans Rd suit une dyna-
mique de type (3.10). Les coefficients b, σ sont homogènes en temps et Z est un processus α-stable symétrique,
i.e.
Z t Z t
Xt = x + b(Xs− )ds + σ(Xs− )dZs ,
0 0

et pour t ≥ 0, u ∈ Rd ,
Z
E[exp(ihu, Zt i)] = exp(ithγ, ui − t |hs, ui|α λ(ds)), (3.21)
S d−1

où λ est une mesure symétrique sur la sphère unité S d−1 , α ∈ (0, 2). On exclut le cas Brownien pour lequel on
peut se reporter à [KM02] pour des résultats similaires.
Dans ce cadre, le schéma d’Euler associé à (3.11) est bien défini et explicitement simulable, cf. [PT97] ou
[ST94]. Pour un pas h > 0 et la grille de discrétisation {(ti := ih)i∈N∗ } associée, on pose X0h = x et pour i ≥ 0,

Xthi+1 = Xthi + b(Xthi )h + σ(Xthi )(Zti+1 − Zti ).

En revanche, les hypothèses d’intégrabilité sur la mesure de Lévy requises dans Protter et Talay [PT97] pour
l’analyse de l’erreur faible ne sont pas satisfaites par les processus stables. Ces derniers ont en effet des queues
lourdes. Dans le cas monodimensionnel et pour des fonctions tests bornées, le contrôle de l’erreur faible avait
été obtenu par Hausenblas [Hau02] à l’aide du calcul de Malliavin. L’approche parametrix permet de contrôler
l’erreur entre les densités.
Sous de bonnes hypothèses spécifiées ci-après, Kolokoltsov a donné une représentation en série (3.6) de la
densité p de X. Il a également établi des contrôles à la Aronson avec une échelle de temps caractéristique en
t1/α généralisant le cas Gaussien, cf. Théorème 3.2 de [Kol00]. Nous avons de façon similaire au procédé présenté
en Section 3.1.3 obtenu le développement (3.13) pour la densité du schéma d’Euler. En suivant le programme
détaillé en fin de Section 3.1.3 nous avons pu établir un développement d’erreur sur les densités à l’ordre h. Pour
ce faire il a fallu dans l’analyse de stabilité raffiner les arguments de [Kol00] et utiliser des découpages minutieux
dans les parties d’analyse de Fourier utilisées pour représenter densités gelées et noyau de convolution.
Soit q ∈ N∗ . Introduisons les hypothèses
(A1-q) Pour d ≥ 2, la mesure sphérique λ a une densité de surface C q (S d−1 ) et pour d ≥ 1, il existe 0 < C1 ≤
C2 < +∞, ∀p ∈ Rd , Z
α
C1 |p|α ≤ |hp, si| λ (ds) ≤ C2 |p|α .
S d−1

(A2-q) Les coefficients b, σ et leurs dérivées jusqu’ à l’ordre q sont uniformément bornées. Ainsi, pour 1 < α < 2,
B(x) := b(x) + σ(x)γ est uniformément borné. On suppose pour 0 < α ≤ 1, que B(x) = 0 pour tout x ∈ Rd .
(A3) Il existe Λ ≥ 1 tel que pour tout x ∈ Rd , ξ ∈ Rd , Λ−1 |ξ|2 ≤ hσ(x)ξ, ξi ≤ Λ|ξ|2 .
Les hypothèses (A1-q) et (A3) sont des hypothèses classiques de non-dégénérescence, de la mesure sphérique
λ et du coefficient σ.
La condition de dérive nulle dans (A2-q) provient du fait que pour α ∈ (0, 1] l’ajout d’une dérive d’ordre t
ne correspond pas en temps petit à un terme négligeable par rapport à l’échelle caractéristique en t1/α . Nous
renvoyons à l’Appendice B de [9] pour des détails.

Théorème 3.3.1 Supposons q > d+ 4 et que (A1-q), (A2-q)  vérifiées. Soit 0 < M ≤ q − (d+ 4).
 et (A3) sont
1
Il existe une fonction RM (T, x, y), où |RM (T, x, y)| ≤ CM (T ) 1+|y−x|d+α := ρα,M (T, y−x) pour une constante
positive CM (T ), telle que
M−1
X hl
(p − ph )(T, x, y) = π h (T, x, y) + hM RM (T, x, y) ,
(l + 1)! l
l=1

PM−1
où l=1 πlh (T, x, y) ≤ ρα,M (T, y − x), i.e. les contrôles sur les (πlh )l∈[[1,M−1]] sont uniformes en h.

28
 Les termes dans le développement précédent dépendent de h. Il est néanmoins possible de faire des dévelop-
pements supplémentaires des quantités πlh de sorte à obtenir des termes indépendants de h dans l’expression
ci-dessus, ce que l’uniformité en h du contrôle laisse d’ailleurs à penser.
Le théorème précédent permet également de relâcher l’hypothèse de bornitude de la fonction test de [Hau02]
et d’obtenir un développement de l’erreur faible EhD := Ex [g(XTh )] − Ex [g(XT )] pour des fonctions g mesurables
vérifiant la condition de croissance ∃C > 0, ∀x ∈ Rd , |g(x)| ≤ C(1 + |x|β ), β < α, qui est optimale pour le cas
considéré. En particulier aucune hypothèse de régularité sur g n’est ici requise, contrairement aux approches
utilisées dans [TT90], [PT97]. Rappelons que d’un point de vue numérique, le développement de l’erreur faible
permet d’utiliser la méthode d’extrapolation de Romberg Richardson pour accélérer la convergence, cf. Chapitre
1.3.1.
Enfin, le développement du Théorème 3.3.1 est utile pour étudier la sensibilité de l’erreur faible EhD par
rapport à x sans conditions supplémentaires sur g. Ce point est fondamental pour les applications en finance
(couverture). Nous renvoyons à Guyon [Guy06] pour des détails das le cas diffusif.

Extensions et perspectives
Nous avons dans les articles [8],[9] utilisé la méthode parametrix de McKean et Singer qui dans des cas non
dégénérés ou faiblement dégénérés est très robuste et se prête naturellement à la discrétisation.
Pour l’instant les théorèmes limites locaux sur les densités ont été obtenus à partir des développements
parametrix (3.6), (3.13) et du théorème limite local en régime Gaussien. Il serait intéressant d’étendre ce type
de résultats à des processus dont la densité n’a pas de comportement de type Gaussien. A titre d’exemple, si
l’on remplace en Section 3.3 les accroissements stables par des variables i.i.d. renormalisées qui sont dans le
bassin d’attraction de la loi stable initiale, quel type de théorème limite local peut-on espérer établir ? Quelle
serait la vitesse en fonction de l’indice α du stable ?
Nous avons également jusqu’à présent toujours énoncé des résultats à horizon T > 0 fini. Dans le cas où
le processus admet une mesure invariante la question de l’existence du théorème limite local et de la vitesse
associée pour une estimation par chaı̂ne de Markov, éventuellement à pas décroissant, reste ouverte.

29
 Chapitre 4

Approximation d’Equations
Différentielles Rétrogrades (EDSR)

Nous présentons dans ce chapitre les travaux [2], [5] (en collaboration avec F. Delarue de l’Université
de Nice), [6] (en collaboration avec B. Bouchard de l’Université Paris Dauphine) relatifs à l’approximation
d’équations différentielles stochastiques (progressives) rétrogrades (EDS(P)R). Introduites par Pardoux et Peng
[PP90], les EDSR ont permis d’étendre les réprésentations de Feynman-Kac à des EDP paraboliques ou ellip-
tiques semi ou quasi-linéaires. Ce sont par ailleurs des objets qui apparaissent de façon naturelle en finance, cf.
El Karoui et al. [KPQ97], dans la mesure où les problèmes de réplication d’options se formulent directement en
terme d’EDSR.
C’est un sujet qui a fait l’objet d’une active recherche tant du point de vue théorique que numérique.
Nous présenterons en Section 4.1 un bref aperçu de la théorie des EDSR et de leurs liens avec les EDP non
linéaires. En terme numérique l’approximation d’EDSR fournit un procédé alternatif aux méthodes déterministes
de différences/éléments finis pour les problèmes non linéaires. Nous verrons ci-après que la complexité des
algorithmes proposés dépend de façon hautement non linéaire de la dimension du problème sous-jacent. On
ne retrouve donc pas la complexité linéaire en la dimension spécifique des méthodes de Monte Carlo pour
l’approximation probabiliste de solutions d’EDP linéaires et l’on se heurte à la “fléau de la dimension” bien
connue dans les approximations numériques déterministes d’EDP.
Cette compléxité provient de l’estimation d’espérances conditionnelles, qui est une étape inévitable pour
l’approximation d’EDSR. Trois grands types de techniques ont été introduites à cet effet :
- Les méthodes de Monte Carlo associées à des estimateurs dérivant du calcul de Malliavin. Cette approche a
été développée par Bouchard et Touzi [BT04] pour des EDSR faiblement couplées.
- Les méthodes de régression linéaire, utilisées par Gobet et al. [GLW05] toujours dans le cadre des EDSR
faiblement couplées.
- Les méthodes de quantification, introduites par Bally et Pagès dans le cadre des EDSR réfléchies, problème
d’arrêt optimal pour les options américaines [BP03]. Le premier algorithme de discrétisation d’EDSR proposé
par Chevance [Che97] dans un cadre plus restrictif utilise pour l’estimation des espérances conditionnelles une
méthode qui s’apparente à cette approche. C’est par ailleurs celle que nous avons utilisée dans [2] et [5].
Nous présenterons brièvement en Section 4.2 les algorithmes de [BT04] et [GLW05]. En retenant l’approche
par quantification pour l’approximation des espérances conditionnelles, nous avons pu dans [2], [5] proposer
des schémas de discrétisation d’EDSR fortement couplées associées à des EDP quasi-linéaires. Nous détaillerons
ces algorithmes en Section 4.3. Enfin, nous développerons en Section 4.4, les résultats de [6] où nous avons
analysé l’erreur forte d’un schéma de semi-discrétisation, i.e. on ne propose pas de schéma d’approximation des
espérances conditionnelles, pour une EDSR dans un domaine.

4.1 Quelques résultats sur les EDSR et lien avec les EDP
4.1.1 Introduction à la solvabilité des EDSR et couplage faible
La dynamique d’une EDSR s’écrit de façon générique
Z T Z T
Yt = ξ + f (s, Ys , Zs )ds − Zs dWs , t ∈ [0, T ], (4.1)
t t

31
où (Wt )t≥0 est un mouvement Brownien standard d-dimensionnel sur un espace de probabilité (Ω, F , P) de
filtration naturelle (Ft )t≥0 et T > 0 est un instant arbitraire. La condition terminale ξ est FT mesurable et
généralement supposée de carré intégrable. Le coefficient f est appelé générateur de l’EDSR. Ce dernier peut être
aléatoire à condition de supposer que le processus (f (t, y, z))t≥0 est pour tout (y, z) progressivement mesurable.
Un exemple typique de dépendance aléatoire du générateur est fourni par le couplage faible avec une EDS
progressive, i.e. f s’écrit f (s, Xs , Ys , Zs ), ξ = g(XT ) où X est la solution d’une EDS progressive. L’EDSR
faiblement couplée s’écrit alors :
Z T Z T
Yt = g(XT ) + f (s, Xs , Ys , Zs )ds − Zs dWs ,
t t
Z t Z t
Xt = x + b(s, Xs )ds + σ(s, Xs )dWs , t ∈ [0, T ]. (4.2)
0 0

Ce sont les EDSR de type (4.2) qui permettent d’établir une connexion avec les EDP semi-linéaires.

Principes de résolution
Résoudre une EDSR de type (4.1) consiste à trouver un couple (Yt , Zt )t∈[0,T ] progressivement mesurable
tel que (4.1) soit vérifiée pour tout t ∈ [0, T ], P-presque sûrement. L’hypothèse de progressive mesurabilité est
essentielle. Le terme de droite dans (4.1) est a priori FT mesurable là où l’on cherche à ce que celui de gauche
soit Ft mesurable. C’est le rôle du processus (Zt )t∈[0,T ] que de ramener condition terminale et générateur dans
la bonne tribu.
Ce point s’illustre directement dans le cas d’un générateur nul. En effet, si f = 0 et que l’on considère une
condition terminale ξ de carré intégrable, le théorème de représentation des martingales donne l’existence d’un
RT
unique processus (Zt )t∈[0,T ] progressivement mesurable et de carré intégrable tel que ξ = E[ξ] + 0 Zs dWs .
Rt
Ainsi Yt := E[ξ|Ft ] = E[ξ] + 0 Zs dWs est bien solution de (4.1). L’unicité découle de celle du processus Z de
représentation.
En bref, c’est la condition de mesurabilité sur le couple (Yt , Zt )t∈[0,T ] qui permet de “fermer” l’équation
(4.1) à deux inconnues : à la relation différentielle on ajoute la condition de mesurabilité.
Sous des conditions de structure du générateur f , l’équation (4.1) est alors bien posée. Lorsque f est Lipschit-
RT
zien en (y, z), on peut construire une suite (Y n , Z n )n∈N définie par (Y 0 , Z 0 ) = 0, Ytn+1 = ξ+ t f (s, Ysn , Zsn )ds−
R T n+1
t Zs dWs , t ∈ [0, T ], n ≥ 0. La condition de Lipschitz continuité de f permet d’identifier le processus
n+1
RT
Z au processus de représentation associé à la variable aléatoire ξ + t f (s, Ysn , Zsn )ds qui est bien de carré
intégrable. Les arguments usuels d’analyse différentielle, lemme de Gronwall, et stochastique, inégalités BDG,
permettent de déduire que cette suite admet un unique point fixe solution forte de (4.1). La même démarche
reste possible si l’on remplace l’hypothèse de Lipschitz continuité du générateur en y par une condition de mo-
notonie associée à une croissance au plus linéaire. Ces résultats de solvabilité sont initialement dûs à Pardoux et
Peng [PP90]. Par la suite de nombreux auteurs ont affaibli ces hypothèses pour la solvabilité de (4.1), indiquons
Briand et al. [BDH+ 03] pour le cas d’une condition terminale Lp , p ∈ (1, 2) ou Kobylanski [Kob00] pour le cas
d’un générateur à croissance quadratique en Z, hypothèse par ailleurs naturelle du point de vue des EDP, cf.
Ladyzhenskaya et al. [LSU68].

Lien avec les EDP semi-linéaires

Nous considérerons pour l’illustration du lien entre EDSR et EDP le cas parabolique, le cas elliptique a été
abordé par Darling et Pardoux [DP97].

De l’EDP vers l’EDSR. Comme pour la formule de Feynman-Kac dans le cadre linéaire c’est l’application
de la formule d’Itô qui permet de relier EDP et EDSR. Considérons l’EDP parabolique semi-linéaire
(
(∂t + Lt )u(t, x) + f (t, x, u(t, x), ((σ∇u)(t, x)) = 0, (t, x) ∈ [0, T ) × Rd ,
(4.3)
u(T, x) = g(x), x ∈ Rd ,


où Lt ϕ(x) = hb(t, x), ∇ϕ(x)i+ 12 tr a(t, x)Dx2 ϕ , a(t, x) = σσ ∗ (t, x) est le générateur infinitésimal d’une équation
Rt Rt
de diffusion Xt = x + 0 b(s, Xs )ds + 0 σ(s, Xs )dWs soluble au moins au sens faible. S’il existe une solution
“classique” à (4.3), i.e. bornée et à dérivées d’ordres un en temps et deux en espace bornées, il est possible
d’appliquer sur l’espace de probabilité sous-jacent la formule d’Itô au processus Yt := u(t, Xt ). On identifie alors
Zt = (∇uσ)(t, Xt ) qui permet de retrouver la formulation couplée de l’équation (4.2).

32
 Ainsi, la donnée d’une solution régulière à l’EDP a permis de construire par un argument de vérification
une solution à l’EDSR. Cette méthode peut s’étendre à des solutions de (4.3) ayant des dérivées temps espace
au sens de Sobolev à l’aide de la formule d’Itô-Krylov, cf. [Kry96].

De l’EDSR vers l’EDP. Réciproquement, une question naturelle consiste à essayer de retrouver l’EDP à
partir de la solution de l’EDSR. Supposons que l’EDSR (4.2) admette une solution unique et définissons u :
(t, x) ∈ [0, T ] × Rd 7→ u(t, x) := Ytt,x qui désigne la valeur de l’EDSR (4.2) en t pour la condition initiale Xt = x.
De cette définition de u il est aisé de voir que u(0, x) s’écrit comme la solution en 0 d’une EDSR sur [0, t] avec
condition terminale u(t, Xt0,x ). Ce principe de programmation dynamique satisfait par u permet d’établir qu’elle
est solution de viscosité de (4.3). Nous renvoyons à l’article de Crandall et al [CIL92] pour la définition.
Lorsque les coefficients de l’EDSR sont réguliers il est possible de la dériver par rapport à la condition
initiale à l’aide de techniques usuelles de flots stochastiques sur les composantes progressives et rétrogrades,
voir Kunita [Kun90] ou Ikeda et Watanabe [IW89]. Cette opération permet d’obtenir la différentiabilité de u à
l’ordre un en temps et deux en espace. Le calcul de Malliavin permet ensuite l’identification de (Zt )t∈[0,T ] avec
((σ∇u)(t, Xt ))t∈[0,T ] dont on déduit ensuite grâce à (4.2) que u est solution classique de (4.3). Cette démarche
est celle de l’article de Pardoux et Peng [PP92]. Notons que dans cette approche aucune condition de non
dégénérescence n’est a priori requise, seule la régularité des coefficients intervient. Si l’équation de départ est
en revanche non dégénérée, cette approche n’exploite pas l’effet régularisant du noyau de la chaleur.

4.1.2 EDSR fortement couplées

De nombreuses applications, équations des milieux poreux, problème de gros investisseur en finance, font
intervenir des équations aux dérivées partielles pour lesquelles la solution et/ou son gradient apparaissent dans
l’opérateur différentiel. On parle alors d’EDP quasi-linéaires. De façon générale celles-ci s’écrivent :


1 2
∂t u(t, x) + hb (t, x, u(t, x), (σ∇u)(t, x)) , ∇u(t, x)i + 2 tr a (t, u(t, x)) Dx u(t, x)

+f (t, x, u(t, x), (σ∇u)(t, x)) = 0, (t, x) ∈ [0, T ) × Rd , (4.4)

 d
u(T, x) = g(x), x ∈ R .

Ce sont des équations de ce type pour lesquelles nous avons proposé dans [2], [5] des algorithmes de
discrétisation que nous détaillerons en Section 4.3. En gardant à l’esprit les relations Yt = u(t, Xt ), Zt =
(σ∇u)(t, Xt ) qui lient l’EDSR à l’EDP il apparaı̂t naturel que l’objet probabiliste naturellement associé à (4.4)
sera une EDSR “fortement couplée” de dynamique
Z t Z t
Xt = x + b(s, Xs , Ys , Zs )ds + σ(s, Xs , Ys )dWs
0 0
Z T Z T
Yt = g(XT ) + f (s, Xs , Ys , Zs )ds − Zs dWs . (4.5)
t t

Le couplage est dit fort dans la mesure où les coefficients de la composante progressive dépendent de
l’équation rétrograde. Il n’est ainsi plus possible de résoudre l’équation progressive de façon autonome ce qui
induit une difficulté supplémentaire par rapport à la Section précédente.

Solvabilité en temps court. Le cadre naturel de résolution de (4.5) est celui de coefficients Lipschitziens. On
cherche alors à construire un triplet (Xt , Yt , Zt )t∈[0,T ] adapté tel que E[supt∈[0,T ] |Xt |2 ] + E[supt∈[0,T ] |Yt |2 ] <
RT
+∞, et E[ 0 |Zs |2 ds] < +∞ avec Z progressivement mesurable. L’idée naturelle pour obtenir une telle so-
lution consiste à construire, en suivant la démarche du travail de [PP90], une contraction. Les espaces fonc-
tionnels sont ceux naturellement associés aux conditions ci-avant. Si l’on se donne X0 on définit pour n ≥
RT RT Rt
0, Ytn+1 = g(XTn ) + t f (s, Xsn , Ysn+1 , Zsn+1 )ds − t Zsn+1 dWs , Xtn+1 = x + 0 b(s, Xsn+1 , Ysn+1 , Zsn+1 )ds +
Rt
0
σ(s, Xsn+1 , Ysn+1 )dWs . Sous les hypothèses de Lipschitz continuité uniforme des coefficients b, σ, f, g, Anto-
nelli [Ant93] a montré que le procédé précédent définit une contraction et que son unique point fixe est l’unique
solution forte de (4.5) si T ≤ cK où K désigne la constante de Lipschitz des coefficients. Sans autres hy-
pothèses ce résultat est optimal dans la mesure où il existe des contre-exemples. Considérons le cas particulier
déterministe X0 = x, YT = −XT , dXt = Yt dt, dYt = −Xt dt, t ∈ [0, T ], pour lequel Z = 0. La forme générale
des solutions s’écrit Xt = A cos(t) + B sin(t), Yt = −A sin(t) + B cos(t), t ∈ [0, T ] ce qui impose A = x et
A cos(T ) + B sin(T ) = A sin(T ) − B cos(T ) pour vérifier les conditions de bord. Pour T = 3π/4 cette condition
conduit à A = x et A = 0, ce qui met la solvabilité en défaut.

33
 Schéma en quatre temps. Cette approche est due à Ma, Protter et Yong [MPY94]. L’idée sous-jacente est
de “décorréler” les composantes progressives et rétrogrades de (4.5) en utilisant des propriétés de la solution
supposée connue de l’EDP (4.4) sous-jacente. Précisément, si u solution de (4.4) est une solution “classique”,
l’identification formelle Yt = u(t, Xt ), Zt = (σ∇u)(t, Xt ) précédente conduit à considérer l’équation progressive
autonome
Z t Z t
Xt = x + b(s, Xs , u(s, Xs ), (σ∇u)(s, Xs ))ds + σ(s, Xs , u(s, Xs ))dWs .
0 0

Si les coefficients, b, σ ainsi que u et son gradient sont Lipschitziens, alors l’équation ci-dessus est fortement
soluble. Sous des hypothèses usuelles de croissance du générateur f , l’application de la formule d’Itô donne que
le triplet (Xt , Yt , Zt )t∈[0,T ] = (Xt , u(t, Xt ), (σ∇u)(t, Xt ))t∈[0,T ] est bien solution de (4.5).
C’est encore là une méthode de type vérification. On construit l’objet probabiliste à partir des propriétés
de l’EDP. Pour l’unicité, partant d’une solution arbitraire (X̄, Ȳ , Z̄) de (4.5) il faut l’identifier avec (X, Y, Z)
construite précédemment. L’idée est ici d’utiliser de nouveau la fonction u solution de l’EDP (supposée bornée
ainsi que ses dérivées d’ordre un et deux en espace) pour identifier Ȳt = u(t, X̄t ), Z̄t = (σ∇u)(t, X̄t ). Pour cela,
on applique la formule d’Itô à (u(t, X̄t ) − Ȳt )2 , t ∈ [0, T ]. Les outils usuels en équations rétrogrades : inégalités
BDG, inégalités de Young et Lemme de Gronwall, permettent de conclure.
Une fois effectuée l’identification pour (Ȳ , Z̄), les hypothèses précédentes de Lipschitzianité sur les coefficients
donnent que X = X̄ car les deux composantes progressives vérifient la même EDS. Dans le cadre d’un générateur
Lipschitzien l’identification (Ȳ , Z̄) = (Y, Z) provient ensuite de l’unique solubilité de l’équation rétrograde.
Pour pouvoir appliquer ce programme il faut donc disposer de coefficients b, σ, f Lipschitziens et d’une
solution classique de (4.4) bornée ainsi que ses dérivées d’ordre un et deux en espace. Un jeu typique d’hypothèses
qui garantit ces propriétés est le suivant :
(a) Les coefficients b, f, σ et g sont réguliers et à dérivées d’ordre un en (x, y, z) bornées.
(b) ∃Λ > 0, |(b, f )(t, x, y, z)| ≤ Λ(1 + |y| + |z|), |g(x)| ≤ Λ.
(c) ∃λ > 0, λId ≤ a(t, x, y) ≤ ΛId .
(d) g ∈ C 2 (Rd ) a des dérivées d’ordre un et deux bornées et uniformément η-Hölderiennes, η ∈ (0, 1].
Théorème 4.1.1 (cf. Ladyzhenskaya et al. [LSU68]) Supposons que (a),(b),(c),(d) soient satisfaites. Alors,
l’équation admet une unique solution classique bornée ainsi que ses dérivées d’ordre un et deux en espace.
C’est sous des hypothèses de ce type que nous nous placerons en Section 4.3. Nous reprendrons pour étudier
la convergence de schémas d’approximation d’EDP/EDSR la démarche présentée ci-dessus pour l’unicité en
exploitant les propriétés de u.

Solvabilité en temps arbitraire dans le cadre d’hypothèses “usuelles”.

Sous l’hypothèse d’uniforme ellipticité (c) et de coefficients Lipschitziens, l’unique solvabilité de (4.5) pour
T > 0 arbitraire a été établie par Delarue [Del02]. Pour passer du temps court au temps long, l’idée consiste à
itérer le résultat en temps court de [Ant93] en contrôlant que la solution de viscosité de (4.4) associée à l’EDSR
ainsi construite, u(t, x) := Ytt,x , reste Lipschitzienne en espace (au sens où sa constante de Lipschitz n’explose
pas) sur chacun des intervalles de temps “courts”. Ceci permet en effet d’itérer la procédure en considérant
comme nouvelle condition terminale la solution précédemment obtenue et Lipschitzienne en espace. Celle-ci sert
alors d’entrée à une nouvelle EDSR fortement couplée fortement uniquement soluble en temps court.
L’hypothèse de non-dégénérescence est ici cruciale pour pouvoir construire une solution pour T arbitraire,
comme l’on pouvait l’envisager à partir du contre-exemple ci-avant.
Pour la solvabilité au sens faible d’EDSR fortement couplées dans un cadre non dégénéré nous renvoyons au
travail de Delarue et Guatteri [DG06].

4.2 Algorithmes de discrétisation d’EDSR faiblement couplées

Au regard des parties précédentes le cas faiblement couplé apparaı̂t plus favorable à la discrétisation dans
la mesure où l’approximation de la composante progressive peut être effectuée de façon autonome. Par ailleurs
les méthodes numériques pour celle-ci sont bien comprises, voir Chapitre 1. Nous présentons tout d’abord
un premier schéma de semi-discrétisation pour l’EDSR (4.2). On entend par là un schéma théorique qui ne
propose pas encore d’approximation des espérances conditionnelles en présence. Nous présenterons ensuite deux
techniques numériques pour l’approximation des espérances conditionnelles. Nous analyserons précisément leur
complexité algorithmique.

34
4.2.1 Un premier schéma de semi-discrétisation
Soit {(ti := ih)i∈[[0,N ]]}, N h = T , une grille en temps de pas h donnée. Pour discrétiser (4.2), l’idée naturelle
consiste à partir de l’écriture “locale” en temps :
Z ti+1
Yti = E[Yti+1 + f (s, Xs , Ys , Zs )|Fti ], i ∈ [[0, N − 1]]. (4.6)
ti

Pour i = N − 1, Yti+1 = YT = g(XT ), condition terminale connue. Pour déduire un schéma d’approximation de
l’écriture (4.6), si l’on pense par exemple à une méthode des rectangles pour la discrétisation de l’intégrale, il est
fondamental d’avoir une approximation du processus de représentation. A ce propos, partant de la formulation
(4.2) entre ti et ti+1 , i ∈ [[0, N − 1]], si l’on multiplie les deux membres de la composante rétrograde par
l’accroissement Brownien Wti+1 − Wti et que l’on passe à l’espérance conditionnelle en Fti il vient :
Z ti+1 Z ti+1
E[ Zs ds|Fti ] = E[(Yti+1 − Yti )(Wti+1 − Wti )|Fti ] + E[ f (s, Xs , Ys , Zs )ds(Wti+1 − Wti )|Fti ]
ti ti

= E[Yti+1 (Wti+1 − Wti )|Fti ] + O(h3/2 ), (4.7)

sous les hypothèses typiques de solvabilité de EDSR. Par ailleurs il est naturel de chercher une approximation
de (Zs )s∈[ti ,ti+1 ] qui soit Fti mesurable. Si l’on néglige la contribution associée au générateur dans (4.7), la
meilleure approximation Fti mesurable de (Zs )s∈[ti ,ti+1 ] dans L2 ([ti , ti+1 ] × Ω, dt ⊗ dP) est donnée par :

Ẑti = h−1 E[Yti+1 (Wti+1 − Wti )|Fti ].

Cette technique d’“incrément de martingale” pour l’estimation du Z apparaı̂t initialement chez Bally et Pagès
[BP03]. Supposons que l’on dispose désormais d’une approximation (Xthi )i∈[[0,N ]] de la composante progressive
de (4.2), par exemple le schéma d’Euler. L’itération dérivant des écritures précédentes s’écrit :

Ythi = E[Ythi+1 |Fti ] + hE[f (ti , Xthi , Ythi+1 , Ẑthi )|Fti ],

Ẑthi = h−1 E[Ythi+1 (Wti+1 − Wti )|Fti ], i ∈ [[0, N − 1]]. (4.8)

La vitesse de convergence de (4.8) a été obtenue pas Zhang [Zha04] sous des hypothèses usuelles de solvabilité.
Précisément, lorsque X h est le schéma d’Euler, il établit qu’il existe C dépendant de T et des constantes de
Lipschitz des coefficients dans (4.2) tel que :
N
X −1 Z ti+1
max E[|Ythi − Yti |2 ] + E[ |Ẑthi − Zs |2 ds] ≤ Ch.
i∈[[0,N ]] ti
i=0

On retrouve que l’erreur obtenue sur la composante rétrograde est du même ordre que celle de l’équation
progressive, voir également Gobet et Labart [GL07] pour un développement d’erreur qui met en évidence cette
propriété.
Pour passer à partir de (4.8) à un schéma “implémentable” il est nécessaire de préciser une méthode d’ap-
proximation des espérances conditionnelles.

4.2.2 Méthodes de Monte Carlo et calcul de Malliavin

Cette approche a été utilisée par Bouchard et Touzi [BT04] à partir d’un schéma légèrement différent de
(4.8), i.e. l’itération sur la composante rétrograde est implicite Ythi = E[Ythi+1 |Fti ] + hf (ti , Xthi , Ythi , Ẑthi ). Sous les
hypothèses usuelles et pour h suffisamment petit, ce schéma est également bien défini.
Lorsque le coefficient de diffusion σ est uniformément elliptique, et que les coefficients b, σ sont réguliers, le
schéma d’Euler X h possède une densité. L’estimation des espérances conditionnelles se réduit ainsi à un calcul
de lois conditionnelles. Les auteurs utilisent ensuite des formules d’intégration par partie pour représenter la
densité et des fonctions du schéma d’Euler.
Soit une variable aléatoire R := ρ(Xthi+1 , ξ)A(Wti+1 − Wti ) où ρ et une fonction à valeurs réelles, A une
application affine et ξ une variable aléatoire indépendante de σ{(Xthi ), i ∈ [[1, N ]]}. Les auteurs établissent
l’expression
E[Hx (Xthi )ρ(Xthi+1 , ξ)S(A(Wti+1 − Wti ), x, Xthi )]
E[R|Xthi = x] := , (4.9)
E[Hx (Xthi )S(1, x, Xthi )]

35
 Qd
où Hx (y) = j=1 Ixj ≤yj est un produit de fonctions de Heaviside résultant de l’intégration par parties. L’ex-
pression de S fait intervenir des intégrales de Skorohod et un terme de localisation qui permet de contrôler la
variance pour l’estimation de Monte Carlo associée.
Ainsi, pour une trajectoire du schéma d’Euler, l’estimation par Monte Carlo de la quantité ci-avant va
nécessiter la simulation de M + 1 trajectoires i.i.d. du schéma d’Euler (en incluant la trajectoire initiale notée
X h,0 := X h ). Dans
 la pratique M = KN et à chaque instant (ti )i∈[[0,N −1]] de discrétisation K trajectoires
des (Xth,k
i
)i∈[[0,N ]] sont utilisées pour l’approximation de l’espérance conditionnelle. Ces K variables
k∈[[0,M]]
indépendantes interviennent implicitement dans l’expression de ξ de l’équation (4.9) lorsque R est associée à
l’approximation de l’EDSR. Précisément, si l’on note Ni := {(i − 1)K + 1, ..., iK}, i ∈ [[1, N ]] le schéma s’écrit :
∀l ∈ NN , YTh,l = g(XTh,l ), ∀j ∈ {0} ∪ Ni , i ∈ [[0, N − 1]],
P
1
K l∈Ni+1 HX h,j (Xth,l
i
)Ŷth,l S (Wtli+1 − Wtli , Xth,j
i+1 i i
, Xth,l
i
)
Ẑth,j
ti
i
= 1 P h,l h,j h,l
,
K l∈Ni+1 HX h,j (Xti )Si (1, Xti , Xti )
ti
P
1
K l∈Ni+1 HX h,j (Xth,l
i
)Ŷth,l S (1, Xth,j
i+1 i i
, Xth,l
i
)
Ŷth,j + hf (ti , Xth,j , Yth,j , Ẑth,j
t i
i
= 1 P h,l h,j h,l i i i
).
K l∈Ni+1 HX h,j (Xti )Si (1, Xti , Xti )
ti

La fonction A dans (4.9) est soit l’accroissement Brownien pour la mise à jour du processus de représentation
approché, soit 1 pour la composante rétrograde. L’analyse de l’erreur associée, cf. Théorème 6.2 de [BT04] donne
que pour retrouver l’erreur attendue d’ordre h1/2 en norme L2 , qui était celle du schéma de semi-discrétisation,
il faut prendre K ≥ N 6+d/2 , ce qui conduit à considérer N 7+d/2 trajectoires du schéma d’Euler.
La complexité apparente d’une telle procédure est de l’ordre de K 2 N . Néanmoins les fonctions de Heaviside
intervenant dans l’estimateur permettent de passer à K log(K)(d−1)∨1 N en triant les réalisations. Ce type
d’estimateur reste délicat à calibrer dans la pratique, en particulier le paramètre de troncation implicite dans
S. Enfin, il faut également approcher les intégrales de Skorohod dans ce même terme.

4.2.3 Techniques de régression

Une autre approche est basée sur des méthodes de régression. Elle a été développée par Gobet, Lemor
et Warin [GLW05, GLW06] selon une démarche voisine de celle proposée par Longstaff  et Schwarz [LS01]
dans le cadre des options américaines. L’idée sous-jacente est d’utiliser M trajectoires (Xth,j
i
)i∈[[0,N ]]
j∈[[1,M]]
indépendantes du schéma d’Euler et de “régresser” les approximations naturelles associées à l’EDSR, voir (4.8),
sur une base de fonctions de dimension finie. La phase d’initialisation de l’algorithme reste la même et consiste
à poser YTh,j = g(XTh,j ), j ∈ [[1, M ]]. Pour X à valeurs dans Rd , si l’on se donne (pl )l∈[[1,d]] bases de fonctions
de cardinaux respectifs (Kl )l∈[[1,d]], la première étape pour la mise à jour de l’approximation du Z à l’instant
(ti )i∈[[0,N −1]] consiste à résoudre le problème de moindre carré :
M
1 X −1 h,j 2
inf h Yti+1 (Wti+1 − Wti )l − hα, pl (Xth,j )i , l ∈ [[1, d]].
α∈RKl M j=1 i

Si α̂lti désigne le α optimal, on pose ensuite (Ẑth,j

i
)l = hα̂lti , pl (Xth,j
i
)i, l ∈ [[1, d]]. L’estimation de Yth,j
i
se fait
suivant le même schéma. Si
 
 1 X
M 2
α̂0ti = arginf α∈RK0 Yth,j + hf (ti , Xth,j , Yth,j , Ẑthi ) − hα, p0 (Xth,j )i ,
M i+1 i i+1 i

j=1

où p0 désigne une base de fonctions de cardinal K0 , on pose Yth,j i

= hα̂0ti , p0 (Xth,j
i
)i. L’approximation est ici
comme dans (4.8) explicite.
Nous renvoyons aux articles sus-mentionnés pour le détail de l’analyse de la convergence. La complexité
globale de l’algorithme est de l’ordre de N M d log(K) si l’on choisit K = Kl , l ∈ [[0, d]] et des fonctions
de base qui sont des indicatrices de cellules. Pour retrouver une erreur de l’ordre de h1/2 en norme L2 sous
de bonnes hypothèses de Lipschitz continuité des coefficients, l’équilibre des paramètres N, M, K conduit à
une complexité algorithmique l’ordre de N 2d+4 log(N ) à mettre en regard avec l’ordre de N 7+d/2 log(N ) de
l’algorithme précédent pour obtenir cette vitesse. En revanche, l’approche par régression ne requiert aucune
non dégénérescence des coefficients et peut s’étendre à des EDSR faiblement couplées qui font intervenir dans

36
 leur dynamique rétrograde les variations d’une martingale orthogonale au Brownien ainsi que d’éventuelles
composantes de sauts pour la partie progressive [GLW06].
Nous voyons sur ces exemples que même dans le cas faiblement couplé les complexités algorithmiques as-
sociées à ces modèles non-linéaires sont hautement explosives avec la dimension.

4.3 Discrétisation d’EDSR fortement couplées

La discrétisation d’une EDSR fortement couplée présente les mêmes difficultés que celles rencontrées pour la
construction des objets continus, cf. Section 4.1.2. Il n’est pas possible de séparer estimations des composantes
progressives et rétrogrades à moins de les “décorréler” complètement, c’est à dire de reproduire numériquement
la procédure de point fixe qui a permis de construire l’objet continu en résolvant une suite d’EDSR faiblement
couplées. Dans ce cadre, il devient possible numériquement d’utiliser les algorithmes de la section précédente.
Nous renvoyons pour cette approche aux travaux de Rivière [Riv05] et Bender et Zhang [BZ08].
Nous avons choisi pour les articles [2], [5] d’utiliser les liens forts entre l’EDSR et l’EDP quasi-linéaire et
d’exploiter que, sous les hypothèses du Théorème 4.1.1, la solution u de (4.4) est régulière. Cette approche
est également celle retenue par Milstein et Tretyakov [MT04]. L’idée commune à ces travaux est de procéder
à une itération rétrograde et à chaque étape en temps d’injecter dans la dynamique de l’équation progressive
les estimées de la solution et de son gradient obtenues à l’étape précédente. La régularité de la vraie solution
permet de contrôler la propagation de l’erreur due à cette approximation. De fait, cette procédure fournit un
algorithme probabiliste d’approximation de la solution de l’EDP en exploitant la dynamique locale en temps de
l’EDSR (4.5). L’approximation de l’EDSR elle même se fait dans un second temps après la résolution approchée
de l’EDP de façon similaire à Douglas et al. [DMP96] qui proposaient eux un algorithme déterministe pour (4.4)
avant d’injecter les résultats de l’approximation dans l’EDSR.

4.3.1 Etapes de la discrétisation

Posons pour la suite pour tout (t, x) ∈ [0, T ] × Rd, v(t, x) := σ∇u (t, x). Pour t ∈ [0, T − h], l’équation (4.5)
et le lien avec u solution de (4.4) donnent le principe de programmation dynamique :
Z t+h
u(t, Xt ) := E[u(t + h, Xt+h )|Ft ] + E[ f (s, Xs , Ys , Zs )ds|Ft ].
t

Cette écriture ne peut pas directement s’exploiter numériquement dans la mesure où la dynamique de X entre t
et t+h fait intervenir (Y
s , Zs )s∈[t,t+h] soit la solution

u de (4.4) et son gradient v. Si l’on dispose d’approximations
ũ(t + h, .), ṽ(t + h, .) de u(t + h, .), v(t + h, .) sur tout l’espace, on définit une transition “approchée” entre
t et t + h :

∀x ∈ Rd , T (t, x) := b(t, x, ũ(t + h, x), ṽ(t + h, x))h + σ(t, x, ũ(t + h, x))(Wt+h − Wt ).

On utilise ensuite cette transition dans la dynamique du processus progressif pour mettre à jour l’estimation
du gradient par incrément de martingales, i.e. pour tout x ∈ Rd ,

ṽ(t, x) := h−1 E[ũ(t, x + T (t, x))(Wt+h − Wt )]. (4.10)

On définit enfin la nouvelle approximation ũ(t, .) par

ũ(t, x) := E[ũ(t + h, x + T (t, x))] + hf (t, x, ũ(t + h, x), ṽ(t, x)). (4.11)

Partant de la condition initiale ũ(T, .) = g(.), ṽ(T, .) = (σ∇g)(T, .), ce schéma est bien défini le long d’une grille
de discrétisation en temps {(ti := ih)i∈[[0,N ]]}.
Deux obstacles subsistent toutefois pour sa mise en oeuvre :
(a) Le schéma ne peut être défini en pratique sur tout l’espace.
(b) Les espérances doivent être estimées numériquement.
Grilles spatiales et prolongement. Pour le premier point nous nous sommes restreints à une grille spatiale ce
qui nous rapproche des algorithmes déterministes. Dans le cas découplé, des simulations autonomes permettent
d’avoir une idée a priori de la “géométrie” associée à la composante progressive, i.e. zones plus visitées que
d’autres. Le cas couplé n’autorise pas une telle procédure et sans autre connaissance a priori nous avons retenu
des grilles Cartésiennes tronquées. C’est également le choix de Milstein et Tretyakov.

37
 La restriction à une grille spatiale amène à préciser le sens des relations (4.10),(4.11) pour t = ti , i ∈
[[0, N − 1]]. En effet, si l’on se restreint en ti+1 à l’estimation de u, v sur une grille spatiale Ci+1 , la transition
x + T (ti , x), x ∈ Ci n’appartiendra pas à Ci+1 . Il est alors nécessaire de comprendre ũ(ti+1 , .) dans (4.10), (4.11)
comme un prolongement à tout l’espace de la fonction définie aux points d’une grille spatiale. Dans [2] nous
avons considéré un prolongement au plus proche voisin, i.e. ∀(x, i) ∈ Rd × [[0, N ]], ũ(ti , x) = ũ(ti , ΠCi (x)) où ΠCi
désigne la projection sur la grille Cartésienne Ci (avec une convention pour le centre des mailles).
Dans [5] nous avons en revanche utilisé une interpolation multilinéaire qui peut par ailleurs s’interpréter de
façon naturelle en termes probabilistes. Cette approche a permis de tirer parti de la régularité a priori de u
sous les hypothèses du Théorème 4.1.1 et d’améliorer la vitesse de convergence de notre approximation.

Quantification. Pour le point (b), nous avons choisi d’utiliser une méthode de quantification plutôt que
l’approximation usuelle de Monte Carlo. L’idée de la quantification de variables aléatoires est d’approcher une
variable donnée par une loi discrète à support fini. Cette procédure permet de ramener le calcul d’espérances à
celui d’une somme finie. Nous renvoyons à l’ouvrage de Graf et Luschgy [GL00] pour une référence théorique à
cette technique.
Des équations ci-dessus nous devons quantifier l’accroissement Brownien soit à un changement d’échelle
près la loi normale. Une première approximation consisterait à l’approcher par un produit de lois de Bernoulli
symétriques correctement normalisées, c’est le choix effectué par Milstein et Tretyakov [MT04]. Nous avons
choisi d’utiliser une quantification “optimale” de la loi gaussienne. Si nous notons G la loi gaussienne standard
de Rd et Ĝ une quantifiée optimale ĜM à M points (non nécessairement unique) vérifie (pour M suffisamment
grand) l’inégalité de “distortion” :

∃C := C(G, d), kG − ĜM kL2 ≤ CM −1/d . (4.12)

Les quantifieurs optimaux vérifient également la propriété projective (ou de stationnarité) E[G|ĜM ] = ĜM dont
on déduit que E[ĜM ] = 0 et que pour une fonction f de classe C 2 à dérivées bornées,

E[f (G) − f (ĜM )] = E[h∇f (ĜM ), E[G − ĜM |ĜM ]i]

| {z }
R1 =0
(4.13)
+ 0 dλ(1 − λ)E[hD2 f (ĜM + λ(G − ĜM ))(G − ĜM ), G − ĜM i]
2
:= RM , |RM | ≤ |D 2f |∞ E[|G − ĜM |2 ] ≤ CM −2/d ,

en exploitant la propriété projective pour la première égalité et (4.12) pour la dernière inégalité.
Ainsi, contrairement à la méthode de Monte Carlo, la quantification permet d’exploiter la régularité des
fonctions dont on veut estimer l’espérance. Le comportement de cette méthode se rapproche des techniques de
quadrature ce qui se retrouve sur la vitesse de convergence associée qui dépend fortement de la dimension d
de l’espace. Intuitivement, pour M points, la vitesse de M −1/d provient d’une répartition du même ordre dans
chaque direction. Nous avons tiré tout l’avantage de cette approche dans [5].

4.3.2 Algorithmes
Les deux algorithmes proposés dans [2] et [5] dérivent de (4.10), (4.11) avec les procédures d’extension
précédentes dues aux restrictions à des grilles spatiales. Dans [2], on remplace simplement l’accroissement
Brownien par la quantifiée optimale gaussienne renormalisée en considérant la projection au plus proche voisin.
Il vient

ū(T, x) = g(x), x ∈ CN ,
T (ti , x) := b(ti , x, ū(ti+1 , x), v̄(ti+1 , x))h + σ(ti , x, ū(ti+1 , x))h1/2 ĜM ,
v̄(ti , x) = h−1 E[ū(ti+1 , x + T (ti , x))h1/2 ĜM ],
ū(ti , x) = E[ū(ti , x + T (ti , x))] + hf (ti , x, ū(ti+1 , x), v̄(ti , x)), x ∈ Ci , (4.14)

pour tout i ∈ [[0, N − 1]].

Dans [5], pour l’analyse de l’erreur, il est nécessaire que la matrice de covariance de la variable aléatoire qui
approche la loi normale soit l’identité. Or ce n’est pas forcément le cas pour une quantification ĜM optimale de
G. Pour M assez grand, ΓM := E[ĜM Ĝ∗M ] est inversible et la technique naturelle consiste à préfixer l’incrément
h1/2 ĜM par γM∗
où γM est la matrice triangulaire inférieure de la factorisation de Cholesky de Γ−1 M = γM γM .
∗

38
 On a ainsi l’algorithme

ū(T, x) = g(x), x ∈ CN , (4.15)

∗ 1/2
T (ti , x) := hb(ti , x, ū(ti+1 , x), v̄(ti+1 , x))h + σ(ti , x, ū(ti+1 , x))γM h ĜM ,
v̄(ti , x) = h−1 E[ū(ti+1 , x + T (ti , x))h1/2 γM ∗
ĜM ],
ū(ti , x) = E[ū(ti , x + T (ti , x))] + hf (ti , x, ū(ti+1 , x), v̄(ti , x)), x ∈ Ci , i ∈ [[0, N − 1]],

en considérant l’interpolation multilinéaire comme procédure de prolongement à tout l’espace.

Indiquons que plus que l’aspect optimal c’est l’aspect stationnaire qui est crucial pour l’analyse de conver-
gence dans [5]. Or, il est possible de construire de telles grilles à partir de lois de Bernoulli. Ce point de vue est
celui adopté par Milstein et Tretyakov dans le livre [MT04] et les articles [MT06b, MT07]. Nous verrons par
ailleurs au paragraphe suivant que lorsque l’on exploite la stationnarité, pour une solution régulière de (4.4)
l’erreur de quantification est négligeable devant les autres composantes de l’erreur.

4.3.3 Théorèmes de convergence

Nous présentons les résultats de convergence lorsque le coefficient b dans les algorithmes précédents n’ont pas
de dépendance en gradient. Des solutions alternatives pour traiter ce cas ont été proposées dans [2] et [5]. Dans
les deux cas elles font intervenir une transition intermédiaire additionnelle par rapport à T et induisent une
dépendance déraisonnable entre M et h pour pouvoir analyser la convergence de cette méthode. Par ailleurs,
cette contrainte n’est nullement justifiée par les observations numériques effectuées dans ces cas là. Nous ren-
voyons aux sections 7 de [2] et 7 de [5] pour une description détaillée.

Algorithme avec projection “au plus proche voisin”.

Théorème 4.3.1 (Delarue et Menozzi [2]) Supposons les coefficients b, σ, f dans (4.4) homogènes et uni-
formément Lipschitziens en x, y, z, uniformément bornés en espace et à croissance au plus linéaire en (y, z), b
indépendant de z, σ uniformément elliptique et la fonction terminale g ∈ C 2+α , α ∈ (0, 1] bornée et à dérivées
bornées, sur le modèle de l’énoncé du Théorème 4.1.1. Si l’on prend des grilles spatiales Cartésiennes de pas
d’espace δ > 0 avec (Ci )i∈[[1,N ]] de rayon R + ρ, R, ρ > 0 et C0 de rayon R on a le contrôle suivant pour l’erreur
associée à l’algorithme (4.14) :

sup |u(0, x) − ū(0, x)| ≤ C (E(temps) + E(espace) + E(quantification) + E(troncature)) ,

x∈C0

où C dépend des hypothèses sur les coefficients et est uniforme par rapport aux paramètres de discrétisation. Les
ordres des erreurs sont : E(temps) = h1/2 , E(espace) = h−1 δ, E(quantification) = h−1/2 M −1/d et E(troncature) =
ρ/(ρ + R).
Ces ordres s’expliquent de la façon suivante :
- E(temps). L’ordre 1/2 provient de la régularité en temps de l’accroissement Brownien et de la propriété de
1/2-Hölder continuité en temps de ∇x u sous les hypothèses du théorème.
- E(espace). Pour une fonction régulière l’erreur associée à la projection au plus proche voisin sur une grille
spatiale Cartésienne de pas δ est d’ordre δ. Comme ce type d’erreur apparaı̂t à chaque itération en temps on
obtient bien h−1 δ.
- E(quantification). L’erreur de quantification de l’accroissement Brownien est en norme L1 d’ordre h1/2 M −1/2
d’après (4.12). La vitesse indiquée correspond de nouveau à la sommation de N = T /h erreurs le long de la
grille en temps.
- E(troncature). Ce terme est obtenu par des contrôles de type Tchebychev. Nous n’avons pour l’instant pas su
tirer profit de l’échelle caractéristique de temps en t1/2 de la diffusion pour mieux calibrer nos grilles spatiales.
Précisons que le résultat est présenté pour des coefficients homogènes en temps par souci de cohérence avec
les énoncés de [2]. Ce choix avait été motivé par le souci d’alléger les notations. Le résultat ci-dessus reste-
rait vrai sans difficultés techniques additionnelles pour des coefficients uniformément 1/2-Hölder en temps, cf.
Théorème 4.1.1.

39
 Algorithme interpolé.

Théorème 4.3.2 (Delarue et Menozzi [5]) Supposons les coefficients b, σ, f dans (4.4) homogènes et uni-
formément α-Hölderiens, α ∈ (0, 1], en x, uniformément Lipschitziens en y, z, uniformément bornés en espace
et à croissance au plus linéaire en (y, z), b indépendant de z, σ uniformément elliptique et la fonction terminale
g ∈ C 2+α , bornée et à dérivées bornées. Si l’on prend des grilles spatiales Cartésiennes infinies de pas d’espace
δ > 0 on a le contrôle suivant pour l’erreur associée à l’algorithme (4.15) :

sup |u(0, x) − ū(0, x)| ≤ C (E(temps) + E(espace) + E(quantification)) ,

x∈C0

où C dépend des hypothèses sur les coefficients et est uniforme par rapport aux paramètres de discrétisation.
Les ordres des erreurs sont : E(temps) = hα/2 , E(espace) = h−1 δ 2 , E(quantification) = hα/2 M −2/d .

Le propos du travail [5] était d’estimer le gain découlant de l’interpolation et de l’utilisation intensive de la
stationnarité du quantifieur. Nous n’avons donc pas pris en compte les erreurs associées à la troncation des
grilles Cartésiennes qui pourraient être analysées de façon similaire au travail précédent. Sous les hypothèses
du théorème nous avons toujours que u est une solution classique de (4.4) 1 , cf. Théorème 7.1. du Chapitre 7 de
[LSU68] et [MPY94], à une procédure de régularisation près des coefficients. L’EDSR fortement couplée associée
est alors à comprendre au sens faible, cf. Delarue et Guatteri [DG06].
Les ordres précédents se justifient comme suit :
- E(temps). L’ordre α/2 provient de la composition de la régularité α-Hölder des coefficients en x avec la régularité
en temps de l’accroissement Brownien.
- E(espace). Pour une fonction régulière l’erreur associée l’interpolation linéaire sur une grille spatiale Cartésienne
de pas δ est d’ordre δ 2 . L’itération le long de la grille en temps conduit à l’ordre indiqué.
- E(quantification). L’espérance d’une fonction “régulière” de la quantification Gaussienne optimale renormalisée
par h1/2 est d’après (4.13) d’ordre hM −2/d . La fonction régulière n’est autre dans l’analyse de convergence que
la solution u de (4.4). Ainsi l’itération de la procédure devrait conduire à un contrôle de M −2/d . Nous obtenons
une vitesse de hα/2 M −2/d en utilisant la α-Hölder continuité des dérivées secondes en espace.
Comme nous l’avons indiqué en Section 4.1.2, la démarche pour l’analyse de l’erreur est similaire à celle
utilisée dans [MPY94] pour établir l’unicité de la solution à l’EDSR fortement couplée. Cela nous permet
d’étendre les résultats de Milstein et Tretyakov qui supposent que la solution u de (4.4) est quatre fois dérivable
et utilisent un produit tensoriel de lois de Bernoulli symétriques pour approcher l’incrément gaussien. Par
ailleurs si l’on prend formellement α = 2, M = 2d , δ = h dans le théorème précédent on retrouve bien l’ordre h
obtenu dans [MT04, MT06b, MT07].
Nous conclurons cette section en insistant sur le faible impact de la quantification sur l’erreur globale, ce
dont atteste la borne sur E(quantification).

4.3.4 Résultats numériques

Plusieurs exemples ont été développés dans [2] et [5] : équations de Burgers et KPZ déterministe (voir
ci-après) et des milieux poreux. Par souci de concision nous ne présentons que deux exemples issus de [5] dans
le cadre multidimensionnel.
Equation de Burgers multi-dimensionnelle. Les équations de Burgers sont une forme simplifiée des équations
de Navier-Stokes. Les parties convectives et dissipatives sont les mêmes mais elles ne font intervenir ni le gradient
de pression, ni la contrainte d’incompressibilité. Les équations s’écrivent :

ε2
∂t u − (u.∇)u + ∆u = 0, (t, x) ∈ [0, T ) × Rd , ε > 0,
2 (4.16)
u(T, x) = g(x), x ∈ Rd ,

où ∀i ∈ [[1, d]], ((u.∇)u)i = ∇ui × u. Bien que le Théorème 4.3.2 soit établi pour des fonctions à valeurs réelles,
l’analyse précédente s’étendrait sans difficultés à des systèmes d’équations.
En dimension un il est bien connu que l’équation de Burgers “visqueuse” possède une solution explicite
donnée par une factorisation de type Cole-Hopf, voir e.g. Whitham [Whi73]. Cette approche s’étend au cas
1. i.e. bornée et à dérivées spatiales d’ordre un et deux bornées.

40
multi-dimensionnel lorsque la condition terminale dérive d’un potentiel, i.e. g = ∇g0 où g0 est une fonction à
valeurs réelles. On a dans ce cas : ∀(t, x) ∈ [0, T ] × Rd ,
E[∇g0 (x + εWT −t ) exp(−ε−2 g0 (x + εWT −t ))]
u(t, x) = .
E[exp(−ε−2 g0 (x + εWT −t ))]
Nous avons pour les résultats ci-après considéré l’interprétation “couplée” de l’équation (4.16), i.e. la non linéarité
apparaı̂t en injectant l’approximation de la solution dans la dérive de l’équation progressive. Ce choix semble
plus robuste numériquement et dans le cas particulier de (4.16) permet même d’éviter l’étape d’estimation du
gradient. Q
Nous considérons le cas d = 2 pour la fonction H0 (x) = i=1,2 sin2 (πxi ) qui est symétrique et périodique
permettant ainsi d’éviter les problèmes de troncation. Une des difficultés est de construire une grille spatiale
qui permette de “voir” les transitions de la diffusion approchée qui sont de l’ordre de εh1/2 . Ceci impose que δ
soit plus petit que εh1/2 . Nous prenons le jeu suivant de paramètres : T = 3/8, h = 2.5 × 10−2 , δ = .01, ε2 = .4.
Comme la solution est symétrique nous ne présentons que les résultats associés à la première composante u1 de
u = (u1 , u2 )∗ .
Nous présentons les profils de la solution, ainsi que la différence entre la valeur de référence, obtenue à partir
de la représentation explicite pour 600 points de quantification, et l’approximation de l’algorithme (4.15). Un
∗
autre graphique représente ensuite les résultats obtenus sans facteur de normalisation γM pour la quantification
de l’accroissement Brownien et illustre l’importance de cette procédure. Nous étudions enfin l’influence du
nombre de points M utilisés pour la quantification en prenant successivement M = 4 et M = 150. Les résultats
corroborent la conclusion de la section précédente, à savoir que le nombre de points de quantification n’a qu’un
faible impact.
Burgers equations d=2: profiles Burgers equations d=2: Algorithm transition, pointwise absolute error

True Solution t=0 t=0

0.06 0.008
0.04 0.007
0.006
0.02 0.005
0 0.004
-0.02 0.003
0.002
-0.04 0.001
-0.06 0

1 1
0.8 0.8
0 0.6 0 0.6
0.2 0.4 0.2 0.4
0.4 0.4
0.6 0.2 0.6 0.2
0.8 0.8
1 0 1 0

Burgers equations d=2: "Natural" transition, pointwise absolute error Burgers equations d=2: pointwise absolute difference between M=4 and M=150

t=0 t=0

0.12 0.0007
0.1 0.0006
0.08 0.0005
0.0004
0.06
0.0003
0.04 0.0002
0.02 0.0001
0 0

1 1
0.8 0.8
0 0.6 0 0.6
0.2 0.4 0.2 0.4
0.4 0.4
0.6 0.2 0.6 0.2
0.8 0.8
1 0 1 0

Equation de KPZ déterministe.

L’équation KPZ 2 peut se voir comme la primitive de l’équation de Burgers. Elle s’écrit :
1 1
∂t u + |∇u|2 + ∆u = 0, (t, x) ∈ [0, T ) × Rd ,
2 2
u(T, x) = g(x), x ∈ Rd . (4.17)
Cette équation possède également une représentation explicite obtenue par factorisation Cole-Hopf : u(t, x) =
log(E[exp(g(x + WT −t ))]). C’est une équation quadratique en gradient pour laquelle les approches proposées ci-
avant ne sont pas tout à fait satisfaisantes quant à l’analyse de l’erreur. En revanche, dans l’exemple suivant, le
comportement numérique des algorithmes “couplés”, où l’on voit le gradient comme la dérive dans la composante
progressive, et découplés, où le carré du gradient apparaı̂t au générateur, est très différent.
Soit d = 2, g(x) = 10 cos(5|x|2 ), T = .1, h = .02, δ = .02, M = 4. On tronque la grille pour |x| ≥ 1. La valeur
de référence est obtenue par 200 points de quantification. Il vient :
2. Cette équation apparaı̂t en mécanique statistique où elle est de plus forcée par un bruit blanc, cf. [KPZ86]. Elle intervient
dans des modèles de croissance ou d’interface.

41
 KPZ equation d=2: coupled case and true solution KPZ equation d=2: backward case

Reference value, t=0 t=0

Coupled algorithm, t=0

8.8 2.5e+19
8.6
8.4 2e+19
8.2
8 1.5e+19
7.8
7.6 1e+19
7.4
7.2 5e+18
7
6.8 0

0.4 0.4
0.2 0.2
-0.4 0 -0.4 0
-0.2 -0.2
0 -0.2 0 -0.2
0.2 -0.4 0.2 -0.4
0.4 0.4

Dans la pratique l’intégration numérique de grands gradients provoque des overflow, et incite dans ce cas de
figure à préférer l’interprétation couplée plus “robuste”. Indiquons qu’un exemple tridimensionnel pour (4.17)
est présenté dans [5]. Pour la grille Cartésienne de pas δ, le nombre de noeuds pour lesquels il faut mettre à jour
la solution et son gradient est de l’ordre de δ −d donnant une complexité globale de h−1 M δ −d ce qui interdit
dans la pratique les grandes dimensions.

4.3.5 Extensions et perspectives

1. L’analyse précédente se base sur les propriétés de régularité des solutions d’EDP quasi-linéaires à coefficients
réguliers et terme de diffusion elliptique. Ce contexte permet de mener à bien une analyse de stabilité “usuelle”.
Un enjeu serait de proposer une extension de la méthode précédente pour l’approximation de solutions faibles
ou variationnelles de l’équation (4.4). Un point de départ pourrait être le travail de Barles et Lesigne [BL97] qui
établit dans le cadre de problèmes semi-linéaires sous forme divergence une équivalence entre formulation varia-
tionnelle usuelle et formulation variationnelle associant l’EDSR et des fonctions tests aléatoires, cf. également
Bally et Matoussi pour un cadre infini-dimensionnel [BM01].
2. Le problème stationnaire pour les EDP quasi-linéaires non dégénérées est bien compris, cf. Ladyzenskaya et
Uralts’eva [LU68] ou Gilbarg et Trudinger [GT77]. En revanche l’approximation numérique reste ouverte. En
effet comment procéder autrement que par la résolution approchée d’une succession de problèmes “découplés”
avec temps d’atteinte ? Dans ce cas là, que dire du nombre d’itérations nécessaire pour obtenir une convergence ?
Un autre problème lié à ce contexte consiste d’un point de vue théorique à s’interroger sur l’extension des
résultats d’existence et d’unicité à l’EDSR fortement couplée avec temps d’arrêt sous les hypothèses naturelles
de Lipschitz continuité et non dégénérescence des coefficients. Là aussi l’approche itérative par succession de
problèmes découplés semble la plus naturelle.
3. Nous avons abordé dans la partie numérique le cas de l’équation de Burgers. Belopolskaya et Milstein [BM03]
ont proposé pour les équations de Navier-Stokes dans le cas périodique tridimensionnel en temps court un
algorithme d’approximation basé sur une itération rétrograde. L’idée est d’exploiter la périodicité pour faire
un développement en série de Fourier de la solution approchée et d’utiliser la décomposition d’Helmoltz pour
conserver pour lePschéma la contrainte d’incompressibilité. Le gradient de pression est mis à jour en estimant
l’équation ∆p = (i,j)∈[[1,d]]2 ∂xi uj ∂xj ui qui relie le Laplacien de la pression aux gradients de la solution supposés
réguliers, i.e. au moins continus. La difficulté de cette approche est de passer des contrôles L2 naturels à ce
contexte aux contrôles ponctuels nécessaires pour avoir la convergence.
Dans l’optique de mieux appréhender ce problème, il serait intéressant d’envisager une représentation de
la solution en terme d’EDSR. Une idée serait de procéder par pénalisation pour avoir l’incompressibilité. La
gestion du gradient de pression est plus délicate dans la mesure où ce terme est formellement lié à des dérivées
secondes que l’on retrouveraient en terme source. Notons que les formulations vortex semblent mal adaptées
aux interprétations rétrogrades. La structure de convolution rend naturelles l’approche de Mc Kean et Vlasov
et les méthodes numériques particulaires associées, cf. Méléard [Mél01].
Nous discutons avec Bruno Bouchard de Dauphine à propos de problèmes d’existence et d’unicité pour des
EDSR avec contraintes sur le Z.
4. Le cas quadratique est mal compris du point de vue numérique alors qu’il est assez naturel et bien compris tant
du point de vue des EDSR [Kob00] que des EDP [LSU68]. Pour l’approximation ce sont bien sûr les contrôles
sur le gradient qui sont à raffiner.

4.4 Semi-discrétisation d’EDSR dans un domaine

Nous présentons dans cette Section le travail [6], en commun avec B. Bouchard de l’Université Paris Dau-
phine, relatif à l’approximation forte d’EDSR dans un domaine. Il s’agit là de l’analyse d’un schéma de semi-

42
 discrétisation dont la spécificité est de prendre en compte des domaines non-réguliers et des coefficients de
diffusion dégénérés dans un contexte parabolique, i.e. un instant d’observation final T > 0 est fixé.

4.4.1 Cadre de travail

Nous considérons l’EDS faiblement couplée
Z t Z t
Xt = x + b(Xs )ds + σ(Xs )dWs ,
0 0
Z T Z T
Yt = g(τ, Xτ ) + Is<τ f (Xs , Ys , Zs )ds − Zs dWs , (4.18)
t t

où les coefficients b, σ, f, g de (4.18) sont uniformément Lipschitziens en espace et τ est le temps de sortie de
(t, Xt ) d’un domaine cylindrique D = [0, T ) × D où D ⊂ Rd est un domaine régulier par morceaux.
Sous de bonnes hypothèses on peut montrer que u(t, x) := Ytt,x est l’unique solution de viscosité de
(
(∂t + L)u(t, x) + f (t, x, u(t, x), (σ∇u)(t, x)) = 0, (t, x) ∈ [0, T ) × D,
u(t, x) = g(t, x), (t, x) ∈ [0, T ) × ∂D ∪ {T } × D̄,

où L désigne le générateur infinitésimal de X, cf. Darling et Pardoux [DP97] pour le cas elliptique qui peut
s’adapter à (4.18).
Considérons maintenant le schéma de semi-discrétisation suivant le long de la grille en temps (ti := ih)i∈[[1,N ]]
où N h = T :

YTh = g(τ h , Xτhh ),

Ythi = E[Ythi+1 |Fti ] + Iti <τ h hf (Xthi , Ythi , Zthi ), Zthi = h−1 E[Ythi+1 (Wti+1 − Wti )|Fti ], i ∈ [[0, N − 1]],(4.19)

où X h est le schéma d’Euler de l’équation progressive de (4.18), τ h := inf{ti := ih : (ti , Xthi ) 6∈ D} et g désigne
une extension de la condition de bord à tout l’espace. Remarquons que Ythi Iti ≥τ h = g(τ h , Xτhh ) et par conséquent
Zthi Iti ≥τ h = 0.
Le propos est de fournir en fonction du pas de discrétisation un contrôle sur l’erreur :
Z T
Err(h, T )2 := max E[ sup |Yt − Ythi |2 ] + E[ h
|Zt − Zφ(t) |2 dt],
i≤N t∈[ti ,ti+1 ] 0

où φ(t) := sup{ti := ih : ti ≤ t}.

4.4.2 Résultats
Supposons que la condition terminale g vérifie :
(Hg) g ∈ C 1,2 ([0, T ] × Rd ) et il existe C0 > 0 telle que |∂t g|∞ + |∇g|∞ + |D2 g|∞ ≤ C0 .
Nous avons alors le résultat suivant :
Proposition 4.4.1 Il existe C dépendant uniquement des coefficients de (4.18), de la borne dans (Hg) et de
T telle que   Z τ 
2 h h h 2
Err(h, T ) ≤ C h + R(Y ) + R(Z) + E ξ|τ − τ | + Iτ h <τ |Zs | ,
τh

où ξ est une variable aléatoire positive telle que E[ξ p ] ≤ C p , p ≥ 2, et R(Y )h := maxi≤N E[supt∈[ti ,ti+1 ] |Yt −
RT Rt
Yti |2 ], R(Z)h := E[ 0 |Zt − Ẑφ(t) |2 dt], Ẑti := h−1 E[ tii+1 Zs ds|Fti ], i ∈ [[0, N − 1]].

Les termes R(Y )h et R(Z)h sont usuels dans l’analyse des schémas de semi-discrétisation d’EDSR et appa-
raissent déjà dans [Zha04]. A l’aune du lien entre EDSR et EDP, il apparaı̂t néanmoins clairement qu’obtenir
un contrôle à l’ordre h pour R(Y )h , R(Z)h nécessite d’avoir des propriétés de Lipschitz continuité en espace de
la solution de viscosité de l’EDP (4.19). Pour obtenir ces estimées et contrôler l’écart entre les temps d’arrêts il
est nécessaire de rajouter quelques hypothèses de “structure” sur le domaine, cf. Chapitre 2.
Introduisons :

43
 ∗
(D) Le domaine D := ∩m 2
l=1 Dl , m ∈ N et les (Dl )l∈[[1,m]] sont de classe C à frontière compacte. Nous suppo-
serons également que D vérifie une condition de sphère extérieure uniforme et de cône intérieur uniforme (cf.
Section 2.2 de [6] pour une définition précise de ces deux dernières conditions).
Tm
Désignons par C := l6=k=1 ∂Dl ∩ ∂Dk l’ensemble des coins et par B(C, η) un voisinage de rayon η > 0 de
cet ensemble. Nous supposerons également :
(C) Le coefficient de diffusion a vérifie une condition de frontière non caractéristique uniforme hors d’un voisinage
B(C, η), η > 0 des coins, i.e ∃C > 0, ∀x ∈ ∂D\B(C, η), ha(x)n(x), n(x)i ≥ C où n désigne la normale rentrante
au domaine qui est bien définie hors d’un voisinage des coins, cf. Chapitre 2. Le coefficient a est uniformément
elliptique dans un voisinage des coins, i.e. ∃C > 0, ∀x ∈ D̄ ∩ B(C, η), ∀ξ ∈ S d−1 , ha(x), ξ, ξi ≥ C.
Ces conditions permettent d’utiliser des arguments d’EDP de type fonctions barrières pour contrôler le
module de continuité au bord du domaine. A l’intérieur du domaine, où le coefficient de diffusion peut dégénérer,
ce même contrôle est obtenu par des techniques probabilistes plus standards de flots stochastiques. Nous obtenons
ainsi le résultat suivant :
Théorème 4.4.1 Si les coefficients b, σ, f, g sont uniformément Lipschitziens de constante K et que les hy-
pothèses (Hg), (D), (C) sont vérifiées

∃C := C(K, T, (Hg), (D), (C)), R(Y )h + R(Z)h ≤ Ch,

et l’unique solution de viscosité u de (4.19) dans la classe des solutions continues à croissance polynomiale est
uniformément 1/2 Hölder continue en temps et Lipschitz en espace, i.e. pour le même C :

|u(t, x) − u(t′ , x′ )| ≤ C(|t − t′ |1/2 + |x − x′ |), ∀(t, x), (t′ , x′ ) ∈ D.

Nous avons par ailleurs le contrôle suivant pour le temps d’arrêts.

Théorème 4.4.2 Sous les hypothèses du Théorème 4.4.1, pour tout ε ∈ (0, 21 ), ∃C := C(ε, K, T, (Hg), (D), (C)),

E[|τ − τ h |] ≤ Ch1/2−ε .

Cette estimée s’obtient par des techniques de martingales et des arguments de changement de temps assez
similaires à ceux introduits au Chapitre 2. Par ailleurs, à ε près ce contrôle est optimal, cf. Théorème 2.2.2.
C’est ainsi l’approximation du temps d’arrêt qui domine dans l’erreur et l’on déduit de la Proposition 4.4.1
et des Théorèmes 4.4.1, 4.4.2 la borne finale :
1
∀ε ∈ (0, ), ∃C := C(ε, K, T, (Hg), (D), (C)), Err(h, T )2 ≤ Ch1/2−ε .
2

4.4.3 Extensions et perspectives

Nous avons retenu pour la clarté de la présentation dans [6] un domaine cylindrique. Le résultat précédent
s’étendrait sans difficultés supplémentaires à domaine temps-espace vérifiant des hypothèses similaires à (C),(D).
De même, pour un coefficient de diffusion non dégénéré il semble naturel de retrouver les même estimées pour
le problème elliptique, dans la mesure où les moments des temps de sorties sont bien contrôlés pour le schéma
d’Euler, voir Gobet et Maire [GM05] ou [7]. En revanche les situations elliptiques dégénérées restent à ana-
lyser, en particulier le transfert des propriétés des temps d’atteinte associés à l’objet continu au schéma de
discrétisation.

44
Chapitre 5

Estimation de densité : cas dégénérés

et applications

Nous présentons ici les résultats des articles [10], écrit en collaboration avec François Delarue de l’Université
de Nice, et [11], en collaboration avec Vincent Lemaire de l’Université Paris 6.
Nous commençons par une synthèse des principales techniques intervenant dans l’estimation quantitative de
densités de processus stochastiques. On soulignera en particulier les hypothèses de nature différente en fonction
des outils utilisés : hypoellipticité et régularité des coefficients pour le calcul de Malliavin, et du côté analytique,
uniforme ellipticité et mesurabilité pour des opérateurs sous forme divergence et uniforme ellipticité et Hölder
continuité pour des opérateurs sous forme non divergence, voir aussi le Chapitre 3. Dans ce dernier cas, la notion
de problème de martingale de Stroock et Varadhan permet d’unifier l’approche probabiliste et analytique.
Nous détaillons ensuite l’article [10] où nous considérons une EDS vérifiant une hypothèse d’Hörmander
faible, i.e. c’est le coefficient de dérive qui propage le bruit dans le système. Sous de bonnes hypothèses, proches
de celles naturelles du point de vue analytique, nous retrouvons des estimées Gaussiennes de type Aronson en
régime non diffusif. Pour obtenir la minoration, une première approche consiste à utiliser la transformée de
Fleming qui permet d’exprimer le logarithme de la densité de la diffusion comme solution d’un problème de
contrôle stochastique. Nous donnerons également une seconde preuve plus géométrique. La majoration découle
quant à elle de l’approche parametrix.
Les techniques développées dans [10] peuvent s’adapter à des schémas de discrétisation de certains processus
dégénérés ou non. Dans [11] nous avons obtenu des estimées Gaussiennes pour les densités de certains schémas
d’Euler. L’utilisation de techniques de concentration de la mesure : inégalités de Sobolev logarithmique, argu-
ment de Herbst, nous a permis d’en déduire des intervalles de confiance non asymptotiques pour l’approximation
de Monte Carlo.

5.1 Estimation de densité : quelques techniques et résultats

5.1.1 Opérateurs non dégénérés sous forme divergence
Considérons un opérateur sous forme divergence
1
Lϕ(x) = div(a(x)∇ϕ(x)) + hb(x), ∇ϕ(x)i, ϕ ∈ C0∞ (Rd ),
2
où les coefficients a, b sont mesurables, b est borné, a est uniformément elliptique borné, i.e ∃Λ ≥ 1, ∀ξ ∈
Rd , Λ−1 |ξ|2 ≤ ha(x)ξ, ξi ≤ Λ|ξ|2 . Dans ce cadre, des techniques analytiques désormais classiques donnent que
pour tout T > 0, il existe des constantes C := C(T, Λ, |b|∞ ), c := c(Λ), (c, C) ∈ (0, 1] × [1, +∞) telles que si p
désigne la solution fondamentale associée à L,
   
C −1 −1 |x − y|
2
C |x − y|2
exp −c ≤ p(t, x, y) ≤ exp −c , ∀(t, x, y) ∈ (0, T ] × Rd × Rd . (5.1)
td/2 t td/2 t
Si par ailleurs b = 0 les estimées sont uniformes en temps, i.e. les constantes (c, C) := (c, C)(Λ) et (5.1) est
valable pour tout t > 0. Nous appellerons pour la suite de ce chapitre inégalités de Aronson 1 des inégalités de
la forme (5.1), avec éventuellement une échelle diagonale non diffusive.
1. qui dans son travail fondateur [Aro67] a établi les inégalités (5.1)

45
 La borne supérieure s’obtient typiquement à l’aide des inégalités de Nash et d’opérations sur la forme de
Dirichlet associée à L. La borne inférieure s’obtient via la représentation spectrale du noyau de la chaleur tué
de L et un argument de chaı̂nage, voir Chapitre 7 de Bass [Bas97] ou Stroock [Str88].
Dans cette approche seule la non dégénérescence de a compte. Pour b = 0, l’aspect auto-adjoint est également
très important dans les preuves et donne l’uniformité en temps.

5.1.2 Opérateurs sous forme non divergence

Soit maintenant
1
Lt ϕ(x) = tr(a(t, x)D2 ϕ(x)) + hb(t, x), ∇ϕ(x)i, t ≥ 0, ϕ ∈ C0∞ (Rd ).
2
Nous avons vu au Chapitre 3.1.2 que pour a uniformément elliptique borné, η-Hölder continu en espace, η > 0 et
b mesurable borné, on obtenait la borne supérieure de (5.1) à l’aide de la représentation en série parametrix, voir
aussi Friedman [Fri64]. Cette représentation fournit également la borne inférieure de (5.1) pour |x − y|/t1/2 ≤
C0 , t ≤ T0 , ce qui avait été initialement remarqué par Il’in et al. [IKO62]. L’extension à |x − y|/t1/2 > C0 , t ≤ T0
se fait par un argument de chaı̂nage tout à fait similaire à celui effectué pour les opérateurs sous forme divergence.
Nous reviendrons sur la structure du chaı̂nage que l’on peut en fait relier à un problème de contrôle déterministe
en Section 5.2. L’extension à un temps arbitraire s’obtient par convolution ou par des arguments de changement
d’échelle (voir Section 2.3 de [10]).
Par cette approche la constante C de (5.1) dépend du temps même pour b = 0. Nous obtenons toutefois
les bornes Gaussiennes de façon assez directe pour la densité de l’unique solution du problème de martingale
associé à (a, b). Cette analyse s’étend sous les mêmes hypothèses au schéma d’Euler, voir Section 5.3. Même si
l’on ne peut plus ignorer toute régularité spatiale sur a, la condition est assez peu contraignante.

5.1.3 Opérateurs hypoelliptiques et calcul de Malliavin

P
Considérons l’opérateur L = 21 ri=1 A2i + A0 sous forme de Hörmander, où les (Ai )i∈[[1,r]] sont des champs
Rt Rt
de vecteurs réguliers sur Rd , et qui est le générateur de la diffusion Xt = x + 0 A0 (Xs )ds + 0 Ai (Xs ) ◦ dWsi
où ◦dW désigne l’intégrale de Stratonovitch. Si la dimension de l’espace vectoriel engendré par

A1 , · · · , Ar , [Ai , Aj ]0≤i,j≤r , [Ai , [Aj , Ak ]]0≤i,j,k≤r , · · ·

au point initial x est Rd (où [, ] désigne le crochet de Lie entre champs de vecteurs) 2 , alors pour tout t > 0, Xt
admet une densité pt , cf. Norris [Nor86] ou Nualart [Nua98]. La question de l’estimation de cette densité est
dans le cadre hypoelliptique général très difficile. Il est en effet très délicat de mettre en évidence une échelle
globale intrinsèque au système.
Plusieurs résultats ont néanmoins été obtenus l’aide du calcul de Malliavin. Cet outil (voir Malliavin [Mal97],
Stroock [Str83] et Nualart [Nua95]) permet de quantifier la sensibilité du système par rapport au bruit et cette
sensibilité se lit sur sur les dérivées de Malliavin et la matrice de covariance associée. La non-dégénérescence de
cette matrice donne l’existence de la densité. L’utilisation de cette technique pour l’estimation de densité est
due à Kusuoka et Stroock [KS84, KS85, KS87]. Dans le dernier travail, les auteurs établissent des estimées de
type Aronson sous des hypothèses d’Hörmander fortes, i.e. le bruit ne se propage dans le système qu’à travers la
composante diffusive. La normalisation diagonale provient alors du volume des boules hypoelliptiques de rayon
t1/2 centrées au point initial. C’est en effet une normalisation diffusive mais les boules ne sont plus euclidiennes
et le caractère “local” de leur volume apparaı̂t clairement. La décroissance hors diagonale est donnée par la
métrique de Carnot-Carathéodory, ou métrique de contrôle, associée aux champs (Ai )i∈[[1,r]] . Ces résultats ne
sont valables que pour une dérive nulle ou engendrée par les champs de vecteurs de la partie diffusive.
Pour d’autres exemples d’application du calcul de Malliavin nous renvoyons à Bally [Bal90], pour une
discussion sur les liens entre matrice de Malliavin et hypothèses de Hörmander, à Cattiaux [Cat90], pour les
propriétés de résolvantes de diffusions hypoelliptiques, à Cattiaux et Mesnager [CM02] pour l’existence de
la densité dans le cadre non-homogène. Pour des résultats plus spécifiques sur l’estimation de densité nous
renvoyons à Ben Arous et Léandre [BL91], qui ont mis en évidence plusieurs régimes pour le comportement en
temps petit de la densité de certaines diffusions hypoelliptiques en fonction de la forme de leur dérive.
En résumé, dans un cadre dégénéré l’utilisation du calcul de Malliavin pour obtenir des estimées de densités
demande de la régularité sur les coefficients. Cette régularité n’est pas toujours nécessaire. Dans la section
suivante nous arrivons à obtenir pour une diffusion faiblement hypoelliptique des estimées de densités sous
2. Dans l’article initial [Hör67], Hörmander supposait cette condition vérifiée en tout point de l’espace.

46
 des hypothèses proches de celles naturelles pour les opérateurs non dégénérés sous forme non divergence. En
particulier peu de régularité des coefficients est nécessaire.

5.2 Estimations de densités pour un bruit propagé dans une chaı̂ne

d’Equations différentielles
Le modèle que nous avons considéré pour l’article [10] est le suivant :

dXt1 = F1 (t, Xt1 , . . . , Xtn )dt + σ(t, Xt1 , . . . , Xtn )dWt ,

dXt2 = F2 (t, Xt1 , . . . , Xtn )dt,
dXt3 = F3 (t, Xt2 , . . . , Xtn )dt, t ≥ 0, (5.2)
···
dXtn = Fn (t, Xtn−1 , Xtn )dt,

où (Wt )t≥0 désigne un mouvement Brownien d-dimensionnel et pour tout i ∈ [[1, n]], (Xti )t≥0 , est aussi à valeurs
dans Rd .
Un exemple typique pour (5.2) est donné par un système de n oscillateurs couplés, chacun d’entre eux
bougeant verticalement et étant directement relié à son plus proche voisin, lorsque le premier oscillateur est
excité par un bruit aléatoire.

BRUIT ⇒

Précisons que les Fi dans (5.2) ne peuvent pas dépendre des positions Xt1 , . . . , Xti−2 : d’un point de vue physique
le bruit doit traverser les (i − 1) premiers oscillateurs avant d’atteindre le ième . Dans la dynamique différentielle
n’intervient que la composante qui transmet le bruit. Intuitivement cette condition permet la séparation des
échelles. Supposons que le bruit soit en effet transmis par la composante X i−1 à la composante X i , ce qui
sera assuré par la condition de non dégénérescence des Dxi−1 Fi (hypothèse de type Hörmander faible). Alors,
l’échelle caractéristique de X i est ti−1/2 pour tout i ∈ [[1, n]].
Des systèmes de ce type apparaissent dans différents champs d’applications. Pour n = 2, (5.2) décrit la
dynamique de certains systèmes hamiltoniens, voir Soize [Soi94] pour une synthèse ou Talay [Tal02] et Hérau
et Nier [HN04] pour l’aspect convergence à l’équilibre. Toujours pour n = 2, (5.2) correspond à la dynamique
utilisée en finance pour les options asiatiques, cf. [BPV01] pour des questions de régularité du prix associé. Pour
n ≥ 2, cette équation apparaı̂t dans des modèles de conduction de la chaleur, cf. Eckmann et al. [EPRB99] et
Rey-Bellet et Thomas [RBT00] où la chaı̂ne d’oscillateurs est chauffée aux deux extrémités par des bains de
chaleur. Mentionnons également le travail récent de Bodineau et Lefevere [BL08]).
Pour le reste de la section notre jeu d’hypothèses (noté (A)) est le suivant :

(a) Le spectre de la fonction matricielle (t, x) := (t, x1 , · · · , xn ) ∈ R+ × Rnd 7→ a(t, x) = [σσ ∗ ](t, x) est inclus
dans [Λ−1 , Λ]pour Λ ≥ 1.
(b) Les fonctions F1 , . . . , Fn et σ sont respectivement uniformément Lipschitz et η-Hölder continues (η ∈ (0, 1])
par rapport aux variables d’état, pour une constante positive κ.
(c) Pour tout i ∈ [[2, n]], (t, (xi , . . . , xn )) ∈ R+ × R(n−i+1)d , la fonction xi−1 ∈ Rd 7→ Fi (t, xi−1 , . . . , xn ) est
continûment différentiable, la dérivée notée (t, xi−1 , . . . , xn ) ∈ R+ × R(n−i+2)d 7→ Dxi−1 Fi (t, xi−1 , . . . , xn ), est
supposée η-Hölder continue de constante κ en la première variable spatiale xi−1 .
(d) Il existe un ensemble convexe fermé Ei−1 ⊂ GLd (R) (ensemble des matrices inversible d × d) tel que pour
tout t ≥ 0 et (xi−1 , . . . , xn ) ∈ R(n−i+2)d , la matrice Dxi−1 Fi (t, xi−1 , . . . , xn ) appartienne à Ei−1 . Par exemple,
Ei , 1 ≤ i ≤ n − 1, peut-être une boule fermée incluse dans l’ouvert GLd (R).

Le point (a) est une condition classique d’uniforme ellipticité de la matrice de diffusion. La condition (b) est
assez naturelle du point de vue des EDO/EDP. Elle garantit que l’EDO sous-jacente à (5.2), i.e. lorsque l’on

47
supprime le forçage stochastique, est bien posée. Par ailleurs on retrouve sur a l’hypothèse usuelle de régularité
pour les opérateurs sous forme non divergence. La condition (c) est utilisée pour l’analyse de densité : on suppose
un peu plus de régularité sur Fi par rapport à la composante xi−1 qui transmet le bruit.
Le point (d) traduit une condition de type Hörmander faible qui permet la propagation du bruit dans le
système. L’hypothèse de convexité est technique. Nous renvoyons à la Section 3 de [10] pour plus de détails.
En particulier, on retrouve comme au Chapitre 3 que les coefficients, supposés mesurables, peuvent être
irréguliers en temps (F (t, 0) borné).

Théorème 5.2.1 Supposons (A) vérifiée est que l’EDS (5.2) est uniquement soluble en loi. Alors, la solu-
tion (faible) de (5.2) admet pour tout t > 0 une densité p(t, ., .) qui vérifie des estimées de type d’Aronson.
Précisément, pour tout T > 0, ∃CT ≥ 1, CT := CT (Λ, η, κ, (Ei )i∈[[1,n−1]] , n, d), telle que pour tout t ∈ (0, T ],
(x, y) ∈ Rnd × Rnd :
2
2
2
2

CT−1 t−n d/2 exp −CT t Tt−1 θt (x) − y ≤ p(t, x, y) ≤ CT t−n d/2 exp −CT−1 t T−1 t θt (x) − y .

Ci-dessus, Tt est la matrice d’“échelle” du système. C’est une matrice diagonale de taille nd × nd qui se
décompose en n blocs diagonaux de taille d × d où le ième bloc diagonal est donné par ti Id , où Id est la matrice
identité d × d. De plus, θt (x) = (θt1 (x), θt2 (x), . . . , θtn (x)) désigne la valeur dans Rnd au temps t de la solution
de l’EDO déterministe associée à (5.2) :
1
θ˙t = F1 (t, θt1 , . . . , θtn ),
i
θ˙t = Fi (t, θti−1 , . . . , θtn ), 2 ≤ i ≤ n, θ 0 (x) = (θ01 (x), . . . , θn1 (x)) = x. (5.3)

Le transport de la condition initiale x par le flot déterministe θ peut se voir comme une généralisation non linéaire
de l’exemple de Kolmogorov abordé en Section 3.2.1. La matrice d’échelle traduit le fait que Xti , i ∈ [[1, n]] oscille
autour de θti (x) avec des fluctuations d’ordre ti−1/2 . Pour i = 1, l’action du flot (d’ordre t) est négligeable devant
le bruit (d’ordre t1/2 ) en temps petit. A l’opposé, pour i ≥ 2, les fluctuations de Xti peuvent être beaucoup plus
petites que la distance entre θi (x) et xi : le terme de transport joue alors un rôle clé dans l’estimée de densité
pour ces coordonnées.
L’effet multi-échelle dans l’estimation ci-dessus provient de la structure des coefficients (Fi )i∈[[1,n]] et de la
condition de Hörmander (d) qui nous donne que le bruit se propage suivant une seule direction dans le système.
Nous avons supposé ci-avant que F1 était Lipschitz, en autorisant de fait les croissances linéaires. Sous réserve
d’unicité en loi le résultat précédent resterait valable pour F1 mesurable borné.
Pour la suite nous réécrirons la dynamique de (5.2) sous la forme condensée

dXt = F(t, Xt ) + Bσ(t, Xt )dWt ,

où F = (F1 , . . . , Fn ) et B désigne la matrice d’injection de Rd dans Rnd .

5.2.1 Transformée de Fleming

Pour prouver le Théorème précédent nous allons exploiter la formulation de − log p(t, x, y) en terme de
problème de contrôle stochastique. Ce procédé, connu sous le nom de transformée de Fleming, a été utilisé par
Sheu [She91] pour retrouver les estimées d’Aronson pour des opérateurs sous forme non-divergence dans le cadre
du Paragraphe correspondant à la Section 5.1. La représentation Feynman-Kac associée à cette transformation
amène à considérer un processus contrôlé où le paramètre de contrôle intervient comme dérive additionnelle
sur la première composante, non dégénérée, de l’équation initiale. Un choix pertinent de contrôle permettra
d’obtenir une minoration de la densité.
Nous supposerons par la suite, pour justifier l’existence de la densité que les coefficients (Fi )i∈[[1,n]] et σ
sont réguliers en espace et en temps de sorte que le théorème de Hörmander appliqué à l’opérateur parabolique
∂t + Lt donne l’existence de la densité. Tous les résultats que nous obtenons sont uniformes par rapport aux
paramètres de régularisation. Le théorème se déduit ensuite par passage à la limite, sous réserve d’unicité en
loi, par des arguments de type Radon-Nikodym.
Fixons (T, y0 ) ∈ (0, +∞) × Rnd et posons pour ε > 0 et (t, x) ∈ [0, T − ε] × Rnd , uε (t, x) = E[ηε (Xt,x
T −ε )] où
Xt,x désigne la solution de (5.2) pour Xt = x et (ηε )ε>0 est une suite d’approximations régulières strictement
positives de δy0 . La fonction uε est solution du problème de Cauchy
(
∂t uε (t, x) + Lt uε (t, x) = 0, (t, x) ∈ [0, T − ε) × Rnd ,
uε (T − ε, x) = ηε (x), x ∈ Rnd .

48
Par ailleurs limε→0 uε (0, x) = limε→0 E[ηε (X0,x T −ε )] = p(T, x, y0 ).
Introduisons maintenant pour (t, x) ∈ [0, T − ε] × Rnd , Jε (t, x) = − ln[uε (t, x)]. Cette fonction vérifie l’EDP
quadratique suivante
(
∂t Jε (t, x) + Lt Jε (t, x) − 12 ha(t, x)Dx1 Jε (t, x), Dx1 Jε (t, x)i = 0, (t, x) ∈ [0, T − ε) × Rnd ,
(5.4)
Jε (T − ε, x) = − ln(ηε (x)), x ∈ Rnd .

Le point clé consiste à remarquer que la partie différentielle précédente peut se réécrire
 1 
∂t Jε (t, x) + Lt Jε (t, x) + inf hv, Dx1 Jε (t, x)i + ha−1 (t, x)v, vi = 0, (5.5)
v∈R d 2

où l’infimum est atteint en v = v(t, x) := −a(t, x)Dx1 Jε (t, x) et redonne (5.4). Notons P(T − ε) l’ensemble des
R T −ε
processus progressivement mesurables (vt )t∈[0,T −ε] à valeurs dans Rd tels que E[ 0 |vt |2 dt] < +∞. A tout
(vt )0≤t≤T −ε ∈ P(T − ε) on peut associer une version contrôlée de (5.2), i.e.

dχt = [F(t, χt ) + Bvt ] dt + Bσ(t, χt )dWt , χ0 = x. (5.6)

En exploitant (5.5), nous pouvons écrire Jε (0, x) comme la fonction valeur d’un problème de contrôle stochas-
tique : " Z #
1 T −ε −1 0,x 0,x
Jε (0, x) = inf E ha (t, χt )vt , vt idt − ln[ηε (χT −ε )] . (5.7)
(vt )t ∈P(T −ε) 2 0

L’équation (5.7) associée à l’identité − ln(p(T, x, y0 )) = limε→0 Jε (0, x) constitue l’élément clé de notre preuve.
Elle intervient pour les deux bornes du théorème.

5.2.2 Etapes pour la borne inférieure

L’idée principale est la suivante. Pour obtenir une borne inférieure sur p(T, x0 , y0 ), x0 ∈ Rnd , il est suffisant
d’obtenir une borne supérieure pour Jε (0, x0 ), uniformément en ε > 0. D’après la formulation (5.7), tout le
problème consiste à trouver un contrôle (vt )t∈[0,T −ε] pertinent pour tout ε > 0.
Le choix de (vt )0≤t≤T −ε s’effectue en trois étapes :
(i). Nous avons d’abord caractérisé la densité dans le cas linéaire, i.e. lorsque la dérive F est affine en espace,
et la matrice σ ne dépend que du temps. La solution (Xt )t≥0 de (5.2) est alors un processus Gaussien. Sous (A),
la densité de transition existe et possède une forme explicite. De plus le terme hors diagonale de la densité peut
s’identifier à l’énergie du contrôle optimal pour le problème de contrôllabilité linéaire déterministe, cf. Section
3 de [10]. En utilisant la convexité des ensembles du point (d) de (A), on peut alors déduire que la densité est
homogène aux bornes du théorème. Enfin, le contrôle optimal dans (5.6–5.7) s’écrit explicitement.
(ii). Pour se rapprocher le plus possible du cas précédent, l’idée naturelle est ensuite de linéariser l’équation
contrôlée (5.6). Le point crucial est celui du choix de la courbe déterministe (φt )t∈[0,T ] autour de qui opérer
cette linéarisation. On écrit F(t, χt0,x0 ) = F(t, φt ) + Dx F(t, φt )(χt0,x0 − φt ) + o(|χt0,x0 − φt |), où Dx F ∈ Mnd (R)
est la dérivée spatiale de F. De même nous approchons le coefficient de diffusion σ(t, χt0,x0 ) par σ(t, φt ). Comme
(φt )0≤t≤T est déterministe, l’application z ∈ Rnd 7→ F(t, φt ) + Dx F(t, φt )(z − φt ) est bien affine, le coefficient
de diffusion (σ(t, φt ))0≤t≤T est déterministe, et l’on retrouve le cadre de (i). Un choix naturel consiste à prendre
(φt )0≤t≤T comme la solution du problème de contrôle déterministe pour l’ODE naturellement associée à (5.2),
i.e.
φ̇t = F(t, φt ) + Bϕt , t ∈ [0, T ]; φ0 = x0 , (5.8)
où (ϕt )0≤t≤T est un contrôle déterministe de L2 ([0, T ], Rd ) solution du problème de minimisation :
Z T

I(T, x0 , y0 ) := inf |ϕt |2 dt : φ0 = x0 , φT = y0 . (5.9)
0

Des arguments usuels de contrôle, cf. Coron [Cor07] nous ont permis d’établir que I(T, x0 , y0 ) est finie,
i.e. l’équation (5.8) est contrôllable. Nous avons également établi que I(T, x0 , y0 ) est du même ordre que
T |T−1 2
T (θ T (x0 ) − y0 )| . La comparaison avec l’énoncé du Théorème 5.2.1, nous donne donc que l’énergie du
contrôle optimal déterministe va fournir la décroissance hors diagonale de la densité. De cette étape on décide
de choisir le contrôle vt dans (5.6) sous la forme vt := ϕt + vt0 .

49
 (iii). Il reste donc à choisir (vt0 )t∈[0,T −ε] . Pour cela remarquons que d’après (5.6) et (5.8), la dynamique de
(χ0,x
t − φt )0≤t≤T −ε s’écrit :
   
d χ0,x
t − φt = Dx F(t, φt )(χ0,x t − φt ) + B(vt − ϕt ) dt + Bσ(t, φt )dWt + dRt , (5.10)

où (Rt )0≤t≤T −ε désigne un terme de reste qu’une analyse de stabilité technique, car non linéaire, permet
de contrôler. Observons que χ0,x 0 − φ0 = 0 et que χ0,x T −ε − φT −ε est a priori proche de y0 − y0 = 0. Ces
observations heuristiques font qu’un candidat naturel pour v 0 est le contrôle optimal Gaussien qui partant de 0
permet d’atteindre 0 en T . Par l’étape (i), nous savons également calculer son coût (au sens de la formulation
2
(5.7)) : celui-ci est à une constante additive près de l’ordre de − ln(T −n d/2 ) et correspond dans le Théorème
5.2.1 au contrôle de la densité sur la “diagonale”.
En résumé, de la formulation contrôle (5.7) on choisit un contrôle particulier vt = ϕt + vt0 , où (ϕt )t∈[0,T −ε]
permet d’isoler les contributions hors diagonale, (vt0 )t∈[0,T ] donne la vitesse diagonale. C’est l’étape de contrôle
des restes dans (5.10) qui est délicate. En effet, nous devons à ce propos utiliser des arguments de type Gron-
wall non linéaire qui font que l’on établira d’abord les bornes du théorème en temps court. Un argument de
convolution ou de changement d’échelle permet ensuite de les étendre à temps arbitraire fixé.
La fonction I(T, x, y) introduite dans (5.9) est la fonctionnelle d’action des grandes déviations. Elle établit
un lien naturel entre le problème de contrôle déterministe (5.8–5.9) et la densité de transition de (5.2), cf.
Freidlin and Wentzell [FW98] pour le cas non-dégénéré et Ben Arous et Léandre [BL91] pour des résultats
sous condition de Hörmander forte. L’utilisation de la transformée de Fleming, introduite par Sheu [She91] et
développée dans cette section permet une extension de ce lien lorsque le temps ne tend pas forcément vers 0.
Comme indiqué ci-avant, nous conservons à cause de l’étude de stabilité une contrainte de temps court.

5.2.3 Parametrix et majoration

La formulation (5.7) n’est pas la plus adaptée pour obtenir une majoration. En effet, il faut minorer Jε pour
espérer majorer la densité. L’expression en terme d’infimum oblige alors à partir du contrôle optimal.
En revanche la structure du transport intervenant dans le terme exponentiel donne une bonne idée pour
construire une famille de processus Gaussiens de densités homogènes à celle du Théorème 5.2.1 et qui serviront
de parametrix.
Notons par p̃T,y0 la densité Gaussienne qui en temps petit et y0 doit être “proche” de p(T, ·, y0 ). L’exposant
(T, y0 ) de p̃T,y0 indique que moyenne et covariance de p̃ dépendent du temps final et du point d’arrivée où l’on
considère la densité. Comme pour la minoration l’idée est de linéariser la dynamique de la diffusion (5.2) autour
d’une courbe déterministe. Comme c’est la variable spatiale d’arrivée qui est connue on va développer autour
de la solution de l’EDO rétrograde
d

θt,T (y0 ) = F t, θt,T (y0 ) , t ∈ [0, T ]; θT,T (y0 ) = y0 .
dt
Le modèle Gaussien associé s’écrit :


dX̃t = F(t, θ t,T (y0 )) + Dx F(t, θt,T (y0 )) X̃t − θt,T (y0 ) dt + Bσ(t, θ t,T (y0 ))dWt , t ∈ [0, T ].

Désignons par (L̃tT,y0 )0≤t≤T le générateur de X̃ gelé au temps T et au point y0 . On rappelle que la quantité
fondamentale pour le développement parametrix est le noyau
 T,y0  T,y0

H(t, T, x, y0 ) = Lt,x − L̃t,x p̃ (t, T, x, y0 ) ,

cf. [MS67] ou Chapitre 3. Nous obtenons pour tout N ≥ 1 le développement :

N Z
X T Z
p(T, x0 , y0 ) = p̃(0, T, x0 , y0 ) + p̃(0, t, x0 , z)H ⊗k (t, T, z, y0 )dzdt
k=1 0 Rnd
Z Z (5.11)
T
+ p(t, x0 , z)H ⊗(N +1) (t, T, z, y0 )dzdt,
0 Rnd

avec p̃(0, t, x0 , z) := p̃t,z (0, t, x0 , z) et où l’on rappelle

Z T Z
H ⊗(k+1) (t, T, z, y0 ) := H(t, s, z, z′ )H ⊗k (s, T, z′ , y0 )dz′ ds.
t Rnd

50
Malheureusement nous n’arrivons pas à faire tendre N vers l’infini dans (5.11). A cause du transport non
linéaire θ, nous n’avons pas réussi à contrôler les noyaux itérés (H ⊗k )k≥1 uniformément en k. Nous renvoyons
au Chapitre 3 pour le cas linéaire auquel la démarche usuelle s’applique.
Nous devons donc ici tronquer la série. Le point clé est que la formulation (5.7) permet d’établir que, pour
N suffisamment grand, le terme de reste
Z TZ
p(t, x0 , z)H ⊗(N +1) (t, T, z, y0 )dzdt
0 Rnd

dans (5.11) admet une borne Gaussienne dont on déduit l’estimée énoncée. Cette technique est à notre connais-
sance nouvelle.

5.2.4 Une seconde approche pour la minoration : parametrix et chaı̂nage

Nous avons développé dans [10] une deuxième méthode pour obtenir la minoration du Théorème 5.2.1. A
partir du développement (5.11), en exploitant la propriété régularisante en temps du noyau de convolution H et
les contrôles utilisés pour la borne supérieure, nous pouvons déduire dans l’esprit de [IKO62] la borne inférieure
sur les “compacts” de la métrique en temps petit. Il existe T0 > 0 tel que pour tout (T, x, y) ∈ (0, T0 ] × (Rnd )2
−n2 d/2
lorsque la métrique dT (x, y) := T 1/2 |T−1
T (θ T (x) − y)| ≤ C0 , p(0, T, x, y) ≥ C1 T .
Il reste ensuite à l’obtenir pour dT (x, y) > C0 . La propriété de Markov donne
Z N
Y −2
p(T, x, y) := p(0, t1 , x, z1 ) p(ti , ti+1 , zi , zi+1 )p(tN −1 , T, zN −1 , y)dz1 ...dzN −1 , (5.12)
(Rnd )N −1 i=1

où t0 = 0 < t1 < t2 < ... < tN = T, N ∈ N∗ .

L’idée est ensuite de choisir N , les (ti )i∈[[0,N ]] et des ensembles (Bi )i∈[[0,N ]] (avec les conventions B0 =
{x}, BN = {y}) tels que pour (zi , zi+1 ) ∈ Bi × Bi+1 , i ∈ [[0, N − 1]], l’on puisse appliquer la borne sur les
compacts de la métrique pour p(ti , ti+1 , zi , zi+1 ). Si (ti+1 −ti )1/2 |T−1
ti+1 −ti (θ ti+1 ,ti (zi )−zi+1 )| ≤ C0 , où θ̇ s,ti (zi ) =
2
F(s, θs,ti (zi )), s ≥ ti , θ ti ,ti (zi ) = zi , alors p(ti , ti+1 , zi , zi+1 ) ≥ C1 /(ti+1 − ti )n d/2
. De (5.12) on écrit
Z N
Y −2
p(T, x, y) ≥ QN −1
p(0, t1 , x, z1 ) p(ti , ti+1 , zi , zi+1 )p(tN −1 , T, zN −1 , y)dz1 ...dzN −1
i=1 Bi i=1
N
Y −1
C1 C1
≥ n2 d/2
|Bi |.
t1 i=1
(ti+1 − ti )n2 d/2
2
Si maintenant ti+1 − ti ≤ T /N, |Bi | ≥ C2 (T 1/2 dT (x, y)/N )n d , ∀i ∈ [[0, N − 1]], il vient
 n2 d/2  (N −1)n2 d/2
N d2T (x, y)
p(T, x, y) ≥ C1N C2N −1 .
T N
C3N
Pour N = ⌈Kd2T (x, y)⌉ ≥ 1, nous obtenons p(t, x, y) ≥ T n2 d/2
, où C3 = C1 C2 /(K + 1) < 1 pour K suffisamment
grand. On retrouve donc la borne inférieure en temps petit à partir de l’écriture p(T, x, y) ≥ exp(N ln(C3 ))
T −n2 d/2
en
utilisant la définition de N . L’extension à un temps arbitraire se fait de nouveau par convolution ou changement
d’échelle.
La construction des (ti )i∈[[0,N ]] , (Bi )i∈[[0,N ]] vérifiant les conditions précédentes, est présentée en appendice
de [10]. L’idée est d’interpréter la construction usuelle dans le cas elliptique, cf. Chapitre 7 de Bass [Bas97], en
terme de problème de contrôle déterministe de sorte à pouvoir le transposer à notre cadre.
Pour le cas non dégénéré, les ensembles Bi sont des boules dont les centres sont uniformément répartis sur
la ligne droite reliant point de départ et d’arrivée. La ligne droite est ici la solution optimale du problème de
contrôllabilité déterministe, φt = x+ Tt (y−x). Le déplacement y−x T est constant ce qui explique l’équirépartition.
L’idée dans notre cas dégénéré consiste à disposer le centre des ensembles Bi le long de la trajectoire optimale
(φt )t∈[0,T ] solution de (5.8). A une localisation technique près 3 on définit t0 = 0, ∀i ≥ 1, ti := inf{t ≥ ti−1 :
Rt
ti−1
|ϕs |2 ds = I(T, x, y)/N }, i.e. on considère une équirépartition en terme de niveau d’énergie. Les centres des
(Bi )i∈[[1,N −1]] sont associés à (φti )i∈[[1,N −1]], la géométrie de ces ensembles est liée à la métrique naturelle.
3. qui consiste à prendre l’infimum avec T /N dans la définition ci-après en s’assurant que le nombre de pas de temps reste de
l’ordre de N .

51
5.2.5 Extensions et perspectives
Ce travail ouvre de nombreuses perspectives. Tout d’abord le Théorème 5.2.1 est énoncé sous réserve d’uni-
cité en loi. Si les coefficients de (5.2) vérifient (A) avec η > 1/2 nous obtenons cette identité en loi comme
corollaire d’un principe de comparaison de Ishii et Lions [IL90] pour les solutions de viscosité d’équations non
linéaires. Cet argument n’est pas très approprié dans notre cas dans la mesure où il n’utilise pas du tout l’aspect
régularisant du noyau de la chaleur. Par ailleurs, dans les preuves aucune limite de ce type sur η n’est apparue.
Dans un travail en cours 4 , je m’intéresse au problème de martingale pour (5.2), dans le cas où a est uni-
formément elliptique et continu en espace, en gardant inchangées les autres hypothèses dans (A). Comme dans
le cas non dégénéré, cf. [SV79], l’idée est d’établir des inégalités de Calderón-Zygmund pour le modèle Gaussien.
Une des difficultés est dans notre cadre d’introduire une bonne métrique parabolique qui tienne compte de l’as-
pect multi-échelle. De l’unicité du problème de martingale on déduirait immédiatement que le Théorème 5.2.1
est valable pour η ∈ (0, 1]. Les estimées de Calderón-Zygmund permettraient également d’aborder des EDP
dégénérées semi-linéaires. J’essaie également de donner des contrôles ponctuels sur les dérivées de la densité à
partir de l’analyse précédente.

En terme d’applications, on peut penser aux systèmes Hamiltoniens, deux problèmes associés à l’équation
(5.2) considérée apparaissent naturellement : le développement asymptotique de la densité en temps petit,
et l’estimation du régime stationnaire, ainsi que de la vitesse de convergence vers ce dernier, sous de bonnes
conditions de stabilité. Pour l’estimée en temps petit, on s’attend à ce que le contrôle optimal dans la formulation
(5.7) soit proche du contrôle optimal déterministe dans (5.8). Il faudrait raffiner l’analyse de stabilité pour
l’établir rigoureusement. Adapter la démarche précédente à l’estimation en temps long sous conditions de type
Lyapunov reste à faire. Une idée pourrait être de choisir comme contrôle dans (5.7) un contrôle optimal explicite
associé à un processus “canonique” pour le régime limite considéré, i.e. Ornstein-Ulhenbeck si la force de rappel
est à croissance linéaire.
Si l’on rajoute dans (5.2) un bruit sur la dernière composante on retrouve le cadre des modèles microscopiques
de mécanique statistique pour la diffusion de chaleur. Dans un modèle linéaire d’interaction aux plus proches
voisins comme celui de Bodineau et Lefevere [BL08], il serait intéressant de reprendre l’analyse précédente et
d’essayer d’obtenir des bornes de déviations non asymptotiques en temps. Nous avons prévu avec Raphaël Le-
fevere de discuter à ce propos.

Nous avons considéré ici un régime Gaussien qui permettait de faire apparaı̂tre naturellement, et avec de
bonnes estimées, le lien avec le contrôle déterministe. Une extension naturelle consiste à considérer d’autres
régimes, i.e. lorsque des crochets de Lie d’ordre strictement supérieur à 1 interviennent
Rt pour engendrer l’espace.
L’exemple le plus simple consiste à étudier la densité du couple (x + Wt , 0 (x + Ws )k ds)t≥0 où (Wt )t≥0 est le
mouvement Brownien réel standard.
Pour k = 2 des formules “explicites” existent, cf. formule 1.9.8 de Borodin et Salminen [BS03] mais ne sont
pas directement exploitables pour déduire des contrôles à la Aronson sur la densité. Pour x = 0 et k impair
Chiara Cinti et Sergio Polidoro, des Universités de Bologne et Modena, ont obtenu dans [CP08] une borne
inférieure. Les outils sont ceux de la théorie du potentiel pour les opérateurs sous-elliptiques, cf. Bonfiglioli et
al. [BLU07]. Nous avons avec ces auteurs étendu ces résultats pour k quelconque. Ce travail est en cours de
rédaction.
Par ailleurs, avec V. Konakov de l’Université de Moscou, nous avons montré par des techniques de calcul de
Malliavin que ces estimées sont précises, i.e. qu’il existe une borne supérieure du même ordre. Pour x1 = 0 et k
quelconque nous obtenons ∃c ∈ (0, 1], C := C(T ) ≥ 1 telles que pour tout (t, (0, x2 ), (y1 , y2 )) ∈ (0, T ] × ({0} ×
R) × R2 :
     
C −1 |y1 |2 |y2 − x2 |2/k C |y1 |2 |y2 − x2 |2/k
exp −c−1 + ≤ p(t, (0, x2 ), (y ,
1 2y )) ≤ exp −c + .
t(k+3)/2 t t1+2/k t(k+3)/2 t t1+2/k

Ce travail est également en cours. L’idée est de bien comprendre la structure du transport de la condition initiale
pour x 6= 0. Nous y arrivons lorsque |y2 − x2 | ≥ t|x1 |k , i.e. lorsque l’on est en grand régime par rapport au
2 k 2/k
transport déterministe. Le terme dans l’exponentielle précédente s’écrit alors |y1 −x
t
1|
+ (|y2 −x2t1+2/k
|−t|x1 | )
.
Le transport en petits régimes reste à analyser. Une piste peut être d’utiliser des inégalités de Tcheby-
chev exponentielles à l’aide des transformées de Laplace, voir Mansuy et Yor [MY09] pour le cas k = 2. La
caractérisation complète, i.e. pour tous les régimes, est une étape nécessaire en vue d’appliquer une méthode
paramétrix pour remplacer le Brownien par une diffusion.
4. [Link] “Well posedness of the martingale problem for Kolmogorov’s equations”

52
5.3 Estimation de densité et concentration
Nous présentons ici une application des estimées de type Aronson : l’obtention d’intervalles de confiance non
asymptotiques pour la méthode de Monte Carlo. La domination Gaussienne permet en effet d’utiliser les outils
usuels en concentration de la mesure, inégalité de Sobolev logarithmique et argument de Herbst (cf. Ledoux
[Led99]), pour déduire des inégalités de déviation non asymptotiques.
Nous avons abordé le cas des diffusions :
Z t Z t
Xt = x + b(s, Xs )ds + σ(s, Xs )dWs , (5.13)
0 0

où (Xt )t≥0 est à valeurs dans Rd et (Wt )t≥0 est un mouvement Brownien de Rd . Par la suite nous ferons
référence à l’équation précédente en parlant du cas (a). Nous avons également considéré :
Z t Z t
Xt1 = x1 + b(s, Xs )ds + σ(s, Xs )dWs ,
0 0
Z t
2
Xt = x2 + Xs1 ds, (5.14)

où Xt = (Xt1 , Xt2 ) est à valeurs dans Rd × Rd . L’équation (5.14) est un cas particulier de (5.2) pour n = 2 et
F (x1 , x2 ) = (b(x1 , x2 ), x1 )∗ . Nous en parlerons comme du cas (b).
Les coefficients b, σ des équations ci-dessus vérifient les hypothèses de la partie opérateurs sous forme non
divergence de la Section 5.1 dans le cas non dégénéré (a), et (H) de la Section 5.2 en supposant F1 = b borné
pour le cas dégénéré (b).
Dans le cas (a), nous avons considéré le schéma d’Euler de pas de temps h := T /N, N ∈ N∗ . Posons
∀i ∈ N, ti = ih, et pour t ≥ 0, φ(t) := sup{ti : ti ≤ t}. La dynamique s’écrit :
Z t Z t
Xth = x + h
b(φ(s), Xφ(s) )ds + h
σ(φ(s), Xφ(s) )dWs . (5.15)
0 0

Pour (b) nous avons défini :

Z t Z t
Xth,1 = x1 + h
b1 (φ(s), Xφ(s) )ds + h
σ(φ(s), Xφ(s) )dWs ,
0 0
Z t
Xth,2 = x2 + Xsh,1 ds, (5.16)
0

où Xth := (Xth,1 , Xth,2 ). L’équation (5.16) definit un schéma complètement simulable à incréments Gaussiens.
Sur chaque pas de temps, les d dernières composantes sont l’intégrale d’un processus Gaussien.
L’erreur faible pour ces modèles est bien comprise, voir Chapitre 1 et 3. Il reste donc à contrôler, avec les
notations de l’équation (1.5) :
M
S 1 X
EM,h := f (T, (XTh )i ) − Ex [f (T, XTh )], (5.17)
M i=1


où T > 0 est fixé et (XTh )i i∈N∗ est une suite i.i.d de réalisations du schéma d’Euler (5.15) dans le cas (a) et
(5.16) pour le cas (b).
Dans le cadre ergodique, le contrôle non asymptotique de l’erreur a été étudié par Malrieu et Talay [MT06a].
Les auteurs ont réussi à établir que le schéma d’Euler d’un processus de diffusion multi-dimensionnel à coefficient
de diffusion constant vérifie une inégalité de Sobolev logarithmique, qui implique la concentration gaussienne
(voir e.g. Ledoux [Led99]). La borne de concentration non asymptotique en découle aisément. Toujours dans le
contexte ergodique, mentionnons le travail de Joulin et Ollivier pour les chaı̂nes de Markov [JO09].
Notre but est ici différent. Nous souhaitons donner un contrôle pour T > 0 et ne pouvons appliquer di-
rectement la machinerie log-Sobolev, par ailleurs assez rigide pour le passage aux discrétisations, directement.
Sous les hypothèses précédentes nous avons néanmoins pu établir des bornes Gaussiennes dans les cas (a) et
(b). L’idée pour en dériver une borne non asymptotique est d’utiliser ensuite que la densité Gaussienne qui do-
mine vérifie une inégalité de Sobolev logarithmique. On peut alors modifier l’argument de Herbst pour obtenir
l’inégalité de déviation souhaitée. Il faut toutefois gérer la constante qui apparaı̂t devant la densité Gaussienne
dans l’estimée d’Aronson suivante.

53
Théorème 5.3.1 (Estimées d’“Aronson” Gaussiennes pour le schéma d’Euler) Il existe c ∈ (0, 1], C ≥
1, C := C(T, d, Λ, |b|), c := c(d, Λ), tels que pour 0 ≤ j < j ′ ≤ N :

C −1 pc−1 (tj ′ − tj , x, x′ ) ≤ ph (tj , tj ′ , x, x′ ) ≤ Cpc (tj ′ − tj , x, x′ ), (5.18)

où ph désigne la densité du schéma d’Euler. Pour 0 ≤ s < t ≤ T , dans le cas (a), pc (t − s, x, x′ ) :=
 d/2 ′
−x|2

c
2π(t−s) exp −c |x
2(t−s) et dans le cas (b)

√ !d   ′ x +x′ 
′ 3c |x1 − x1 |2 |x′2 − x2 − 1 2 1 (t − s)|2
pc (t − s, x, x ) := exp −c +3 .
2π(t − s)2 4(t − s) (t − s)3

Pour établir ce théorème, les techniques sont les mêmes qu’auparavant : i.e. représentation parametrix discrète
de la densité des schémas pour obtenir la borne supérieure, et méthode de chaı̂nage adaptée au cas semi-
Markovien du schéma d’Euler pour la minoration. Néanmoins, ce résultat est à ma connaissance le premier sous
les hypothèses usuelles précédentes. Il peut dans le cas non dégénéré se déduire de Gobet et Labart [GL08] sous
des hypothèses plus fortes de régularité de b, σ.
Nous obtenons ensuite le résultat de concentration.
Théorème 5.3.2 (Concentration Gaussienne ) Considérons pour le cas (a) une fonction f : Rd → R
Lipschitzienne en espace et mesurable  bornée en temps,
 vérifiant |∇f (T, .)|∞ ≤ 1. Pour le cas (b) on suppose
−1 −1 Id 0
que f (T, x) = g(T, TT x) où TT = et g est une fonction Lipschitzienne en espace et mesurable
0 T −1 Id
bornée en temps, vérifiant |∇g(T, .)|∞ ≤ 1.
Pour les constantes c, C du Théorème 5.3.1, pour tout h = T /N, N ≥ 1,
 
 S  M 2
∀r > 0, ∀M ≥ 1, Px |EM,h | ≥ r + δC,α(T ) ≤ 2 exp − r ,
α(T )

où (
c
1 2T dans le cas (a),
= √ c
α(T ) (4 − 13) 2T dans le cas (b),
p
et δC,α(T ) = 2 α(T ) log C.
Lorsque de plus f tend vers l’infini dans une direction, nous établissons une borne inférieure du même ordre.
Nous renvoyons à [10] pour un énoncé complet du théorème.
Le terme correctif δC,α(T ) correspond à une borne supérieure pour la distance de Wasserstein entre les lois
de densités ph (T, x, .) et pc (T, x, .) lorsque celles-ci vérifient la condition de domination de l’équation (5.18). Si
C = 1 ce terme correctif disparaı̂t. Dans ce cas là on est directement Gaussien. C’est en revanche le prix à payer
dans le cas général. Ce terme est fixé une fois pour toutes. Pour M grand, l’ordre de déviation typique donné par
le théorème central limite est r = √cM 0
pour c0 > 0 dépendant d’un seuil fixé et de la variance de f (T, XTh ). Du
fait de la technique employée, i.e. domination par une densité qui vérifie une inégalité de log-Sobolev, notre borne

1/2
n’est pas pertinente pour M tendant vers l’infini à T fixé. En revanche elle l’est tant que r = √cM 0
≥ 2T c

1/2
et en particulier lorsque r = √cM 0
et 2Tc sont du même ordre. Ainsi, si le temps T est “petit”, la borne
non asymptotique restera valable jusqu’à M assez grand, relié à T par l’inégalité précédente. Un autre cas
intéressant est celui où pour T > 0 fixé il est difficile de simuler un nombre important de trajectoires de Monte
Carlo. On peut penser pour illustrer ce dernier cas à un schéma d’Euler dont l’évaluation des coefficients ferait
appel à des codes de calculs lourds.

54
Annexe A

Unicité en loi pour certains systèmes

dégénérés

Nous indiquons dans cette section les idées principales de la note [12] qui établit un lien entre la méthode
parametrix à la Mc Kean et Singer présentée ci-avant et la technique employée par Bass et Perkins [BP09] pour
prouver l’unicité du problème de martingale pour un opérateur sous forme non divergence non dégénéré.
Cette identification et les contrôles obtenus dans l’article [10] permettent ensuite de prouver l’unicité du
problème de martingale associé à l’équation dégénérée (5.2) sous les hypothèses (A) de la Section 5.2 pour tout
η ∈ (0, 1]. Ceci étend la validité du Théorème 5.2.1, comme nous l’espérions (cf. Section 5.2.5).

A.1 Rappel sur le modèle et la construction du noyau H

Nous rappelons la dynamique de Xt qui sous forme compacte s’écrivait :

dXt = F(t, Xt ) + Bσ(t, Xt )dWt , (A.1)

La linéarisation utilisée dans la méthode parametrix pour la majoration de la densité faisait intervenir pour
T > 0, y ∈ Rnd l’équation linéaire :


dX̃T,y
t = F(t, θt,T (y)) + DF(t, θ t,T (y)) X̃T,y
t − θ t,T (y) dt
(A.2)
+ Bσ(t, θ t,T (y))dWt , 0 ≤ t ≤ T,

 t,T (y) = F(t, θ t,T (y)), t ≥ 0, avec condition terminale

où (θ t,T (y))t≥0 est solution de l’EDO [d/dt]θ  θT,T (y) = y
0 ··· ··· ··· 0
 Dx1 F2 (t, x) 0 ··· ··· 0 
 
 .. 
nd
et ∀(t, x) ∈ [0, T ] × R , DF(t, x) =   0 Dx2 F3 (t, x) 0 0 . 
 est la sous-
 .. . .. 
 . 0 . . . 
0 ··· 0 Dxn−1 Fn (t, x) 0
diagonale de la matrice Jacobienne Dx F. Désignons par p̃T,y (t, T, x, .) la densité de X̃T,y T partant de x en
t.
Pour 0 ≤ s < t ≤ T, z, y ∈ Rnd , définissons le noyau de convolution H comme :
1
H(s, t, z, y) = {hF(s, z), Dz p̃t,y (s, t, z, y)i + Tr(a(s, z)Dz21 p̃t,y (s, t, z, y))} −
2
{hF(s, θs,t (y)) + DF(s, θ s,t (y))(z − θs,t (y)), Dz p̃t,y (s, t, z, y)i
1
+ Tr(a(s, θ s,t (y))Dz21 p̃t,y (s, t, z, y))} = (Ls,z − L̃t,y
s,z )p̃
t,y
(s, t, z, y), (A.3)
2

où L est le générateur de X et L̃t,y

s,z et p̃
t,y
(s, t, z, .) désignent respectivement le générateur à l’instant s et
t,y t,y
la densité à l’instant t de X̃ , X̃s = z avec coefficients “gelés” en t, y. L’indice z précise le paramètre de
différentiation.
Nous avons établi dans [10] (cf. Lemma 5.5)

55
Lemme A.1.1 Soit T0 > 0 fixé. Il existe une constante CA.1.1 > 0, dépendant de (A) et T0 t.q. pour tout
t ∈ [0, T ), T ≤ T0 et z, y ∈ Rnd ,
η

|H(t, T, z, y)| ≤ CA.1.1 (T − t) 2 −1 gCA.1.1 ,T −t z − θt,T (y) ,
2d
où pour tout a > 0, t > 0, ga,t (y) = t−n 2 exp(−a−1 t|T−1 2
t y| ), a, t > 0, y ∈ Rnd .
Par le même type de techniques que celles utilisées pour prouver le lemme précédent, on peut également
établir :

Lemme A.1.2 Soit une fonction h ∈ C 0 (Rnd , R). Définissons pour tout (s, x) ∈ [0, T ] × Rnd ,
Z
∀ε > 0, Ξε (s, x) := dyh(s, y)p̃s+ε,y (s, s + ε, x, y).
Rnd

Alors Ξε (s, x) converge de façon bornée et ponctuellement vers h(s, x) lorsque ε → 0.

A.2 Unicité du problème de martingale sous (A)

A.2.1 Résultats principaux
Le résultat principal de la note est le suivant :
Théorème A.2.1 Sous (A) et pour η ∈ (0, 1], le problème de martingale associé au générateur L dans (A.3)
est bien posé. En particulier, on a l’unicité en loi pour l’EDS (A.1).
Nous en déduisons le corollaire suivant.
Corollaire A.2.1 Sous (A) et pour η ∈ (0, 1], l’unique solution faible de (A.1) admet pour tout t > 0 une
densité qui vérifie les estimées d’Aronson du théorème 5.2.1.

A.2.2 Preuve du théorème A.2.1

Supposons que l’on ait deux solutions P1 , P2 du problème de martingale associé à (Lt )t∈[s,T ] partant de x
en s. Nous pouvons sans perte de généralité supposer T ≤ 1. Pour une fonction f : [0, T ] × Rnd → R borélienne
bornée, définissons : "Z #
T
i
S f := Ei dtf (t, Xt ) , i ∈ [[1, 2]],
s

où (Xt )t∈[s,T ] désigne ici le processus canonique associé à (Pi )i∈[[1,2]] . Précisons (comme il est mentionné dans
[BP09]) que Sf est seulement une fonctionnelle linéaire et non une fonction dans la mesure où Pi ne provient
pas forcément d’un processus de Markov. Posons :

S ∆ f := S 1 f − S 2 f, Θ := sup |S ∆ f |.
kf k∞ ≤1

Remarquons que Θ ≤ T − s.
Si f ∈ C02 ([0, T ) × Rnd , R), par définition du problème de martingale nous obtenons :
"Z #
T
f (s, x) + Ei dt(∂t + Lt )f (t, Xt ) = 0, i ∈ [[1, 2]]. (A.4)
s

Pour y ∈ Rnd fixé et ε ≥ 0, introduisons ∀f ∈ C01,2 ([0, T ) × Rnd , R) la fonction de Green

Z T Z
nd ε,y
∀(s, x) ∈ [0, T ) × R , G f (s, x) := dt dzp̃t+ε,y (s, t, x, z)f (t, z). (A.5)
s Rnd

On vérifie aisément que

(∂s + L̃t+ε,y
s,x )p̃
t+ε,y
(s, t, x, z) = 0, ∀(s, x, z) ∈ [0, t) × (Rnd )2 , p̃t+ε,y (s, t, x, .) −→ δx (.). (A.6)
t↓s

56
Si l’on définit pour toute fonction f ∈ C01,2 ([0, T ) × Rnd , R), (s, x) ∈ [0, T ) × Rnd ,
Z T Z
ε,y
Ms,x f (s, x) := dt dzL̃t+ε,y
s,x p̃
t+ε,y
(s, t, x, z)f (t, z),
s Rnd

nous obtenons des équations (A.5), (A.6) :

(∂s Gε,y f + Ms,x

ε,y
f )(s, x) = −f (s, x), ∀(s, x) ∈ [0, T ) × Rnd . (A.7)

Introduisons désormais pour une fonction régulière h ∈ C01,2 ([0, T ) × Rnd , R) et pour tout (s, x) ∈ [0, T ) × Rnd :

Φε,y (s, x) := p̃s+ε,y (s, s + ε, x, y)h(s, y),

Z
Ψε (s, x) := dyGε,y (Φε,y )(s, x).
Rnd

L’identité (A.5) donne :

Z Z T Z
Ψε (s, x) := dy dt dzp̃t+ε,y (s, t, x, z)Φε,y (t, z)
Rnd s Rnd
Z Z T Z
= dy dt dzp̃t+ε,y (s, t, x, z)p̃t+ε,y (t, t + ε, z, y)h(t, y)
Rnd s Rnd
Z Z T
= dy dtp̃t+ε,y (s, t + ε, x, y)h(t, y), (A.8)
Rnd s

en utilisant la propriété de semi-groupe de la densité gelée. Ecrivons pour tout (s, x) ∈ [0, T ) × Rnd ,
Z
(∂s + Ls )Ψε (s, x) = dy(∂s + Ls )(Gε,y Φε,y )(s, x)
Rnd
Z Z
= dy(∂s Gε,y Φε,y + Ms,x
ε,y ε,y
Φ )(s, x) + dy(Ls Gε,y Φε,y − Ms,x
ε,y ε,y
Φ )(s, x)
R nd R nd
Z Z
= − dyΦε,y (s, x) + dy(Ls Gε,y Φε,y − Ms,xε,y ε,y
Φ )(s, x) := I1ε + I2ε ,
Rnd Rnd

en utilisant (A.7) pour l’avant dernière égalité. Le Lemme A.1.2 donne I1ε −→ −h(s, x). Par ailleurs,
ε→0

Z T Z
I2ε = dt dy(Ls − L̃t+ε,y
s,x )p̃
t+ε,y
(s, t + ε, x, y)h(t, y)
s Rnd
Z T Z
= dt dyH(s, t + ε, x, y)h(t, y).
s Rnd

Du Lemme A.1.1,
Z T Z
|I2ε | ≤ CA.1.1 khk∞ dt(t + ε − s)−1+η/2 dygCA.1.1 ,t+ε−s (x − θs,t+ε (y))
s Rnd
Z T Z
−1+η/2
≤ khk∞ CA.9 dt(t + ε − s) dygCA.9 ,t+ε−s (θ t+ε,s (x) − y) ≤ CA.9 {(T − s) ∨ ε}η/2 khk∞ ,
s Rnd
(A.9)

en utilisant la propriété de bi-Lipschitz continuité du flot pour l’avant dernière inégalité, en modifiant CA.9 dans
la dernière (cette constante modifiée ne dépend toutefois que de (A)). Ainsi, pour T suffisamment petit nous
avons par (A.9)
1
|I2ε | ≤ khk∞ . (A.10)
2
L’équation (A.4) et la définition précédente de S ∆ donnent :

S ∆ ((∂s + Ls )Ψε ) = 0 ⇒ |S ∆ I1ε | = |S ∆ I2ε |.

57
Par la partie convergence bornée du Lemme A.1.2 et (A.10), nous obtenons :

khk∞
|S ∆ h| = lim |S ∆ I1ε | = lim |S ∆ I2ε | ≤ Θ lim sup |I2ε | ≤ Θ .
ε→0 ε→0 ε→0 2

Par un argument de classe monotone, l’inégalité précédente reste vraie pour toute fonction h mesurable bornée,
à support compact dans [0, T ) × Rnd
hR. En prenanti le supremum
hR sur khki∞ ≤ 1, on arrive à Θ ≤ 21 Θ ce qui donne
T T
Θ = 0 comme Θ < +∞. Ainsi, E1 s dth(t, Xt ) = E2 s dth(t, Xt ) ce qui prouve le résultat sur l’intervalle
[0, T ]. Les probabilités conditionnelles régulières permettent d’étendre le résultat sur R+ (voir Chapitre 6.2 de
Stroock et Varadhan [SV79]).

58
Bibliographie

[Als94] G. Alsmeyer. On the Markov renewal theorem. Stoch. Proc. Appl., 50(1) :37–56, 1994.
[Ant93] F. Antonelli. Backward-Forward Stochastic Differential Equations. Ann. Appl. Probab., 3-3 :777–793,
1993.
[Aro67] D. G. Aronson. Bounds for the fundamental solution of a parabolic equation. Bull. Amer. Math.
Soc., 73 :890–896, 1967.
[Bal90] V. Bally. On the connection between the Malliavin covariance matrix and Hörmander’s condition.
J. Funct. Anal., 96–2 :219–255, 1990.
[Bal95] P. Baldi. Exact asymptotics for the probability of exit from a domain and applications to simulation.
Ann. Prob., 23(4) :1644–1670, 1995.
[Bas97] R. F. Bass. Diffusions and Elliptic Operators. Springer, 1997.
[BCI99] P. Baldi, L. Caramellino, and M.G. Iovino. Pricing general barrier options : a numerical approach
using sharp large deviations. Math. Finance, 9(4) :293–322, 1999.
[BDH+ 03] P. Briand, B. Delyon, Y. Hu, E. Pardoux, and L. Stoica. Lp solutions of backward stochastic
differential equations. Stoch. Proc. Appl., 108 :109–129, 2003.
[BGJ87] K. Bichteler, J.B. Gravereaux, and J. Jacod. Malliavin Calculus for processes with jumps. Stochastics
Monographs, volume 2, 1987.
[BGK97] M. Broadie, P. Glasserman, and S. Kou. A continuity correction for discrete barrier options. Ma-
thematical Finance, 7 :325–349, 1997.
[BGK99] M. Broadie, P. Glasserman, and S. Kou. Connecting discrete and continuous path-dependent options.
Fin. and Stoch., 3 :55–82, 1999.
[BL82] A. Bensoussan and J. L. Lions. Contrôle impulsionnel et inéquations quasi variationnelles. Gauthiers
Villars, 1982.
[BL91] G. Ben Arous and R. Léandre. Décroissance exponentielle du noyau de la chaleur sur la diagonale.
II. Prob. Th. Rel. Fields, 90(3) :377–402, 1991.
[BL97] G. Barles and E. Lesigne. SDE, BSDE and PDE. Backward stochastic differential equations (Paris
95-96). Stat. and Prob. Letters, Pitman research, 364 :47–80, 1997.
[BL08] T. Bodineau and L. Lefevere. Large deviations of lattice Hamiltonian dynamics coupled to stochastic
thermostats. Journal of Statistical Physics, 133 :1–27, 2008.
[BLP07] N. Bruti-Liberati and E. Platen. Strong approximations of stochastic differential equations with
jumps. J. Comput. Appl. Math., 205–2 :982–1001, 2007.
[BLU07] A. Bonfiglioli, E. Lanconelli, and F. Uguzzoni. Stratified Lie groups and potential theory for their
sub-Laplacians. Springer, 2007.
[BM01] V. Bally and A. Matoussi. Weak solutions for SPDEs and backward doubly stochastic differential
equations. J. Theoret. Probab., 14–1 :125–164, 2001.
[BM03] Y. Belopolskaya and G. N. Milstein. An approximation method for Navier–Stokes equations based
on probabilistic approach. Stat. and Prob. Letters, 64 :201–211, 2003.

59
 [BP03] V. Bally and G. Pagès. A quantization algorithm for solving discrete time multi-dimensional optimal
stopping problems. Bernoulli, 9–6 :1003–1049, 2003.
[BP09] R.F. Bass and E.A. Perkins. A new technique for proving uniqueness for martingale problems. From
Probability to Geometry (I) : Volume in Honor of the 60th Birthday of Jean-Michel Bismut, pages
47–53, 2009.
[BPR06] A. Beskos, O. Papaspiliopoulos, and G. O. Roberts. Retrospective exact simulation of diffusion
sample paths with applications. Bernoulli, 12–6 :1077–1098, 2006.
[BPV01] E. Barucci, S. Polidoro, and V. Vespri. Some results on partial differential equations and asian
options. Math. Models Methods Appl. Sci, 3 :475–497, 2001.
[BR76] R. Bhattacharya and R. Rao. Normal approximations and asymptotic expansions. Wiley and sons,
1976.
[BS97] R. Bañuelos and R.G. Smits. Brownian Motion in cones. Prob. Th. Rel. Fields, 108 :299–319, 1997.
[BS03] A. N. Borodin and P. Salminen. Handbook of Brownian motion : facts and formulae. Birkhauser,
2003.
[BT96a] V. Bally and D. Talay. The law of the Euler scheme for stochastic differential equations : I. Conver-
gence rate of the distribution function. Prob. Th. Rel. Fields, 104-1 :43–60, 1996.
[BT96b] V. Bally and D. Talay. The law of the Euler scheme for stochastic differential equations, II. Conver-
gence rate of the density. Monte-Carlo methods and Appl., 2 :93–128, 1996.
[BT04] B. Bouchard and N. Touzi. Discrete time approximation and Monte-Carlo simulation of Backward
Stochastic Differential Equations. Stoch. Proc. Appl., 111 :175–206, 2004.
[BZ08] C. Bender and J. Zhang. Time discretization and Markovian Iteration for Coupled FBSDEs. Annals
of Applied Probability, 18 :143–177, 2008.
[Cat90] P. Cattiaux. Calcul stochastique et opérateurs dégénérés du second ordre - I. Résolvantes, théorème
de Hörmander et applications. Bull. Sc. Math., 2ème série, 114 :421–462, 1990.
[Cat91] P. Cattiaux. Calcul stochastique et opérateurs dégénérés du second ordre - II. Problème de Dirichlet.
Bull. Sc. Math., 2ème série, 115 :81–122, 1991.
[Che97] D. Chevance. Numerical methods for Backward Stochastic Differential Equations. Publ. Newton
Inst., Cambridge University Press, pages 232–234, 1997.
[CIL92] M. G. Crandall, I. Ishii, and Lions P. L. User’s guide to viscosity solutions of second order partial
differential equations. Bull. Amer. Math. Soc., 27 :1–67, 1992.
[CJ59] H. Carslaw and J. Jaeger. Conduction of Heat in Solids. Oxford Univ. Press, 1959.
[CKG06] C. Costantini, N. El Karoui, and E. Gobet. Boundary sensitivities for diffusion processes in time
dependent domains. Appl. Math. Optim., 54–2 :159–187, 2006.
[CM02] P. Cattiaux and L. Mesnager. Hypoelliptic non-homogeneous diffusions. Prob. Th. Rel. Fields,
123(4) :453–483, 2002.
[Cor07] J.-M. Coron. Control and nonlinearity. Mathematical Surveys and Monographs, 136, AMS, 2007.
[CP08] C. Cinti and S. Polidoro. Harnack inequalities and lifting procedure for evolution hypoelliptic equa-
tions. Lect. Notes Semin. Interdisc. Mat., 7 :93–105, 2008.
[Del02] F. Delarue. On the existence and uniqueness of solutions to FBSDEs in a non-degenerate case.
Stoch. Proc. and App., 99 :209–286, 2002.
[DG06] F. Delarue and G. Guatteri. Weak Existence and Uniqueness for FBSDEs. Stoc. Proc. and Appl.,
116 :1712–1742, 2006.
[DK06] R. A. Doney and A. E. Kyprianou. Overshoots and undershoots of Lévy processes. Ann. Appl.
Probab., 16–1 :91–106, 2006.

60
 [DL06] M. Deaconu and A. Lejay. A random walk on rectangles algorithm. Methodology And Computing
In Applied Probability, 8(1) :135–151, 2006.
[DMP96] J. Douglas, J. Ma, and P. Protter. Numerical methods for Forward-Backward Stochastic Differential
Equations. Ann. Appl. Prob., 6 :940–968, 1996.
[DP97] R. Darling and E. Pardoux. Backward SDE with random terminal time and applications to semilinear
elliptic PDE. Ann. Probab, 25 :1135–1159, 1997.
[Dyn63] E. B Dynkin. Markov Processes. Springer Verlag, 1963.
[EPRB99] J.-P. Eckmann, C.-A. Pillet, and L. Rey-Bellet. Non-equilibrium statistical mechanics of anharmonic
chains coupled to two heat baths at different temperatures. Comm. Math. Phys., 201–3 :657–697,
1999.
[FL01] C.D. Fuh and T.L. Lai. Asymptotic expansions in multidimensional Markov renewal theory and first
passage times for Markov random walks. Adv. in Appl. Probab., 33(3) :652–673, 2001.
[Fre85] M. Freidlin. Functional integration and Partial differential equations. Annals of Mathematics studies,
Princeton University Press, 1985.

[Fri64] A. Friedman. Partial differential equations of parabolic type. Prentice-Hall, 1964.

[Fri75] A. Friedman. Stochastic differential equations and applications. Vol. 1. Academic Press, 1975.
[FW98] M. Freidlin and A. Wentzell. Random Perturbations of Dynamical Systems. Springer, 1998.
[GL00] S. Graf and H. Lushgy. Foundations of quantization for random vectors. LNM-1730, Springer-Verlag,
2000.
[GL07] E. Gobet and C. Labart. Error expansion for the discretization of Backward Stochastic Differential
Equations. Stochastic Processes and Applications, 117–7 :803–829, 2007.
[GL08] E. Gobet and C. Labart. Sharp estimates for the convergence of the density of the euler scheme in
small time. Electron. Commun. Probab., 13 :352–363, 2008.
[GLW05] E. Gobet, J.P. Lemor, and X. Warin. A regression-based Monte-Carlo method to solve backward
stochastic differential equations. Ann. Appl. Proba., 15 :2172–2202, 2005.
[GLW06] E. Gobet, J.P. Lemor, and X. Warin. Rate of convergence of an empirical regression method for
solving generalized backward stochastic differential equations. Bernoulli, 12 :889–916, 2006.
[GM02] M.G. Garroni and J. L. Menaldi. Elliptic problems for integro-differential operators. Chapman and
Hall, 02.
[GM05] E. Gobet and S. Maire. Sequential Control Variates for Functionals of Markov Processes. Siam.
Journal of Num. Analysis, 43-3 :1256–1275, 2005.
[Gob00] E. Gobet. Euler schemes for the weak approximation of killed diffusion. Stoch. Proc. Appl., 87 :167–
197, 2000.
[Gob01] E. Gobet. Euler schemes and half-space approximation for the simulation of diffusions in a domain.
ESAIM : Prob. Stat., 5 :261–297, 2001.
[Gri92] P. Grisvard. Variational methods in the study of singularities. LMN, Springer, 1992.
[GT77] D. Gilbarg and N.S. Trudinger. Elliptic partial differential equations of second order. Springer Verlag,
1977.
[Guy06] J. Guyon. Euler scheme and tempered distributions. Stochastic Process. Appl., 116–6 :877–904,
2006.
[Hau02] E. Hausenblas. Error analysis for approximation of stochastic differential equations driven by Poisson
random measures. SIAM J. Num. Anal., 40–1 :87–113, 2002.

61
 [HN04] F. Hérau and F. Nier. Isotropic hypoellipticity and trend to equilibrium for the fokker-planck
equation with a high-degree potential. Arch. Ration. Mech. Anal., 171–2 :151–218, 2004.
[Hör67] L. Hörmander. Hypoelliptic second order differential operators. Acta. Math., 119 :147–171, 1967.
[HP57] E. Hille and R. S. Philips. Functional Analysis and semigroups. AMS, 57.
[IKO62] A. M. Il’in, A. S. Kalashnikov, and O. A. Oleinik. Second-order linear equations of parabolic type.
Uspehi Mat. Nauk, 17–3(105) :3–146, 1962.
[IL90] H. Ishii and P.-L. Lions. Viscosity solutions of fully nonlinear second-order elliptic partial differential
equations. J. Diff. Equations, 83 :26–78, 1990.
[IP06] P. Imkeller and I. Pavlyukevich. First exit times of sdes driven by stable lÈvy processes. Stoc. Proc.
Appl, 116–4 :611–642, 2006.
[IW89] N. Ikeda and . Watanabe. Stochastic Differential Equations and Diffusion processes. North Holland,
1989.
[Jac04] J. Jacod. The Euler scheme for lévy driven stochastic differential equations : limit theorems. Ann.
Probab., 5(32) :1830–1872, 2004.
[JKMP05] J. Jacod, T. G. Kurtz, S. Méléard, and P. Protter. The approximate euler method for Lévy driven
stochastic differential equations. Ann. Inst. H. PoincarÈ Probab. Statist, 41–3 :523–558, 2005.
[JMW96] A. Janicki, Z. Michna, and A. Weron. Approximation of stochastic differential equations driven by
α-stable Lévy motion. Appl. Math. (Warsaw), 24–2 :149–168, 1996.
[JO09] A. Joulin and Y. Ollivier. Curvature, concentration, and error estimates for Markov chain Monte
Carlo. http ://[Link]/abs/0904.1312, 2009.
[JP98] J. Jacod and P. Protter. Asymptotic error distributions for the euler method for stochastic differential
equations. Ann. Probab, 26–1 :267–307, 1998.
[KHP94] A. Kohatsu-Higa and P. Protter. The Euler scheme for SDEs driven by semimartingales. in Stochastic
analysis on infinite dimensional spaces, 1994. H. Kunita and [Link] (Eds.), 310 :141–151, 1994.
[KHT10] A. Kohatsu Higa and P. Tankov. Jump-adapted discretization schemes for Lévy-driven sdes. Preprint,
2010.
[KM00] V. Konakov and E. Mammen. Local limit theorems for transition densities of Markov chains conver-
ging to diffusions. Prob. Th. Rel. Fields, 117 :551–587, 2000.
[KM02] V. Konakov and E. Mammen. Edgeworth type expansions for euler schemes for stochastic differential
equations. Monte Carlo Methods Appl., 8–3 :271–285, 2002.
[KM09] V. Konakov and E. Mammen. Small time edgeworth-type expansions for weakly convergent nonho-
mogeneous markov chains. Probab. Theory Related Fields, 143–1-2 :137–176, 2009.
[KMR97] V.A. Kozlov, V.G. Maz’ya, and J. Rossmann. Elliptic boundary values problem in domains with
singularities. Math. Surveys and Monographs, 52. AMS, 1997.
[Kob00] M. Kobylanski. Backward stochastic differential equations and partial differential equations with
quadratic growth. Ann. Prob., 28 :558–602, 2000.
[Kol34] A. N. Kolmogorov. Zufällige Bewegungen (zur Theorie der Brownschen Bewegung). Ann. of Math.,
2-35 :116–117, 1934.
[Kol00] V. Kolokoltsov. Symmetric stable laws and stable-like jump diffusions. Proc. London Math. Soc.,
80 :725–768, 2000.
[KP91] T.G. Kurtz and P. Protter. Wong-Zakai corrections, random evolutions and numerical schemes for
SDEs. Academic Press, 1991.
[KP92] P. Kloeden and E. Platen. Numerical solution of stochastic differential equations. Springer-Verlag
Berlin Heidelberg New-York, 1992.

62
[KPQ97] N El Karoui, S. Peng, and M.C. Quenez. Backward stochastic differential equations in finance. Math.
Finance, 7 :1–71, 1997.
[KPZ86] M. Kardar, G. Parisi, and Y.-C. Zhang. Dynamic scaling of growing interfaces. Phys. Rev. Lett.,
56 :889–892, 1986.
[Kry96] N. V. Krylov. Controlled diffusion processes. Springer, 1996.
[KS84] S. Kusuoka and D. Stroock. Applications of the Malliavin calculus. I. Stochastic analysis (Ka-
tata/Kyoto, 1982), North-Holland Math. Library, 32 :271–306, 1984.
[KS85] S. Kusuoka and D. Stroock. Applications of the Malliavin calculus. II. J. Fac. Sci. Univ. Tokyo
Sect. IA Math, 32 :1–76, 1985.
[KS87] S. Kusuoka and D. Stroock. Applications of the Malliavin calculus. III. J. Fac. Sci. Univ. Tokyo
Sect. IA Math, 34 :391–442, 1987.
[Kun90] H. Kunita. Stochastic flows and stochastic differential equations. Cambridge University Press, 1990.
[Kus01] S. Kusuoka. Approximation of expectation of diffusion process and mathematical finance. Taniguchi
Conference on Mathematics Nara ’98. Adv. Stud. Pure Math., Math. Soc. Japan, Tokyo, 31 :147–165,
2001.
[KV06] S. A. Klokov and A. Yu Veretennikov. On mixing and convergence rates for a family of markov
processes approximating SDEs. Random Oper. Stochastic Equations, 14–2 :103–126, 2006.
[Led99] M. Ledoux. Concentration of measure and logarithmic sobolev inequalities. Séminaire de Probabilités
XXXIII. LNM, 1709, page 120–216, 1999.
[Lem07] V. Lemaire. An adaptive scheme for the approximation of dissipative systems. Stochastic Process.
Appl., 117–10 :1491–1518, 2007.
[Lie96] G.M. Lieberman. Second Order parabolic differential equations, 1st edn. World Scientfic, River Edge,
NJ, 1996.
[LP02] D. Lamberton and G. Pagès. Recursive computation of the invariant distribution of a diffusion.
Bernoulli, 8–3 :367–405, 2002.
[LP03] D. Lamberton and G. Pagès. Recursive computation of the invariant distribution of a diffusion : the
case of a weakly mean reverting drift. Stoch. Dyn., 3–4 :435–451, 2003.
[LS01] F. A. Longstaff and R. S. Schwartz. Valuing American options by simulation : a simple least-square
approach. Rev. Financ. Studies, 14 :113–147, 2001.
[LSU68] O.A. Ladyzenskaja, V.A. Solonnikov, and N.N. Ural’ceva. Linear and quasi-linear equations of
parabolic type. Vol.23 Trans. Math. Monog., AMS, Providence, 1968.
[LU68] O.A. Ladyzhenskaya and N.N. Ural’ceva. Linear and quasilinear equations of elliptic type. Academic
Press, 1968.
[LV04] T. Lyons and N. Victoir. Cubature on wiener space, Stochastic analysis with applications to mathe-
matical finance. Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci, 460 :169–198, 2004.
[Mal97] P. Malliavin. Stochastic analysis. Springer-Verlag, Berlin, 1997.
[Mél01] S. Méléard. Monte-Carlo approximations for 2d Navier-Stokes equations with measure initial data.
Prob. Th. Rel. Fields, 121–3 :367–388, 2001.
[Mil74] G.N. Milshtein. Approximate integration of stochastic differential equations. Theory Probab. Appl.,
19 :557–562, 1974.
[Mil76] G.N. Milshtein. A method of second order accuracy integration of sde. Theory Probab. Appl.,
23 :396–401, 1976.
[Mil85] G.N. Milshtein. Weak approximation of solutions of systems of stochastic differential equations.
Theory Probab. Appl., 30 :750–766, 1985.

63
 [Mil98] G. N. Milstein. Weak approximation of a diffusion process in a bounded domain. Stoch. Stoch.
Reports, 64 :211–233, 1998.
[MPY94] J. Ma, P. Protter, and J. Yong. Solving Forward-Backward Stochastic Differential Equations expli-
citly - a four step scheme. Prob. Th. Rel. Fields, 98 :339–359, 1994.
[MS67] H. P. McKean and I. M. Singer. Curvature and the eigenvalues of the Laplacian. J. Differential
Geometry, 1 :43–69, 1967.
[MSH02] J. C. Mattingly, A. M. Stuart, and D. J. Higham. Ergodicity for sdes and approximations : locally
Lipschitz vector fields and degenerate noise. Stochastic Process. Appl., 101–2 :185–232, 2002.
[MT99] G. N. Milstein and M. V. Tretyakov. Simulation of a space-time bounded diffusion. Ann. Appl.
Prob., 9–3 :732–779, 1999.
[MT04] G.N. Milstein and M.V. Tretyakov. Stochastic numerics for mathematical physics. Springer-Verlag,
Berlin, 2004.
[MT06a] F. Malrieu and D. Talay. Concentration inequalities for euler schemes. Monte Carlo and Quasi-Monte
Carlo Methods 2004, H. Niederreiter and D. Talay (Eds.), pages 355–37, 2006.

[MT06b] G. N. Milstein and M. V. Tretyakov. Numerical algorithms for forward-backward stochastic diffe-
rential equations. SIAM J. of Sci. Comput., 28 :561–582, 2006.
[MT07] G. N. Milstein and M. V. Tretyakov. Discretization of forward-backward stochastic differential
equations and related quasi-linear parabolic equations. SIAM J. of Sci. Comput., 1 :24–44, 2007.
[MY09] R. Mansuy and M. Yor. Aspects of Brownian motion. Springer, 2009.
[NN09] M. Ninomiya and S. Ninomiya. A new higher-order weak approximation scheme for stochastic
differential equations and the runge-kutta method. Finance and Stochastics, 13–3 :415–443, 2009.
[Nor86] J. R. Norris. Simplified Malliavin Calculus. Séminaire de Probabilités, XX :101–130, 1986.
[Nua95] D. Nualart. The Malliavin calculus and related topics. Probability and its Applications. Springer-
Verlag, New York, 1995.
[Nua98] D. Nualart. Analysis on Wiener space and anticipating stochastic calculus. In Lectures on probability
theory and statistics (Saint-Flour, 1995), pages 123–227. LNM 1690. Springer, Berlin, 1998.
[NV08] S. Ninomiya and N. Victoir. Weak approximation of stochastic differential equations and application
to derivative pricing. Appl. Math. Finance, 15–1-2 :107–121, 2008.
[Pan08] F. Panloup. Computation of the invariant measure for a Lévy driven SDE : rate of convergence.
Stochastic Process. Appl., 118–8 :1351–1384, 2008.
[PP90] E. Pardoux and S.G. Peng. Adapted solution of a Backward Stochastic Differential Equation. Systems
Control Lett, 14-1 :55–61, 1990.
[PP92] E. Pardoux and S.G. Peng. Backward Stochastic Differential Equations and quasilinear parabolic
partial differential equations. Stochastic partial differential equations and their applications (Char-
lotte, NC, 1991), 176 :200–217, 1992.
[PT97] P. Protter and D. Talay. The Euler scheme for lévy driven stochastic differential equations. Ann.
Probab., 1(25) :393–423, 1997.
[RBT00] L. Rey-Bellet and L. Thomas. Asymptotic behavior of thermal nonequilibrium steady states for a
driven chain of anharmonic oscillators. Comm. Math. Phys., 215–1 :1–24, 2000.
[Riv05] O. Rivière. Equations Différentielles Stochastiques Rétrogrades couplées : EDPs et discrétisation.
Thèse de doctorat. Université Paris Descartes, 2005.
[Rub03] S. Rubenthaler. Numerical simulation of the solution of a stochastic differential equation driven by
a Lévy process. Stochastic Process. Appl., 103–2 :311–349, 2003.

64
 [RY99] D. Revuz and M. Yor. Continuous martingales and Brownian motion. 3rd ed. Grundlehren der
Mathematischen Wissenschaften. 293. Berlin : Springer, 1999.
[She91] S. J. Sheu. Some estimates of the transition density of a nondegenerate diffusion Markov process.
Ann. Probab., 19–2 :538–561, 1991.
[Sie79] D. Siegmund. Corrected diffusion approximations in certain random walk problems. Adv. in Appl.
Probab., 11(4) :701–719, 1979.
[SK74] B.W. Stuck and B. Z. Kleiner. A statistical analysis of telephone noise. The Bell System Technical
Journal, 53 :1263–1320, 1974.
[Soi94] C. Soize. The Fokker-Planck equation for stochastic dynamical systems and its explicit steady state
solutions. Series on Advances in Mathematics for Applied Sciences, 17. World Scientific Publishing
Co., Inc., River Edge, NJ, 1994.
[ST94] G. Samorodnittsky and M. Taqqu. Stable non Gaussian random processes, stochastic models with
infinite variance. Chapman and Hall, New York, 1994.
[Str83] D.W. Stroock. Some applications of stochastic calculus to partial differential equations. In : Eleventh
Saint Flour probability summer school—1981 (Saint Flour, 1981), Lecture Notes in Math., 976,
Springer, Berlin, pages 267–382, 1983.
[Str88] D. W. Stroock. Diffusion semigroups corresponding to uniformly elliptic divergence form operators.
Séminaire de Probabilités, XXII :316–347, 1988.
[SV71] D.W. Stroock and S.R.S. Varadhan. Diffusion processes with boundary conditions. Comm. Pure
Appl. Math., 24 :147–225, 1971.
[SV79] D.W. Stroock and S.R.S. Varadhan. Multidimensional diffusion processes. Springer-Verlag Berlin
Heidelberg New-York, 1979.
[SY82] D. Siegmund and Y.S. Yuh. Brownian approximations for first passage probabilities. Z. Wahrsch.
verw. Gebiete, 59 :239–248, 1982.
[SZ92] J. Sokolowski and J. P. Zolésio. Introduction to shape optimization. Shape sensitivity analysis,
volume 16. Springer Series in Computational Mathematics, 1992.
[Tal86] D. Talay. Discrétisation d’une eds et calcul approché d’espérances de fonctionnelles de la solution.
Math. Model ling Numer. Anal., 20 :141–179, 1986.
[Tal02] D. Talay. Stochastic Hamiltonian dissipative systems : exponential convergence to the invariant
measure, and discretization by the implicit Euler scheme. Markov Processes and Related Fields,
8–2 :163–198, 2002.
[Tem01] E. Temam. Couverture Approchée d’Options Exotiques, Pricing des Options Asiatiques. Ph. Dis-
sertation. University of Paris VI, 2001.
[TT90] D. Talay and L. Tubaro. Expansion of the global error for numerical schemes solving stochastic
differential equations. Stoch. Anal. and App., 8-4 :94–120, 1990.
[Ver03] A. Yu. Veretennikov. On large deviations for approximations of SDEs. Probab. Theory Related
Fields, 125–1 :135–152, 2003.
[Web51] M. Weber. The fundamental solution of a degenerate partial differential equation of parabolic type.
Trans. Amer. Math. Soc., 71 :24–37, 1951.
[Whi73] N. Whitham. Linear and Non linear waves. Wiley, 1973.
[Woo82] M. Woodroofe. Nonlinear renewal theory in sequential analysis. SIAM, 1982.
[Zha88] C.H. Zhang. A nonlinear renewal theory. Ann. Prob., 16-2 :793–824, 1988.
[Zha04] J. Zhang. A numerical scheme for BSDEs. Annals of Applied Probability, 14–1 :459–488, 2004.

65