Conception d'un pro logiciel interactif sous r pour la simulation de processus de diffusion

N° d'ordre : 26/2012-M/MT

RÉPUBLIQUE ALGÉRIENNE DÉMOCRATIQUE ET POPULAIRE
MINISTÈRE DE L'ENSEIGNEMENT SUPÉRIEUR ET DE LA RECHERCHE SCIENTIFIQUE

UNIVERSITÉ DES SCIENCES ET TECHNOLOGIE DE HOUARI BOUMEDIENNE
FACULTÉ DES MATHÉMATIQUES



MÉMOIRE
PRÉSENTÉ POUR L'OBTENTION DU DIPLÔME DE MAGISTER.

EN MATHÉMATIQUES
SPÉCIALITÉ : PROBABILITÉ & STATISTIQUE

Par : Arsalane Chouaib GUIDOUM
Sujet

Conception d'un Pro Logiciel Interactif Sous

R Pour La Simulation de Processus de

Diffusion

Soutenu publiquement le 25/02/2012, devant le jury composé de :

Mr. Mustapha MOULAI Professeur à l'USTHB Président.

Mr. Kamal BOUKHETALA Professeur à l'USTHB Directeur de Mémoire.

Mr. Rachid OUAFI Maître de Conférences/A à l'USTHB Examinateur.
Mr. Abdelkader TATACHAK Maître de Conférences/A à l'USTHB Examinateur.

Table des matières

Table des matières Liste des figures Introduction Générale 1 Généralité sur Les Processus Stochastiques							iii v 3 6
	1.1		Introduction				7
	1.2		Définitions des processus				7
	1.3		Processus stochastiques particuliers				8
	1.4		Processus stochastiques établi à partir de la distribution gamma				10
	1.5		Processus stochastiques établi à partir de la distribution de Student				11
	1.6		Processus de Markov				12
	1.7		Processus du second ordre				13
	1.8		Processus ergodiques				15
	1.9		Martingales et temps d'arrêt				16
			1.9.1 Temps d'arrêt				16
			1.9.2 Théorème d'arrêt				17
2	Mouvement Brownien						18
	2.1		Introduction				19
	2.2		Construction du mouvement brownien				20
			2.2.1 Construction par un processus gaussien				20
			2.2.2 Construction par une limite d'une marche aléatoire				21
			2.2.3 Construction par le développement de Karhunen-Loève (D.K.L)				22
	2.3		Semi-groupe du mouvement brownien				24
			2.3.1 Propriété de Markov				24
			2.3.2 Mouvement brownien multidimensionnel				25
	2.4		Approximation la dérive d'un mouvement brownien standard par un bruit blanc
			gaussien				27
	2.5		Continuité des trajectoires				27
	2.6		Régularité des trajectoires				28
	2.7		Mouvement brownien arithmétique				32
	2.8		Mouvement brownien géométrique				34
	2.9		Pont brownien				36
		2.9.1 Construction par processus contraint 2.9.2 Construction par le développement de Karhunen-Loève (D.K.L) 2.10 Martingales exponentielles 2.10.1 Caractérisation de Paul Lévy du mouvement brownien 2.10.2 La loi du temps d'atteinte du mouvement brownien 2.10.3 Les temps de passage du mouvement brownien				37 39 39 40 42 43
		2.11 Conclusion				45
3		Processus de Diffusion				46
		3.1 Introduction				47
		3.2 Intégrale stochastique				48
		3.2.1 L'intégrale d'Itô				48
		3.2.2 L'intégrale de Stratonovitch				50
		3.2.3 Processus d'Itô				51
		3.2.4 Formule d'Itô				52
		3.2.5 La règle de multiplication				55
		3.3 Èquations différentielles stochastiques				56
		3.3.1 Introduction et définitions				56
		3.3.2 Existence et unicité des solutions de l'ÉDS				59
		3.3.3 Équation de Langevin				61
		3.3.4 Bruit blanc, bruit coloré				62
		3.3.5 Transformée de Lamperti				64
		3.4 Schémas numériques				65
		3.4.1 Simulation numérique				68
		3.4.2 Relation entre le schéma d'Euler et Milstein				72
		3.5 Les modèles attractives				75
		3.5.1 Modèle d'une diffusion en attraction M ó s=1(Vt)				76
		3.5.2 Modèle de deux diffusion en attraction M ó m>0(V(1) t		) +-' M 0 ó (V(2) t	) . . . .	80
		3.6 Conclusion				83
4		Comportement Asymptotique Des Processus de Diffusion				84
		4.1 Introduction				85
		4.2 Équation de Fokker-Planck				85
		4.2.1 L'origine de l'équation de Fokker-Planck				87
		4.2.2 Modélisation d'une équation physique				89
		4.2.3 Existence d'une solution				94
		4.3 Processus de diffusion stationnaires				95
		4.4 Classification des processus de diffusion linéaire				96
		4.4.1 Processus de diffusion de type N				97
		4.4.2 Processus de diffusion de type G				101
		4.4.3 Processus de diffusion de type B				105
		4.5 Calculs de l'instant de premier passage				108
		4.5.1 IPP d'une diffusion en attraction M ó s=1(Vt)				109

4.5.2 IPP de deux diffusion en attractionMó m>0(V(1)

t)?-M 0 ó(V(2)

t ) 114

4.6 Conclusion 115

Conclusion Générale 116

Annexe A : R Code 117

Annexe B : Packages Sim.DiffProc & Sim.DiffProcGUI 122

Bibliographie 131

Table des figures

1.1 Bruit blanc gaussien avec m = 0 et a² = 1. 9

1.2 Densité spectrale d'un bruit blanc gaussien. 9

1.3 Processus aléatoire établi à partir de la distribution ['(0.5,2). 11

1.4 Processus aléatoire établi à partir de la distribution St(2). 12

2.1 Trajectoire brownienne simulée a partir d'une distribution gaussienne. 21

2.2 Flux de trajectoires brownienne simulées a partir d'une distribution gaussienne. 21

2.3 Trajectoire brownienne comme limite d'une marche aléatoire. 22

2.4 Approximation d'un mouvement brownien par le D.K.L. 24

2.5 Mouvement brownien 2--D simulée a partir d'une distribution gaussienne 26

2.6 Approximation d'un mouvement brownien 2--D par une marche aléatoire 26

2.7 Mouvement brownien 3--D simulée a partir d'une distribution gaussienne 26

2.8 Fonction de covariance empirique d'un mouvement brownien standard. 29

2.9 Le mouvement brownien standard est non différentiables 31

2.10 La limite de mouvement brownien standard par rapport au temps. 31

2.11 Trajectoire d'un mouvement brownien arithmétique avec 9 = 2 et a = 1 34

2.12 Flux d'un mouvement brownien arithmétique avec 9 = 2 et a = 1. 34

2.13 Trajectoire d'un mouvement brownien géométrique avec 9 = 2 et a = 1 35

2.14 Flux d'un mouvement brownien géométrique avec 9 = 2 et a = 1. 35

2.15 Trajectoire d'un pont brownien à partir de Xt0 = --2 et XT = 1. 38

2.16 Flux de 100 trajectoires d'un pont brownien standard Xt0 = XT = 0. 38

2.17 Approximation d'un pont brownien standard par le D.K.L. 39

3.1 Simulation l'intégrale stochastique f 0 t WsdWs vs J 0 t W_s odW_s. 51

3.2 Trajectoire d'un processus d'Ornstein-Uhlenbeck avec r = 2 et a = 1. 59

3.3 Flux d'un processus d'Ornstein-Uhlenbeck avec r = 2 et a = 1. 59

3.4 Trajectoire simulée de l'équation de Langevin avec a = 2 et D = 1. 62

3.5 L'équation de Langevin en deux dimensions avec a = 2 et D = 1. 62

3.6 Simulation une seule trajectoire du modèle Radial Ornstein-Uhlenbeck par le

schéma d'Euler 70
3.7 Simulation un flux de 100 trajectoires du modèle Radial Ornstein-Uhlenbeck par

le schéma d'Euler. 70

3.8 Simulation une seule trajectoire du modèle dXt = (0.03tXt -X3 t )dt +0.2dWt par

le schéma de Milstein. 71
3.9 Simulation un flux de 100 trajectoires du modèle dXt = (0.03tXt - X3 t )dt +

0.2dWt par le schéma de Milstein. 71
3.10 Simulation un flux de 100 trajectoires du modèle dXt = cos(t)dt + sin(t)dWt par

le schéma de Itô-Taylor. 72
3.11 Transformation de modèle Cox-Ingersoll-Ross (CIR) dXt = (0.1 - 0.2Xt)dt +

0.05/XtdWt . 75

3.12 Trajectoire du polluant dans une surface d'eau turbulente en 2-D avec s = 1. . 79

3.13 Trajectoire du polluant dans une surface d'eau turbulente en 3-D avec s = 1. . 79

3.14 Ttrajectoire du polluant dans une surface d'eau turbulente en 2-D avec s > 1. . 80

3.15 Ttrajectoire du polluant dans une surface d'eau turbulente en 3-D avec s > 1. . 80

3.16 L'interaction entre deux insectes en 2-D. 83

3.17 L'interaction entre deux insectes en 3-D. 83

4.1 L'oscillateur de Van Der Pol, régime permanent sinusoïdal a = 0. 94

4.2 L'oscillateur de Van Der Pol, régime permanent non sinusoïdal a > 0. 94

4.3 Simulation un échantillon de taille 100 à partir du modèle VAG. 99

4.4 Ajustement de la distribution stationnaire du modèle VAG par la méthode d'his-

togramme 101 4.5 Ajustement de la distribution stationnaire du modèle VAG par la méthode du noyau.101 4.6 Ajustement de la distribution stationnaire du modèle CIR par la méthode d'his-

togramme 104 4.7 Ajustement de la distribution stationnaire du modèle CIR par la méthode du noyau.104 4.8 Ajustement de la distribution stationnaire de processus de diffusion de Jacobi par

la méthode d'histogramme 107
4.9 Ajustement de la distribution stationnaire de processus de diffusion de Jacobi par

la méthode du noyau 107

4.10 L'instant de premier passage du modèle Mó s=1(Vt) en 2-D 111

4.11 L'instant de premier passage du modèle Mó s=1(Vt) en 3-D 111

4.12 Ajustement de la distribution de 1/'r^s=1

c par la méthode d'histogramme. 112

4.13 Ajustement de la distribution de 1/'rs=1

c par la méthode du noyau. 112

XIV Graphical User Interface for Sim.DiffProc package at start-up 126

Remerciements.

T

out d'abord je tiens à remercier Dieu de m'avoir donné le courage, la morale et la santé pour mener à bien ce travail.

J

e tiens à remercier avec tous mes sentiments de respectueuse gratitude mon promoteur Mr. Kamal BOUKHETALA Professeur à l'USTHB pour sa proposition de sujet ainsi pour son soutien, ses orientations et ses précieux conseils.

J

'exprime aussi ma profonde gratitude à Mr. Mustapha MOULAI Professeur à l'USTHB, pour avoir accepté de présider le jury de soutenance.

J

e remercie également : Mr. Rachid OUAFI Maître de conférences à l'USTHB, et Mr. Abdelkader TATACHAK Maître de conférences à l'USTHB, pour avoir accepter d'examiner cette thèse.

E

nfin, je remercie chaleureusement toutes les personnes qui m'ont aidé, et qui ont contribué de proche ou de loin à la réalisation de ce travail.

Résumé

D

ans ce travail, on propose un nouveau package Sim.DiffProc [9] pour la simulation des processus de diffusion, muni d'une interface graphique (GUT), sous langage R. Le déve-

loppement de l'outil informatique (logiciels et matériels) ces dernières années, nous a motivé de réaliser ce travail. A l'aide de ce package, nous pouvons traiter beaucoup de problèmes théoriques difficiles liée à l'utilisation des processus de diffusion, pour des recherches pratiques, tels que la simulation numérique trajectoires de la solution d'une EDS. Ce qui permet à beaucoup d'utilisateurs dans différents domaines à l'employer comme outil sophistiqué à la modélisation de leurs problèmes pratiques. Le problème de dispersion d'un polluant [2, 4], en présence d'un domaine attractif que nous avons traité dans ce travail en est un bon exemple. Cet exemple montre l'utilité et l'importance pratique des processus de diffusion dans la modélisation simulation de situations réelles complexes. La fonction de densité de la variable aléatoire 'r_c "instant de premier passage" de la frontière de domaine d'attraction peut être utilisée pour déterminer le taux de concentration des particules polluantes à l'intérieur du domaine. Les études de simulation et les analyses statistiques mises en application à l'aide du package Sim.DiffProc, se présentent efficaces et performantes, comparativement aux résultats théoriques explicitement ou approximativement déterminés par les modèles de processus de diffusion considérés.

Introduction générale

L

a nature aléatoire de nombreux phénomènes évolutifs, dans des domaines très divers, tels ceux de la physique, astronomie, biologie, les mathématiques financières, géologie, analyse

génétique, épidémiologie et beaucoup d'autres champs de la science et de l'ingénierie, nécessite fréquemment la description de tels phénomènes par des équations différentielles stochastiques, qui peuvent être un outil puissant pour la modélisation. L'étude des équations stochastiques s'est beaucoup développée ces dernières années, c'est un domaine plein de perspectives et de recherche. Dans beaucoup de literatures, on rencontre les équations stochastiques avec légère variations d'écriture, ainsi que leurs divers applications, voir par exemple [1, 18, 21, 36, 40].

Soit x0 E Rⁿ, Considérons, pour u : Rⁿ --+ Rⁿ, champ vectoriel régulier, l'équation différentielle ordinaire :

Jdx dt (t) = u(t,x(t)), Vt ~ 0 (1)

x(0) = x0

la solution de l'équation (1), si elle est unique, est représentée par une trajectoire. Dans la plupart des applications où une telle équation différentielle ordinaire intervient, les trajectoires mesurées expérimentalement ne sont que rarement conformes aux solutions analytiques de l'équation. Des effets aléatoires viennent se superposer à la trajectoire idéale, et il semble donc raisonnable de modifier l'équation (1) en y introduisant un processus aléatoire perturbant le système. Formellement, la modification s'écrit :

(

Vt ~ 0 (2)

dX dt (t) = u(t,Xt) +a(t,Xt)ît, X0 = x0 où : a : Rⁿ --+ R^nxm, et ît est un bruit blanc gaussien m-dimensionnel. Cette approche soulève

les problèmes suivants :

1. Définir ît de manière rigoureuse.

2. En déduire l'influence de ît dans la résolution de l'équation (1).

3. Montrer que l'équation (1) a une solution, discuter de l'unicité, du comportement asymptotique, du rôle de a, de x0,...

Considérons tout d'abord que le bruit blanc gaussien ît est la dérivée formelle du processus de Wiener. Dans le cas général, l'équation (2) peut alors se réécrire :

(

dXt = u(t,Xt)dt + a(t,Xt)dWt, Vt = 0 (3)

X0 = x0

l'équation obtenue alors est une équation différentielle stochastique.

Dans notre travail on s'intéresse a les solutions des équations différentielles stochastiques qu'on va étudier, sont des processus Markoviens à valeurs dans Rⁿ, qu'on appelle les processus de diffusion qui constituant la plus importante classe. Cependant, le mouvement Brownien Wt (où processus de Wiener) est utilisé comme un modèle de diffusion homogène. Les modèles de diffusions sont représentés par des équations différentielles stochastique de la forme (3), dite équation de type Itô.

Le mathématicien Kiyoshi Itô [27], a donné un sens à ces équations, et a montrer l'existence et l'unicité de la solution de l'équation (3) sous certaines conditions de régularité des fonctions u(t,x) et ó(t,x), appelées respectivement coefficient de dérive et coefficient de diffusion. Xt est solution de l'équation (3) si et seulement si :

Z t Z t

Xt = x0 + ₀ u(s,X_s)ds + ₀ ó(s,X_s)dW_s, t = 0 (4)

Xt est donc défini à l'aide d'intégrales dites intégrales stochastiques d'Itô.

Notre travail est structuré de la façon suivante :

Dans le premier chapitre nous rappelons quelques notions fondamentales et les principaux aspects théorique des processus stochastiques. Les martingales sont un outil essentiel et important, qui nous a permis d'établir de nombreux résultats ainsi le théorème d'arrêt.

Le second chapitre est consacré à l'étude détaillée du mouvement Brownien, sa construction, ces propriétés, ainsi sa simulation unidimensionnel et multidimensionnel. La difficulté de modélisation du mouvement brownien réside dans le fait que ce mouvement est aléatoire et que statistiquement, le déplacement est nul, c'est-à-dire que le coefficient de dérive est nul, il n'y a pas de mouvement d'ensemble. Pour cela il est aussi possible de définir la notion de mouvement brownien avec dérive. Il s'agit d'un mouvement brownien arithmétique (dans la finance modèle de Merton). Le mouvement brownien géométrique (où le modèle de marché de Black et Scholes) et le pont brownien, sont des objets mathématiques de la théorie des probabilités, représentent l'un des plus importants processus stochastique et ont de nombreuses applications dans la finance, et ce en les illustrant par des exemples de simulation, afin de rendre plus facile la compréhension des propriétés qui le caractérise. On termine ce chapitre par la détermination les lois des temps de passage, et la caractérisation de Paul Lévy du mouvement brownien.

Le troisième chapitre est composé de quatre parties. La premier partie sera consacrée à donner un sens à l'intégrale stochastique lorsque celle-ci est prise par rapport à un mouvement brownien, notons juste que la construction de cet intégrale au sens d'Itô [27] fait appel à un changement de la mesure d'intégration habituelle par une nouvelle mesure, ainsi nous parlons de la formule d'Itô dans le cas unidimensionnel et multidimensionnel, qui sera illustrée par quelques exemples

d'applications. Dans la deuxième partie on s'intéresse à définir les équations différentielles stochastiques d'une manière générale, à l'exception les processus de diffusion, naturellement nous annonçons le théorème de l'existence et unicité des solutions des ÉDS, nous donnons la définition et la différence entre un bruit blanc et un bruit coloré, ainsi la transformée de Lamperti [30] qui sera très utile sur le plan numérique pour améliorer la précision. La détermination de la solution exacte de Xt d'une ÉDS est très difficile à déterminer surtout lorsque la partie aléatoire est exprimée en fonction du processus Xt, c'est ainsi qu'on se tourne vers les méthodes numériques de l'ÉDS dont la base est fondée sur la discrétisation du temps. Nous présenterons alors dans la troisième partie, l'approche numérique, en citant les principaux schémas qui la concernent et qui permettent d'approcher la solution exactes de ces équations. Nous illustrerons par des exemples de simulation des différents méthodes. Dans la dernier partie de ce chapitre on s'intéresse à l'utilisation des processus de diffusion dans la modélisation de deux phénomènes, le premier est de modéliser une trajectoire d'un polluant, qui se déplace sur une surface d'eau turbulente en présence d'un mécanisme d'attraction [2, 4], et le deuxième porte sur une modélisation d'un phénomène d'attraction entre deux insectes mâle et femelle [5].

Le dernier chapitre, nous commençons par la présentation de deux équations fondamentales permettant de d'écrire l'évolution des lois de probabilités relatives à un processus de diffusion. Les équations de Fokker-Planck [35] ou les équations de Kolmogorov, qui décrit l'évolution dynamique de la densité de probabilité d'un système hors d'équilibre, ainsi sa distribution stationnaire si elle existe. Nous annonçons deux théorèmes donnant l'existe d'une relation simple entre les deux différentielles Itô et Stratonovitch, ces derniers sont très utiles pour la modélisation d'une équation physique, nous illustrerons cette modélisation par l'exemple de l'oscillateur de Van Der Pol. Par la suite nous parlerons de classification des processus de diffusion linéaire, en distinguant trois types. Ce chapitre se terminer par l'étude et l'analyse statistique de la variable aléatoire l'instant de premier passage "IPP" dans le cas du modèle d'une diffusion en attraction [2, 4, 6, 7], et l'instant de la première rencontre entre deux insectes, c'est-à-dire dans le cas du modèle de deux diffusion en attraction.

On terminera ce travail avec une conclusion générale, et quelques perspectives. Il est à noté que dans tous nos programmes, nous utiliserons le langage R [32], qui sont présentée dans l'annexe A, ainsi pour toutes les exemples. Dans l'annexe B nous donnons quelques règles pour la création des packages sous R, et nous donnons aussi une présentation de deux packages:

(1) Sim.DiffProc : Simulation of Diffusion Processes [8, 9].

(2) Sim.DiffProcGUI : Graphical User Interface for Simulation of Diffusion Processes [10].

1

Généralité sur Les Processus Stochastiques

Sommaire
1.1	Introduction	7
1.2	Définitions des processus	7
1.3	Processus stochastiques particuliers	8
1.4	Processus stochastiques établi à partir de la distribution gamma	10
1.5	Processus stochastiques établi à partir de la distribution de Student . . . .	11
1.6	Processus de Markov	12
1.7	Processus du second ordre	13
1.8	Processus ergodiques	15
1.9	Martingales et temps d'arrêt	16

1.1 Introduction

L

'origine des processus stochastiques remonte aux progrès faits au début du XX^e siècle dans certaines branches appliquées, telles que la mécanique statistique (par Gibbs, Boltzmann,

Poincaré, Smoluchowski et Langevin). Les bases théoriques ont été formulées plus tard par [17, 28] et d'autres (1930-1940). C'est durant cette période que le mot "stochastique", qui provient du grec stokhastikos "conjectural", a commencé à être employé. D'autres progrès ont été faits aussi à partir du mouvement brownien en physique (par Einstein, Lévy et Wiener).

Nous introduisons dans ce chapitre les principaux processus aléatoires à l'exception du mouvement brownien qui fait l'objet d'un chapitre séparé. Nous retrouverons la plupart de ces processus dans les problèmes de calcul stochastique. Les processus du second ordre ont de nombreuses applications en théorie du signal.

1.2 Définitions des processus

Soit (~,.i,P) un espace de probabilité, et T un ensemble d'indices (T = [a,b],T = [0,oo[,...) un processus stochastique X(t,(0) à valeur dans un espace mesurable (E, ) est une application de T ×~ dans E qui est mesurable par rapport à la mesure du produit A·P où A est la mesure de Lebesgue sur T. Il est noté indifféremment Xt((0) ou X(t,(0). La fonction t '-? X(t,(0) est appelée trajectoire ou réalisation de Xt. A t fixé, la fonction (0 i-? X(t,(0) est une variable aléatoire. Xt est adapté à la filtration _t si Xt est t-mesurable. Le théorème de Kolmogorov assure l'existence des processus stochastiques. Xt est un processus centré si son espérance est nulle E(Xt) = 0, et si Xt est dans L² (E|Xt|² < oo), on définit :

La moyenne du processus

Zm_x(t) = E(Xt((0)) = Ù Xt((0)dP((0)

La variance

a2 _x(t) = E[|Xt -E(Xt)|²]

La fonction de covariance

['(s,t) = E(X_s -E(X_s))(Xt -E(Xt))

= E(X_sXt) - E(X_s)E(Xt)

La régularité des trajectoires est déterminer par le théorème de Kolmogorov.

Théorème 1.1 (Kolmogorov) Soit (Xt)t=0 un processus stochastique tel que pour tout t, t + h dans [a,b], il existe des constantes p > 0, c > 0 et r > 0 vérifiant

E[|Xt+h -Xt|^p] = c|h|^1+r

alors presque toutes les trajectoires sont continues.

Les théorèmes suivants fondent la représentation spectral des processus stationnaires.

Théorème 1.2 (Herglotz) Soit c une fonction semi-définie positive de Z dans C. Il existe une unique mesure positive u sur ] - ð,+ð] telle que pour tout entier n E Z,

Z +ð

c(n) = -ð einëdu(ë)

Théorème 1.3 (Bochner) Soit c une fonction continue et semi-définie positive de R dans C. Il existe une unique mesure bornée u telle que pour tout t E R,

Z +8

c(t) = -8 eitëdu(ë)

Définition 1.1 (Filtration) Soit (Ù,.i,P) un espace de probabilité, la filtration est une famille croissante de sous tribus de .i, noté par ( t,t = 0). La tribu _t est une description mathématique de toute l'information dont on dispose à l'instant t. Cette information nous permet d'attribuer des probabilités cohérentes aux événements pouvant intervenir.

Définition 1.2 (Processus adapté) Un processus {Xt,t = 0} est dit adapté à la filtration ( t,t = 0) si pour chaque t, Xt est t-mesurable. Un processus adapté est celui pour lequel une description probabiliste est réalisable.

1.3 Processus stochastiques particuliers

Définition 1.3 (Processus strictement stationnaire) Un processus Xt est strictement stationnaire si pour tout entier n et pour tous réels t1,t2,...,t_n et pour tout h, les variables aléatoires (Xt1,Xt2,...,Xt_n) ^et (Xt1+h,Xt2+h,...,Xt_n+h) ont même loi.

Définition 1.4 (Processus stationnaire) Un processus Xt est stationnaire (ou faiblement stationnaire) si :

- son espérance E(Xt) est une constante indépendante du temps t.

- sa fonction de corrélation R(s,t) ne dépend que de la différence ô = t - s.

- R(ô) est continue (à l'origine).

Si Xt est un processus stationnaire, sa fonction de corrélation est continue en 0 et définie positive. La fonction

R(ô)

Ö(ô) = R(0)

est une fonction caractéristique d'après le théorème 1.3 de Bochner. Par conséquent, il existe une fonction Z(ù) appelée fonction de répartition spectrale telle que Z(ù) s'annule au voisinage de -8, et reste fini lorsque ù tend vers +8 telle que

Z +8

R(ô) = -8 eiùôdZ(ù)

De plus, si Z est absolument continue, il existe une fonction S(ù) appelée densité spectrale du processus Xt, telle que S(ù) = dZ(ù)/dù et

+8 .

R(ô) = f e^iùôS(ù)dù

Par application de la transformée de Fourier inverse, si f |R(ô)|dô < 8 alors la densité spectrale est donnée par

S(ù) = 1 f +8 _e- i

^ùôR(ù)dù

2ð _8

Exemple 1.1 Un bruit blanc est un processus stationnaire dont la densité spectrale est constante SX(ù) = a. Sa fonction de corrélation est donc de la forme RX(ô) = aä(ô). Un tel processus contient toutes les fréquences (d'oil son nom de bruit blanc). Sa variance n'est pas bornée, il n'est donc pas réalisable en théorie.



Les figures 1.1 et 1.2 montre une seule trajectoire simulée d'un bruit blanc gaussien, avec sa densité spectrale estimée.



FIGURE 1.1 - Bruit blanc gaussien avec m = 0 et ó² = 1.

FIGURE 1.2 - Densité spectrale d'un bruit blanc gaussien.

Définition 1.5 (Processus à accroissements indépendants) Un processus Xt est un processus à accroissements indépendants si pour tous réels t1 < t2 <
·
·
· < t_n, les variables aléatoires Xt0,Xt1 - Xt0,...,Xt_n - _Xtn-1 sont indépendantes. Si Xt est un processus à accroissements indépendants, alors pour tout s et t tels que 0 < s < t, Xt -X_s est indépendant de as = ó(X_u,u < s). La loi de Xt est entièrement déterminée par la loi de Xt - X_s.

Définition 1.6 (Processus à accroissements indépendants stationnaires) Un processus Xt est un processus à accroissements indépendants stationnaires si Xt est un processus à accroissements indépendants et si Xt -X_s a même loi que _Xt-s. Si Xt est un processus à accroissements indépendants, alors sa fonction caractéristique s'écrit

Öt1,...,t_n(x1,...,xn) = Öt1(x1 + ··· + xn)Öt1,t2(x2 + ··· + xn)...Öt_n-1,t_n(xn)

avec

Öti,ti+1(x) = E[expix(Xti+1 - Xt_i)]

si Xt est un processus à accroissement indépendants stationnaires alors sa loi est indéfiniment divisible. La loi du processus Xt est entièrement déterminée par la loi de X1.

Théorème 1.4 Soit Xt un processus à accroissement indépendants stationnaires, (at) la filtra-
tion engendrée par le processus Xt et T un temps d'arrêt. Alors le processus Yt = XT+t - Xt est
un processus à accroissements indépendants stationnaires de même loi que Xt et indépendant de

T.

Définition 1.7 (Processus gaussien) Un processus Xt est gaussien si pour tout entier n et pour tous réels t1,t2,...,t_n les variables aléatoires (Xt1,Xt2,...,Xt_n) ont une distribution gaussienne. Si on note m(t) l'espérance de Xt et R(s,t) la fonction de corrélation du processus, sa fonction caractéristique s'écrit

(Öt1 ,...,tn (x1 , . . . , xn) = E expi		n ? j=1		xjXtj (ù))
(= exp i	n ? j=1		^n 1 xjm(tj) - ? 2 j=1		n ? k=1	)xjxkR(tj,tk)

un processus gaussien est entièrement déterminé par la donnée de sa valeur moyenne m(t) et de sa fonction de corrélation R(s,t).

1.4 Processus stochastiques établi à partir de la distribution gamma

La distribution gamma, caractérisée par les deux paramètres á et â, se définit comme suit :

f(x|á,â) = â á á-1 p ( -;)

(á)^{x ex} où (á) est la fonction gamma définie comme :

(á) =

_yá-1

0 exp(-y)dy

A noter que si á = n et que n est un entier, on a alors le résultat suivant :

(n) = (n--1)!

La moyenne d'une variable aléatoire X qui obéit à une distribution gamma est de :

E(X) = áâ

Et sa variance est de :

var(X) = áâ²

La figure 1.3 montre l'évolution d'un tel processus. Comme on peut le constater à la lecture de cette figure, de nombreux sauts se manifestent dans un tel processus stochastique, en ce sens que ces sauts s'éloignent de beaucoup plus d'écart-types de la moyenne que ne l'autorise la distribution normale.



FIGURE 1.3 - Processus aléatoire établi à partir de la distribution (0.5,2).

1.5 Processus stochastiques établi à partir de la distribution de Student

La distribution de Student est une autre distribution qui est de plus en plus utilisée dans la modélisation des prix des produits dérivés.

Soit Z une N(0,1) et ÷²_r) une variable aléatoire Chi-deux avec r degrés de liberté. Alors (

T = Z a une distribution Student qui est notée _St(r). Elle est symétrique autour de 0, mais

\14₎/r

tend vers la distribution normale lorsque le nombre de degrés de liberté se dirige vers l'infini.

La figure 1.4 prend acte de l'évolution d'un tel processus. Encore une fois, des sauts très apparents se font jour lors du déroulement du processus stochastique de la série, ils sont d'autant plus amples que l'on réduit le nombre de degrés de liberté de la distribution de Student.



FIGURE 1.4 - Processus aléatoire établi à partir de la distribution St(2).

1.6 Processus de Markov

Soit (S ,A,P) un espace de probabilité et (E,B) l'espace des états, B désigne la tribu des boréliens de E, Un processus Xt est un processus de Markov si pour tout u et t = 0 et pour tout A ? B, on a

P(Xt+u ? A|X_s,s = t) = P(Xt+u ? A|Xt)

Ce qui signifie que le processus ne dépend que du dernier instant et non de toute son histoire

P(Xt_n < xn|Xtn-1= xtn-1,...,Xt1 = xt1) = P(Xt_n < xtn|Xtn-1 = xtn-1)

Pour démontrer qu'un processus est un processus de Markov, il suffit de montrer par le théorème
de classe monotone que pour tout n, pour tous réels t1 = t2 = ··· = t_n et pour toute fonction f

borélienne bornée

E(f (Xt_n)) = E(f (Xtn)|Xtn-1)

La probabilité de transition pour passer de l'état x au temps s à un état appartenant à A à l'instant t est notée pour s < t

Ps,t(x,A) = p(s,x;t,A) = P(Xt ? A|X_s = x)

La fonction A 7? p(s,x;t,A) est une probabilité sur B. La probabilité de transition vérifie l'équation de Chapman-Kolmogorov qui s'écrit sous les formes suivantes. Soit s < u < t tel que X_u = y, on a

Ps,t(x,A) =LPs,u(x,dy)Pu,t(y, A)

et dans le cas d'un espace E dénombrable

Ps,t(x,z) = L Ps,u(x,y)Pu,t(y,z)

y?E

Si on note pt₀(A) la distribution de Xt0

pt₀(A) = P(Xt0 ? A) alors pour tous réels t0 < t1 < ··· < tn,

n-1

(Xt1 ? B1,...,Xt_n ? B_n) = f ;tk

Pt_o(dX0)

P

nf p(tk_i,xk_,,dxk)p(tn-1,xn-1;tn,Bn)

k=1 Bk

En particulier,

P(Xt ? B) = f p(t0,x;t,B)pt₀(dx)

Le processus Xt est un processus de Markov homogène si pour tout s et t de T, la transition

p(s,x;t,B) = p(x,T,B)

ne dépend que de
·r = t -s. Les opérateurs Ps,t = Pt-s sont notés simplement Pt. Ils forment dans le cas homogène un semi-groupe. On a PtP_s = Pt+_s, pour tout s,t > 0.

1.7 Processus du second ordre

Un processus Xt est un processus du second ordre si E|Xt|² < 00. Une série chronologique est un processus du second ordre à temps discret. On dit que le processus du second ordre Xt converge en moyenne quadratique (i.e. dans L²) vers une variable aléatoire Y quand t tend vers t0 si et seulement si sa fonction de corrélation converge quand s et t tendent vers t0 et dans ce cas

R(s,t) ---? s,t?t0

E|Y|^2

On dit que Xt un processus du second ordre est continu en moyenne quadratique si

h?0

lim E|Xt+h - Xt|² = 0

Proposition 1.1 Un processus du second ordre Xt est continu en moyenne quadratique si et seulement si sa fonction de corrélation R(t,t) est continue en t.

La continuité en moyenne quadratique n'entraîne pas la continuité des réalisations du processus. Si Xt est un processus de Poisson, sa fonction de corrélation R(s,t) = ëmin(s,t) est continue, mais presque toutes les réalisations ont des discontinuités sur un intervalle de temps fini.

2

On dit que Xt un processus du second ordre est dérivable en moyenne quadratique si la limite du taux d'accroissement (Xt+h - Xt)/h converge dans L² vers une variable notée X⁰_t .

Jiⁱ?0 fh X: = 0

Xt+h - Xt

Proposition 1.2 Un processus du second ordre Xt est dérivable en moyenne quadratique si et seulement si sa fonction de corrélation R(s,t) est dérivable en (s,t)

La dérivation en moyenne quadratique est linéaire. Si Xt est dérivable en moyenne quadratique, alors Xt est continue en moyenne quadratique.

Proposition 1.3 Si Xt est un processus du second ordre centré et dérivable en moyenne quadratique et si les dérivées partielles existent, alors

?k+1

RX(k)X(l) (s, t) = E[X^(k) (s)X^(l) (t)] = ?ks?ltRX(s,t)

En particulier, un processus stationnaire est différentiable en t si et seulement si sa fonction de corrélation R(u) admet une dérivée du second ordre en u et dans ce cas

ce⁺¹

RX(k)X(l)(u) = (-1 )k duk+1RX(u)

Si Xt est dérivable en moyenne quadratique et si Xt admet la densité spectrale SX(ù), alors le processus dérivé admet une fonction spectrale SX/(ù) donnée par

SX^,(ù) = ù²SX(ù)

Plus généralement,

= (i ù)^k (i ù)^l S _X (ù)

SX(k) X(l) (ù)

En particulier, la fonction de corrélation de la dérivée d'ordre k d'un processus Xt est

d2k

RX(k)(u) = (-1)k _du2kRX (u)

Sa fonction spectrale est donnée par la formule

SX(k)(ù) = (-1)^k(ù)^2kSX(ù)

Théorème 1.5 (Karhunen-Loève) Soit Xt un processus du second ordre pour t ? [a,b], continu en moyenne quadratique, centré. Il existe un et un seul développement, appelé développement de Karhunen-Loève de la forme

Xt(ù) =

8 ? n=1

Yn(ù)Ön(t)

ot les Y_n sont des variables aléatoires du second ordre telles que E(YiYj) = 0 si i =6 j et E(|Yi|²) = ëi ot les ëi sont les valeurs propres et les Öi sont les vecteurs propres de l'opérateur R : L² ? L² symétrique compact

b

Rf(s) = Z R(s,t)f(t)dt

vérifiant

Za b R(s,t)Ö_n(t)dt = ë_nÖ_n(t)

Exemple 1.2 Pour un mouvement brownien {Wt,0 = t = 1} (Chapitre 2 Section 2.2.3), le développement s'écrit pour une suite de variables aléatoires Z_n de loi normale centrée réduite N(0,1) telles que E|Z_n|² = 1,

Wt = v2

8 ? n=1

sin(n + 1/2)ðt Zn (n+ 1/2)ð

Et pour un pont brownien standard {Xt,0 = t = 1} (Chapitre 2 Section 2.9.2), nous avons le D.K.L

8

?

n=1

Xt = v2

sin(ðnt)

ðn

Zn

1.8 Processus ergodiques

Un processus du second ordre Xt de moyenne m(t) et de fonction de corrélation R(s,t) est ergodique si

R(s,t)dsdt-^mq? 0

1

T

_I

0

T²

T

_I

0

1

T

T

_I

0

ce qui équivaut à

(Xt - m(t))dt-^mq? 0

si Xt est stationnaire, Xt est ergodique si et seulement si

0

R(ô) ---?

|ô|?8

ou :

lim

T?8

T 1 - T) R(ô)ch = 0

1 fT ( ô

Ainsi la connaissance de la fonction de corrélation et de la moyenne permet de déterminer si le processus est ergodique.

1.9 Martingales et temps d'arrêt

Les martingales à temps discret ou à temps continu sont l'outil essentiel du probabiliste. Elles permettent d'établir de nombreux résultats.

Définition 1.8 Soit (Ù,A,P) un espace de probabilité, muni d'une filtration (Ft)t=0. Un processus à valeurs réelles (Mt)t=0 est une Ft-martingale si :

- il est adapté à la filtration (Ft)t=0, ce qui veut dire que pour tout t, Mt est Ft-mesurable. - chaque variable Mt est intégrable, et :

s = t E(Mt|F_s) = M_s

on dit que Mt est une Ft-sur-martingale (resp. Ft-sous-martingale) si l'égalité ci-dessus est remplacée par:

E(Mt|F_s) = M_s (resp.E(Mt|F_s) = M_s)

En particulier, l'espérance E(Mt) d'une martingale, (resp. d'une sur-martingale, sous-martingale), est une fonction constante du temps (resp. décroissante, croissante). De manière évidente, une martingale est un processus qui est à la fois une sur-martingale et une sous-martingale et si Mt est une martingale, alors -Mt est une sous-martingale.

Propriété 1.1 (1). (Mt)t=0 est une martingale si et seulement si (Mt)t=0 est à la fois une surmartingale et une sous-martingale.

(2). (Mt)t=0 est une sous-martingale si et seulement si (-Mt)t=0 est une sur-martingale.

(3). La somme de deux martingales (resp. sous-martingale, sur-martingale) est une martingale (resp. sous-martingale, sur-martingale).

1.9.1 Temps d'arrêt

Cette notion joue un rôle très important en théorie des probabilités.

Définition 1.9 Soit (Ft)t=0 une filtration. Une application T : Ù 7? [0,8[ est un temps d'arrêt si {T = t} ? Ft pour tout t = 0.

Un temps d'arrêt est donc un temps aléatoire, tel que sur chaque ensemble {ù : T(ù) = t}, l'application ù 7? T(ù) dépend seulement de ce qui s'est passé avant le temps t.

Un joueur honnête, qui ne peut pas anticiper sur les événements futurs, peu décider d'arrêter le jeu au temps aléatoire T uniquement si T est un temps d'arrêt. Un exemple trivial de temps d'arrêt est donné par T(ù) = t pour tout ù.

En dehors des temps constants, l'exemple fondamental de temps d'arrêt est le temps d'atteinte d'un ensemble borélien A par un processus Xt à trajectoires continues à droite et adapté à la filtration Ft. On définit plus précisément :

T = inf(t = 0;Xt ? A)

1.9.2 Théorème d'arrêt

Soit Mt une martingale. La propriété des martingale peut facilement être étendue aux temps d'arrêt bornés.

Théorème 1.6 Si S et T sont deux temps d'arrêt et si a E R, alors :

E(MT |FS) = MS sur l'ensemble{S < T < a} (1.1)

En particulier, si T est un temps d'arrêt qui est borné, on a :

E(MT) = E(M0) (1.2)

Quand Mt désigne de nouveau le gain d'un joueur au temps t, la propriété (1.2) peut être interprétée comme suit. Quelle que soit la stratégie non anticipante que le joueur choisit pour arrêter le jeu, et s'il doit finir de jouer avant un temps déterministe donné (aussi grand que soit ce temps), alors la valeur espérée de son gain est constante et égale à son capital initial.

Observons que la relation (1.1) est, en général, fausse sur l'ensemble {S < T}, et de même (1.2) est fausse si T n'est pas borné. Par exemple si Mt = Bt ou Bt est un mouvement brownien et si T = inf(t : Mt = 1), alors E(M0) = 0 < E(MT) = 1. Dans ce cas, le temps aléatoire T est presque sûrement fini, mais n'est pas borné et a même une espérance infinie.

2

Mouvement Brownien

Sommaire

2.1 Introduction 19

2.2 Construction du mouvement brownien 20

2.3 Semi-groupe du mouvement brownien 24

2.4 Approximation la dérive d'un mouvement brownien standard par un bruit blanc gaussien 27

2.5 Continuité des trajectoires 27

2.6 Régularité des trajectoires 28

2.7 Mouvement brownien arithmétique 32

2.8 Mouvement brownien géométrique 34

2.9 Pont brownien 36

2.10 Martingales exponentielles 39

2.11 Conclusion 45

2.1 Introduction

L

e mouvement brownien est une description mathématique du mouvement aléatoire d'une grosse particule immergée dans un fluide et qui n'est soumise à aucune autre interaction que des chocs avec les petites molécules du fluide environnant. Il en résulte un mouvement très irrégulier de la grosse particule, qui a été décrit pour la première fois en 1827 par le biologiste Robert Brown [11] alors qu'il observait du pollen de Clarkia pulchella (une espèce de fleur sauvage nord-américaine), puis de diverses autres plantes, en suspension dans l'eau.

La description physique la plus élémentaire du phénomène est la suivante

1. Entre deux chocs, la grosse particule se déplace en ligne droite avec une vitesse constante.

2. La grosse particule est accélérée lorsqu'elle rencontre une molécule de fluide ou une paroi.

Cette description permet de décrire avec succès le comportement thermodynamique des gaz (théorie cinétique des gaz), ainsi que le phénomène de diffusion. Brown n'est pas exactement le premier à avoir fait cette observation. Il signale lui-même que plusieurs auteurs avaient suggéré l'existence d'un tel mouvement (en lien avec les théories vitalistes de l'époque). Parmi ceux-ci, certains l'avaient effectivement décrit, on peut mentionner en particulier l'abbé John Turberville Needham (1713-1781), célèbre à son époque pour sa grande maîtrise du microscope.

La réalité des observations de Brown a été discutée tout au long du XXe siècle. Compte tenu de la médiocre qualité de l'optique dont il disposait, certains ont contesté qu'il ait pu voir véritablement le mouvement brownien, qui intéresse des particules de quelques micromètres au plus. Les expériences ont été refaites par l'Anglais Brian Ford [12] au début des années 1990, avec le matériel employé par Brown et dans les conditions les plus semblables possibles. Le mouvement a bien été observé dans ces conditions, ce qui valide les observations de Brown.

Quelques années plus trad en 1905, Albert Einstein mit en évidence les étranges relations que le processus entretenait avec l'équation de la chaleur. Vers 1909, Jean Perrin entreprit son étude expérimentale et Paul Langevin posa la première équation. Mais il faudra attendre 1925 et les travaux de Norbert Wiener pour que le mouvement brownien ait véritablement un sens mathématique comme modèle d'un bruit blanc. À partir des années 1950, Kiyoshi Itô [27] l'utilisa pour définir l'intégrale qui porte son nom et jeta les bases du calcul stochastique.

Nous nous intéressons dans ce chapitre principalement a les principaux processus aléatoires à l'exception du mouvement brownien, ces propriétés, ces caractéristiques, ainsi sa construction mathématique et sa simulation par trois approches, nous retrouverons la plupart de ces processus dans les problèmes de calcul stochastique. Nous introduisons la notion de martingale exponentielle qui joue un rôle très importante dans les calcules des lois des temps de passage du mouvement brownien, ainsi la théorème de caractérisation de Paul Lévy.

2.2 Construction du mouvement brownien

Ce processus de diffusion peut être construit par differentes approches. Les definitions les plus usuelles du mouvement brownien sont les suivantes.

2.2.1 Construction par un processus gaussien

Un processus stochastique Wt est un mouvement brownien ou un processus de Wiener si W0 = 0 (on dit que Wt est issu de 0) et si pour tous reels 0 < t1 < t2 < ··· < t_n, les variables aleatoire Wt1 -Wt0, . . . ,Wt_n -Wtn-1 sont independants et suivent une distribution gaussienne centree reduite (on dit que le mouvement brownien est standard si m = 0 et ó = 1) telle que :

{

E(Wt+h -Wt) = ⁰ E(Wt+h -Wt)² = h

Dans le cas general, lorsque le mouvement brownien n'est pas centre reduit, on a :

E(Wtk - Wtk-1) = m(tk - tk-1)

{

E((Wt+h -Wt - m(tk -tk-1))2 = ó²(tk -tk-1)

le vecteur (Wt0,Wt1,...,Wt_n) est un vecteur gaussien. Le processus Wt suit une loi gaussienne de moyenne mt et de variance ó²t. On peut facilement simuler une trajectoire de mouvement brownien dans un intervalle de temps [0,T], il suffit de fixee un pas de temps Ät > 0 et d'ecrire

v

W(Ät) = W(Ät) -W(0) ~ N(0,Ät) ~ ÄtN(0,1)

Les accroissements (WnÄt -W(n-1)Ät) etant independants et gaussiens, il suffit donc de simuler une loi gaussienne

Wt+Ät - Wt ~ N(0, Ät) ~ vÄtN(0, 1)

Ainsi, nous pouvons simuler facilement une seule trajectoire brownienne de la façon suivante. On considère la subdivision de l'intervalle de temps [0,T] suivante 0 = t1 < t2 < ··· < tN < tN+1 = T, avec _ti+1 -ti = Ät, pour i = 1 on a W(0) = W(t1) = 0. On donne l'algorithme suivant :

1. Generee un nouveau variable aleatoire Z de la distribution gaussienne N(0,1).

2. i = i+1.

3. W(ti) = W(ti-1)+ZvÄt.

4. Si i = N + 1, reiterez a l'etape 1.

La fonction BMN permet de simuler un mouvement brownien standard {Wt,t = 0} dans l'intervalle de temps [t0,T] avec un pas Ät = (T -t0)/N, et la fonction BMNF permet de simuler un flux brownienne standard (C = ó²)

R> BMN(N = 1000, t0 = 0, T = 1, C = 1)

R> BMNF(N = 1000, M = 100, t0 = 0, T = 1, C = 1)



FIGURE 2.1 - Trajectoire brownienne simulée a partir d'une distribution gaussienne.

FIGURE 2.2 - Flux de trajectoires brownienne simulées a partir d'une distribution gaussienne.

2.2.2 Construction par une limite d'une marche aléatoire

Une caractérisation du mouvement brownien indique qu'il peut voir en tant que limite d'une marche aléatoire dans le sens suivant. Considérons une suite de variables aléatoires indépendants Xi centrées de variance ó² et la marche aléatoire S_n = X1 +X2 + ··· +X_n, où

(

+1 si p = 1/2

Xi =

-1 si p = 1/2

On définit une suite de variables Y_n par la formule suivante :

	Yn(t) =	S[nt] + (nt - [nt])X[nt]+1
où [.] est la partie entière.		óvn

Ce résultat fondamental est donné par le théorème de Donsker (1951) et est, en fait, au niveau des processus, une version du théorème usuel de la limite centrale

Théorème 2.1 (Principe d'invariance de Donsker) Soit (Xn)n=1 une suite de variables aléatoires réelles indépendantes, identiquement distribuées, avec E(X_n) = 0 et E(X² _n) = 1. Soit S_n = Y.1=i=nXi avec S0 = 0.

Les processus des sommes normalisées Yⁿ_t = 1 vnS[nt] (oil [nt] désigne la partie entière de nt) convergent en loi, en tant que processus, vers le mouvement brownien.

Cette convergence donne une définition du mouvement brownien comme l'unique limite (en loi) de marches aléatoires.

Le code 1, permettre de simulée un mouvement brownien standard comme l'unique limite de marche aléatoire. La figure 2.3 donne une représentation de l'approximation d'un mouvement brownien par une marche aléatoire pour n = 10,n = 100 et n = 1000.



FIGURE 2.3 - Trajectoire brownienne comme limite d'une marche aléatoire.

2.2.3 Construction par le développement de Karhunen-Loève (D.K.L)

Pour un mouvement brownien {Wt,0 = t = 1}, le développement de Karhunen-Loève [15, 16] s'écrit pour une suite de variables aléatoires Z_n de loi normale centrée réduite telles que E|Z² _n| = 1,

v Wt = 2

8 ? n=1

sin(n + 1/2)ðt Zn ?t ? [0,1] (2.1) (n+1/2)ð

Preuve A partir du théorème de Karhunen-Loève 1.5, on a l'équation aux valeurs propres

Z 1

0 (s?t)ö_n(s)ds = ë_nö_n(t)

s'écrit

Z t Z 1

0 sö_n(s)ds + _t tö_n(s)ds = ë_nö_n(t)

soit en dérivant

Z 1

ë_nö⁰ _n(t) = t ön(s)ds

deux fois

ë_nö⁰⁰ _n(t) = -ö_n(t)

avec ö⁰_n(1) = 0. La résolution de cette équation conduit à

\/

ö_n(t) = csin(t/ ë_n)

La constante c est déterminée par le fait que les fonctions propos sont orthonormées

Z0 1 ö2 _n(s)ds = 1

v 2. La condition limite ö⁰_n(1) = 0 donne l'équation cos(1/vë_n) = 0 qui fournit les

d'où c =

valeurs propres

1

ë_n = (n+1/2)ð2

d'où le développement

v Wt = 2

8 ? n=1

Zn

sin(n+ 1/2)ðt (n+ 1/2)ð

Le code 2 (Annexe A), permettre de simulée un mouvement brownien standard a partir de D.K.L. La figure 2.4 donne une représentation graphique du une approximation d'un mouvement brownien par le D.K.L pour n = 10,n = 100 et n = 1000.



FIGURE 2.4 - Approximation d'un mouvement brownien par le D.K.L.

2.3 Semi-groupe du mouvement brownien

2.3.1 Propriété de Markov

Considérons un mouvement brownien Wt sur (1,A,P), et la filtration (Ft)t=0 qu'il engendre. Puisqu'il est à accroissements indépendants, la variable Y := Wt+s -Wt est indépendante de la tribu Ft. On a pour chaque fonction f borélienne bornée sur R :

Z 1

E( f (Wt+s)|Ft) = E( f (Wt +Y)|Ft) = _v2ðse^-x2/^2s f (Wt + x)dx (2.2)

R

Cette formule montre que conditionnellement à Ft, la loi de _Wt+s ne dépend pas de tout le passé (c'est-à-dire de toutes les variables W_r pour r = t), mais seulement de la valeur "présente" Wt du processus. On dira que le mouvement brownien est un processus de Markov.

Définition 2.1 Un processus Xt est un processus de Markov si, étant donné la filtration Ft engendrée par le processus, celui-ci vérifie la propriété de Markov, à savoir que pour tous s,t = 0 et pour toute fonction f borélienne bornée sur R :

E(f(Xt+s)|Ft) = E(f(Xt+s)|Xt) (2.3)

dy (2.5)

Dans le cas du mouvement brownien, les variables W_r pour r = t, sont independantes, conditionnellement à la valeur de Wt. De plus, la loi de _Wt+s sachantt depend bien sûr de s, mais pas de t. On dit que le mouvement brownien est un processus de Markov homogène.

Définition 2.2 Soit Xt un processus de Markov homogène. On appelle semi-groupe de transition de Xt la famille (Pt)t=0 d' operateurs positifs lineaires :

Pt : Ö ? L⁸ (R^d) PtÖ ? L⁸ (R^d)

PtÖ(x) = E(Ö(Xt)|X0 = x) = f_d Ö(y)Pt(x,dy)

R

qui satisfait Pt1 = 1 et la propriete de semi-groupe :

P0 = Id , Pt+s = P_s ? Pt ?s, t = 0 (2.4)

Dans le cas du mouvement brownien, qui est un processus de markov homogène, le semigroupe (Pt)t=0 est donne par :



~ ¹Pt (x, dy) = exp _-(y - x)²

2ðt 2t

comme le montre immediatement la formule (2.2).

En utilisant les relations (2.2) et (2.5), on obtient facilement les proprietes :

E(Wt|F_s) = W_s ; E(W2 _t |F_s) = W2 s +t - s , ?s < t (2.6)

2.3.2 Mouvement brownien multidimensionnel

Le mouvement brownien multidimensionnel est très utilise dans les modèles de marche en temps continu. Par exemple, lors de la modelisation simultanee des prix de plusieurs actifs risques.

Définition 2.3 Un mouvement brownien d-dimensionnel est une collection W = (Wⁱ)1=i=d de d mouvements browniens à valeurs reelles Wi = (W_t ⁱ)t=0, qui sont indépendants entre eux.

Ce processus est encore un processus de Markov homogène (et même un processus à accroissements independants). Son semi-groupe vaut alors :

~ 1 Pt(x,dy) = (2ðt)d/2 exp -Hy 2tx ||2)

dy (2.7)

où x et y appartiennent à R^d, ||.|| designe la norme euclidienne sur R^d, et dy la mesure de Lebesgue sur R^d.

Les deux fonctions BMN2D et BMRW2D permettre de simulee un mouvement brownien 2 dimensions, respectivement par la distribution gaussienne et par approximation une marche aleatoire.

R> BMN2D(N = 10000, t0 = 0, T = 1, x0 = 0, y0 = 0, Sigma = 1) R> BMRW2D(N = 10000, t0 = 0, T = 1, x0 = 0, y0 = 0, Sigma = 1)



FIGURE 2.5 - Mouvement brownien 2--D simulée a partir d'une distribution gaussienne.

FIGURE 2.6 - Approximation d'un mouvement brownien 2--D par une marche aléatoire.

La fonction BMN3D permettre de simulée un mouvement brownien 3-dimensions.

R> BMN3D(N = 10000, t0 = 0, T = 1, X0 = 0, Y0 = 0, Z0 = 0, Sigma = 1)



FIGURE 2.7 - Mouvement brownien 3--D simulée a partir d'une distribution gaussienne.

2.4 Approximation la dérive d'un mouvement brownien standard par un bruit blanc gaussien

Un bruit blanc est une réalisation d'un processus aléatoire dans lequel la densité spectrale est la même pour toutes les fréquences, on parle souvent de bruit blanc gaussien, il s'agit d'un bruit blanc qui suit une loi normale de moyenne et variance données.

Définition 2.4 Un processus åt est qualifié de bruit blanc gaussien si : - E(åt) = 0 et E(å2 _t ) = ó2

- åt et å_s sont indépendants ?t =6 s

- åt ~ N(0,ó²)

Par analogie avec le bruit blanc en temps discret, défini comme une suite de variables aléatoires, centrées, du second ordre et indépendants, on cherche à définir {åt}t=0 comme un processus stochastique vérifiant ?t > 0 et ?h > 0 :

- E(åt) = 0

- E(åtåt+h) = ä0(h)

où ä0 est la mesure de Dirac en 0.

Un tel processus n'existe pas. Son idéalisation est la dérivée d'un mouvement brownien standard. Si pour Ät > 0 fixé, on considère le processus

åt =

Wt+Ät -Wt

Ät

Il est facile de montrer grâce aux propriétés données par la définition de mouvement brownien

que :

( ) _(|h| )

1 1 - |h|

E(åt) = 0 et E(åtåt+h) = ¹[0,1]

Ät Ät Ät

Quand Ät -? 0,E(åtåt+h) converge vers ä0(h). Il est donc clair que la dérivée formelle dWt

dt a

les propriétés d'un bruit blanc gaussien. Ce qui justifie l'affirmation concernant l'idéalisation du

bruit blanc. Donc on peut écrire formellement :	åt =	dWt
		dt

2.5 Continuité des trajectoires

Dire qu'un processus aléatoire {Xt,t = 0} est continu c'est, par définition dire que

lim |Xt+h -Xt| = 0 h?0

Selon le type de convergence de cette variable aléatoire, on obtient une continuité plus ou moins forte. La plus faible des notions de continuité est liée à la convergence en loi. Elle est évidement vérifiée. Nous allons démontrer une continuité en probabilité pour le mouvement brownien standard.

Proposition 2.1 Soit å > 0 et {Wt,t = 0} un mouvement brownien standard. On a

1

lim

h?0 h

P(|Wt+h -Wt| > å) = 0

Preuve Soit h > 0, par définition, l'accroissement Wt+h -Wt admet pour loi N(0,h). Donc

1

2 8 1 x2

_hP(|Wt+h - Wr| > å) = h , 2ðhe- 2h dx

8 1 1

å

x

< 2

J e-

E .V22.c h3/2 2h

ie v

1

2h

-

_å2
2h

e

å

h³/²

v

2

=

_vð

Le dernier terme converge vers 0 lorsque h ? 0.

2.6 Régularité des trajectoires

Le mouvement brownien a de nombreuse propriétés dont certaines peuvent être prise comme définition.

Proposition 2.2 Le processus Wt est un processus à accroissements indépendants de fonction de covariance

(s,t) = E(W_sWt) = ó²min(s,t)

Preuve Le mouvement brownien est un processus centré, les accroissements étant indépendants, on a pour 0 = s < t,

E(W_sWt) = E(W_s(Wt -W_s)) +E(W²_s )

= E(W_s)E(Wt -W_s) + var(W_s) = 0 + ó²s = ó²s

et pour 0 = t < s on a,

E(WtW_s) = E(Wt(W_s -Wt)) + E(W²

t )

= E(Wt)E(W_s -Wt)+var(Wt) = 0 + ó²t = ó²t

d'ou : E(W_sWt) = ó²min(s,t).

dx

La fonction BMcov donne une représentation graphique (Figure 2.8) de la fonction de covariance empirique d'un mouvement brownien standard (pour M = 500 trajectoires simulée de taille N = 1000, C = ó2).

R> BMcov(N = 1000, M = 500, T = 1, C = 1)



FIGURE 2.8 - Fonction de covariance empirique d'un mouvement brownien standard.

Proposition 2.3 La densité de ÄW = (Wt1 -Wt0,...,Wt_n -Wtn-1) est donnée par

fÄW(x1,...,xn) =	n ? j=1	1 p2ð(tj -tj-1) exp	-x2 j
			2(tj -tj-1)

Proposition 2.4 La densité de W = (Wt1,...,Wt_n) est

f(x1,...,xn) =

n ? j=1

1 p2ð(tj -tj-1) exp

n ? j=1

-(xj -xj-1)^2 2(tj -tj-1)

Proposition 2.5 Pour un mouvement brownien standard (m = 0,ó = 1) et pour tout n, l'espérance

Z +8

E(Wt+h -Wt)²ⁿ = 1

v

-8 x²ⁿe^-x2/^2hdx = 1.3...(2n - 1)hⁿ 2ðh Presque toutes les réalisations du mouvement brownien sont continues (appliquer le théorème 1.1 de Kolmogorov avec c = 3,p = 4 et r = 1).

Proposition 2.6 Soit Wt un mouvement brownien standard. On a presque sûrement,

lim t?+8

sup

Wt vt

= +8 , lim t?0

sup

Wt vt

= +8,

inf ^Wt vt

, lim

t?0

inf ^Wt vt

= -8

= -8,

lim Wt = 0.

t?+8 t

lim

t?+8

et

Preuve Comme pour tout s > 0,

U = lim t?+8	sup	Wt	= lim t?+8	sup	Wt-s-Wt
		vt			vt

cette limite est indépendante de la tribu ó(W_u,u = s) et donc de ó(W_u,u = 0) et par conséquent
d'elle-même. On a soit P(U = +8) = 1 soit P(U = a) = 1. Supposons que la limite sup soit

)atteinte en a. Pour tout b > a, P (^Wt_vt > b ) tend vers 0, mais P (^Wt_vt > b= P(W1 > 0) = 1, Donc

a ne peut être qu'infini. Pour la limite inférieure, il suffit de considérer -Wt. Pour les limites au voisinage de zéro, on considérera le processus Xt = tW1/t, pour établir la dernière formule, on pose s = 1/t et on considère le mouvement brownien Xt = tW1/t, on a Wt/t = sW1/s = X_s ? 0 p.s. puisque Xt est un mouvement brownien standard.

Les deux codes 3 et 4, permettes de vérifier par simulation la proposition 2.6, le mouvement brownien standard est simulée a partir de développement de Karhunen-Loève. La figure 2.9 montre clairement que le mouvement brownien standard est non différentiables, et la figure 2.10 montre que la limite de mouvement brownien standard par rapport au temps tend vers 0 quand t?+8.



FIGURE 2.9 - Le mouvement brownien standard est non différentiables.

FIGURE 2.10 - La limite de mouvement brownien standard par rapport au temps.

Proposition 2.7 Soit 0 = t0 < t1 < ··· < t_n < tn+1 = t une subdivision de l'intervalle [0,t] dont le pas tend vers 0, la variation quadratique converge dans L² vers t

2 L2

V_n = Eⁿk=0(Wtk+1 -Wtk) ---?t Preuve Le calcul de la norme donne

||Vn -t||^2L2 = E		n k=0		(Wtk+1 -Wtk)^2 - (tk+1 -tk)	!2
= E	n k=0		((Wtk+1 -Wtk)^2 - (tk+1 -tk))^2 +E i6=j			E(XiXj)

avec

Xi = ((Wti+1 -Wt_i)² - (ti+1 -ti))

Le produit des termes croisés est nul, car on a pour i < j E(XiXj) = E(E(XiXj)|Ft_j) = E(E(Xi)E(Xj|Ft_j)) = 0 et puisque Wti+1 -Wti est Ft_j-mesurable. Il reste donc

n

||V_n -t||²_L2 = E((Wtk+1 -Wt_k)² _-(tk+1 -tk))²

k=0

Mais les accroissements (Wtk+1 ^--Wtk)²/(tk+1 -- tk) ont même loi que Z², où Z est une variable aléatoire centrée réduite. D'où

n

||V_n --t||²_L2 = E(Z² -- 1)² k=0 (tk+1 --tk)²

< E(Z² -- 1)²t _sup(tk+1 --tk) -+ 0

Proposition 2.8 Si Wt est un mouvement brownien, alors il en est de même pour les processus suivante

(1). Xt = 1aWa2t pour a constante non nulle (invariance par changement d'échelle).

(2). Xt = tW1/t pour t > 0 et X0 = 0 (invariance par inversion de temps).

(3). Xt = WT--t --WT avec T > 0 et t E [0,T] (invariance par retournement du temps). Preuve Il suffit de vérifier que le processus est gaussien et de même covariance que Wt(t A s). Vérifions (1), on a

1

E(X_sXt) = _a2 E(Wa2sWa2t) = a2 min(a²s,a²t) = min(s,t)

1

De mêmee pour (2), on a

E(X_sXt) = tsE(W1/sW1/t)) = tsmin(1/s,1/t)) = min(s,t)

et pour (3), on a

E(X_sXt) = E((WT--ss ^--WT)(WT--tt --WT)))

= E(WT--sWT--t)^--E(WT--sWT)) -- E(WTWT--t)+E(W²T ))

=min( T--s, T--t)-- T+s+t

Sit > s, E(X_sXt ) = T-- t-- T+s+ t =s Sit <s, E(X_sXt ) =T --s--T +s+t =t d'ouu : E(X_sXt) = min(s,t).

2.7 Mouvement brownien arithmétiquee

Le mouvement brownien standard ou processus de wiener que nous avons étudier comporte certaines lacunes. D'abordd sa dérivee est nul, or plusieurs processus stochastiques comportent une tendance. Par exemple les indices boursiers font montre d'unee tendance àa la hausse àa long terme. De plus, la variance d'unn processus stochastique est égale au pas At, comme cette variance ne peut accepter qu'unn nombre trèss limité de processus stochastiques, ilt y a lieu de choisir la partie aléatoire d'unn processus stochastique par la variance observée de la série. Le mouvement

brownien arithmétique où mouvement brownien avec dérive (Dans la finance modèle de Merton (1973)) corrige ces deux déficiences du processus de wiener. Il s'écrit comme suit

dXt = èdt + ódWt (2.8) où è est la dérive de processus et ó son écart-type, Wt est un mouvement brownien standard. A partir de l'équation différentielle stochastique (2.8) en remarque que c'est la partie stochastique du processus qui domine à court terme et sa tendance à long terme. 1

Pour simuler une seule trajectoire du mouvement brownien arithmétique sur un intervalle de temps [t0,T] avec un pas Ät = (T -t0)/N, nous avons utilisé la fonction ABM. Et pour un flux de trajectoires utilisant la fonction ABMF.

On considère la subdivision de l'intervalle de temps [t0,T] suivante t0 < ··· < tN < tN+1 = T, avec ti+1 - ti = Ät, pour i = 0 on a W(0) = W(t0) = 0 et X(0) = X(t0) = x0, on a l'algorithme suivant :

1. Générer un nouveau variable aléatoire Z de la distribution gaussienne N(0,1).

2. i = i+1. v

3. W(ti) = W(ti-1) + Z Ät.

4. X(ti) = Xti-1 + èÄt + ó(Wti -Wti-1)

5. Si i = N + 1, réitérez a l'étape 1.

Remarque 2.1 Si è = 0 on a un mouvement brownien.

Comme en prend acte la figure 2.11, le mouvement brownien arithmétique est caractérisé par un dérive à long terme ponctué de déviations qui dépendent de l'écart-type du processus stochastique.

R> ABM(N = 1000, t0 = 0, T = 1, x0 = 0, theta = 2, sigma = 1)

R> ABMF(N = 1000, M = 100, t0 = 0, T = 1, x0 = 0, theta = 2, sigma = 1)

v

1. En effet, è est multiplié par dt et ó par dt.



FIGURE 2.11 - Trajectoire d'un mouvement brownien arithmétique avec 9 = 2 et a = 1.

FIGURE 2.12 - Flux d'un mouvement brownien arithmétique avec 9 = 2 eta = 1.

2.8 Mouvement brownien géométrique

Un mouvement brownien arithmétique est inapproprié pour décrire l'évolution du prix d'une action, étant donné la croissance espérée du prix de cette action, désignée par 9, et l'écart-type du taux de rendement de l'action, représenté par a. En effet, cela supposerait que le rendement total de l'action soit dSt

St , aurait tendance à diminuer au cours du temps, ce qui est contraire aux données observées sur les rendements des actions. On fait donc l'hypothèse que le prix d'une action obéir à un mouvement brownien géométrique où le modèle de marché de Black et Scholes, c'est-à-dire

dSt = 9Stdt + aStdWt (2.9)

La dérive et l'écart-type sont donc multipliés par St, soit le niveau du prix de l'action. Il s'ensuit le taux de rendement de l'action suit un mouvement brownien arithmétique

dSt
St

= 9dt +adWt (2.10)

Utilisant le lemme d'Itô que nous examinerons plus en détail dans le chapitre suivant, l'équation différentielle stochastique (2.9) admet pour solution

~~ ~ ~

9 - a2

St = S0 exp t + aWt , S0 > 0 (2.11)

2

( ~

Comme le mouvement brownien standard Wt est de loi N(0,t), è - ^ó2 t + óWt est de loi

2

(( ) )

è - ó2

N t,ó²tet St est de loi lognormale.

2

Pour simuler une seule trajectoire du mouvement brownien géométrique (utilisant l'équation (2.11)) sur un intervalle de temps [t0,T] avec un pas Ät = (T - t0)/N, nous avons utilisé la fonction GBM. Et pour un flux de trajectoires utilisant la fonction GBMF.

R> GBM(N = 1000, t0 = 0, T = 1, x0 = 1, theta = 2, sigma = 1)

R> GBMF(N = 1000, M = 100, t0 = 0, T = 1, x0 = 1, theta = 2, sigma = 1)



FIGURE 2.13 - Trajectoire d'un mouvement brownien géométrique avec è = 2 et ó = 1.

FIGURE 2.14 - Flux d'un mouvement brownien géométrique avec è = 2 et ó = 1.

Si l'équation différentielle stochastique du prix de l'action se conforme à un mouvement brownien géométrique, alors le prix de l'action suit une loi lognormale. Et si tel le cas le logarithme de St suit une loi normale. Pour le montrer soit l'équation (2.9) du prix de l'action, et soit la fonction g qui est égale à ln(St). Cette fonction dépend donc de la variable aléatoire St. Selon le lemme d'Itô, l'équation différentielle de la fonction g s'écrit

?g 1 ?²g

dgt = ?SdSt + ?S2 dS² (2.12)

t

2

Le lemme d'Itô s'apparente donc à une expansion de Taylor du second degré. Or,

et

1

=

S

?g
?S

1
S2

?²g

?S2 =

En substituant ces dérivées dans l'équation (2.12) et en remplaçant dSt par l'équation (2.9), on obtient

1 ó2S2 t

dgt = (èStdt + óStdWt) - (2.13)

St 2S2 t

Avec dS2 t = ó2S2 _t dt, nous le justifierons dans le chapitre suivant. Après simplification de l'équation (2.13), on obtient

( )

è - 1

dgt = ₂ó²dt +ódWt (2.14)

Le logarithme de St, soit la fonction g suit bien une loi normale puisque dWt obéit à une loi normale. On peut alors écrire le prix de l'action comme suit

[( ) ]

è - 1

St = St-1 exp ₂ó² dt + ódWt (2.15)

2.9 Pont brownien

Un pont brownien est un objet mathématique de la théorie des probabilités, également appelé mouvement brownien attaché ("tied down brownian motion" en anglais), mouvement brownien attaché en a ("brownian motion tied down at a" en anglais) ou mouvement brownien épinglé ("pinned brownian motion" en anglais).

Considérons un mouvement brownien (Wt)0=t=1 sur l'intervalle de temps [0,1]. On veut calculer la loi conditionnelle de Wt par rapport à la tribu engendrée par la variable aléatoire W1. Notons £(a,.) la loi conditionnelle de Wt sachant W1 = a. Si l'on modifie la loi du processus Wt, celui-ci changer de nom et de propriétés.

Définition 2.5 Le processus Wt sous la loi £(a,.), est appelé le pont brownien, et pont brownien standard si a = 0.

Définition 2.6 Un pont brownien Xt est un processus stochastique à temps continu dont la loi celle d'un processus de Wiener sachant l'événement W0 = W1 = a. Il s'agit d'un processus aléatoire gaussien, c'est à dire que la loi de probabilité de tout vecteur (Wt1,...,Wt_n) conditionnellement à W1 = a est gaussienne. Il est alors caractérisé par sa moyenne et sa fonction de covariance, qui sont données par:

E(Wt|W1 = a) = at ?0 = t = 1

cov(W_s,Wt|W1 = a) = s(1-t) ?0 = s < t = 1

Observons que la moyenne dépend de a, mais pas la covariance. Et si a = 0 on a un pont brownien standard dont

E(Wt|W1 = 0) = 0 et cov(W_s,Wt|W1 = 0) = s(1-t) ?0 = s < t = 1

Son semi-groupe de transition est donné par la formule suivante, pour 0 = s < t < 1

( ~y - 1

1-s(x(1 -t) + a(t - s)))2 )

_P(a) 1

s,t (x,dy) = q2ð(t-s)(1-t) exp dy

- 2(t-s)(1-t)

1-s 1-s

Le pont brownien peut être vu comme la définition de la loi d'un point Xè à la date è ? [s,t], intermédiaire sur la trajectoire d'un mouvement brownien dont X_s et Xt sont connus, un tel point suit une loi normale

( ~

X_s + è - s

t - s (Xt - Xs), (t - è)(è - s)

N t - s

Proposition 2.9 Si (Wt)0=t=1 est un mouvement brownien standard, alors on a

(1). Xt = Wt -tW1 est un pont brownien.

(2). Xt = (1 -t)Wt/1-t est un pont brownien.

Preuve Il suffit de vérifier que le processus est gaussien et de covariance égale à s(1-t) ?0 = s < t = 1. Vérifions (1), on a

E(X_sXt) = E((W_s -sW1)(Wt -tW1))

= E(W_sWt) -tE(W_sW1) -sE(WtW1)+tsE(W²1 )

= min(s,t) -ts = s(1 -t)

De même pour (2), on a

E(X_sXt) = E((1 -s)(1 -t)Ws/1-sWt/1-t)

( s )

= (1 - s)(1 -t)min 1 - s, t

1 -t

= min(s(1 -t),t(1 -s)) = s(1 -t)

2.9.1 Construction par processus contraint

Soit (Wt)t=0 un mouvement brownien standard, a et b des réels quelconques. On définit le processus (Xè)s=è=t, comme la déformation de Wt, à partir de l'instant è = s forcée de passer par a à la date è = t. A tout instant, utilisant la nullité de l'espérance de dWt pour tout t, on écrit donc le drift du processus Xè comme la différence entre a et Xè, sur le temps restant avant t. On obtient donc une description du processus Xè comme :

a - Xè

dXè = t - è dè + dWè, avec X_s = b et Xt = a. (2.16)

L'équation différentielle stochastique (2.16) (voir le chapitre correspondant) admet pour solution:

è - s

_Xt,a

s,b(è) = b +Wè-s - (Wt-s -a+b) (2.17)

t -s

Nous pouvons simuler une trajectoire d'un pont brownien directement à partir de l'équation (2.17). On considère la subdivision de l'intervalle de temps [0,T] suivante 0 = t1 <
·
·
· < tN < tN+1 = T, avec _ti+1 --ti = Ät, pour i = 1 et W(0) = W(t1) = 0, on a l'algorithme suivant :

1. Générer une nouvelle variable aléatoire Z de distribution gaussienne N(0,1).

2. i = i+1.

_V'

3. W(ti) = W(ti--1) + Z Ät.

4.X(ti)=b+Wt_i -- ti T (WT--a+b)

5. Si i N + 1, réitérez a l'étape 1.

Remarque 2.2 Si a = b = 0 on a un pont brownien standard.

La fonction BB (Brownian Bridge) permet de simuler un pont brownien sur l'intervalle de temps [t0,T] avec un pas Ät = (T --t0)/N, et la fonction BBF permet de simuler un flux de pont brownien.

R> b <- - 2; a <- 1;

R> BB(N = 1000, t0 = 0 , T = 1 , x0 = b, y = a)

R> BBF(N = 1000, M = 100, t0 = 0, T = 1, x0 = 0, y = 0)





FIGURE 2.15 - Trajectoire d'un pont brownien à partir de Xt0 = --2 et XT = 1.

FIGURE 2.16 - Flux de 100 trajectoires d'un pont brownien standard Xt0 = XT = 0.

2.9.2 Construction par le développement de Karhunen-Loève (D.K.L)

Pour un pont brownien standard {Xt,0 =t = 1} est un processus gaussien centré, à trajectoires continues (p.s), tel que X(0) = X(1) = 0. Le développement de Karhunen-Loève [15, 16] s'écrit pour une suite de variables aléatoires Z_n de loi normale centrée réduite telles que E|Z² _n| = 1,

v Xt = 2

8 ? n=1

sin(ðnt) Zn ðn

?t ? [0,1] (2.18)

Utilisant le code 5 pour simulée un pont brownien standard a partir de D.K.L. La figure 2.17 donne une représentation graphique d'une approximation d'un pont brownien standard par le D.K.L pour n = 10,n = 100 et n = 1000.



FIGURE 2.17 - Approximation d'un pont brownien standard par le D.K.L.

2.10 Martingales exponentielles

D'autres martingales sont liées au mouvement brownien. Les considérer permet d'obtenir d'intéressants résultats. La martingale exponentielle ou martingale de Wald, permet en particulier d'obtenir des résultats concernant les processus à accroissements indépendants ayant une dérive.

On note

( )

óWt - ó2

Mt = exp 2 t , ?t = 0

Mt est une martingale par rapport au mouvement brownien, i.e., pour s < t,

E(Mt|W_s,s < t) = Ms

2.10.1 Caractérisation de Paul Lévy du mouvement brownien

Soit {Wt,t > 0} un mouvement brownien. Les processus suivantes sont des martingales continues relativement à la filtration (Ft)t>0

Xt = Wt (2.19)

Yt =Wt ² --t (2.20)

_ft

Zt = J f (s)dW_s oùf E L² (version continue de l'intégrale stochastique) (2.21)

0

Lt 2

ó

= exp ( óWt -- 2 t ) (2.22)

Preuve Pour (2.19), si s < t alors Xt -- X_s est indépendante de F_s donc on a E(Xt --X_s|F_s) = E(Wt --W_s|F_s) = E(Wt ^--Ws) = 0

d'où

E(Wt --W_s|F_s) = E(Wt|F_s) --W_s = ⁰ Pour (2.20), de même on a

E(W²_t --W²_s |F_s) = E((Wt --W_s)² + 2W_s(Wt --W_s)|F_s)

= E((Wt --W_s)²|F_s)+2E(W_s(Wt --W_s)|F_s) = E((Wt --W_s)²|F_s)+ 2W_sE((Wt --W_s)|F_s) = E(W²t--s) = t -- s

d'où

E(W²_t --t|F_s) = W_s² -- s Vs < t

Pour (2.21), Si f est une fonction étagée,

f = ? ëi¹]ti,ti+1] i

le processus Zt

Zt = f f (s)dWs = ?_iëi(Wti+1nt ^--Wtint)

est une martingale comme somme de deux martingales arrêtées. De plus,

2

E [(f_o^t f(s)dWs ) = E (f^t f²(s)ds)

0

Si f est dans L², il existe une suite de fonctions f_n qui converge vers f . La somme

f_o f_n(s)dW_s

est une suite de martingales. En passant à la limite, on voit que Zt est une martingale. Pour (2.22), on a

E(Lt|F_s) = E (_eóWt-^ó2 ²t+óW_s-óW_s

= e-
= e-

= e-
= e-

_ó2

^{2 t}E(e^{óWt+óWs-óWs})|F_s)

_ó2

2 tE ( eóWseó(Wt-Ws) |Fs)
2 teóWsE ( eó(Wt-Ws) |Fs)

ó2 t óW _E aWt_ ) 2 es e s

et comme on Wt est un processus gaussien à accroissements indépendants et stationnaires, donc Wt-s ~ N(0,t - s) = vt - sN(0,1), d'où

E(Lt|F_s) = e- ^ó2 t eóWsE (eóvt-sN(0,1))

= e- ó2 2 teóWseó2 2 (t-s)

_ó2

= exp ( óW_s - 2 s)

D'une manière plus générale, le processus Mt est une martingale

_rt

Mt = exp f (s)dW_s - f f² (s)ds)oùf ? L2

0 2 f

En prenant pour fonction ëf on construit une infinité de martingales, le processus Mt devient

ë2 ft 2

Mt = exp (ë f f(s)dW_s - 2 0 f (s)ds) 0

En dérivant selon le paramètre ë, on obtient de nouvelles martingales. Pour f = 1, on trouve que Mt = exp ( ëWt - 22 t) est une martingale. En dérivant, on obtient une suite de martingales ?Mt/?ë|ë=0 = Wt, ?2Mt/?ë2|ë=0 = Wt ² -t, ?³Mt/?ë³|ë=0 = Wt 3 - ³tWt,...

Théorème 2.2 (Caractérisation de Paul Lévy du mouvement brownien) Soit Xt un processus aléatoire réel continu. Pour que Xt soit un mouvement brownien il faut et il suffit que Xt et X2 t -t soient des martingales pour la filtration propre de Xt.

Autrement dit, toute martingale continue est un mouvement brownien. L'hypothèse continue est essentielle. Le processus de Poisson Nt est une martingale, même que N2 t - t, mais il n'est pas continu et donc n'est pas un mouvement brownien.

2.10.2 La loi du temps d'atteinte du mouvement brownien

Soit Wt un mouvement brownien, les martingales exponentielles permettent de calculer les lois des temps de passage. Soit a > 0, on définit T_a le temps d'arrêt (Chapitre 1 Section 1.9.1), avec T0 = 0

T_a = inf{t > 0 : Wt = a}

représentant le premier instant où le mouvement brownien atteint la valeur a. On pose T_a =
inf(/0) = 8. On peut calculer la transformée de Laplace de la variable aléatoire T_a, en considérant

la martingale

( )

óWt - ó2

Ut = exp 2 t

En arrêtant Ut à l'aide de T_a, on introduit la martingale

( )

óWt?T_a - ó2

Xt = Ut?T_a = exp 2 (t ? T_a)

qui converge vers

~ ~

óa - ó2

exp 2 T_a ¹{T_a<8}

Pour ó > 0, Xt est une martingale bornée puisque _Wt?Ta est majorée par a. On peut donc appliquer
la relation (1.2) à la martingale Xt et au temps d'arrêt T_a, ce qui donne E(XT_a) = E(X0) = 1.

Puisque

( )

óa - ó2

XT_a = exp 2 Ta

Si l'on pose è = ó²/2, on obtient :

v

E(exp(-èT_a)) = exp(-a 2è) (2.23)

On obtient ainsi la transformée de Laplace de la variable aléatoire T_a, cette transformée de Laplace peut être inversée, et cela montre que T_a admet une densité sur R⁺ donnée par:

/ ~

a _-a²

f_a(x) = v exp ; a > 0, t > 0 (2.24)

2ðx³ 2x

La loi de T_a s'appelle une loi stable d'indice 1/2 où distribution de Lévy. En particulier, on en déduit que

P(T_a < 8) = 1 et E(T_a) = 8

Des formules similaires existent avec -a au lieu de a si a est négatif, par symétrie du mouvement brownien. Le processus -Wt est encore un mouvement brownien, et le temps d'atteinte de a par -Wt est égale au temps d'atteint de -a par Wt.

Exemple 2.1 Soit Wt un mouvement brownien standard. On note

T_x = inf{t = 0,Wt = x}

et on considère une variable exponentielle X de paramètre 1, alors on a

Z 8 v

2x

P(T_x < X) = 0 P(T_x < t)e^-tdt = E(e^-Tx) = e-

2.10.3 Les temps de passage du mouvement brownien

Considérons maintenant a > 0 et b > 0, et soit T = min(T_a,T-b). La martingale Mt =Wmin(T,t) est bornée, donc uniformément intégrable, et on peut appliquer la relation (1.2) au temps T. cela entraîne que :

0 = E(M0) = E(MT) = aP(T = T_a) - bP(T = T-b)

Mais {T = T_a} = {T_a < T-b} et {T = T-b} = {T-b < T_a}, de telle sorte que finalement on obtient

b P(Ta < T-b) =	P(T-b passage = = =	a < Ta) =	(2.25) (2.26) (2.27)
a + b et Proposition 2.10 Les temps de E(exp(-èT-b)1{T-b<Ta} ) E(exp(-èTa)1{Ta<T-b} ) E(exp(-è(Ta ?T-a)))		a+b sont donnés par sinh(a 2è) v
		sinh((b + a) 2è) v sinh(b 2è) v
		sinh((b + a) 2è) v 1
		cosh(a 2è) v

( ) ( )

Preuve Pour (2.25) on a exp óWt - ó2 2 t et exp -óWt - ó2 2 tsont des martingales, le processus

Xt = (áexp(óWt) + âexp(-óWt))exp(-ó²t/2)

est une martingale. Soit T = inf(T_a,T-b) un temps d'arrêt. Le processus _Xt?T est une martingale bornée (puisque -b = Wt?T = a). Lorsque t tend vers l'infini, le processus _Xt?T tend vers

(áexp(óWT) + âexp(-óWT))exp(-ó²T/2)

oil WT vaut

WT = a¹{T_a<T-b} - b¹{T-b<T_a}

Choisissons a et p de sorte que aexp(aa) + pexp(-aa) = 0. Par exemple en prenant, a = exp(-aa)/2 et p = -exp(aa)/2. Le processus

Xt = sinh(a(Wt - a))exp(-a²t/2)

est une martingale. Donc le processus arrêté

Xt?T = sinh(a(Wt?T - a))exp(-a²(t ? T)/2)

est aussi une martingale. Par conséquent, d'après le théorème d'arrêt (Chapitre 1 Section 1.9.2) E(XT-b) = E(X0) = sinh(-aa) ? sinh(a(-b - a))E(exp(-a²T-b/²)¹{T-b<T_a}) = sinh(-aa) et comme la fonction sinh est une fonction impair (sinh(-x) = -sinh(x)), donc on a

E(exp(-a2T-b/2)1{T sinh(aa)

_b<T_a}) = .

smh(a(b + a))

Si l'on pose 0 = a²/2, on obtient

E(exp(-⁰T-b)¹{T-b<T_a}) =

sinh(av20)

sinh((b + a) v20)

On démontre la formule (2.26) de la même manière, choisissons a et p de sorte que aexp(-ab)+ pexp(ab) = 0, en prenant a = exp(ab)/2 et p = -exp(-ab)/2. Dons le processus

Xt = sinh(a(Wt + b))exp(-a²t/2) est une martingale. Donc on le processus arrêté

Xt?T = sinh(a(Wt?T + b))exp(-a²(t ? T)/2) est aussi une martingale. Appliquons le théorème d'arrêt

E(XT_a) = E(X0) = sinh(ab) ? sinh(a(a + b))E(exp(-a²Ta/²)¹{T_a<T-b}) = sinh(ab)

2 sinh(ab)

? E(exp(-aTa/2)1 {T_a<T b}) = sinh(a(a + b))

On a 0 = a²/2, on obtient

sinh(bv20)

E (exp(-0Ta)1{Ta<T-b}) -- sinh((b + aW20)

La troisième formule (2.27) est la somme des deux premières (2.25) et (2.26) dans laquelle on prendra a = b, donc on

2 sinh(av20)

E(exp(-0(T_a ?T-_a))) =

sinh(2av20)

On sait que sinh(2x) = 2 sinh(x) cosh(x) (x = av20), d'oil on trouver

1

E(exp(-0(T_a ? T-a))) =

cosh(av20)

2.11 Conclusion

Le mouvement brownien, ou processus de Wiener, joue un rôle fondamental dans de nombreux domaines. A l'heure actuelle, un rôle important en mathématiques financières. Comme on a pu le constater dans ce chapitre, il y a plusieurs façons de simulée une processus stochastique d'une variable aléatoire. Le processus de Wiener est à la base de toutes ces représentations stochastiques. Il entre justement dans la formulation de la partie stochastique de ces divers modèles. Le mouvement brownien arithmétique est le plus simple des processus stochastiques. Mais, pour modéliser le prix d'une action, son principal désavantage est que le rendement de cette action est une fonction décroissante du prix de l'action. On recourt par conséquent au mouvement brownien géométrique pour représenter le mouvement stochastique du prix d'une action. Dans ce chapitre, nous avons présenté le théorème de Paul Lévy qui caractérise un mouvement brownien, et on a montre l'importance du théorème d'arrêt dans les calculs des temps de passages, à cause de la propriété de martingale.

3

Processus de Diffusion

Sommaire
3.1	Introduction	47
3.2	Intégrale stochastique	48
3.3	Equations différentielles stochastiques	56
3.4	Schémas numériques	65
3.5	Les modèles attractives	75
3.6	Conclusion	83

3.1 Introduction

A

fin de prendre en compte les modélisations de force aléatoires qui interviennent dans les équations de la dynamique comme des termes complémentaires, on a cherché à donner un sens à l'intégrale lorsque celle-ci est prise par rapport à un processus stochastique. Ce chapitre est consacré dans un premier temps à étudier l'intégrale

Za b f(t,ù)dW(t,ù)

dans la quelle Wt est un mouvement brownien standard et f un processus stochastique. Si f est une fonction déterministe de classe C¹, l'intégrale est une intégrale de Stieltjes classique. Mais si f dépend aléatoirement de la variable ù, comme le mouvement brownien standard n'est nulle part différentiable et que presque toutes ces réalisations n'ont pas de variations bornées, l'intégrale n'a pas de sens. [27] a été le premier à donner un sens à cette intégrale. À la suite des travaux d'Itô, Stratonovitch a proposé une autre définition de l'intégrale stochastique, la construction de l'intégrale d'Itô est assez technique. Dans le contexte de l'intégrale stochastique, nous employons indifféremment la notion de processus de Wiener et celle de mouvement brownien.

Les systèmes apparaissant dans les applications sont souvent soumis à des perturbations que l'on peut considérer comme aléatoires. Dans le cas où le bruit åt est un bruit blanc gaussien, on peut modéliser l'équation par un mouvement brownien et la traiter comme une équation d'Itô. On note X(t) = (X1(t),...,X_n(t)) le vecteur décrivant l'état du système à l'instant t. On va considérer des processus Xt qui sont solution d'une forme perturbée de l'équation différentielle :

dXt
dt

= u(t,Xt) +ó(t,Xt)åt (3.1)

où le vecteur X(t) = (X1(t),...,X_n(t)) décrit l'état du système à l'instant t. L'équation (3.1) peut s'écrire formellement :

dXt
dt

= u(t,Xt) + ó(t,Xt)^dWt (3.2)

dt

où dWt

dt est la dérivée formelle par rapport au temps d'un mouvement brownien Wt (Chapitre 2 Section 2.4). En fait, la dérivée dWt

dt n'existe pas (Dans la littérature d'ingénieur, la dérivée de Wiener est souvent appelée bruit blanc). L'équation (3.2) doit être réécrite avec la notation différentielle

dXt = u(t,Xt)dt + ó(t,Xt)dWt (3.3)

avec u(t,Xt) et ó(t,Xt) sont des fonctions mesurables localement bornées sur Rⁿ. Dans le cas où le processus Xt est markovien, on a ce qu'on appelle une diffusion.

Nous décrivons dans ce chapitre des méthodes numériques pour la résolution d'équation différentielles stochastiques (ÉDS). Pour approcher numériquement l'équation (3.3), on utilise

des schémas aux différences classiques et le fait que pour un pas h donné, les variables W(n+1)h - Wnh suivent des lois gaussiennes indépendantes de variance h. On note xt le processus approché et on considère la subdivision

0 = t0 < t1 < ··· < tN-1 < tN = T

de pas régulier

h = Ät = tn+1 -t_n

Dans le cas multidimensionnel u = (u1,...,ud) et Xt = (X1(t),...,Xd(t)) sont des vecteurs de R^d. Le mouvement brownien a p composantes Wt = (W1(t),...,Wp(t)) et ój = (ó¹ _j,ó² _j,...,ó^d_j) pour j = 1,...,d. L'équation (3.3) s'écrit

p

Xt = u(t,Xt)dt + ? ój(t,Xt)dWj(t)

i=1

et se traite de la même manière.

La dernière partie de ce chapitre concerne l'utilisation des processus de diffusion dans la modélisation et l'étude analytique et par simulation de deux problèmes :

· Le premier problème est de modéliser une trajectoire d'un polluant, qui se déplace sur une surface d'eau turbulente en présence d'un mécanisme d'attraction, qui est donnée par (Boukhetala 1994,1996 [2, 4]).

· Le deuxième problème porte sur une modélisation d'un phénomène d'attraction entre deux insectes mâle et femelle, qui est donnée par (Boukhetala 1998 [5]).

3.2 Intégrale stochastique

3.2.1 L'intégrale d'Itô

On considère une espace de probabilité (Ù,.i,P) muni d'une filtration 3t, i.e. une suite de tribus croissantes pour l'inclusion. On appelle tribu des prévisibles sur Ùx [0,8[ la plus petite tribu rendant mesurable tous les processus continus adaptés à la filtration 3t. Un processus ou un ensemble est prévisibles s'il est mesurable par rapport à cette tribu. Supposons donné un processus de Wiener standard Wt, adapté à la filtration 3t et tel que pour tout 0 s t l'accroissement Wt(ù) -W_s(ù) soit indépendant de 3t. Sur un intervalle de temps [a,b], on note H² l'ensemble des processus f(t,ù) définis pour t E [a,b], 3t-mesurable et de carré intégrable presque sûrement. Dans ces conditions, si f est dans H² et si a = t0 < t1 < ··· < t_n < tn+1 = b est subdivision de l'intervalle [a,b], alors f est indépendant des incrementsWt_j+1 -Wtj, en d'autre termes f est prévisible. Pour toute fonction f de H², on définit l'intégrale stochastique d'Itô comme la limite dans L² des accroissements ci-dessous (on notera que toutes les limites de cette section sont des limites quadratiques, i.e. dans L²). On définit ainsi l'intégrale stochastique comme la limite des sommes Riemann.

Za

b n

f(t,ù)dW(t,ù) = lim ?f(ti,ù)(W(ti+1) -W(ti))

^n?8i=0

Cette définition est cohérente avec les propriétés usuelles de l'intégrale. On a de plus quelques propriétés complémentaire liées à la dépendance aléatoire de f

(1) Si f ? H² et si f_a^b E(f²(t))dt < 8, alors

Za b E(f²(t))dWt = 0

(2) Si f ,g ? H² et si R b ~E~ f ²(t) + Eg²(t)dt < 8, alors

a

~Z b ~~Z b ~ Z b

E _a f (t)dWt a g(t)dWt = a E(f (t)g(t))dt

En particulier, on a

~Z b _~2 Z b

a E f ²(t)dt

E _a f (t)dWt =

Le résultat suivante montre que l'intégrale stochastique est une martingale. Si f (t,ù) ? H² et si pour tout t ? [a,b], la somme f_a^b E(f²(t))dt < 8, alors l'intégrale stochastique

Z b

a

f (t,ù)dW(t,ù)

est une martingale (la démonstration dans la page 40) est ces trajectoires sont p.s. continues.

Exemple 3.1 Calcul de I0 = ft₀^tWsdWs

Considérons la subdivision de l'intervalle [t0,t] : t0 < t1 < ··· < t_n < tn+1 = t et appliquons la définition. On rappelle que toutes les limites sont des limites prises dans L²

I0 = lim

n ? i=0

wti [Wti+1 -wti]

En considérant la somme

I1 = lim

n ? i=0

Wti+1 [Wti+1 -wti]

on peut former la différence et appliquer les propriétés du mouvement brownien

I1 -I0 = lim

n ? i=0

[Wti+1 -Wti]^2 = t -t0

De manière générale, considérons

Ië = (1 -ë)I0 +ëI1

On a Ië = lim	n ? i=0			[ëwti+1 + (1 - ë)Wti][Wti+1 -^wti]
= lim	?n h ti + ë~Wti+1 -Wti ~2i WtiWti+1 -W2 i=0
= lim	n ? i=0	h i ~ ~ 1 W2 ti+1 -W^2 ë - 1 + (t -t0) ti 2 2

~ ~

1 ~W2 ~ +

= t -W² ë - 1 (t -t0)

t0

2 2

d'où la valeur de l'intégrale d'Itô obtenue pour ë = 0

rt 1t - t0

I0 = J2 iW2 -W0 J

2i

2

t0

Remarquons que ce type de calcul ne respecte pas les lois du calcul différentiel ordinaire puisqu'on voit apparaître un terme supplémentaire. Pour t0 = 0 l'équation

ft_o

WsdWs =¹ 2W2 t - t 2

peut s'écrire

Wt 2 = Z ds + 2 f W_sdW_s

soit sous forme différentielle

dW_t² = dt +2WtdWt

ce qui montre clairement le terme supplémentaire (dt). Dans les calculs stochastiques, la dérivation doit être poursuivie à l'ordre 2, les termes croisés sont nuls dWt · dt = dt · dWt = 0, de même que le terme dt · dt = 0. En revanche, le terme dWt · dWt = dt donne naissance à un terme supplémentaire.

3.2.2 L'intégrale de Stratonovitch

Qui est définie par la limite quadratique

Za	lim? ^nf(ti,ù )+ f (ti+1,ù) (Wti+1 -Wt-i) b f(t,ù)?dW(t,ù) = 2 i=0

satisfait les règles du calcul différentiel ordinaire. En reprenant le calcul précédent pour ë = 1/2, on trouve

_ft

I1/2 = J

r0

ws ? dW_s = 21 [Wt 2 Wt02]

Dans ce cas, il n'y a pas de terme supplémentaire. On préfère toutefois développer le calcul stochastique en utilisant la définition d'Itô.

La fonction ItovsStra permet de simuler l'intégrale stochastique I0 par les deux types Itô et Stratonovitch, sur l'intervalle de temps [0,T] avec un pas Ät = T/N.

Nous avons observé lors de la construction de l'intégrale stochastique I0 au sens d'Itô dans la figure 3.1, que les trajectoires du processus {J 0 t W_sdW_s,0 < t < T} pouvaient être négatives à certains instants.

R> ItovsStra(N = 1000, T = 1)



FIGURE 3.1 - Simulation l'intégrale stochastique f 0 t WsdWs vs J 0 t W_s ?dW_s.

3.2.3 Processus d'Itô

On appelle processus d'Itô, un processus {Xt,0 < t < T} à valeur dans R tel que :

Z t Z t

Xt = X0 + ₀ p_sds + 0 asdWs (3.4)

avec

(1) p = {pt,0 <t < T} eta = {at,0 <t < T} sont des processus adaptés à la filtration {Ft}t=0

[i T i

(2) P ₀ |p_s|ds < 8 = 1

[i T ]

(3) P ₀ |ó_s|²ds < 8 = 1

Écrit sous sa forme différentielle, le processus d'Itô devient

dXt = utdt + ótdWt (3.5)

Remarque 3.1 Par conséquent, si le coefficient de dérive est nul, c'est-à-dire que ut = 0 pour tout 0 = t = T, alors le processus d'Itô Xt

Z t

Xt = X0 + 0 ósdWs (3.6)

[i 0 T |ós|2]

est une t-martingale si et seulement si E ds < 8 est vérifier.

3.2.4 Formule d'Itô

Pour un processus stochastique Xt vérifiant

Z t2 Z t2

X(t2) - X(t1) = u(t)dt + ó(t)dWt

t1 t1

on note sous forme différentielle

dXt = u(t,Xt)dt + ó(t,Xt)dWt

Si f (t,x) est une fonction de classe C², alors f (t,Xt) admet une intégrale stochastique par rapport au même processus de Wiener donnée par la formule d'Itô

(?f(t,Xt) ~

?t + u(t)? f (t,Xt)

d f (t,Xt) = ?x + ¹ 2ó2(t)?2 f (t,Xt) dt + ó(t)? f (t,Xt)

?x dWt (3.7)

?x²

Cette formule permet de trouver les solutions de l'équation différentielle stochastique (3.3), et le calcul d'espérance d'un processus donné. Elle permet aussi de montrer que certains processus sont des martingales (puisque l'intégrale stochastique est une martingale).

Par exemple, si f (x) = x², alors

f (Wt) = W2 t et f (W0) = W2 0 = 0

et

? f ?² f

_?w(Ws) = 2W_s et ?w2 = 2.

En remplaçant dans la formule (3.7), nous obtenons

dW²

t = dt + 2WtdWt

sous forme intégrale

Wt 2 = 2 f WsdWs + f ds = 2 f WsdWs +t

0 0 0

ce qui implique

L.

t W2 t -t

WsdWs = 2

D'une manière générale la dérivé stochastique du processus Xt = Wn t . Posons f(t,x) = xⁿ et appliquons la formule d'Itô. Il vient

1

dWt n = 2 n(n-1)W_t ^n-2dt +nW_t ^n-1dwt

sous forme intégrale

1 n-1 ft

W^n-2ds

on 2

Wt = n(n --1) f W^n-2ds+n f

, 2 0 s _wn- 1 dWs ^f W^n-1dW_s = 1 Wtn

0

Pour n = n+ 1, on obtient

fo

1 ⁿ f t

WⁿdW = Wⁿ⁺ⁱ - W_s ^n-1ds

s s n1 t 2 0

Exemple 3.2 Soit le mouvement brownien géométrique oil le modèle de marché de Black et Scholes (Chapitre 2 Section 2.8)

dSt = èStdt +óStdWt (3.8)

on chercher une solution unique explicite à ce modèle, nous appliquons la formule d'Itô (3.7) à l'équation (3.8), posant f (t,x) = ln(x), on a

dln(St) = ₍0 +èSt ¹ + ¹ ó²S² dt +óSt 1 dWt

St 2 t S2 St

t

~ = - 12 ó²) dt + ódWt sous forme intégrale

~ln(St) = ln(S0) + è -

0

1 ^a2) f ^{t d s} + ^a I dW s

0 rt

2

= ln (S0) + (è - ₂ó²) t + óWt

1

d'oil la solution de l'équation (3.8), pour t = 0 et S0 > 0 est

~~ ~ ~

è - ó2

St = S0 exp t + óWt

2

Exemple 3.3 Soit à calculer la dérivée stochastique du processus Yt = e^Wt, oil Wt est un mouvement brownien standard. Posons Xt = Wt, f (t,x) = e^x et appliquons la formule (3.7). il vient

1

d (eWt) = ₂e^Wtdt + e^WtdWt

soit sous forme intégrale

eWt = 1 + 12 o ft e^Wsds + f te^WsdW_s

o

En prenant l'espérance de chaque membre et en appliquant les règles de calcul du paragraphe précédent (l'espérance de l'intégrale en dW_s est nulle).

^t

E (eWt) = 1 + 2 JO E (eWs) ds

1

En posant

y(t) = E(eWt)

l'équation précédente est une équation différentielle du premier ordre déterministe en y(t) avec comme condition initiale y(0) = 1

y (t) = 21 y(t)

qui admet comme solution y(t) = E(eWt) = e^t/²

Remarquons que l'application de la formule d'Itô à la fonction f (t,x) = e^iux permet le calcul de la fonction caractéristique de Wt.

=

E (_eiuWt) _e-tu²/2

Formule d'Itô Multidimensionnel

Soit Xt un processus de dimension m vérifiant l'équation

dXt = A(t)dt + B(t)dWt (3.9)

oil A(t) = (u1(t),...,u_m(t)) et B(t) = (óij(t)) est une matrice m x n. Le mouvement brownien étant de dimension n Wt = (W1(t),...,W_n(t)). Soit f (t,x) une fonction définie pour t E [a,b] et x un vecteur de R^m à valeurs dans Rp de classe C², alors le processus f (t,Xt) admet une dérivée stochastique donnée par

df (t ,Xt) = + 2, ui(t) + L ? ói,k(t)ó k(t) i, ?x2

?t i=1 ?x 2 k=1i, j=1 ?2 f(t,Xt)) dt

(?f (t ,Xt) _`--,^m ?f (t ,Xt) 1 ,-ⁿ -, m

(3.10)

n

+ ?

k=1

m

?

i=1

? f(t,Xt)
ói,k(t) dWk(t)
?x

Dans le cas, oil la fonction f est indépendante du temps et Xt = Wt est une mouvement brownien sur Rⁿ, la formule d'Itô se simplifie en

t

f(Wt) = f(W0)+ _f ?f(W,0dW_s + 1 f ^t Äf (WOds (3.11)

0 2 0

3.2.5 La règle de multiplication

La règle de multiplication ou intégration par parties est utile lorsque nous voulons étudier, par exemple le comportement du prix actualisé d'un actif lorsque nous connaissons les processus pour l'évolution du prix de l'actif et celui du facteur d'actualisation.

Théorème 3.1 (Intégration par parties) Soit Xt et Yt deux processus d'Itô tels que

dXt = utdt + ótdWt et dYt = ^-utdt +^-ótdWt (3.12)
la formule d'Itô (3.10) conduit à la formule d'intégration par parties

dXtYt = XtdYt +YtdXt + d(X,Y)t

(3.13)

= (^-utXt + uYt + ót^-ót)dt + (^-ótXt + óYt)dWt

où la variation quadratique est caractérisée par

d(X,Y)t = ót^-ótdt

Preuve En particulier lorsque les browniens sont les mêmes et que la matrice B(t) se réduit à un vecteur, la dérivée s'écrit

df = ( ?t f + _ut?x f +^-ut?_y f + 2 t

(ó²?_xx f +2ót^-ót?xy f + ó?yy f)) dt + (ót?x f + ^-ót?y f)dWt appliquons cette formule a la fonction f(t,x,y) = xy, on trouve directement le résultat

dXtYt = (utYt +^-utXt + ót^-ót)dt + (ótYt +^-ótXt)dWt

= utYtdt +^-utXtdt + ót^-ótdt + ótYtdWt +^-ótXtdWt

= Yt(utdt + ótdWt) + Xt(^-utdt + ^-ótdWt) + ót^-ótdt

= XtdYt +YtdXt +d(X,Y)t

Exemple 3.4 (La valeur actualisée d'un actif risqué) Supposons que le processus stochastique St satisfaisant l'équation (3.8), est l'évolution du prix d'un titre risqué. Le processus {Yt = _e--rt ,t > 0} est notre facteur d'actualisation. Notons que

d dtYt = --re^--rt = --rYt

c'est-à-dire que le processus Yt est un processus d'Itô satisfaisant l'équation différentielle

dYt = --rYtdt

le processus Zt = YtSt represente l'evolution de la valeur actualisee du prix du titre. La règle de multiplication (3.13) entraîne que

dZt = dYtSt

= YtdSt + StdYt + d(S,Y)t

= Yt(èStdt + óStdWt)+ St(--rYtdt)

= (è -- r)StYtdt + óStYtdWt

= (è -- r)Ztdt + óZtdWt

sous sa forme integrale, cette dernière equation s'ecrit

t t

f

Zt = Z0 + (è -- r) f Z_sds

+ ó ZsdWs 0 0

3.3 Equations différentielles stochastiques

3.3.1 Introduction et définitions

De manière informelle, on appelle equation differentielle stochastique une equation differentielle ordinaire perturbee par un terme stochastique. Plus precisement, c'est une equation du type suivant :

dXt = u(t,Xt)dt +ó(t,Xt)dWt, X0 = x0 (3.14)

Dans cette equation, dWt est la differentielle d'un mouvement brownien standard Wt, et u, ó sont les coefficients de l'equation (ce sont des fonctions de R⁺ x R dans R), et x0 E R est la valeur initiale. Tous ces termes sont donnes. La notation (3.14) est la plus usuelle.

Définition 3.1 Rechercher une solution de l'equation (3.14) consistera à rechercher un processus {Xt,t > 0f satisfaisant l'equation integrale :

t t

Xt = x0 + f u(s,X_s)ds + f ó(s,X_s)dW_s (3.15)

0 0

oil la seconde integrale est une integrale stochastique.

L'equation (3.14) ou l'equation (3.15) etaient jusqu'à present unidimensionnelles. On peut egalement definir une equation d-dimensionnelle de la manière suivante. Le processus inconnu

s _tX = (Xi)1 <1 .<d est une famille de processus a valeurs reelles Xi = (Xi)t>0, la condition initiale

_ _ _

x0 = (xi0 ) 1<i<d appartient à R^d, le mouvement brownien W = _(W /) est q-dimensionnel,

1<j<q

et les coefficients ont les dimensions appropriees, soit u = (ui)1<i<d et ó = (ói,j)1<i oil

<d,1<j<_q'

les coefficients uⁱ et ó^i,j sont des fonctions de R⁺ x R dans R. On ecrit encore l'equation sous les formes (3.14) ou (3.15), mais cela signifie maintenant que l'on a :

t

i

Xⁱ_t =xi + f p(s x)ds+ ? a^i,j(s,X_s)dW_s^j;i= 1, .. . ,d (3.16)

j=1

0 , s q f^t --

o 0

Définition 3.2 Quand les coefficients p et a ne dependent pas du temps et sont seulement des fonctions définies sur R^d, on dit que l'équation est homogène.

Le coefficient p est appelé le coefficient de dérive, tandis que a est le coefficient de diffusion. Un processus qui résout l'équation (3.14), ou de manière équivalente (3.16), est appelé processus de diffusion ou, plus simplement, une diffusion.

Définition 3.3 (Processus de diffusion) Un processus de diffusion est un processus de Markov à trajectoires continues vérifiant l'équation d'Itô (3.14). Soit (Xt)t=0 un processus aléatoire défini sur l'espace de probabilité (Ù,A,P) à valeurs réels, muni d'un filtration (Ft,t = 0). On dit que Xt est une processus de diffusion caractérisée par

(1) la limite donnant la dérive

lim h?0	E(Xt+h -Xt|Xt = x)	= p(x,t)
	h

(2) la limite donnant la diffusion

lim h?0	E((Xt+h -Xt)^2|Xt = x)	= a^2(x,t)
	h

(3) la condition de Dyukin

?å > 0, lim h?0	P(|Xt+h -Xt| > å|Xt = x)	= 0
	h

Remarque 3.2 Le mouvement brownien standard est un processus de diffusion avec coefficient de dérive nulle et coefficient de diffusion égale à 1.

Preuve

1

1

lim h?0 h

E(Wt+h -Wt|Wt = x) = lim

h?0 h

E(Wt+h -Wt|Wt -W0 = x)

= lim

1 _hE(Wt+h -Wt)

h?0

1

= lim _h(E(Wt+h) - E(Wt)) = 0

h?0

lim 1
h?0 h

E((Wt+h -Wt)²|Wt = x) = lim _hE((Wt+h -Wt)²|Wt -W0 = x)

1

h?0

1

= lim_hE~(Wt+h -Wt)²~ h?0

1

= lim _hh = 1

h?0

Exemple 3.5 Le processus d'Ornstein-Uhlenbeck est la diffusion solution de l'équation (de langevin Section 3.3.3)

dXt = -rXtdt +ódWt

où r et ó sont des constantes positives et Wt un processus de Wiener standard. On suppose de plus qu'à l'instant initial X0 = x0. En posant

Yt = Xte^rt

et en appliquant la formule d'Itô (3.7) à la fonction f(t,x) = xe^rt, on obtient

dYt = óe^rtdWt

sous forme intégrale

Z t

Yt = Y0 + ó e^rsdW_s

t0

d'où la solution pour t = t0,

Z t

Xt = Xt₀e^-r(t-t0) + ó e^-r(t-s)dW_s

t0

En notant x0 = E(Xt₀), le processus d'Ornstein-Uhlenbeck a pour moyenne

E(Xt) = x0e-r(t-t0)

En appliquant la formule d'Itô au processus Zt = X2 t , on a

dZt = -2rZtdt + ó²dt + 2ódWt

En prenant la moyenne, on trouve en résolvant une équation différentielle ordinaire

(

E(Zt) = E(X2 _t ) = ó2 1 - e-2r(t-t0)) + E(X2 t0) 2r

D'où la variance du processus d'Ornstein-Uhlenbeck

(

var(Xt) = E(X2 _t ) - [E(Xt)]² = ó2 1 - e-2r(t-t0))

2r

La fonction de corrélation du processus d'Ornstein-Uhlenbeck vaut

ó2 e-r|t-s|

R(t,s) =

2r

Nous pouvons simuler une trajectoire du processus d'Ornstein-Uhlenbeck comme suit. On considère la subdivision de l'intervalle de temps [t0,T] suivante t0 < t1 < ··· < tN < tN+1 = T, avec _ti+1 -ti = Ät, pour i = 1 on a W(t0) = W(t1) = 0 et X(t0) = x0, on a l'algorithme suivant :

1. Générée un nouveau variable aléatoire Z de la distribution gaussienne N(0,1).

2. i = i+1. v

3. W(ti) = W(ti-1) + Z Ät.

4. Ii = e-^r(ti+¹-^ti)(Wi+1 -Wi)

5. X(ti) = X(t0)e-^r(ti-¹-^t0) + a?i _j=1 Ij.

6. Si i = N + 1, réitérez a l'étape 1.

La fonction OU permet de simuler une seule trajectoire de Xt dans l'intervalle [t0,T] avec un pas Ät = (T -t0)/N, et la fonction OUF permet de simuler un flux de Xt.

R> OU(N = 1000, t0 = 0, T = 10, x0 = 10, r = 2, sigma = 1)

R> OUF(N = 1000, M = 100, t0 = 0, T = 10, x0 = 10, r = 2, sigma = 1)



FIGURE 3.2 - Trajectoire d'un processus d'Ornstein-Uhlenbeck avec r = 2 et a = 1.

FIGURE 3.3 - Flux d'un processus d'OrnsteinUhlenbeck avec r = 2 eta = 1.

3.3.2 Existence et unicité des solutions de l'ÉDS

Soit T > 0, u(t,x) une fonction mesurable de [0,T] × Rⁿ dans Rⁿ et a(t,x) une fonction mesurable de [0,T] × Rⁿ dans Rⁿ vérifiant :

(1) Condition de croissance : il existe une constante C telle que

|u(t,x)| + |a(t,x)| = C(1 + |x|)

(2) Condition de Lipschitz : il existe une constante K telle que

|u(t,x) -u(t,y)| + |a(t,x) -a(t,y)| = K|x-y|

< 0°,

(3) X0 est une variable indépendante de la tribu a(W_s,s = 0) et E X2 0

Alors l'équation différentielle stochastique d'Itô (3.14) admet une solution unique Xt dont presque toutes les réalisations sont continues, qui est un processus de Markov de loi initiale X0 et de probabilité de transition

p(s,x,t,A) = P(Xt ? A|X_s = x)

De plus, si les fonctions u et a sont continues, alors Xt est un processus de diffusion de dérive u(t,x) et de matrice de diffusion aa^T (symétrique et semi-définie positive). En particulier lorsque u et a sont lipschitziennes, l'équation (3.14) admet une solution unique.

Remarque 3.3 (1) Une équation peut admettre une solution locale sans admettre de solution globale. Par exemple sur [0,1], l'équation

dXt
dt

= X2 t , X0 = 1

admet une solution unique locale sur [0,1[

Xt = 1 -t , t ? [0,1[

1

mais n'admet pas de solution globale sur l'intervalle [0,1].

(2) La condition de croissance évite que les solutions explosent, i.e. que les solutions |Xt| tendent vers l'infini en un temps fini. L'équation

dXt = ₂aX³

1 t dt, Xt = y 1 et y > 0

admet une solution

Xt =

1

\/

y2 - at

qui diverge dans la direction y² = at. Si la condition de croissance n'est pas satisfaite, l'équation peut quand même avoir une solution.

(3) La condition de Lipschitz garantit l'unicité. L'équation

dXt = 3X^2/3

t dt, X0 = 0

admet plus d'une solution

Xt = (t-a)³1t>_a

car la fonction u(t,x) = 3x²/³ ne vérifie pas la condition de Lipschitz.

3.3.3 Equation de Langevin

Des particules dans un fluide sont soumises à une force de frottement proportionnelle à la vitesse de ces particules et à une force aléatoire. La vitesse des particules est processus aléatoire. L'équation fondamentale de la dynamique s'écrit

dv
dt

= -av + c(t)

avec a = 6mir/m, la masse d'une particule étant m et i désigne la viscosité du fluide. Les particules sont assimilées à des petites sphères de rayon r. La force aléatoire c(t) est un bruit blanc, son espérance est nulle et sa fonction de corrélation est constant, D est appelé coefficient de diffusion d'Einstein.

E(c(t)) = 0 R(s,t) = 2D8(t -s)

En termes mathématiques, ce problème s'écrit sous la forme d'une équation d'Itô

dXt = -aXtdt +v2DdWt (3.17)

t

Xt = X0e^-at + v2D f e^a(t-s)dW_s

o

Wt est un processus de Wiener standard. Appliquons la formule d'Itô (3.7) à Yt = Xte^at, en posant f (t,x) = xe^at, on a at f = _axeat_, a_x f =e^at et a_xx f = 0. On en déduit que la solution est le processus de diffusion, markovien, gaussien Xt donné par

De cette expression, on posant x0 = E(X0), on trouve l'espérance

E(Xt) = x0e^-at

et la variance du processus

var(Xt) = D(1 -e-2at)

a

La fonction de corrélation du processus est

R(s,t) = D a

e-a|t-s|

En supposant que var(X0) < 00, il s'en suit que lorsque t ? 00, on a

E(Xt) ---? 0

et

var(Xt) ---? D/a

La distribution de Xt approche de la loi N (0, Da ) lorsque t ? 00. Ainsi, quelque soit la distribu-
tion initiale, la solution de l'équation différentielle stochastique pour des temps très grands suit
approximativement une loi gaussienne centrée et de variance D/a (Chapitre 4 exemple 4.2), qui

représente un équilibre entre la force aléatoire de perturbation dWt (que nous avons appelé dans l'introduction un bruit blanc) et la force de frottement -aXt.

Le code 6 permet de simulée une seule trajectoire de l'équation de Langevin Xt avec un pas Ät = 10-⁴ (peut être fixé par l'utilisateur),et les deux paramètres a = 2 et D = 1, qui est présentée par la figure 3.4. Et le code 7 permet de simulée l'équation de Langevin en deux dimensions avec les conditions initiales x0 = 5 et y0 = 6, comme c'est illustré par la figure 3.5. C'est-à-dire que on a système d'équation différentielle stochastique de la forme suivante :

( v

dXt = -aXtdt + 2DdW¹

t , X0 = x0

v

dYt = -aYtdt + 2DdW²

t , Y0 = y0

Avec W1

t et W2

t deux mouvements browniens standards indépendants, a = 2 et D = 1.



FIGURE 3.4 - Trajectoire simulée de l'équation de Langevin avec a = 2 et D = 1.

FIGURE 3.5 - L'équation de Langevin en deux dimensions avec a = 2 et D = 1.

3.3.4 Bruit blanc, bruit coloré

L'interprétation en sciences physiques des équations stochastiques repose sur la notion de bruit blanc ou de bruit coloré. Un processus de Markov continu a des increments X(t +dt)-X(t) qui ne dépendant que du temps t, dt et de X(t) et de manière différentiable pour X(t) = x de t et de x. La continuité de X(t) traduit le fait que l'increment X(t + dt) - X(t) sachant que X(t) = x tend vers 0 quand dt tend vers 0. Cet increment suit une équation de la forme

v

X(t + dt) - X(t) = u(t,x)dt + ó(t,x)N(0,1) dt (3.18)

où u(t,x) et ó(t,x) sont des fonctions arbitraires et N(0,1) est une variable aléatoire gaussienne centrée réduite. En remplaçant x par Xt, on a

X(t +dt) -X(t)	ó(t,Xt)N(0,1) = u(t,Xt)dt + vdt
dt

Le processus aléatoire N(0, 1)/vdt suit une loi gaussienne N(0,1/dt). La limite de ce processus définit le bruit blanc gaussien

î(t) = lim dt?0

1 N ( 0, dt)

Ce bruit conduit à l'équation de Langevin (3.17)

dXt
dt

= u(t,Xt)dt + ó(t,Xt)î(t)

Il vérifie les propriétés

(î(t)) = E(î(t)) = 0

(î(s),î(t)) = R(s,t) = ä(t - s)

La fonction de corrélation ne dépend que de la différence ô = t - s. La densité spectrale du

processus vaut

+8 +8

S(ù) = 2J R(ô)e^-iùôdô = 2/ -8 ä(ô)e^-iùôdô = 2

Le spectre du bruit blanc est constant et indépendant des fréquences.

Un bruit coloré est un bruit dont la densité spectrale dépend des fréquences. La couleur est la sélectivité ¹ du spectre de fréquences. Donnons deux exemples de bruit coloré

(1) Le processus d'Ornstein-Uhlenbeck stable (l'exemple 3.5) est un bruit coloré. Il est obtenu à partir du processus d'Ornstein-Uhlenbeck Xt dépendant d'un instant t0 et en faisant tendre le temps initial vers -8

åt = lim t?-8

Xt = N ( 0, 2) 2r

Ce processus est indépendant de la condition initiale x0 et de t. Il vérifie les propriétés

(å(t)) = E(å(t)) = ⁰

2

Ws), å(t)) = R(s,t) = ó

_2r-r|t-s|

Lorsque le paramètre r tend vers l'infini, ce bruit coloré tend vers un bruit blanc. Sa densité spectrale est

S(ù) =

ó2 +^82a^2 | e -iùô _7 2a2 uô = r _ ù2 + r2 . I 8

1. Télécommunications capacité d'un récepteur à capter les signaux d'un émetteur tout en éliminant les signaux émis par des fréquences voisines.

óx (Xt)) dt + dWt (3.21)

Preuve Appliquant la formule d'Itô (3.7) à la fonction Yt = f (t,x) = fj ól:), avec

(2) Un autre exemple de bruit coloré est le bruit des télégraphistes. Le processus stochastique est un processus Xt = #177;a, qui ne prend que deux valeurs. L'intervalle de temps entre les changements de signe est distribue exponentiellement. Si on note u le nombre moyen de changement de signe par unité de temps, le processus a pour fonction de corrélation

(X(s),X(t)) = R(s,t) = _a2_e-2u|t-s|

Sa densité spectrale vaut

8u

S(ù) = a2 ù2 + 4u²

3.3.5 Transformée de Lamperti

Avant d'aborder les méthodes de résolution numérique, nous allons introduire la transformée de Lamperti [30] qui est très utile pour traiter de nombreux processus stochastiques markoviens. Sous la forme différentielle suivante

dXt = u(t,Xt)dt + ó(Xt)dWt (3.19)

On appelle la transformée de Lamperti de X, la fonction Y définie de la manière suivante :

Xt du

Yt = F(Xt) = fo ó(u) (3.20)

Avec ce changement de variable, l'équation stochastique (3.19) devient

dYt = uY (t,Yt)dt + dWt

avec

u (t,F-^1(y)) uY (t, Yt) = ó (F-^1(y))

1 ? 2 ?x ó (F^-1 (y))



d'où, on a

dYt =

( 1

p(t

,

X_t)
ó(Xt)₂

?f(t,x)
?t

?f(t,x) 1 ?^2 f(t,x) óx(x) = 0 , = et = ?x ó(x) ?x^2 -ó^2(x)

on obtient directement le résultat

dYt = df(t,x) = (0 +u(t,Xt)

1 1 2 óx PO 1 ó (Xt) 2 ) dt + ó(Xt) dWt ó (Xt) 2 ó (Xt) óW(t)

(u(t,Xt) 1

=

ó(Xt) ₂óx(Xt)) dt + dWt

L'intérêt de cette transformation est qu'elle supprime le bruit multiplicatifs de l'équation différentielle stochastique initiale au profit d'une équation de Langevin non linéaire avec un bruit blanc simple. Nous allons voir par la suite que sur le plan numérique ce changement de variable est très utile et peut améliorer la précision numérique. Cette transformation a été l'objet de peu de travaux jusqu'à un passé récent, elle est étudié en détail ses propriétés et son intérêt dans [19]. En fait, la transformation de Lamperti peut être envisagée d'au moins deux points de vue. D'une part, elle permet une manipulation indirecte des processus auto-similaires (filtrage, prédiction, etc.) en opérant sur leurs contreparties stationnaires. D'autre part, des outils stationnaires classiques peuvent être "Lampertisés", offrant ainsi de nouvelles façons de manipuler les processus autosimilaires².

3.4 Schémas numériques

De manière similaire aux équations différentielles ordinaires où la résolution numérique passe par une discrétisation du temps et un schéma d'approximation concernant l`intervalle de temps élémentaire sur lequel l'intégration est faite, il est nécessaire de procéder de manière similaire avec les équations différentielles stochastiques, à quelques différences prés :

(i) Pour une équation ordinaire, la trajectoire étant déterministe (en tout cas pour les exemples simples à une particule, oscillateur harmonique ou anharmonique), on peut contrôler avec la solution exacte la qualité de l'approximation. Avec un processus de Wiener, nous avons vu que deux trajectoires sont très différentes (Figure 2.2), cette différence s'accroissant avec le temps. Si par contre, on cherche à résoudre un processus d'Ornstein-Uhlenbeck (Figure 3.3), on sait que le système évolue vers une distribution stationnaire que la variance des trajectoires est finie.

(ii) Les schémas d'approximation des méthodes de résolution des équations différentielles sont basés sur le calcul différentiel usuel. Dans le cas des équations différentielles stochastiques, ces schémas reposent sur le calcul différentiel stochastique de nature assez différente.

(iii) Des questions fondamentales sur l'existence et l'unicité d'une solution pour une équation existent de manière similaire aux équations ordinaires. Sans bien évidemment rentrer dans le détail qui relève de travaux très techniques, il y a deux critères pour prouver l'existence et l'unicité d'une solution sur un intervalle donné (Section 3.3.2)

Dans beaucoup de literatures peut être vue dans les références, les schémas numériques pour la résolution numérique de l'équation (3.22) ont été proposés dans [22, 29, 31, 37, 40].

Soit à nouveau la forme différentielle de l'équation différentielle stochastique

dXt = f(Xt)dt +g(Xt)dWt, X0 = x0, (3.22)

avec Wt est un processus de Wiener et x0 la valeur initiale. On utilise des schémas aux diffé-
rences classiques et le fait que pour un pas h donné, les variables W(n+1)h ^-Wnh suivent des lois

2. On appelle processus auto-similaires les processus dont la loi est invariante par changement d'échelle des temps.

gaussiennes indépendantes de variance h. On note xt le processus approché et on considère la subdivision

0 = t0 < t1 < ··· < tN-1 < tN = T

(T -t0)

de pas régulier

h = Ät = tn+1 -t_n =

N

oil

t_n = nÄt, n ? {1,2,...,N}.

Posant la notation suivante X_n = Xt_n, les trois variables aléatoires suivantes seront employées dans les schémas



ÄW_n = Wtn+1 -Wt_n,

tn+1 f_ns

nÄZ_n = dW

Ils sont obtenus comme des variables aléatoires normales, utilisant la transformation suivante ÄWn = în,1(Ät)1/2,

1

n1 + ^în,2) (Ät)3/2,

ÄZ_n = 2 ( , _v3

( ÄZ_n = 21 n1 -

.\n/7;) (Ät)3/2.

Ainsi qu'eux nous employons plus loin Ä l7Vn ¹/², avec P

= în,2 (Ät) -,n,1, în,2 deux variables indé-

pendantes de loi N(0,1).

Schéma d'Euler

Xn+1 = Xn + f_nÄt + gnÄWn (3.23)

Schéma de Milstein

Xn+1 = X_n + f_nÄt + _gnÄWn + ₂g⁰

1 _ng_n((ÄW_n)² - Ät) (3.24)

sds,

ÄZn = tn+1 fs

dsdWs.

tn n

Schéma de Milstein Seconde

Xn+1 = X_n + f_nÄt + gnÄWn + ₂g⁰

1 _ng_n((ÄW_n)² - Ät)+ f⁰_ngnÄZn

~ ~

+ g0 n f_n + ¹ ₂g⁰⁰ _ng² ÄZ_n + ¹ 6(g02 _n g_n + g⁰⁰ _ng² _n)((ÄW_n)³ - 3ÄtÄW_n)

n

(3.25)

Schéma de Itô-Taylor D'ordre 1.5 Xn+1 = Xn + fnÄt + gnÄWn + 2g^0 1 ngn((ÄWn)^2 - Ät) + f^0ngnÄZn + (gn^t fn + 21 gngn) ÄZn +^61 (gn^t2gn + gn^00gn^2)((ÄWn)^3 - 3ÄtÄWn) ~ ~ 1 + f n 0 fn + 1 2 f n 00 g2 (Ät)2 n 2 Schéma de Heun D'ordre 2

(3.26)

1 Xn+1 = X_n + ₂ [F1 + F2]Ät + ¹² [G1 + G2]ÄW_n (3.27)

Avec

F1 = F(Xn), G1 = g(X_n),

F2 = F(X_n +F1Ät + G1ÄW_n), G2 = g(Xn +F1Ät + G1ÄW_n),

~ ~

f - 1

F_x = ₂g⁰g (x).

Schéma de Runge-Kutta D'ordre 3

Xn+1 = Xn +

~ ~ 1 1 1 4[F1 + 3F3]Ät + 4[G1 + 3G3]ÄWn + f 0 2v3 ngn - g0 n fn - 1 2g^00 ng^2 ÄtÄ ^eWn (3.28) n

Avec

F1 = F(X_n), G1 = g(X_n),

~ ~ ( )

X_n + 1

F2 = F ₃F1Ät + 1 X_n + 1

₃G1ÄWn , G2 = g ₃F1Ät + 1 ₃G1ÄWn ,

( ~ ( ~

X_n + 2

F3 = F ₃F2Ät + 2 X_n + 2

₃G2ÄWn , G3 = g ₃F2Ät + 2 ₃G2ÄWn ,

[ ]

f -- 1

F_x = ₂g'g(x).

3.4.1 Simulation numérique

L'intérêt pratique de la simulation d'équations différentielles stochastiques est très important, car ce n'est pas toujours facile à déterminer la solution analytiquement. Cela rend difficile l'étude de l'évolution dynamique d'un phénomène, où par exemple l'analyse statistique de la variable aléatoire l'instant de premier passage (IPP) correspondant à la solution de l'équation, qui sera illustré dans le chapitre suivant. Aujourd'hui, le développement de l'outil informatique motive les scientifiques pour mettre au point des schémas numériques pour la résolution approchée des ÉDS.

L'approcher est simple, simuler numériquement Xt jusqu'au temps T, puis approcher l'espérance E(Xt) à l'aide de la loi des grands nombres, c'est-à-dire la moyenne de M trajectoires indépendantes de Xt. Cette dernière technique présente un certain intérêt, surtout lorsque la dimension de Xt est grande. En effet, les méthodes où les schémas numériques deviennent dans ce cas vite lourdes à mettre en oeuvre. La fonction snssde ³ permettre de simulée numériquement la solution approchée des ÉDS par les schémas numériques qui sont en détail dans la section précédente.

R> help("snssde")

R> example("snssde")

R> snssde(N, M, T = 1, t0, x0, Dt, drift, diffusion, Output = FALSE, + Methods = c("SchEuler", "SchMilstein", "SchMilsteinS",

+ "SchTaylor", "SchHeun", "SchRK3"), ...)

3. simulation numerical solution of stochastic differential equations.

Détails:

N La taille de processus.

M Le nombre de trajectoire à simulée.

T L'instant final (par défaut égale à 1).

t0 L'instant initial.

x0 La valeur initial.

Dt La discrétisation où le pas (Si Dt est fixée alors par défaut

T = t0 + Dt * N).

drift⁴ Coefficient de dérive (une expression qui dépende de x et de t).

diffusion⁵ Coefficient de diffusion (une expression qui dépende de x et de t). Output Pour sauvegardée les résultats de simulation sous forme Excel
(par défaut c'est FALSE).

Methods Différents méthodes de simulation (par défaut schéma d'Euler),

avec SchEuler(3.23), SchMilstein(3.24),SchMilsteinS(3.25), SchTaylor(3.26), SchHeun(3.27), SchRK3(3.28).

Exemples illustratifs :

Radial Ornstein-Uhlenbeck (ROU)

Ce modèle est défini par l'équation différentielle stochastique suivante :
( è )

dXt = - Xt dt + ódWt, x0 > 0.

Xt

avec è,ó > 0, en simuler numériquement par la méthode d'Euler (3.23), une seule trajectoire de la solution approché ce modèle qui est présentée dans la figure 3.6. En suite un flux de M = 100 trajectoires indépendantes avec leurs moyenne qui est illustré par la figure 3.7. Posant par exemple è = 1 et ó = 0.25, pour un pas Ät = 0.01 et la valeur initiale x0 = 3.

R> f <- expression( (1/x) - x )

R> g <- expression( 0.25 )

R> snssde(N = 1000, M = 1, t0 = 0, x0 = 3, Dt = 0.01, drift = f,

+ diffusion = g)

R> snssde(N = 1000, M = 100, t0 = 0, x0 = 3, Dt = 0.01, drift = f,

+ diffusion = g)

4. u(t,Xt) peut être constant où une fonction qui dépende de x ou t.

5. ó(t,Xt) peut être constant où une fonction qui dépende de x ou t.



FIGURE 3.6 - Simulation une seule trajectoire du modèle Radial Ornstein-Uhlenbeck par le schéma d'Euler.

FIGURE 3.7 - Simulation un flux de 100 trajectoires du modèle Radial Ornstein-Uhlenbeck par le schéma d'Euler.

Coefficient de dérive en fonction de x et t

On consider l'équation différentielle stochastique suivante :

dXt = (átXt -X3 _t )dt +ódWt, x0 = 0.

Avec á,ó > 0, cette équation n'admet pas une solution de forme explicite. En utilise le schéma de Milstein (3.24), pour simulée numériquement une solution approché de Xt. En remarque que cette équation différentielle est spéciale car les trajectoires simulées de ce processus sont comprises entre #177;vát qui est illustré dans la figure 3.8, c'est-à-dire que la variance de Xt dépende de temps t. La simulation d'un flux de 100 trajectoires de Xt montre que son espérance tend vers 0 quand t tend vers l'infini, qui est donné par la figure 3.9. Posant á = 0.03 et ó = 0.2, pour un pas Ät = 0.1 et la valeur initiale x0 = 0.

R> f <- expression( (0.03 * t * x - x^{^}3)) R> g <- expression( 0.2 )

R> snssde(N = 1000, M = 1, t0 = 0, x0 = 0, Dt = 0.1, drift = f,

+ diffusion = g,Methods="SchMilstein")

R> snssde(N = 1000, M = 100, t0 = 0, x0 = 0, Dt = 0.1, drift = f,

+ diffusion = g,Methods="SchMilstein")

R> G <- function(t) sqrt(0.03 * t)

R> t <- seq(0, 100, by = 0.1)

R> points(t, +G(t), type="l", col="blue", lwd=2)
R> points(t, -G(t), type="l", col="blue", lwd=2)



FIGURE 3.8 - Simulation une seule trajectoire du modèle dXt = (0.03tXt -X3 _t )dt +0.2dWt par le schéma de Milstein.

FIGURE 3.9 - Simulation un flux de 100 trajectoires du modèle dXt = (0.03tXt - X3 t )dt + 0.2dWt par le schéma de Milstein.

Coefficients de dérive et de diffusion en fonction de t

Soit l'équation différentielle stochastique suivante :

dXt = cos(t)dt + sin(t)dWt, x0 = 0.

Cette équation n'admet pas une solution de forme explicite. En utilise le schéma de Itô-Taylor D'ordre 1.5 (3.26), pour simulée numériquement une solution approché de Xt. D'après la figure 3.10 en remarque que l'espérance de Xt égale à x0 + sin(t).

R> f <- expression( cos(t) ) R> g <- expression( sin(t) ) R> snssde(N = 1000, M = 100, t0 = 0, x0 = 0, Dt = 0.1, drift = f,

+ diffusion = g,Methods="SchTaylor")

R> G <- function(t) sin(t)

R> t <- seq(0, 100, by = 0.1)

R> points(t, G(t), type="l", col="blue", lwd=2)



FIGURE 3.10 - Simulation un flux de 100 trajectoires du modèle dXt = cos(t)dt + sin(t)dWt par le schéma de Itô-Taylor.

Remarque 3.4 (Temps d'exécution) Le langage R n'étant pas un langage compilé mais interprété ligne à ligne, les instructions de ce type sont exécutées lentement et on préférera, lorsque c'est possible, utiliser les outils de calcul sur les vecteurs et les matrices (Optimisés sous R). La fonction system.time permet de mesurer le temps d'exécution d'une fonction et de comparer la vitesse d'exécution de plusieurs programmes. On peut gagner un temps considérable (10 à 100 fois plus rapide) en repérant les fonctions qui prennent du temps et en les programmant dans un langage compilé (C en particulier). L'avantage de langage R c'est les calculs formelles (pas besoin de donnée les dérivées de coefficient de dérive et de diffusion dans les schémas numérique).

3.4.2 Relation entre le schéma d'Euler et Milstein

L'approximation de schéma de Milstein (3.24) a un ordre fort de convergence égal à 1. Cette méthode améliore donc les instabilité numériques par rapport à la méthode d'Euler (3.23). Toutefois, il y a un lien entre les deux méthodes dans le cas où on peut réaliser une transformée de Lamperti (3.21) de l'équation différentielle stochastique de départ. En effet, dans le cas où l'équation stochastique de départ n'a pas de bruit multiplicatif, comme par exemple dans l'exemple du processus d'Ornstein-Uhlenbeck, la méthode d'Euler a un ordre de convergence fort qui devient égal à 1. Or, avec une transformée de Lamperti (si g(Xt) est indépendante du temps), on peut transformer l'équation stochastique en une autre équation sans bruit multiplicatif. Ainsi, on peut montrer que le schéma d'Euler de l'équation transformée est identique au schéma de Milstein

sur l`équation originale. Dans le cas où la transformée de Lamperti est difficile à obtenir analytiquement, il est utile d'utiliser le schéma de Milstein qui est plus précis.

Preuve Soit le schéma de Milstein de l'équation différentielle stochastique (3.22)

1

Xn+1 = X_n + f_nÄt +g_nÄW_n +

₂gn⁰gn((ÄWn)² - Ät)

Soit la fonction y = F(x) et l'inverse x = G(y), appliquons la transformée de Lamperti (3.21) à Yt = F(Xt), on trouver

0

dYt = (fn gn - 1 2gn) dt +dWt

En suite appliquons le schéma d'Euler (3.23) a l'équation Yt

~ fn ~

- 1

ÄY = Yn+1 -Y_n = ₂g⁰ Ät + ÄW_n

n

g_n

Et en appliqué le développement de Taylor à l'inverse de la transformation de Lamperti, i.e. à Xt = G(Yt), on trouver

G(Y_n +ÄY) = G(Y_n)+ G⁰(Y_n)ÄY + 21 Of (Y_n) (Ä11² + O (ÄY³)

Notons

d

G⁰(Y_n) = F^-1(y) = 1

F0(G(y)) = g_n

dy

et

Gⁿ(Yn) = G⁰(Y_n)g'_n = gngn^'

Donc finalement, on a

G(Y_n +ÄY) - G(Y_n) = g_nÄY + 21 gngn^' (ÄY)² + O(ÄY³)

~~ fn ~ ~ ~~ fn ~ _~2

- 1 - 1

= g_{n 2}g⁰ + 1

Ät + ÄW_{n 2}gng⁰ ₂g⁰ Ät + ÄW_n

n n n

g_n g_n

+ O(ÄY³)

~ ~~ fn ~ ~~~~ fn ~ ~

g_n + 1

= ₂gng⁰ - 1 ₂g⁰ - 1

Ät + ÄW_{n 2}g⁰ Ät + ÄW_n

n n n

g_n g_n

+ O(ÄY³) ~ ~ ~ ~

g_n + 1

= ₂gng⁰ f_n - 1

_nÄWn ÄW_n + ₂gng⁰ Ät + O(Ät³/²)

n

= f_nÄt + gnÄWn +

₂g_n⁰gn((ÄWn)² - Ät)+ O(Ät³/²)

Puisque on a ÄW_n = î_n,1(Ät)¹/² et (ÄW_n)² = î2 n,1Ät, d'où on trouver le schema de Milstein

G(Y_n +ÄY) -G(Y_n) = Xn+1 -X_n = f_nÄt + _gnÄWn + 21 g_n⁰ g_n((ÄW_n)² - Ät) + O(Ät³/²)

Exemple 3.6 (Transformation de modèle Cox-Ingersoll-Ross (CIR)) Le modèle Cox-IngersollRoss (CIR) est utilisé en mathématiques financières pour modéliser l'évolution des taux d'intérêt court terme. Il s'agit de la solution de l'équation différentielle stochastique

dXt = (è1 - è2Xt)dt + è3vXtdWt, X0 = x0 > 0. (3.29)

oil sous cette forme

dXt = ê(è - Xt)dt + óvXtdWt, X0 = x0 > 0. (3.30)

Notons que la solution de cette de l'équation (3.30) reste strictement positive sous la condition 2êè > ó². Le paramètre è donne la moyenne à long terme, et ê > 0 donne la vitesse à laquelle le processus va converger vers cet équilibre.

Appliquons le schéma de Milstein (3.24) à l'équation (3.29), on a :

ÄX = ( (è1 - è2Xn) - 41 èi) Ät + è3 vX_nÄW_n + ⁴¹ è3 (ÄW_n)² (3.31)

Maintenant appliquons la formule d'Itô (3.7) à Yt = vXt pour l'équation (3.29), en posant y = f(t,x) = vx, on a ?t f = 0, ?_x f = 1/2vx et ?xx f = -1/4x³/². On obtient :

dYt =

2Yt , 1 4 2 1 ((è1 - è2Y`t^,) - è^,1) dt + 2 è3dWt -

(3.32)

Appliquons le schéma d'Euler (3.23) à Yt, on a

ÄY =	2Ynn 1 1 ( (è1 - è2Y^2) - 4 è32 ) Ät + 2 è3ÄWn

Et en appliquant le développement de Taylor à l'inverse de y = vx, c'est-à-dire à G(y) = x², on obtient :

G(Y_n +ÄY) - G(Y_n) = (Yn + ÄY)² -Y²_n

= (ÄY)² + 2Y_nÄY

= [ 2Yn ((è1 - è2Y 4

2) -- 1 è32 ) Ät + ²¹ è3ÄWn1 2

+ ( (è1 - è2Y2 ) - 41 èi) Ät + Ynè3ÄWn

= ( (è1 - è2Y2 ) - 14 è3²) Ät + _è3YnÄWn + ₄ è3 (ÄW_n)² + O (Ät²) On remplace Yt = vXt, on trouve la même résultat (3.31).

Exemple de simulation

Le code 8 nous permettre de comparée entre le processus origine Xt (3.29) avec son transformation Yt (3.32), qui est illustré par le graphie 3.11. Posons (è1,è2,è3) = (0.1,0.2,0.05), avec un pas Ät = 0.1 et x0 = 1. (2è1 > è2 3 est vérifie)



FIGURE 3.11 - Transformation de modèle Cox-Ingersoll-Ross (CIR) dXt = (0.1 - 0.2Xt)dt + 0.05vXtdWt .

3.5 Les modèles attractives

Le phénomène de dispersion est un problème complexe, souvent caractérisé par sa diffusion. Beaucoup de situations réelles se comportent naturellement a ce phénomène de dispersion, notamment en environnement, biologie, chimie, etc... La modélisation mathématique ou par simulation du comportement dynamique de tels phénomènes est très difficile, et nécessite souvent un outil puissant, pour décrire le phénomène observé.

Différentes approches, déterministes et non déterministes ont été proposées par plusieurs auteurs, pour décrire l'évolution des phénomène de dispersion, voir par exemple [23, 26, 24, 2, 4, 5]. Les modèles mathématiques de diffusion dont l'origine est un phénomène biologique (voir le cas particulier du mouvement Brownien), peuvent être un outil de modélisation intéressent pour beaucoup de phénomènes de dispersion. Dans un premier cas, on propose un modèle de dispersion d'un polluant qui se déplace sur une surface d'eau turbulente, en présence d'un mécanisme d'attraction (une cible), qui oriente le polluant vers un point de localisation. Le deuxième modèle décrit le comportement d'une population d'insectes, où deux insecte mâle-femelle s'attirent l'un vers l'autre. On fait les hypothèses suivantes :

H1 Le mouvement des insectes (les polluantes) est Markovien.

H2 Le domaine D est compact.

H3 Le mécanisme d'attraction est suffisamment fort pour attirer les insectes (les polluantes).

H4 Les insectes (les polluantes) se déplacent, indépendamment les uns des autres.

H5 Les forces extérieures opposées déplacement de l'insecte (les polluantes) sont négligeables.

H6 Les taux d'immigration et de mortalité sont négligeables.

Sous ces hypothèses, la famille de modèles que nous proposons est une idéalisation mathématique d'un système réel complexe, c'est aussi une sorte de simulation par des processus de diffusion du phénomène observé.

3.5.1 Modèle d'une diffusion en attraction 914ós=1(Vt)

La pollution est un problème qui menace aujourd'hui la nature des humains et des espèces. En particulier, la pollution de l'eau est un sujet important, auquel beaucoup de scientifiques s'intéressent actuellement. Généralement on utilise les équations aux dérivées partielles déterministes, pour décrire la densité du polluant. Les processus de diffusion sont utilisés comme modèles mathématiques pour décrire le comportement dynamique des particules polluantes [24]. Le modèle que nous proposons ici décrit le comportement dynamique d'une particule en mouvement sur une surface d'eau turbulente [2, 4].

On suppose qu'on observe à l'instant t, la position Vt = (Xt,Yt) ? D ? 1R² d'un particule. Le point V0 = (X0,Y0) représente une position initiale du polluant observé. La nature du phénomène observé implique que pour tout t ? [0,T], Vt peut être considéré comme processus de diffusion. On désigne par L(x,y,t) la profondeur de la surface d'eau au point (x,y) et à l'instant t, U_w(x,y,t) et V_w(x,y,t) sont respectivement les champs de vitesses, dus au mouvement d'eau, suivant les directions x et y. U_a(x,y,t) et V_a(x,y,t) sont respectivement les champs de vitesses, dus au mécanisme d'attraction, suivant les directions x et y. Le coefficient de dispersion est notée par D(x,y,t). À chaque instant t la position de la particule est supposée être un processus de diffusion Vt = (Xt,Yt), solution de système d'équation différentielle stochastique suivante [4]

dVt =	? ? ?	dXt = (-Ua +Uw + ??Lx Dx ) dt+v2DdWt1 , t ? [0,T] (3.33) dYt = - Va +Vw + ay L D +^?D) dt +v2DdWt^2 a y
Avec Ua =			Kx V= (V x2 + y2 s+1 a	Ky ~p x2 + y2s+1 , D(x,y,t) = 2ó^2 1

ó et K sont des constantes positives, s = 1 et 2K > ó², W_t¹ et W²_t deux mouvements browniens standards indépendants. U_w(x,y,t) et V_w(x,y,t) sont supposés négligeable et la profondeur L(x,y,t) est constante.

Pour chaque particule, on suppose qu'il y a conservation de la masse, ce qui est traduit par l'équation suivante [24]

?L ?t +	?UL + ?x	?VL		= 0 (3.34)
		?y

Avec

U(x,y,t) = U_w(x,y,t)+U_a(x,y,t) et V (x,y,t) = V_w(x,y,t)+V_a(x,y,t)

D' où le système (3.33) devient

? ???

???

dVt =

dXt =

V -KX iq+Yt2Y+1dt + ódW_t¹

, t ? [0,T] (3.35)

dYt = dt +ódW_t²

(0Q+Yt²y⁺¹

Appliquons la formule d'Itô multidimensionnel (3.10) au processus Rt = ||(Xt,Yt)||, avec Rt est la distance depuis l' attraction (la cible) a l' instant t. Posant f(x,y,t) = Vx² +y², on a ?t f = 0

?_x f = x

?xx f = _y2 , ?_yf = y ? = x2

Vx2 + _y2 (x² + _y2)3/2 Vx2 + y2 yy f (_x2 + _y2)3/2

d'où, la formule d'Itô entraîne

_~q ~ !

+ D X2

d X2 t +Y2 t

-Ua Xt -Va Yt t +Y2

= _{p p} dt

t X2 t +Y2 X2 t +Y² (X2 _t +Y2

t )3/2

t t

v2D ^Xt

+

Yt

VXt 2 Y2 dWt ¹ + v2D Xt ² +Yt2dWt2

+ D

(VXt ² +Yt 2 )s+2 Xt ² +Yt2)s+2

?
? ?

=

KY²

t

KX_t²

(X²_t +Y2t )3/2

X²_t +Y²_t

?
? ?

dt

v 2D

~XtdW1 ~

+ p t +YtdW2 t

X2 t +Y2

t

=	? ? ?	K~X2 ~ t +Y2 t (VXt ^2+ Yt2)s+2		D	? v ?2D dt + (XtdWt ^1 + YtdWt^2) VXt 2+Yt2
			p X2 t +Y2t

Posant le changement de variable en coordonnées polaire suivant

{

Xt = Rt cos(èt) Yt = Rt sin(èt) Ce qui impliquera

~ K D

dRt = - R.; + _Rt) dt + v2D (cos(èt)dWt ¹ + sin(èt)dW_t²)

Proposition 3.1 Soit Wt ¹ et W²_t deux mouvements browniens standards indépendants, alors on a
d ^-Wt = cos (èt)dWt ¹ + sin(èt)dW²_t

est un mouvement brownien standard.

Preuve d ^-Wt est un mouvement brownien standard, vérifiant :

E ( d Wt) = cos(èt)E(dWt ¹) + sin(èt)E(dWt ²) = 0

et

(- 2

E ( dWt) = cos²(èt)dt + sin²(èt)dt = (cos²(èt) + sin²(èt)) dt = dt

Donc d'après cette proposition le processus Rt est

KD

dRt = ( - + dt + v2Dd liit

R^s_t Rt

₍₇2

= ( - K + - ) dt + ód ^-Wt R^s_t 2Rt

On obtient, après simplification, le processus de diffusion Rt solution de l'équation différentielle stochastique suivante

{

ó2 Rs-1-v dRt = ( 2 "t s " Rt dt + ód Wt^- , R0 = a

a > 0, t ? [0,T] (3.36)

On s'intéresse particulièrement a la simulations aux deux modèles _Mós=1(Vt) ^et _Mós>1(Vt) en 2- et 3- dimensions.

Simulation du modèle Mós=1(Vt) L'équation de ce modèle, est de la forme

dXt = -KXt

t +Y2 t dt + ódW¹

_X2 t

dVt = , t ? [0,T] (3.37)

dYt = -KYt

X2 t +Y2 t dt + ódW²

t

La fonction RadialP2D_11 ⁶ permet de simuléee le phénomènee du polluant qui se déplacee sur une surface d'eauu turbulente en 2-Dimensionss (figure 3.12). Et la fonction RadialP3D_11 simuler le phénomènee en 3-Dimensionss (figure 3.13). att1in3DD

6. Le constant v est un rayon d'unee sphèree oùn un cercle qui nous permettre de donnéee lapremièree instant d'attendree la cible par une polluant.

0.001,	T =	1,	X0	=	2,	Y0	=	1,
0.2)
0.001,	T =	1,	X0	=	1,	Y0	=	1.25,

R> RadialP2D_1(N = 10000, t0 = 0, Dt = + v = 0.3, K = 2, Sigma = R> RadialP3D_1(N = 10000, t0 = 0, Dt = + Z0 = 0.5, v = 0.2, K = 2, Sigma = 0.2)



FIGURE 3.12 - Trajectoire du polluant dans une surface d'eau turbulente en 2--D avec s = 1.

FIGURE 3.13 - Trajectoire du polluant dans une surface d'eau turbulente en 3--D avec s = 1.

Simulation du modèleMó s>1(Vt)

L'équation de ce modèle est donnée par le système d'équation (3.35). La fonction RadialP2D_2 permet de simulée le phénomène du polluant qui se déplace sur une surface d'eau turbulente en 2--Dimensions (figure 3.14). Et la fonction RadialP3D_2 simuler le phénomène en 3--Dimensions (figure 3.15).

1,	X0	=	2,	Y0	=	1,
1,	X0	=	1,	Y0	=	1.25,

R> RadialP2D_2(N = 10000, t0 = 0, Dt = 0.001, T = + v = 0.3, K = 2, s = 2, Sigma = 0.2) R> RadialP3D_2(N = 10000, t0 = 0, Dt = 0.001, T = + Z0 = 0.5, v = 0.2, K = 2, s = 2, Sigma=0.2)



FIGURE 3.14 - Ttrajectoire du polluant dans
une surface d'eau turbulente en 2-D avec s > 1.

FIGURE 3.15 - Ttrajectoire du polluant dans
une surface d'eau turbulente en 3-D avec s > 1.

3.5.2 Modèle de deux diffusion en attraction 91,4_m

_>0(0¹⁾) ?-91,4₀^ó (V⁽²⁾

t )

Dans le paragraphe suivante, nous allons proposer, un modèle de deux diffusion en attraction [5], l'une de dérive nulle 91q(V_t⁽²⁾) et l'autre 91,4'm'>0(Vt ⁽¹⁾), de dérive non nulle. Ce qui peut modéliser le comportement de deux insectes, l'un attire l'autre.

On considère V_t⁽¹⁾ = (Xt,1,Xt,2) et Vt⁽²⁾ = (Yt,1,Yt,2) deux processus aléatoire de diffusion, qu'on suppose représentant respectivement les positions d'un insecte mâle et d'un insecte femelle, et que le mâle est attiré par la femelle. Le comportement de la femelle est supposé être un processus de mouvement brownien, défini par l'équation suivante

dV_t⁽²⁾ = óI2×2dWt (3.38)

alors que le comportement du mâle est supposé être un processus de diffusion, dont la dérive est
orientée, à chaque instant t, vers la position de la femelle, et qui est donné par l'équation suivante

dVt (1) = dVt ⁽²⁾ + _um+1 (||Dt||)Dtdt _+óI2×2d ^eWt (3.39)

où

i Dt = Vt (1) - V_t⁽²⁾ tum (||d||) = - ||dKmt||

K et m sont des constantes positives, et K > ó².

Le modèle proposé est de la forme suivante

dVt (1) =

? ???? ????

K(Xt,1-Yt,1) dXt,1 = dYt,1 - ~v(Xt,1-Yt,1)2+(Xt,2-Yt,2)2~m+1 dt + ód K(Xt,2-Yt,2) dXt,2 = dYt,2 - ~v(Xt,1-Yt,1)2+(Xt,2-Yt,2)2~m+1 dt + ód

Cf^7t^1 eW^2 t

,t ? [0,T] (3.40)

(

dYt,1 = ódW¹

dV(2) t

t = ,t ? [0,T] (3.41)

dYt,2 = ódW2 t

avec eW1

t , eW2 t , W1

t et W2

t sont des mouvements Browniens indépendants.

Appliquons la formule d'Itô multidimensionnel (3.10) au processus Rt = ||(V_t ⁽¹⁾,V_t ⁽²⁾)||, avec Rt représente la distance entre les deux insectes a l'instant t. Posant f (x1,x2,y1,y2,t) = _p

(x1 -y1)² + (x2 -y2)², on a ?t f = 0

x1 - y1 (x2 - y2)²

?x1 f = p (x1 - y1)² + (x2 - y2)2 , ?x1x1 f = ((x1 -y1)² + (x2 -y2)²)^3/2

x2 -y2 (x1 -y1)²

?x2 f =p (x1 - y1)² + (x2 - y2)2 , ?x2x2 f = ((x1 -y1)² + (_x2 -y2)²)^3/2

-

?y1 f =

(x1 - y1) (x2 - y2)²
p (x1 - y1)² + (x2 - y2)2 , ?y1y1 f = ((x1 -y1)² + (x2 -y2)²)^3/2

- (x2 -y2) (x1 -y1)²
?y2 f = p (x1 - y1)² + (x2 - y2)2 , ?y2y2 f = ((x1 -y1)² + (x2 -y2)²)^3/2 d'oil, la formule d'Itô entraîne

~ ~

?t f + _A1?x1 f + _A2?x2 f + ó2

d f = 2 (?x1x1 f + ?x2x2 f + ?y1y1 f + ?y2y2 f ) dt

~ ~

+ ó (?x1 f + ?y1 f )dW¹

t + (?x2 f + ?y2 f )dW2 _t + ?x1 f d eW¹

t + ?x2 f d eW2 t

= A1	Xt,1 -Yt,1 +A2 \/(Xt,1 -Yt,1)^2+ (Xt,2 - Yt,2)^2		Xt,2 -Yt,2
		p (Xt,1-Yt,1)^2 + (Xt,2 -Yt,2)^2

+ _ó2 (Xt,1 -Yt,1)² + A-,2 -Yt,2)² Xt,1 -Yt,1 d W1

((Xt,1 - Yt,1)² + (Xt,2 - Yt,2)2)3/2)dt + ó -Yt,1)² + (Xt,2 - Yt,2)2

+ (Xt,1 -Yt,1)² + (Xt,2 -Yt,2)2 d eW²

Xt,2 -Yt,2 _pt

Avec
A1 = D'oil	-K (Xt,1 -Yt,1) et A2 = ( 1/ (Xt,1 - Yt 1)^2 + (Xt2 - Yt,2)2) m+1	-K (Xt,2 -Yt,2)
		(1/(Xt,1 -Yt,1)^2+(Xt,2 -Yt,2)^2)^m+1

? ?

df = ?

?

-K ((Xt,1 -Yt,1)² + (Xt,2 -Yt,2)²) ?

₊2 + ó2 ((Xt,1 -Yt 1)² + P(t 2 -Yt 2)²)^3/2

(Xt,1 -Yt,1)² +(Xt,2 -Yt,2)²

dt

(1/ (Xt ,1 - Yt _,1)2 + (Xt ,2 - Y _{t ,}2)²)^m ' -

ó ( (Xt 1 - Yt

V A,1 - Yt,1)² + (Xt,2 -Yt,2)²

,,1)d Wt ¹ + P(t,2 - Yt,2)d W_t²)

+ //

Posant le changement de variable en coordonnées polaire suivant

{

Xt,1 -Yt,1 = Rt cos(èt)

Xt,2 -Yt,2 = Rt sin(èt)

Ce qui impliquera que

df = d (\ I (Xt,1 - Yt,1)² + (Xt,2 -Yt,2)²) = dRt

Donc

-K ó2

dRt = ( Rtm + _Rt) dt + ó ( cos(èt)d iV+ sin(èt)d W^-_t²)

2

( -K + ^ó) =dt + ód ^eWt

Rm Rt

t

On obtient, après simplification, le processus de diffusion Rt solution de l'équation différentielle stochastique suivante

dRt =

(ó²R_t;_m¹-

{ K) dt + ód Wt,

t a > 0, t ? [0,T] (3.42)

R0 = a

Ce modèle permet de réaliser sur ordinateur une simulation dynamique du phénomène réel. Simulation le modèle Móm>0(V(1) _t ) ?-M 0 ó (V(2)

t ) en 2- et 3- Dimensions

La fonction TwoDiffAtra2D ⁷ permet de simulée le phénomène en 2-Dimensions, c'est-àdire simulé les deux systèmes d'équations différentielles stochastiques (3.40) et (3.41), qui est donnée par la figure 3.16. Et la fonction TwoDiffAtra3D simuler le phénomène en 3-Dimensions, la figure 3.17 donne une illustration de l'interaction entre deux insectes.

7. Le constant v est un seuil qui nous permettre de donnée l'instant de la première rencontre entre les deux insectes (voir les détails help("TwoDiffAtra2D")).

R> TwoDiffAtra3D(N = 5000,	t0 =
+		X3_0 =		2,	Y1_0 =
+		K = 3,		m =	0.2,

R> TwoDiffAtra2D(N = 2000, t0 = 0, Dt = 0.001, T = 1, X1_0 = 0.5, X2_0 = 1,

+ Y1_0 = -0.5, Y2_0 = -1,v =

0, Dt -0.5,

0.01, K

	= 2,	m = 0.25, Sigma = 0.1)
= 0.001,	T =	1,	X1_0= 1,	X2_0 = 1,
Y2_0 =	-1,	Y3_0	= 0.25,	v = 0.01,

Sigma = 0.2)



FIGURE 3.16 - L'interaction entre deux insectes en 2--D.

FIGURE 3.17 - L'interaction entre deux insectes en 3--D.

3.6 Conclusion

Dans ce chapitre, nous avons donner les définitions et les règles essentielles du calcul différentiel et intégral stochastique, et en particulier on a montrer comment on intègre une équation différentielle stochastique à coefficients lipschitziens. L'intégrale stochastique aborde définie globalement sur l'espace de toutes les trajectoires par convergence en L². La formule de Itô nous a permet de résoudre des ÉDS, ce qui implique son importance dans les calculs analytiques, et on a montre l'intérêt pratique de la simulation d'équations différentielles stochastiques qui est très important, car ce n'est pas toujours facile à déterminer la solution analytiquement, cela rend difficile l'étude de l'évolution dynamique d'un phénomène. Nous avons vues la puissance de la modélisation par les processus de diffusion en présence d'un mécanisme attractive, qui peut represent certain phénomène réelle telle que la trajectoire d'un polluant dans une surface d'eau turbulente.

4

Comportement Asymptotique Des Processus

de Diffusion

Sommaire
4.1	Introduction	85
4.2	Equation de Fokker-Planck	85
4.3	Processus de diffusion stationnaires	95
4.4	Classification des processus de diffusion linéaire	96
4.5	Calculs de l'instant de premier passage	108
4.6	Conclusion	115

4.1 Introduction

C

e chapitre est consacré à la présentation de deux équations fondamentales permettant de d'écrire l'évolution des lois de probabilités relatives à un processus de diffusion. Il s'agit

maintenant d'obtenir l'information maximale possible dans un cadre probabiliste, c'est-à-dire de déterminer les lois de probabilités elles-mêmes, pour décrire l'évolution temporelle d'un système hors d'équilibre obéissant à une dynamique markovienne.

L'équation de Fokker-Planck ou équation de Kolmogorov progressive est une équation aux dérivées partielles linéaire que doit satisfaire la densité de probabilité de transition p(s,x;t,y) d'un processus de Markov. A l'origine, une forme simplifiée de cette équation a permis d'étudier le mouvement brownien (Chapitre 2), équation qui est assimilable dans ce cas à l'équation de la chaleur. Comme la plupart des équations aux dérivées partielles, elle ne donne des solutions explicites que dans des cas bien particuliers portant à la fois sur la forme de l'équation, sur la forme du domaine où elle est étudiée. À chaque équation d'Itô (3.4) est associée une équation de Fokker-Planck qui décrit l'évolution dynamique de la densité du processus considéré, ainsi sa distribution stationnaire si il existe. Ces équations sont étudiés en détail ses propriétés, méthodes de solutions et simulation dans [18, 21, 25, 35].

À titre d'illustration, ce chapitre se termine par l'étude et l'analyse statistique de la variable aléatoire l'instant de premier passage "IPP" ("first passage time" en anglais). Dans le modèle d'une diffusion en attraction M ó s=1(Vt) (Chapitre 3 section 3.5.1), le taux ô(s)

c du polluant qui passe la frontière d'un voisinage d'une cible est analytiquement étudié pour le cas s = 1, le deuxième cas de s > 1 est étudier par simulation [2, 4, 6, 7]. Dans le deuxième modèle de deux diffusion en attraction (Chapitre 3 section 3.5.2), une étude par simulation est effectuée pour la

( )

_V(1)

densité de probabilité de l'instant de la première rencontre ô _t ,V⁽²⁾ entre deux insectes.

t

4.2 Équation de Fokker-Planck

À chaque équation d'Itô

dXt = u(t,Xt)dt +ó(t,Xt)dWt (4.1)

correspond une équation aux dérivées partielles que vérifie la probabilité de transition p(s,x;t,y) = p. Cette équation est appelée équation de Fokker-Planck ou équation de Kolmogorov progressive, qui est sous cette forme

?p
?t

1?2 (ó²(t,y) · p - ?

= ?y (u(t,y) · p), (s,x) fixé (4.2)

2 ?y²

Il existe une autre équation vérifiée par la densité de probabilité p qui est appelée équation de Kolmogorov rétrograde. Elle porte sur la variable x de départ

?p
?s

= 2ó2(s,x) ?2

1

(4.3)

?x2 p + u(s,x) ? ?x p, (t,y) fixé

L'équation de Fokker-Planck suppose, comme nous allons le voir, que le processus est markovien et que le bruit est un processus gaussien. Elle n'est donc pas valable pour tous les types de bruits.

Sous forme multidimensionnel, l'équation de Fokker-Planck ne change pas d'allure, elle est donnée par:

?p ?t

= L * p (4.4)

où L* est un opérateur connu sous le nom d'opérateur infinitésimal ou opérateur de Dynkin et prend la forme suivante, en dimension n :

L = ? i

ui(x) ? + 1 2 ? ói j(x) ?2 ?xi ?xi?xj i j

Prenons un exemple simple. Supposons un processus de Wiener standard avec un coefficient de diffusion constant égale à ó, c'est-à-dire :

dXt = ódWt

L'équation de Fokker-Planck de ce processus est :

?p ?t (s,x;t,y) = 2ó2 ?2

1 ?x2 p(s,x;t,y) (4.5)

On peut vérifiée facilement que la solution de cette dernière est une densité gaussienne :

1

p(s,x;t,y) =

( )

_-1 (x - y)²

_J exp , t > s et x,y ? R.

ó 2ð(t - s) 2 ó²(t - s)

Ce qui est donnée par l'équation (2.5) (Chapitre 2 section 2.3). Ce résultat est souvent cité dans les livres récents en finance quantitative. Pour ce qui concerne le cas d'un mouvement brownien géométrique (Chapitre 2 section 2.8) la solution de l'équation de Fokker-Planck est la densité lognormale p(t,S;t',S'). En effet, supposons le mouvement brownien géométrique pour l'action S :

dSt
St

= èdt +ódWt, S0 > 0

alors l'équation de Fokker-Planck est donnée par:

?p
?t'

1 ?2

= ?S'2 (ó²S⁰²p) - ? ?S0 (èS⁰p), (t,S) fixé (4.6)

2

où p définit la densité de probabilité de transition d'un état à un autre, t le temps présent et t' le temps futur, S le prix de l'action au temps présent t et S' le prix de l'action à une période future. La solution de l'équation (4.6) est donnée par:

" [log(S/S⁰) + (è - (1/2)ó²)(t -t)]2 #

1

p(t,S;t ,S⁰) = _óS'J exp -

2ð(t' -t) 2ó²(t' -t)

Remarque 4.1 On fait les observations suivantes.

(1) L'équation de Fokker-Planck est une équation aux dérivées partielles linéaire, dont la solution unique sera fixée par la donnée d'une condition initiale p(x,t0) = p0(x) étant une fonction donnée à l'avance. D'autre part, la solution doit être une fonction positive et intégrable, non seulement localement, mais encore sur tout le domaine accessible à la variable x. Ces conditions définissent des conditions aux limites ¹ assurant l'appartenance de l'ensemble des solutions particulières obtenues en tant que modes propres.

(2) L'équation de Fokker-Planck est déterminée par la fonction de dérive u(x,t) qui caractérise le déplacement du mouvement, et la fonction ó(x,t) > 0 qui, elle, caractérise la diffusion.

(3) L'équation de Fokker-Planck est dite linéaire² si :

u(x,t) = a1 +a2x et ó(x,t) = a3

et quasi-linéaire si u(x,t) est non linéaire et ó(x,t) = a3.

(4) Si l'équation différentielle stochastique est linéaire, alors la solution de l'équation de Fokker-Planck est gaussienne.

4.2.1 L'origine de l'équation de Fokker-Planck

Soit Xt un processus de Markov de probabilité de transition P(x,t + h|x0,t). On souhaite calculer la dérivée temporelle de la densité de probabilité p(x,t + h) de la variable aléatoire X à l'instant t + h.

fp(x,t + h) = P(x,t + h|x0,t)p(x0,t)dx0

a partir des relations

fP(x,t + h|x0,t) = ä(z - x)P(z,t + h|x0,t)dz

et

ä(z-x) = ä(z-x+x0-x0)

=	8 ? n=0	(z -x0)^n n!	~ ?)n ä(x - x0) ?x0
=	8 ? n=0	(z -x0)^n n!	~ )n - ? ä(x0 - x) ?x

1. Il convient de noter que certaines conditions aux limites s'expriment physiquement.

2. Le terme linéaire se rapporte ici aux propriétés de u(x,t) et ó(x,t). L'équation de Fokker-Planck est, elle, toujours linéaire pour p.

on établit que

P(x,t +h|x0,t) =

8 ? n=0

1 ? ( axy ä(x0 - x) f (z- x0)^nP(z,t + h|x0,t)dz n!

en introduisant les moments

Mn(x0,t,h) = E((Xt+h -Xt)ⁿ|Xt = x0)

= f (y - x0)ⁿP(y,t +h|x0,t)dy

on a

P(x,t +h|x0,t) =	8 ? n=0	1 ( ? ?x )n ^n! ä(x0 - x)Mn (x0, t , h )
= (1+ = (1+			8 ? n=1 8 ? n=1	1 ( ?x ) n n! Mn (x0 , t, h )) ä(x0 -x) 1 ( ?x) n n! Mn(x,t,h)) ä(x - x0)

Mais en considérant que h est petit, on peut écrire

p(x,t + h) - p(x,t) = h?p(x?t ,t) + O(h²)

= f (P(x,t + h|x0,t) - 1)p(x0,t)dx0

=

8

?

n=1

, , ! ( - ?xy f ä(x-x0)Mn(x,t,h)p(x0,t)dx0

= to ( ? n! Mn(x,t,h) p(x,t)

n=1- ?x)

Et en utilisant un développement de Taylor pour M_n(x,t,h), on a

d

M_n(x,t,h) = f (x - x0)nP(x,t|x0,t)dx+ dh J

1 (x - x0)nP(x,t|x0,t)dX _h=0h+ O(h2)

= hD⁽ⁿ⁾(x,t) + O(h²)

car

P(x,t|x0,t) = ä(x -x0)

et avec

nD(n)(x,t) dh

= M (x,t ,h)h=0

d'oil

(x t)

p(x,t + h)- p(x,t) = hapat , ' +O(h²)

=

co

E

n=1

1 ( _a )n D(n)(x,t)p(x,t) n! ax

On en déduite les équation de Kramers-Moyal [25]

ap(x,t) at

=

co E n=1

n1! ( aa x) n D(n) (x,t)p(x,t)

(4.7)

dont l'équation de Fokker-Planck (4.2) est un cas particulier de l'équation de Kramers-Moyal (4.7). Lorsque le bruit est un processus gaussien, comme un bruit blanc, on démontre que les coefficients de Kramers-Moyal sont nuls pour n = 3, c'est-à-dire

D⁽ⁿ⁾(x,t) = 0, si n = 3

Dans ces conditions, l'équation de Fokker-Planck s'écrit

avec	ap(x,t) at		a1 a2 ax = - D(1)(x" 2 ax^2 t)p(x t) + D^(2) (x,t)p(x,t)
D^(1)(x,t) = lim 1 h?0 h		M1(x,t,h) = lim h?0		1 hE(Xt+h - Xt|Xt = x)

et

_hE((Xt+h - Xt)²|Xt = x) 1

M2(x,t,h) = lim

h?0

D⁽²⁾(x,t) = lim 1

h?0 h

4.2.2 Modélisation d'une équation physique

Avant d'aborder la modélisation, nous annonçons deux théorèmes, donnant l'existe d'une relation simple entre les deux différentielles Itô et Stratonovitch, qui permet de passer de l'une à l'autre pour la résolution d'une équation différentielle stochastique

Théorème 4.1 (Différentielle d'Itô convertit en différentielle de Stratonovitch) Si un processus stochastique Xt satisfait l'équation d'Itô :

dXt = u(t,Xt)dt +a(t,Xt)dWt

alors il satisfait également l'équation de Stratonovitch :

dXt = u(t,Xt)dt + a(t,Xt)?dWt

oil le coefficient de dérive modifié u(t,Xt) est défini par :

u(t,x) = u(t,x) - 1 2 a(t,x)aa a(t,x)

x

Théorème 4.2 (Différentielle de Stratonovitch convertit en différentielle d'Itô) Si un processus stochastique Xt satisfait l'équation de Stratonovitch :

dXt = u(t,Xt)dt + ó(t,Xt) ?dWt

alors il satisfait également l'équation d'Itô :

dXt = u(t,Xt)dt +ó(t,Xt)dWt

oil le coefficient de dérive modifié u(t,Xt) est défini par:

u(t,x) = u(t,x) + ₂ó(t,x)^?ó(t,x)

1

?x

On considère une équation physique déterministe unidimensionnelle soumise à un bruit blanc î(t), sous la forme suivante

ÿx(t) = u(x,t)+ó(x,t)î(t) (4.8)

On suppose que le bruit est un processus centré E(î(t)) = 0 et de fonction de corrélation R(s,t) = E(î(s)î(t)) = ä(t - s). Le bruit blanc est modélisé par un processus de Wiener standard Wt. L'équation physique (4.8) se transcrit directement en équation de Stratonovitch

dXt = u(t,Xt)dt +ó(t,Xt)?dWt (4.9)

et se convertit en équation d'Itô sous la forme (appliquons le théorème 4.2)

dXt = A(t,Xt)dt +B(t,Xt)dWt (4.10)

Avec

A(x,t) = u(x,t) + ₂ó(x,t)^?ó(x,t)

1

?x

et

B(x,t) = ó(x,t)

D'où l'équation physique (4.8) se traduit en équation d'Itô sous la forme suivante

( )

dXt = u(t,Xt) + ¹ ₂ó(t,Xt)^?ó(t,Xt)dt+ó(t,Xt)dWt (4.11)

?x

L'équation de Fokker-Planck pour l'équation physique (4.8), s'écrit alors

?p ?t

t ~ ? u(x,t) · p + 1 ? = - 2ó(x,t) ? + 1 ?x (ó(x,t) · p) ?x2 (ó(x,t) · p) ?x 2

Dans le cas multidimensionnel, on note x = (x1,...,x_n). Le vecteur î(t) = (î1(t),...,î_m(t)) représente des bruits blancs centrés et de fonctions de corrélation normées

Ri j(s,t) = E(îi(t)îj(t)) = ä(t - s)

Les quantités u,ó1,...,ó_m sont des champs de vecteurs sur Rⁿ. L'équation physique multidimensionnel, s'écrit sous la forme

ÿx(t) = u(x,t) + ó1(x,t)î1(t) +
·
·
· + ó_m(x,t)î_m(t) (4.12)

se traduit mathématiquement par l'équation de Stratonovitch

dXt = u(t,Xt)dt +

m ? i=1

ói(t,Xt) 0 dWi(t)

oil W1(t),...,W_m(t) sont m processus de Wiener standard à valeurs dans Rⁿ. Cette équation correspond à l'équation d'Itô pour i = 1,...,n

1 ,---,^m dXl = u^i (t ,Xt) + , L z j=1 Si l'on pose

n ? k=1

j ?ai.(t,Xt)) m · (t x) ^i dt + ? ó^i(t ,Xt)dWj^i (t) , t ?xk j=1

n

?

i=1

L=

n

Ai(x) ? + 1 ? Bi j(x) ?2 ?xi 2 ?xi?xj

i, j=1

avec

?óⁱ _j(t,Xt)

ók j(t,Xt) ?xk

n

?

k=1

1 m

Aⁱ(x) = uⁱ(t,Xt)+ , ?

z j=1

et

n

B^ij(x) = ? óⁱk(x)ó^j_k(x)

k=1

alors le processus de diffusion Xt vérifie l'équation de Fokker-Planck

?p ?t

= Lp

Exemple 4.1 (L'oscillateur de Van Der Pol) On considère l'équation de Van Der Pol pour l'oscillateur ³ électrique, soumis à une fonce d'excitation aléatoire F, qui est centrée et de fonction de corrélation E(FtFt+h) = 2óä(h),

_X2

it + a (_b2 --11 ^Xÿ + ù²0X = F(t) (4.13)

L'équation de Van Der Pol (4.13) est utilisée pour modéliser des oscillateurs entretenus. Elle n'est pas linéaire et n'a pas de solution explicite. Les paramètres caractéristiques de cette équation sont : la pulsation ⁴ propre ù0, le paramètre de réaction ⁵ a et le paramètre de contrôle b.

3. Électricité : dispositif électrique qui sert à produire un courant alternatif périodique de fréquence déterminée.

4. Physique : augmentation momentanée et périodique de l'intensité d'une onde.

5. Physique : processus de modification de la structure d'un noyau atomique avec libération d'énergie.

On écrit cette équation du second degré sous la forme d'un système d'équations du premier degré

(X^ÿ = Y(X2 ) (4.14)
Yÿ = -a b2 - 1 Y - ù² ₀X + F(t) se traduit mathématiquement par l'équation de Stratonovitch

(

dXt = Ytdt ( _(X2 ) )

dYt = - a b2 - 1 Yt + ù2 0Xt dt + 2ó ? dWt

l'équation (4.13) se traduit à un système d'équations d'Itô sous la forme suivante

IdXt = Ytdt( (X2 ) ) (4.15)
dYt = - a b2 - 1 Yt + ù2 0Xt dt + 2ódWt L'équation de Fokker-Planck pour l'équation de Van Der Pol (4.13), s'écrit alors

?p
?t

~ _'x2 ~ ~

? ? + ? 0xp) + 1 ?2

= - _?x(yp) + a b2 - 1 yp ?y^(ù2 ?y2 (2óp)

?y 2

= -y	~x2 ~ ? ? ?y p + ó ?2 ?x p + a b2 - 1 ?y(yp) + ù^2 0x ? ?y2 p

Simulation numérique de l'équation de Van Der Pol

La fonction snssde2D⁶ permettre de simulée numériquement la solution approchée des systèmes d'équations différentielles stochastiques deux dimensions par des schémas numériques classiques.

R> help("snssde2D")

R> example("snssde2D")

R> snssde2D(N, T = 1, t0, x0, y0, Dt, driftx, drifty, diffx, diffy,

+ Step = FALSE, Output = FALSE, Methods = c("SchEuler",

+ "SchMilstein", "SchMilsteinS", "SchTaylor", "SchHeun",

+ "SchRK3"), ...)

6. simulation numerical solution of stochastic differential equations two-dimensional.

Details:

N La taille de processus.

M Le nombre de trajectoire à simulée.

T L'instant final (par défaut égale à 1).

t0 L'instant initial.

x0 & y0 Les valeurs initiaux.

Dt La discrétisation où le pas (Si Dt est fixée alors par défaut

T = t0 + Dt * N).

driftx & drifty Coefficient de dérive (une expression qui dépende de x, y et de t).

diffx & diffy Coefficient de diffusion (une expression qui dépende de x, y et de t).

Output Pour sauvegardée les résultats de simulation sous forme Excel

(par défaut c'est FALSE).

Step Pour l'animation, simulation étape par étape (par défaut c'est FALSE).

Methods Différents méthodes de simulation (par défaut schéma d'Euler),

avec SchEuler(3.23), SchMilstein(3.24),SchMilsteinS(3.25), SchTaylor(3.26), SchHeun(3.27), SchRK3(3.28).

Posant par exemple (b,ù0,ó) = (4,2,0.5) et x0 = y0 = 0. On remarque deux états du ce système physique très différents pour le paramètre a = 0 et a > 0,donc on a l'équation physique

f X· +4X = F(t), a = 0

~X2 ~

X· + 2 16 _ 1 ^Xÿ + 4X = F(t), a > 0

se convertit mathématique en équation d'Itô

(

dXt = Ytdt , a = 0

dYt = _4Xtdt + dWt

et

(

dXt = Ytdt ( _(X2 ) ) , a > 0

dYt = _ 2 16 _ 1 Yt + 4Xt dt + dWt

R> fx <- expression( y ) R> gx <- expression( 0 )

R> fy1<- expression( -4*x ) R> gy <- expression( 1 )

R> snssde2D(N = 10000, T = 1, t0 = 0, x0 = 0, y0 = 0, Dt = 0.1,

+ driftx = fx, drifty = fy1, diffx = gx, diffy = gy)

R> fy2<- expression( -(2*((x/4)^{^}2 -1)*y + 4*x) )

R> snssde2D(N = 10000, T = 1, t0 = 0, x0 = 0, y0 = 0, Dt = 0.1,

+ driftx = fx, drifty = fy2, diffx = gx, diffy = gy)



FIGURE 4.1 - L'oscillateur de Van Der Pol, régime permanent sinusoïdal a = 0.

FIGURE 4.2 - L'oscillateur de Van Der Pol, régime permanent non sinusoïdal a > 0.

4.2.3 Existence d'une solution

Soit £ l'opérateur différentiel

1 ói j(x,t) ?2

£ = 2 ? +?_i ui(x,t) ?

?xi?xj ?xi

i j

On considère l'equation de Fokker-Planck

?p ?t

= L * p

On etablit que lorsque les coefficients u(x,t) et ó(x,t) sont lipschitziens (Chapitre 3 section 3.3.2) et verifient une condition supplementaire d'uniforme ellipticity, l'equation de Fokker-Planck a une solution unique. Plus precisement, on a le resultat suivante.

Théorème 4.3 On suppose que les coefficients de l'équation de Fokker-Planck sont lipschitziens et qu'il existe une constante c > 0 telle que les coefficients uij(x,t) vérifient la condition

? óij(x,t)yiyj = c|y|², ?x,y ?

ij

où |y| est la norme de y, alors il existe une unique solution de l'équation de Fokker-Planck p(s,x;t,y) continue et à dérivées partielles ?xip ^et ?xi?xj p continues.

4.3 Processus de diffusion stationnaires

Un processus aleatoire de diffusion est dite stationnaire si sa densite de probabilite de transition p(s,x;t,y), verifier l'equation suivante

ð(y) = f p(s,x;t,y)ð(x)dx, ?y ?

oil ð(y) est appelee distribution stationnaire du processus si il existe. Pour laquelle on a :

?p

?t = 0

Peut être obtenue explicitement ou numeriquement en resolvant l'une ou l'autre des equations de Fokker-Planck qui ne depend pas du temps (progressive (4.16) oil retrograde (4.17))

d 1 d²

dy (u(y)ð(y))- 2 dy2 (ó2 (y)ð(y)) = 0 (4.16)

oil

2

d

u(x) ^d ð(x) + 2 dx21 ó²(x) ð(x) = 0 (4.17)

Exemple 4.2 (Distribution stationnaire de l'équation de Langevin) Soit à nouveau l'equation de Langevin (Chapitre 3 section 3.3.3),

dXt = -aXtdt +v2DdWt

On a E(Xt) = 0 et la fonction de corrélation R(s,t) = (D/a)e-^a|^t-^s| dépende de h = t -s, donc le processus de Langevin est stationnaire au sens faible, c'est-à-dire qui il existe une distributions stationnaire qui vérifier

?p ?t

= 0

d'après l'équation (4.16) ,on a

d 1 d²

d (yð(y)) + D d2

(-ayð(y)) - dy2 (2Dð(y)) = 0 ? a dy2 ð(y) = 0

dy 2 dy

? ayð(y) +D^d_dyð(y) = 0

? ayð(y) = -D ^d_dyð(y)

a _Dy

= -

ð¹(y)

?

ð(y)

? ln|ð(y)| = - 2aDy2 +C

? ð(y) = K exp (- ₂^a_Dy²)

Avec K = ¹donc on trouve la distributions stationnaire ð(y) de de Langevin,

v2ð(D/a) , l'

qui suit une loi gaussienne N(0,D/a)

1

ð(y) = exp ( - 1 a y2)

\I 2ðD 2D
a 4.4 Classification des processus de diffusion linéaire

Soit {Xt,0 = t = T} un processus de diffusion, solution de l'équation différentielle stochastique suivante

dXt = A(t,Xt)dt +B(t,Xt)dWt, X0 = x0 (4.18)

L'équation (4.18) est dite linéaire si les coefficients de dérive A(t,x) et de diffusion B(t,x) sont de la forme suivante

A(t,x) = C(t) +D(t)x, avec C : [0,T] ? Rⁿ et D : [0,T] ? 9n,n(R)

B(t,x) = E(t)+F(t)x, avec E : [0,T] ? 9n,m(R) et F : [0,T] ? L(Rⁿ,9_n,m(R)) où L(Rⁿ,9_n,m(R)) est l'espace des fonctions linéaires continues de Rⁿ dans 9n,m(R).

Définition 4.1 Une équation différentielle stochastique linéaire est dite homogène si C(t) = E(t) = 0. Et elle est linéaire au sens faible si F(t) = 0.

Maintenant, on considère l'équation différentielle stochastique (4.18) de coefficient de dérive A(t,x) = r(è - x) avec r > 0, et coefficient de diffusion B(t,x) qui peut être constant, linéaire dépend de x , où du type polynômial. Ces trois types de coefficient de diffusion conduite, respectivement, aux distributions stationnaires Gaussiennes, Gamma, et Bêta⁷. Nous ressemblant certaines de leurs propriétés dans le suivant.

4.4.1 Processus de diffusion de type N

Pour ce modèle, nous avons :

Coefficient de dérive : A(t,Xt) = r(è - Xt), r > 0.

Coefficient de diffusion : B(t,Xt) = v2ó, ó > 0.

ÉDS : dXt = r(è - Xt)dt + v2ódWt.

( (x - ²

Distribution stationnaire : ð(x) = /1

V2ðä^exp 2ä0) ) , ä = ó

r .

Statistique : moyenne, mode = è,variance = ä.

Preuve D'après l'équation (4.16) ,on a

d 1 d2d d2

dx (r(è - x)ð(x)) - dx2 (2óð(x)) = 0 ? r _dx ((è - x)ð(x)) - ó _dx2 ð(x) = 0
2

? r(è - x)ð(x) - ó ^d ð(x) = 0

dx

? r(è - x)ð(x) = ó ^d ð(x)

dx

ð^1(x) ? ð(x)

= ór (è -x)

? ln|ð(x)| = - 2r _ó(è - x)² + C ? ð (x) = K exp (- 2ró (x - è)²)

, donc on trouve la distributions stationnaire ð(x) de ce type d'équation, qui

Avec K = v2ð1(ó/r)

suit une loi gaussienne N(0,ó/r)

1 (x- è)²) _s ó

ð(x) _v2ðä exp ( 2ä ) , u = r

7. La loi bêta standard B[0,1]

Simulation numérique de la distribution stationnaire de modèle de type N

Dans cette simulation il s'agit pas de résoudre numériquement l'équation différentielle ordinaire d'ordre deux (4.16) oil (4.17) pour trouver la distribution stationnaire ð(x) de processus Xt. L'idée est de simulée un flux de M trajectoires de ce processus et en coupant à un instant donnée t, donc on obtient un échantillon de taille M de ce modèle à l'instant t = v que l'on note par X_v = (x¹ _v,x² v,...,xM _v ). A partir de sa, en fait une analyse statistique a X_v, ainsi l'ajustement de la distribution de la variable aléatoire X_v utilisant deux méthodes, la méthode de histogramme et la méthode du noyau [38]. Avant d'aborder les simulations, nous donnons une présentation de la méthode du noyau.

Définition 4.2 (Méthode du noyau) La méthode du noyau est une méthode non paramétrique, qui consiste à estimer la densité de probabilité inconnue fX(x) d'une certaine variable aléatoire X, définie sur un domaine D, et dont on connaît qu'un échantillon E = (X1,X2,...,X_n) de sa réalisation. Un estimateur _àfn,h(x), dite du noyau, de la variable X, au vu de l'observation de E, est donné par :

, 1 ^àfn,h(x) = nh

n ? i=1

K (x hXi)

K est une fonction, dit noyau, et l'estimateur de cette fonction est caractérisée par le paramètre de lissage h. Cette fonction est supposé bornée, et elle vérifie :

K(x) = 0 , f K(x)dx = 1

D

Les noyaux les plus utilisées sont Gaussien, cosinus, rectangulaire, triangulaire. Le problème fondamental dans l'utilisation de cette méthode est comment choisir le paramètre h, pour obtenir un estimateur optimal de fX(x).

On considère par exemple le modèle de Vasicek généralisé (VAG) décrit par l'équation différentielle stochastique

dXt = r(è - Xt)dt + ædWt,X0 = x0 (4.19)

avec æ = v2ó, la forme explicite de la solution de l'équation (4.19) est :

t , ,

Xt = è + (X0 -è)e^-rt + æ ^f e-rv-s)dWs

0

La fonction AnaSimX permettre de simulée numériquement un échantillon X_v de taille M à partir d'une EDS.

R> help("AnaSimX")

R> AnaSimX(N, M, t0, Dt, T = 1, X0, v, drift, diff, Output = FALSE,

+ Methods = c("Euler", "Milstein", "MilsteinS", "Ito-Taylor",

+ "Heun", "RK3"), ...)

Posant (r,9,ó) = (2,-2,1) et x0 = 4, donc on ale modèle

dXt = 2(-2 - Xt)dt + 2dWt,X0 = 4 / R> f <- expression( 2*(-2-x) ) R> g <- expression( sqrt(2) )
R> AnaSimX(N = 1000, M = 100, t0 = 0,				Dt = 0.01,	X0 = 4, v =	9,
+	drift = f, diff = g)
R> X
[1]	-2.1943293	-2.2719547	-3.1151950	-0.6626404	-2.5430828	0.3602131
[7]	-3.2372479	-2.0337318	-1.9326694	-1.4953072	-2.9416240	-1.9643427
[13]	-2.0874472	-2.4104522	-2.6901473	-2.1066782	-1.6935978	-2.7531990
.	.	.	.	.	.	.
.	.	.	.	.	.	.
.	.	.	.	.	.	.
[85]	-2.3721165	-2.0848300	-2.5641703	-1.6872738	-0.6147365	-2.0119820
[91]	-2.5983331	-2.0200177	-3.7215718	-1.4605711	-2.3753241	-2.3620768
[97]	-1.9748765	-2.2800181	-2.7551030	-2.1111137	-2.3917035	-1.0311714



FIGURE 4.3 - Simulation un échantillon de taille 100 à partir du modèle VAG.

Statistique descriptive:

R> summary(X)

Min. 1st Qu. Median Mean 3rd Qu. Max.

-3.7570 -2.4160 -2.0520 -2.0080 -1.6070 0.3602 R> var(X)

[1] 0.5546319

La fonction Ajdnorm permettre d'estimée les paramètres de la loi normale par la méthode de maximum de vraisemblance, ainsi l'intervalle de confiance pour á = 0.95.

R> Ajdnorm(X, starts = list(mean = 1, sd = 1), leve = 0.95) Profiling...

$summary

Maximum likelihood estimation

Call:

mle(minuslogl = lik, start = starts)

Coefficients:

Estimate Std. Error mean -2.0063018 0.07482636 sd 0.7445129 0.05291131

-2 log L: 222.5311

$coef

mean sd

-2.0063018 0.7445129

$AIC

[1] 226.5311

$vcov

mean sd

mean 5.598984e-03 -3.358252e-08
sd -3.358252e-08 2.799606e-03

$confint

2.5 % 97.5 %

mean -2.1543897 -1.858205
sd 0.6517166 0.861606

La fonction hist_general permettre d'ajusté la distribution de l'échantillon X_v par la méthode de l'histogramme. Et la fonction Kern_general c'est l'ajustement par la méthode du noyau.

R> help(hist_general) R> help(Kern_general) R> hist_general(Data = X, Breaks = 'Sturges', Law = "Norm")

R> Kern_general(Data = X, bw='Bcv', k = "gaussian", Law = "Norm")



FIGURE 4.4 - Ajustement de la distribution stationnaire du modèle VAG par la méthode d'histogramme.

Test de Kolmogorov-Smirnov :

FIGURE 4.5 - Ajustement de la distribution stationnaire du modèle VAG par la méthode du noyau.

R> ks.test(X, "pnorm", mean = -2, sd = sqrt(1/2) )
One-sample Kolmogorov-Smirnov test

data: X

D = 0.0621, p-value = 0.8357

alternative hypothesis: two-sided

4.4.2 Processus de diffusion de type G

Pour ce modèle, nous avons :

Coefficient de dérive: A(t,Xt) = r(è -Xt), r > 0.

B(t,Xt) = vóXt, ó > 0.

Coefficient de diffusion :

dXt = r(è - Xt)dt + vóXtdWt.

ÉDS :

(x )-1+è ä e- x ä

Distribution stationnaire : ð(x) = _(è ), ä = ó 2r.

ä ä

Statistique : moyenne = è, mode = è- ä, variance = äè.

Preuve D'après l'équation (4.16) ,on a

d 1 d²

d ó d2

dx (r(è - x)ð(x)) - dx2 (óxð(x)) = 0 ? r _dx ((è - x)ð(x)) - dx2 (xð(x)) = 0

2 2

? r(è - x)ð(x) - ó d

2 _dx(xð(x)) = 0

ó

? r(è - x)ð(x) = ₂(ð(x) + xð'(x))

ó ó

? r(è - x)ð(x) - 2 ð(x) = 2 xð0(x)

ð^e(x) ? ð(x)

= 2 (r(è-x)- ó)óx 2

ð^e(x)

?

ð(x)

(2rè ₁ )1 2r ó xó



(2rè

? ln|ð(x)| = - ó 1 ) ln(x)- ó 2r x+C

~ ~

2rè

? ð(x) = Kxó ^-1 exp -2r ó x

Posant ä = ó2r avec K = ä¹-^èä / (2), d'où on trouve la distributions stationnaire ð(x) d' une équation de type G, qui suit une loi gamma ã (èä,1) ä

~x ~-1+è ä e- x ä

ð(x) = _~è ~

ä ä

Simulation numérique de la distribution stationnaire de modèle de type G

On considère par exemple le modèle Cox-Ingersoll-Ross (CIR) (Chapitre 3 exemple 3.6) décrit par l'équation différentielle stochastique

dXt = r(è - Xt)dt + óvXtdWt, X0 = x0 > 0. (4.20)

Avec la condition 2rè > ó². Posant par exemple (r,è,ó) = (0.5,3,1) et x0 = 100,Ät = 0.1, donc

on a le modèle

1

dXt = 2 (3 - Xt)dt +vXtdWt, x0 = 100

simulons numériquement un échantillon X_v=90 de taille 100 à partir de cette équation

R> f <- expression( 0.5*(3-x) ) R> g <- expression( sqrt(x) )

R> AnaSimX(N = 1000, M = 100, t0 = 0, Dt = 0.1, X0 = 100, v = 90,

+ drift = f, diff = g)

R> X

[1]	1.9781677	2.7698461	1.9426347	6.8094888	3.2716436	11.0541449
[7]	1.6495788	2.7213673	5.6279096	4.1147230	3.6330197	1.8322262
.	.	.	.	.	.	.
.	.	.	.	.	.	.
.	.	.	.	.	.	.
[89]	2.1697131	6.5484018	0.8467619	2.5333814	1.1763158	3.3743554
[95]	0.4432751	3.7510871	1.6463437	2.0691975	1.8540492	3.8536741

R> summary(X)

Min. 1st Qu. Median Mean 3rd Qu. Max.

0.5483 1.7760 2.7510 3.1900 4.2450 10.0500
R> var(X)

[1] 3.492991

Estimations des paramètres :

R> Ajdgamma(X, starts = list(shape = 1, rate = 1), leve = 0.95) Profiling...

$summary

Maximum likelihood estimation

Call:

mle(minuslogl = lik, start = starts)

Coefficients:

Estimate Std. Error shape 3.079444 0.4161674 rate 0.962066 0.1412105

-2 log L: 376.7305

$coef

shape rate

3.079444 0.962066

$AIC

[1] 380.7305

$vcov

shape rate

shape 0.17319531 0.05410881
rate 0.05410881 0.01994040

$confint

2.5 % 97.5 %

shape 2.3372820 3.972203
rate 0.7104027 1.265122

Ajustement de la distribution stationnaire ð(x) :

R> hist_general(Data = X, Breaks='Sturges', Law = "GAmma")

R> Kern_general(Data = X, bw='SJ', k = "gaussian", Law = "GAmma")



FIGURE 4.6 - Ajustement de la distribution stationnaire du modèle CIR par la méthode d'histogramme.

Test de Kolmogorov-Smirnov :

FIGURE 4.7 - Ajustement de la distribution stationnaire du modèle CIR par la méthode du noyau.

R> ks.test(X, "pgamma", shape = 3, rate = 1 ) One-sample Kolmogorov-Smirnov test data: X

D = 0.0722, p-value = 0.6741

alternative hypothesis: two-sided

4.4.3 Processus de diffusion de type B

Pour ce modèle, nous avons :

Coefficient de dérive : A(t,Xt) = r(è - Xt), r > 0, è ?]0,1[

Coefficient de diffusion : B(t,Xt) = VóXt(1 -Xt), ó > 0.

ÉDS : dXt = r(è - Xt)dt + VóXt(1 -Xt)dWt, X0 ? [0,1].

Distribution stationnaire : ð(x) =

~1 ~ ~~1-è ä ~è ~x ä ä

-1+è ä (1 - x)-1+1-è ä ,ä = ó . 2r

è - ä è(1 - è)

Statistique : moyenne = è, mode = 1 - 2ä, variance = 1+ä .

Preuve D'après l'équation (4.16) ,on a

d 1 d²

dx (r(è-x)ð(x)) - 2 dx2 (óx(1 - x)ð(x)) = 0 ? (1)

2

d ó d

(1) ? r _dx ((è - x)ð(x)) - 2 dx2 (x(1 - x)ð(x)) = 0

ó d

? r(è - x)ð(x) - 2 _dx(x(1 - x)ð(x)) = 0

2r

? (è - x)ð(x) = ((1 - 2x)ð(x) + x(1 - x)ð⁰(x))

~2r

ó ~

? ó (è - x) - (1 - 2x) ð(x) = x(1 - x)ð'(x)

ð^e(x)

?

ð(x)

1 (2r(è-x)-ó(1 -2x)₎

=

ó x(1 -x)

ð^e(x)

?

ð(x)

=

1 (2rè-ó 2r(1 - è) -ó)
ó x 1 - x

2rè- ó 2r(1 -

?ln|ð(x)| = ó ln |x|+ _ó) - ^óln |1- x| + C

2rè

? ð(x) = Kxó ^-1(1 - x)2r(1-è)

ó -1

Posant ä = ó2r avec K égale à

K=

~1 ~

ä

~è ~~1-è ~ ä ä

d'où on trouve la distributions stationnaire ð(x) d'une équation de type B, qui suit une loi bêta standard B _~è ~

ä, 1-è

ä

ð(x) =

~1 ~ ~~1-è ä ~è ~x ä ä

-^1+^èä(1 - x)-1+1-è ä

Simulation numérique de la distribution stationnaire de modèle de type B

On considère le processus de diffusion de Jacobi Xt, solution de l'équation différentielle stochastique

\/

dXt = r(è - Xt)dt + óXt(1 - Xt)dWt, X0 = x0 ? [0,1]. (4.21)

Avec è ?]0,1[. Posant par exemple (r,è,ó) = (2,0.5,1) et x0 = 0,Ät = 0.01, donc on ale modèle

₍₁ ) \/

dXt = 2 2 - Xt dt + Xt(1 - Xt)dWt, X0 = 0.

Simulons numériquement un échantillon _Xv=9 de taille 100 à partir de cette équation

R> f <- expression( 2*(0.5-x) )

R> g <- expression( sqrt(x*(1-x)) )

R> AnaSimX(N = 1000, M = 100, t0 = 0, + drift = f, diff = g) R> X				Dt = 0.01,	X0 = 0, v =	9,
[1]	0.61568827	0.22363297	0.23333832	0.34456437	0.86391343	0.67870198
[7]	0.82226081	0.06663198	0.60156924	0.58966279	0.57626142	0.63231225
.	.	.	.	.	.	.
.	.	.	.	.	.	.
.	.	.	.	.	.	.
[89]	0.43607965	0.61750493	0.19505968	0.87196355	0.79428433	0.61846946
[95]	0.72321640	0.45797652	0.88082305	0.53119498	0.49339106	0.48472011

R> summary(X)

Min. 1st Qu. Median Mean 3rd Qu. Max.

0.01721 0.35020 0.53090 0.51870 0.71230 0.94530

R> var(X)

[1] 0.0502775

Estimations des paramètres :

R> Ajdbeta(X, starts = list(shape1 = 1, shape2 = 1), leve = 0.95) Profiling...

$summary

Maximum likelihood estimation

Call:

mle(minuslogl = lik, start = starts)

Coefficients:

Estimate Std. Error shape1 2.020047 0.2739796

shape2 1.909258 0.2571918 -2 log L: -23.83837

$coef

shape1 shape2

2.020047 1.909258 $AIC

[1] -19.83837

$vcov
	shape1	shape2
shape1	0.07506481	0.05515151
shape2	0.05515151	0.06614760
$confint
	2.5 %	97.5 %

shape1 1.531554 2.607888
shape2 1.450570 2.460968

Ajustement de la distribution stationnaire ð(x) :

R> hist_general(Data = X, Breaks='Sturges', Law = "Beta")

R> Kern_general(Data = X, bw='Bcv', k = "gaussian", Law = "Beta")



FIGURE 4.8 - Ajustement de la distribution stationnaire de processus de diffusion de Jacobi par la méthode d'histogramme.

FIGURE 4.9 - Ajustement de la distribution stationnaire de processus de diffusion de Jacobi par la méthode du noyau.

Test de Kolmogorov-Smirnov :

R> ks.test(X, "pbeta", shape1 = 2, shape2 = 2 ) One-sample Kolmogorov-Smirnov test data: X

D = 0.0803, p-value = 0.5385

alternative hypothesis: two-sided

Propriété 4.1 (Diffusion de Pearson) Une diffusion de Pearson⁸ {Xt,t = 0} [20], est une solution stationnaire à une équation différentielle stochastique de la forme suivante :

\1 dXt = r(9 - Xt)dt + 2a(aX2 _t + bXt + c)dWt (4.22)

Avec r > 0 c'est la vitesse de convergence vers l'équilibre, et 9 c'est la moyenne à l'équilibre, a > 0. Les paramètres qui caractérise la diffusion a, b et c sont des constants dans R. Remarquons que ce modèle est stationnaire pour certains conditions sur a,b et c, on a

· Si a = b = 0 et c = 1 ? Processus de diffusion linéaire de type N.

· Si a = c = 0 et b = 1 2 ?Processus de diffusion linéaire de type G.

· Si -a = b = 1 2 et c = 0, avec 9 ?]0,1[ ? Processus de diffusion linéaire de type B.

4.5 Calculs de l'instant de premier passage

Les instants de premier passage "IPP" jouent un rôle très important en pratique. L'objet de cette étude est de déterminer la densité de probabilité de la variable aléatoire t l'instant de premier passage d'un processus de diffusion Xt. Dans le chapitre 2 la technique des martingales et théorème d'arrêt qui nous a permet de déterminée la loi du temps d'atteinte du mouvement brownien a un point donnée a (le mouvement brownien est un martingale). Dans cette section, on peut utiliser les techniques de la transformée de Laplace, pour déterminer la loi de probabilité de t.

Soit {Xt,t = 0} un processus de diffusion solution de l'équation différentielle stochastique suivante :

dXt = u(t,Xt)dt +a(t,Xt)dWt, X0 = x0 (4.23)

On admet que les coefficients u(x,t) et a(x,t) remplissent les conditions d'existence et unicité de la solution Xt (Chapitre 3 section 3.3.2). Soit D une partie de R que l'on va supposer ouverte et bornée et soit a un élément de D. On défini la variable aléatoire t_a :

t_a = inf{t = 0 : Xt ? aD,Xt = a}

t_a représentant le premier instant où le processus Xt frappe la frontière du domaine D.

8. La fonction PDP permettre de simulée Pearson Diffusions Process. voir help(PDP) pour plus de détails (Package Sim.DiffProc).

4.5.1 IPP d'une diffusion en attraction M ó s=1(Vt)

Dans une première partie, pour le cas s = 1, nous allons effectuer une etude analytique,

afin de determiner explicitement la densite de la variable ô(s=1)

c , cette loi de probabilite permet

d'approximer pour ce modèle, le taux des polluants qui depasse la frontière ?D(0,c). Pour le cas
de s > 1, la simulation est utilisee pour determiner approximativement la densite de la variable

(s>1)

ô_c [2, 4, 6, 7].

Modèle Mós=1(Vt)

L'equation de ce modèle, pour le processus radial Rt, est de la forme :

R^dR0 = a

t = ó22(1R-t Wt

^k)dt + ad ,

{

a > 0, t ? [0, T] (4.24)

Avec k = 2K ó2 pour des raisons de calcul. On utilise l'equation (4.24) pour determiner explicitement la loi de probabilite de ô(1)

c.

La probabilite de transition ft(r|a) represents la densite de Rt avec R0 = a, qui est donne par l'equation de Fokker-Planck retrograde (4.3) :

? ó² ? _ó2 ?2

?t ft(r|a) = (1-k) 2a ?a ft(r|a) + 2 ?a2 ft(r|a) (4.25)

avec la condition initiale f0(r|a) = ä(a - r), où ä(.) fonction de Dirac.

La forme analytique de la densite ft(r|a) peut être explicitement determinee en inversant son transforment de Laplace fë(r;a), qui est la solution de l'equation suivante :

_ó2 ?2 ó2 ?

2 ?a2 fë(r; a) + (1-k) 2a ?a fë(r; a) - ëfë(r; a) = -ä(a - r) (4.26)

Puisque :

a?t ft (rla) =

t

^e

L ( ) ët ?t ft(r|a)dt = - f0(r; a) + ëfë(r; a)

f

0

Une technique standard [39, 5], la solution de l'equation (4.26) pouvoir être formellement ex-prime comme :

-²g1,ë(r)g2,ë(a)

a > r (4.27)

fë(r;a) =

ó2 (g1,ë(r)g⁰2_,ë(a) - g⁰1_,ë(r)g2,ë(a)

),

avec :

)

g1ë(a) = amIm av2ë )

et g2_,ë(a) = a^mH_m

a v2ë

ó ó

avec m = k - ¹₂, I_m et H_m sont respectivement les fonctions de Bessel modifiée [39] de première et seconde ordre, définies par :

Im(z) =

8 ? i=1

, I-m(z) - Im(z) (-1)^i ( z ) 2i+m n et Hm(z) = i!(i+m+ 2 sin(mð) rn+1) 2 I \

On remplace _g1,ë et _g2,ë dans l'équation (4.27), on trouver :

₂2 (a )m+1 _r2³ Hm ^a2ë _Im rv2A a > r (4.28)

,

fë(r; a) =

ó va r ) ó ó

l'inverse de l'équation (4.28), est donnée par :

~ 1 ~_~a ~k _r ~ ~

r3 _-a² + _r2 ~ _ar ~

ft(r|a) = a exp I_m (4.29)

ó²t r 2ó²t ó²t

Maintenant nous définissons la variable aléatoire ô⁽¹⁾ _cpar :

_ô(1)

c = inf{t = 0 : Rt = c, R0 = a}

On a,

f (1) (ë) = E(e^-ë41) |R0 = a) = g2,ë(a)

a = c

ôc g2 _ë (c) ,

Pour des petites valeurs de c, nous obtenons ([5, 6])



(av_ó2ë)^m

fe) (ë) =

2m-1(m)Hm

av2ë)

en inversant (4.30), on trouve l'expression analytique de la densité de probabilité de la variable aléatoire ô(1) c , donnée par :

1 i a2\ m-1 exp ( a2

fôc (1)(ô) = 2(m) 2ó²ô ) 2ó2ô) (4.31)

Le changement de variable :

donne

X=

a^2 ó^2ô

(1 )m
fX(x) = 2 _xm-1_e-2

(m)

d'où, on a :

a2 1)F,41)(ô) = 1 - F(m1/2) (ó2 ô

1 _ 2K 1- 9

avec m = k - 2 2 et (m, 1/2) est la loi gamma.

-- a²

Remarque 4.2 La condition 2K > ó², doit être vérifier (Chapitre 3 section 3.5.1)

(4.30)

ó

Simulation la variable ô(1)

c

La fonction tho_M1 permettre de simulée un échantillonne de taille M = 50 de la variable

aléatoire ô(1)

c . Avec K = 2 etó= 1, le pas Ät = (T-t0)/N, et le c = v = 0.5.

R> tho_M1(N = 1000, M = 50, t0 = 0, T = 1, R0 = 1, v = 0.5, K = 2,

+ R> FPT		sigma =	1)										[1]
0.245	0.037	0.113	0.050	0.060	0.046	0.107	0.182	0.232	0.056	0.032		[12]
0.219	0.085	0.065	0.090	0.032	0.083	0.147	0.065	0.053	0.072	0.073		[23]
0.074	0.040	0.099	0.101	0.177	0.072	0.089	0.198	0.108	0.088	0.058		[34]
0.071	0.181	0.209	0.117	0.082	0.057	0.136	0.078	0.133	0.026	0.121		[45]
0.079	0.331	0.158	0.090	0.091	0.074

Les deux figures 4.10 et 4.11 donne une illustration de l'instant où le processus Rt frappe la frontière du domaine D. On défini le domaine D par un cercle de rayon c = v = 0.3 (en deux dimensions), et par une sphère de rayon c = v = 0.3 (en trois dimensions).



1

_ô(1)

c

FIGURE 4.10 - L'instant de premier passage du modèle Mó s=1(Vt) en 2-D.

Estimation de la distribution de Y = a2 ó2

FIGURE 4.11 - L'instant de premier passage du modèle Mó s=1(Vt) en 3-D.

R> Y = 1/FPT

R> Ajdgamma(Y, starts = list(shape = 1, rate = 1), leve = 0.95) Profiling...

$summary

Maximum likelihood estimation Call:

mle(minuslogl = lik, start = starts)

Coefficients:

Estimate Std. Error shape 3.5297529 0.68206463 rate 0.2684711 0.05574823 -2 log L: 319.8036

$coef

shape rate

3.5297529 0.2684711

$AIC

[1] 323.8036

$vcov
shape	rate		shape
0.46521216	0.035382963		rate
0.03538296	0.003107865

$confint

2.5 % 97.5 %

shape 2.3624031 5.0474887 rate 0.1731289 0.3925789 R> hist_general(Data = Y, Breaks='Sturges', Law = "GAmma")

R> Kern_general(Data = Y, bw='SJ', k = "gaussian", Law = "GAmma")



FIGURE 4.12 - Ajustement de la distribution de

1/'r^s=1

c par la méthode d'histogramme. FIGURE 4.13 - Ajustement de la distribution de

1/'r^s=1

c par la méthode du noyau.

Modèle M ó s>1(Vt)

Le processus de diffusion radial, qui décrit ce modèle est donné par l'équation différentielle stochastique suivante :

? ? ?

( ó2 ) 2 Rs-1 t -K dRt = dt + ód ^-Wt, R^st R0 = a

a > 0, t E [0,T] (4.32)

L'estimation de la densité de probabilité de la variable aléatoire ô(s>1)

c sera effectuée sur la base

de la simulation d'un flux de trajectoires. Les observations simulées de la variableô^(s>1)

c , seront

traitées statistiquement par deux méthodes, la méthode de l'histogramme et la méthode du noyau, pour estimer sa densité. Ce modèle est analytiquement difficile à résoudre.

Simulation la variable ô(s>1)

c

La fonction tho_M2 permettre de simulée un échantillonne de taille M = 50 de la variable

aléatoire ô(s>1)

c . Avec s = 2, K = 3 et ó = 1, le pas Ät = (T -t0)/N, et le c = v = 0.5.

R> tho_M2(N = 1000, M = + Sigma = 0.3) R> FPTT [1] 0.856 0.827 0.763 [12] 0.840 0.848 0.918 [23] 0.805 0.840 0.557 [34] 0.993 0.599 0.867 [45] 0.880 0.880 0.739

50, 0.795 0.952 0.821 0.790 0.602

t0 = 0, 0.710 0.662 0.833 0.935 0.794

T = 0.963 0.775 0.707 0.639 0.932

1, R0 0.948 0.894 0.661 0.868

= 2, v = 0.5, 0.620 0.895 0.769 0.599 0.597 0.751 0.973 0.579

K = 3, s = 0.788 0.831 0.676 0.971 0.694 0.732 0.832 0.904

2,

On fait l'ajustement de la variable Y = 1/ô^s=2

c par les lois : gamma, exponentiel, lognormale et

weibull. Le meilleur modèle est choisi par le critère AIC (minimum AIC).

R> Mod1 <- Ajdgamma(Y, starts = list(shape = 1, rate = 1), leve = 0.95)

$AIC

[1] 197.0462

R> Mod2 <- Ajdweibull(Y, starts = list(shape = 1, scale = 1), leve = 0.95) $AIC

[1] 203.6206

R> Mod3 <- Ajdexp(Y, starts = list(lambda = 1), leve = 0.95)

$AIC

[1] 294.9443

R> Mod4 <- Ajdlognorm(Y, starts = list(meanlog = 1, sdlog = 1), leve = 0.95) $AIC

[1] 199.8586

En remarque que la loi gamma ajuste mieux la distribution de la variable aléatoire Y = 1/ô^s=2

c ,

selon le critère AIC.

4.5.2 IPP de deux diffusion en attraction M ó m>0(V(1)

t ) ?- M 0 ó (V(2)

t )

Le processus de diffusion Xt qui décrit ce modèle est donné par l'équation différentielle stochastique suivante :

{ (ó2Xm-1 ~

X0 = a

t -K

dXt = dt + ód ^eWt,

Xm t a > 0, t ? [0,T] (4.33)

K et m sont des constantes positives, et K > ó².

L'estimation de la densité de probabilité de la variable aléatoire ô_c(V⁽¹⁾

t ,V(2)

t ) sera effectuée sur la base de la simulation. Cette variable représente l'instant de la première rencontre entre les deux insectes, défini par:

ô_c(V⁽¹⁾

t ,V(2)

t ) = inf{t = 0 : Xt = c}

La fonction tho_02diff permettre de simulée un échantillonne de taille M = 50 de la variable aléatoire ô_c. Avec K = 4, m = 0.5 et ó = 0.2, le pas Ät = 0.01, et le c = v = 0.05.

R> tho_02diff(N = 1000, M = 50, t0 = 0, Dt = 0.001, T = 1, X1_0 = 1,

+ X2_0 = 1, Y1_0 = 0.5, Y2_0 = 0.5, v = 0.05, K = 4, m = 0.5,

+ Sigma = 0.2)

R> FPT

[1]	0.140	0.085	0.104	0.177	0.112	0.098	0.067	0.085	0.142	0.128		0.086	[12]
0.085	0.100	0.080	0.132	0.120	0.108	0.083	0.089	0.074	0.073	0.085		[23]
0.057	0.163	0.114	0.076	0.106	0.167	0.076	0.110	0.105	0.095	0.098		[34]
0.102	0.112	0.096	0.120	0.066	0.097	0.098	0.097	0.091	0.083	0.091		[45]
0.109	0.071	0.181	0.157	0.083	0.093

De même on fait l'ajustement de la variable Y = 1/ô_c par les lois : gamma, exponentiel, lognormale et weibull. Le meilleur modèle est choisi par le critère AIC (minimum AIC).

R> Mod1 <- Ajdgamma(Y, starts = list(shape = 1, rate = 1), leve = 0.95)

$coef

shape rate

16.079041 1.544863 $AIC

[1] 234.4629

R> Mod2 <- Ajdweibull(Y, starts = list(shape = 1, scale = 1), leve = 0.95) $coef

shape scale

4.403975 11.394287 $AIC

[1] 235.836

R> Mod3 <- Ajdexp(Y, starts = list(lambda = 1), leve = 0.95)

$coef

lambda

0.09607808

$AIC

[1] 329.5738

R> Mod4 <- Ajdlognorm(Y, starts = list(meanlog = 1, sdlog = 1), leve = 0.95) $coef

meanlog sdlog

2.3111742 0.2560689 $AIC

[1] 236.0461

En remarque après cette simulation que la loi gamma ajuste mieux la distribution de la variable aléatoire Y = 1/'r_c, mais si on répète les simulations on remarque que les lois weibull et lognormale ajuste aussi mieux la distribution de 1/'r_c.

4.6 Conclusion

Dans ce chapitre nous avons vues l'importance de décrire l'évolution dynamique de la densité de transition du processus considéré, cela repose sur la résolution d'une équation aux dérivée partielle, qui ce n'est pas toujours facile à déterminée sa solution, car les conditions aux limites sont généralement exprimées physiquement. La simulation des processus de diffusion est un objet très puissant pour simulée l'instant de premier passage "IPP" dans les modèles attractives. Les méthodes d'ajustements de la distribution de probabilité, histogramme et méthode de noyau son efficace.

Conclusion générale

I

l y a plus d'un demi siècle que la théorie des processus de diffusion a été introduite, sous certaines conditions, les équations différentielles stochastiques de type Itô, sont des proces-

sus de diffusion. Ce travail a été consacré à l'utilisation de l'aspect théorique des EDS, comme un outil mathématique, dans la modélisation de certains phénomènes réels, présentant un intérêt pratique, où on démontre, à travers ces exemples pratiques, l'importance et l'intérêt de la simulation. Après une analyses statistique des résultats obtenus par simulation, on peut les utilisés pour une aide à la décision. Généralement, il est à remarquer que une solution analytique exacte d'une EDS n'est pas toujours facile à obtenir, car les règles classiques de différentiabilité ne sont pas toujours applicables, à cause de la propriété de martingale que doit vérifier l'intégrale stochastique d'Itô. La formule d'Itô, permet de donne une transformation simplifiée de l'équation initiale, souvent utilisée pour la résolution des EDS. Dans ce travail nous avons fourni au lecteur un package Sim.DiffProc, muni d'une interface graphique sous langage R, de manière à ce qu'on puisse simulé des modèles de diffusion d'EDS et analysé statistiquement les résultats de simulation.

Cette étude nous à permet de tracer quelques lignes perspectives, tel que le problème de la détermination analytique de la fonction de densité de transition p(s,x;t,y) ou de la variable de sortie du processus solution, nécessite la résolution d'équations aux dérivées partielles, qui ne sont pas toujours faciles à résoudre, la méthode des différence finies est parfois utilisée pour déterminer une solution numérique approchée du problème, d'ou la nécessité de développement et l'utilisation des algorithmes générales pour déterminée les densités de transition pour n'importe quel processus de diffusion, ainsi des méthodes de validation, pour les résultats obtenus par les simulations.

L

'objective de l'annexe A est de montre l'efficacité et la puissance de la simulation sous lan gage R, et donnée une idée sur la programmation mathématique, qui n'est pas très diffèrent celle de Matlab.

Code 1

R> n <- 10 ;T <- 1;

R> t <- seq (0,T, length =100)

R> S <- cumsum (2*( runif (n ) >0.5) -1)

R> W <- sapply (t, function (x) ifelse (n*x >0,S[n*x] ,0)) R> W <- as.numeric (W)/ sqrt (n)

R> plot (t,W, type ="l",xlab=expression(W[t]),ylim=c(-1,1),las=1) R> n <- 100

R> S <- cumsum (2*( runif (n ) >0.5) -1)

R> W <- sapply (t, function (x) ifelse (n*x >0,S[n*x] ,0)) R> W <- as.numeric (W)/ sqrt (n)

R> lines (t,W, lty =2)

R> n <- 1000

R> S <- cumsum (2*( runif (n ) >0.5) -1)

R> W <- sapply (t, function (x) ifelse (n*x >0,S[n*x] ,0)) R> W <- as.numeric (W)/ sqrt (n)

R> lines (t,W, lty =3)

R> legend("topleft",c("n=10","n=100","n=1000"),lty=c(1,2,3), + lwd=1,cex=0.9)

Code 2

R> phi <- function (i,t){

+ (sqrt(2)/ ((i + 0.5) *pi)) * sin ((i + 0.5) *pi*t) }

R> N <- 1000

R> t <- seq (0,1, length =N +1)

R> W <- numeric (N +1)

R> n <- 10

R> Z <- rnorm(n)

R> for (i in 2:( N +1)) W[i] <- sum (Z* sapply (1:n,

+ function(x) phi(x,t[i])))

R> plot (t,W, type ="l",ylim =c( -1 ,1),xlab=expression(W[t]),las =1) R> n <- 100

R> Z <- rnorm(n)

R> for (i in 2:( N +1)) W[i] <- sum (Z* sapply (1:n,

+ function(x) phi(x,t[i])))
R> lines(t,W, lty =2)

R> n <- 1000

R> Z <- rnorm(n)

R> for (i in 2:( N +1)) W[i] <- sum (Z* sapply (1:n,

+ function (x) phi(x,t[i])))
R> lines(t,W, lty =3)

R> legend("topleft",c("n=10","n=100","n=1000"),lty=c(1,2,3), + lwd=1,cex=0.9)

Code 3 R> phi <- function (i,t,T){					+ (2* sqrt (2*T))/ ((2 *i +1) *pi)
* sin (((2 *i	+1)	*pi*t)/(2*T))		+ }
				R> s <- 0.1; n <- 100; T <- 1;
				R> Z <- rnorm (n)
				R> Delta <- seq(0, 0.01, length =50)
				R> W <- sum (Z* sapply (1:n, function
(x) phi (x ,s	,T	)))		R> for (i in Delta ){W_h <- sum (Z* sapply
(1:n,

+ function (x) phi (x ,s+i,T )))}

R> lim_W <- abs(W_h - W)/ Delta

R> plot(Delta , lim_W , type ="l",xlab = expression ( Delta *t),las=1, + ylab = expression ( abs (W (s+ Delta *t)- W (s)) / Delta *t))
R> max(lim_W ,na.rm=T)

[1] Inf

Code 4

R> phi <- function (i,t,T){

+ (2* sqrt (2*T))/ ((2 *i +1) *pi) * sin (((2 *i +1) *pi*t)/(2*T))

+ }

R> T <- 100; N <- 1000;

R> t <- seq (0,T, length =N +1)

R> W <- numeric (N +1)

R> n <- 100

R> Z <- rnorm (n)

R> for (i in 2:( N +1)) W[i] <- sum (Z* sapply (1:n,

+ function (x) phi (x,t[i],T )))

R> plot (t,W, type ="l",ylab=expression(W[t]))

R> lines(t,W/t,col="red")

R> legend("topleft",c(expression(W[t]),expression(lim(frac(w[t],t),

+ t%->%+infinity) %~~%0)),lty=c(1),col=c("black","red"),

+ lwd=1,cex=0.9)

Code 5

R> phi <- function (i,t){

+ (sqrt(2)/ (pi * i)) * sin (pi*i*t) }

R> N <- 1000

R> t <- seq (0,1, length =N +1)

R> X <- numeric (N +1)

R> n <- 10

R> Z <- rnorm(n)

R> for (i in 2:( N +1)) X[i] <- sum (Z* sapply (1:n,

+ function(x) phi(x,t[i])))
R> plot (t,X, type ="l",ylim=c(-1,1), las =1)

R> n <- 100

R> Z <- rnorm(n)

R> for (i in 2:( N +1)) X[i] <- sum (Z* sapply (1:n,

+ function(x) phi(x,t[i])))
R> lines(t,X, lty =2)

R> n <- 1000

R> Z <- rnorm(n)

R> for (i in 2:( N +1)) X[i] <- sum (Z* sapply (1:n,

+ function (x) phi(x,t[i])))
R> lines(t,X, lty =3)

R> legend("topleft",c("n=10","n=100","n=1000"),lty=c(1,2,3), + lwd=1,cex=0.9)

Code 6

R> N = 10000; t0 = 0; Dt = 0.0001; x0 = 6; a = 2; D = 1; R> time <- c(t0 ,t0+ cumsum(rep(Dt,N)))

R> u = runif(N,0,1)

R> for (i in 1:length(u)){

+ if ( u[i] >= 0.5)

+ u[i] = +1

+ else u[i] = -1 }

R> w = cumsum(c(0,u))*sqrt(Dt)

R>

R> Ito.sum <- c(0,sapply(1:(N+1),function(x){

+ exp(-a*(time[x+1]-time[x]))*(w[x+1]-w[x])}))

R> X <- sapply(1:(N+1),function(x){

+ x0*exp(-a*time[x])+sqrt(2*D)*sum(Ito.sum[1:x])})

R> plot(time,X,type="l",las=1,xlab="time",ylab=expression(X[t])) R> mtext("Langevin Process",line=2.5,cex=1.2 )

R> mtext(bquote(dX[t]==-.(a)*X[t]*dt+sqrt(.(2*D))*dW[t]),

+ line=0.25,cex=1.2,col="red")

R> mtext(bquote(x[.(0)]==.(x0)),line=0.1,cex=0.9,adj=1,col="red") R> mtext(bquote(t[0]==.(t0)),line=0.9,cex=0.9,adj=1,col="red")

R> mtext(bquote(Delta*t==.(delta.time)),line=0.4,cex=1,adj=0,col="red") R> data.frame(time,X)

Code 7

R> N = 10000; t0 = 0; Dt = 0.0001; a = 2; D = 1;

R> x0 = 5; y0 = 6;

R> time <- c(t0 ,t0+ cumsum(rep(Dt,N)))

R> wx = cumsum(rnorm(N+1,mean=0,sd=sqrt(Dt)))

R> wy = cumsum(rnorm(N+1,mean=0,sd=sqrt(Dt)))

R> Ito.sumx <- c(0,sapply(1:(N+1),function(x){

+ exp(-a*(time[x+1]-time[x]))*(wx[x+1]-wx[x])}))

R> Ito.sumy <- c(0,sapply(1:(N+1),function(x){

+ exp(-a*(time[x+1]-time[x]))*(wy[x+1]-wy[x])}))

R> X <- sapply(1:(N+1),function(x){

+ x0*exp(-a*time[x])+sqrt(2*D)*sum(Ito.sumx[1:x])})

R> Y <- sapply(1:(N+1),function(x){

+ y0*exp(-a*time[x])+sqrt(2*D)*sum(Ito.sumy[1:x])})

R> plot(X,Y,type="l",las=1,xlab=expression(X[t]),ylab=expression(Y[t])) R> mtext("Langevin Process In 2D",line=2.5,cex=1.2 )

R> data.frame(X,Y)

Code 8

R> N = 1000; t0 = 0; Dt = 0.001; theta1 = 0.1; theta2 = 0.2;theta3=0.05; R> x0 = y0 = 1;

R> Error1 <- (2*theta1 > (theta3)^{^}2)

R> time <- c(t0 ,t0+ cumsum(rep(Dt,N)))

R> w = cumsum(rnorm(N+1,mean=0,sd=sqrt(Dt)))

R> Dw <- diff(w)

R> X <- Y <- numeric()

R> X[1] = Y[1] <- x0

R> for (i in 2:(N+1)){

+ X[i] = X[i-1] + (( theta1 - theta2 * X[i-1])-0.25* (theta3)^{^}2) * Dt +

+ theta3 * sqrt(X[i-1])*Dw[i-1] +0.25 * theta3 *(Dw[i-1])^{^}2

+ Y[i] = Y[i-1] + ((theta1 - theta2 * (Y[i-1])^{^}2)-0.25*(theta3)^{^}2)*Dt +

+ theta3*Y[i-1]*Dw[i-1]+0.25*theta3*(Dw[i-1])^{^}2

+ }

R> plot(time,X,type="l",col="blue")

R> lines(time,Y,type="l",col="red")

R> mtext(bquote(dX[t]==(.(theta1)-.(theta2)*X[t])*dt+.(theta3)*sqrt(X[t])* + dW[t]),line=2.5,cex=1,adj=0)

R> mtext(bquote(dY[t]==frac(1,2*Y[t])*bgroup("(",(.(theta1)-.(theta2)*

+ Y[t]^{^}2)-frac(1,4)*.(theta3)^{^}2 ,")") *dt+frac(1,2)*.(theta3)*dW[t]),

+ line=0.1,cex=1,adj=0)

R> legend("topright",border="gray",c(expression(X[t]),expression(Y[t])),

+ lty=c(1,1),col=c("blue","red"),lwd=2)

R> Error2 <- sum(abs(X-Y))/N

R> Error1 [1] TRUE R> Error2 [1] 0.0009837406

Annexe B : Packages Sim.DiffProc &

Sim.DiffProcGUI

L

'objective de l'annexe B est de données quelques règles pour la création des packages sous langage R. Et aussi une présentation des deux packages Sim.DiffProc et Sim.DiffProcGUI. Création d'un package

Pour créer un package sous langage R [32], il vous faut installer sur votre ordinateur un certain nombre de logiciels, tous sont disponibles gratuitement sur le web, puis les configurer.

· Perl : est un langage optimisée pour l'extraction d'informations de fichiers textes et la génération de rapports.

http://www.activestate.com/Products/ActivePerl/Download.html

· Rtools : Les Rtools sont des outils Unix qui fournissent une couche d'émulation pour le système Windows. Ils rendent possible l'exécution de programmes Unix sous Windows.

http://www.murdoch-sutherland.com/Rtools/tools.zip

· MinGW : permet de compiler du code C, C++ et FORTRAN. Si vous n'incluez pas de langage autre dans votre code, vous n'avez pas besoin de MinGW.

Sinon : http://prdownloads.sourceforge.net/mingw/MinGW-5.0.0.exe

· HTML Help Workshop: Ce logiciel permet de produire des aides au format .chm, le format d'aide propriétaire de Windows. http://msdn.microsoft.com/library/en-us/ htmlhelp/html/hwmicrosofthtmlhelpdownloads.asp

· Un compilateur LATEX : permet de produire une aide au format pdf. Plusieurs compilateurs sont disponibles, MiKTeX est assez stable. http://www.miktex.org/

· FTP : quand vous aurez terminé votre package, il vous faudra le poster sur le site du CRAN. Cela se fait grâce à un logiciel gérant le FTP (File Transfert Protocol). Là encore, plusieurs choix sont possibles. http://filezilla-project.org/

La création d'un package se fait via une fenêtre de commande DOS ou une fenêtre terminal sous Linux. Tous les nouveaux programmes qui viennent d'être installés doivent être accessibles depuis cette fenêtre. Pour cela, il faut préciser l'endroit où ils ont été installés sur votre ordinateur. Cela se fait en modifiant la variable PATH.

Le plus simple est de créer un fichier Rpath.bat contenant la définition de la variable PATH. Á noter que les espaces ne sont pas acceptés dans les noms de fichier. Sous Windows, Programme File doit donc être abrégé en PROGRA~1. La séparation des différents chemins peut se faire grâce à (;). Pour plus de lisibilité, il est possible de définir le PATH sur plusieurs lignes :

		Annexe B : Packages Sim.DiffProc & Sim.DiffProcGUI
SET PATH =C:\ Rtools \bin SET PATH =% PATH %;C:\ Perl \bin
SET PATH =% PATH	%;C:\	Rtools \ MinGW \bin
SET PATH =% PATH	%;C:\	PROGRA~1\R\R -2.13.1\	bin
SET PATH =% PATH	%;C:\	PROGRA~1\R\R -2.13.1\	include
SET PATH =% PATH	%;C:\	PROGRA~1\ MIKTEX~ 2.9\	miktex	\bin
SET PATH =% PATH	%;C:\	PROGRA~1\ HTMLHE~1
SET PATH =% PATH	%;C:\	WINDOWS
SET PATH =% PATH	%;C:\	WINDOWS \ system32

Si vous sauvegardez Rpath.bat dans le répertoire racine C:/, il vous suffit ensuite de taper C:/Rpath dans la fenêtre système que vous venez d'ouvrir et votre variable PATH est modifiée comme il convient. Il est également possible de modifier la variable PATH en allant explorer les variables d'environnement. Mais ces modifications sont permanentes jusqu'à ce qu'elles soient rectifiées. Cela veut dire qu'à chaque fois que vous allumerez votre ordinateur, le PATH sera modifiée même si vous n'avez pas de package à compiler ce jour-là. Votre ordinateur sera un peu moins rapide. Aussi, il est plus intéressant de créer un fichier Rpath.bat que l'on exécutera les jours où c'est nécessaire.

Un package est un ensemble de plusieurs fichiers et répertoires, tous réunis dans un répertoire racine. Le répertoire racine a pour nom le futur nom du package (par exemple Monpackage). Il contient un fichier nommé DESCRIPTION, un fichier nommé NAMESPACE, plus les répertoires /R/, /man/, /data/ et /tests/.

DESCRIPTION : contient une description du package.

NAMESPACE : sert à définir la visibilité de nos fonctions et des fonctions des autres packages, via import et export.

/R/ : contient le code des programmes.

/man/ : contient les fichiers d'aide.

/data/ : contient les jeux de données.

/tests/ : contient les fichiers permettant de tester notre programme (tests selon la R Core Team¹)

package.skeleton est une fonction qui crée pour vous l'arborescence des fichiers (répertoire racine, répertoires /R/, /man/ et /data/). Elle crée aussi les fichiers d'aide (dans /man/), les fichiers codes (dans /R/), éventuellement les données (dans /data/), le fichier DESCRIPTION et éventuellement le fichier NAMESPACE. L'idée est d'utiliser package.skeleton pour la création initiale de votre package, puis de modifier ensuite à la main les fichiers qu'il a crées.

R> help(package.skeleton)

Lorsque les fichiers sources (programme, aides et données) sont prêts, il reste a les vérifier, à construire le package et éventuellement à construire des documentations au format pdf. Ces trois étapes se font dans une fenêtre de commandes DOS sur Windows ou une fenêtre terminal sur Linux. Pour plus de détails [32, 33].

1. La R Core Team a mis au point un système de gestion de tests automatiques très performant.

Présentation du package

Le package Sim.DiffProc² [9] sous la version de Windows de langage R à été développée pour simulée et traité statistiquement des équations différentielles stochastiques d'une façon générale à l'exception les processus de diffusion [8]. Le langage R n'incorpore pas une interface graphique GUI ³, mais il inclus des outils pour construire des interfaces graphiques. Basé sur le package tcltk⁴ [13, 14, 41], qui fournit la possibilité de crée une interface à la boîte à outils de Tcl/Tk. Le package Sim.DiffProcGUI⁵ [10] fournit une interface graphique pour les fonctions qui sont dans le package Sim.DiffProc.

Installation du package

L'installation⁶ du package Sim.DiffProcGUI sous R est simple, il suffira de écrire dans R Console la commande suivante :

R> install.packages("Sim.DiffProcGUI")

Pour le chargement de package après l'installation et pour des informations sur le package, est décrit par les commandes suivantes :

R> library("Sim.DiffProcGUI")

Le chargement a nécessité le package : Sim.DiffProc Le chargement a nécessité le package : tcltk Chargement de Tcl/Tk... terminé

Le chargement a nécessité le package : tcltk2 Le chargement a nécessité le package : stats4 Le chargement a nécessité le package : rgl

Le chargement a nécessité le package : xlsx

Le chargement a nécessité le package : xlsxjars Le chargement a nécessité le package : rJava

Sim.DiffProcGUI version 2.0-Sat Feb 26 15:58:10 2011. BOUKHETALA Kamal < [e-mail]>.

GUIDOUM Arsalane < [e-mail]>.

University of Sciences and Technology Houari Boumediene (USTHB)

Faculty of Mathematics

Department of Probabilities and Statistics. Copyright (C) 2011 Algeria.

User Interface :: Sim.DiffGUI().

R> library(help="Sim.DiffProc")

R> library(help="Sim.DiffProcGUI")

2. Sim.DiffProc : Simulation of Diffusion Processes.

3. GUI: Graphical User Interface.

4. tcltk : Tool Command Language.

5. Sim.DiffProcGUI : Graphical User Interface for Simulation of Diffusion Processes.

6. Il faut être connecté à Internet pour effectuer ces opérations, où téléchargé manuellement les deux packages [9, 10].

Information sur le package Sim.DiffProc :

Package: Sim.DiffProc

Type: Package

Title: Simulation of Diffusion Processes

Version: 2.0

Date: 2011-02-09

Author: BOUKHETALA Kamal < [e-mail]>,

GUIDOUM Arsalane < [e-mail]>

Maintainer: BOUKHETALA Kamal < [e-mail]>

Depends: R (>= 2.11.0), tcltk, tcltk2, stats4, rgl, xlsx

License: GPL (>= 2)

URL: http://www.r-project.org

Repository: CRAN

LazyLoad: yes

Packaged: 2011-02-12 17:35:17 UTC;

Date/Publication: 2011-02-13 16:11:13

Built: R 2.13.1; ; 2011-09-29 07:46:05 UTC; windows

Information sur le package Sim.DiffProcGUI :

Package: Sim.DiffProcGUI

Type: Package

Title: Graphical User Interface for Simulation of

Diffusion Processes

Version: 2.0

Date: 2011-02-26

Author: BOUKHETALA Kamal < [e-mail]>,

GUIDOUM Arsalane < [e-mail]>

Maintainer: BOUKHETALA Kamal < [e-mail]>

Depends: Sim.DiffProc (>= 2.0)

License: GPL (>= 2)

LazyLoad: yes

Packaged: 2011-02-26 17:46:40 UTC;

Repository: CRAN

Date/Publication: 2011-02-27 16:28:40

Built: R 2.13.1; ; 2011-09-29 07:48:20 UTC; windows

Pour quelques démonstrations graphique de package Sim.DiffProc, est décrit par la commande suivante :

R> demo(Sim.DiffProc)

Les objectifs de conception de l'interface graphique est : renforcer le package Sim.DiffProc, pour une utilisation facile, une interface simple et familière de menu/dialogue-boîte. Le menu

principale est : File, Edit, Brownian Motion, Stochastic Integral, Stochastic Models, Parametric Estimation, Numerical Solution of SDE, Statistical Analysis, Help "?". Après chargement du package, la fenêtre GUI devrait apparaître plus ou moins comme dans la figure XIV.



FIGURE XIV - Graphical User Interface for Sim.DiffProc package at start-up.

Cette interface graphique d'écran de GUI à été crée sous Windows Seven. En cas de l'utilisation une autre version Windows, ou une autre plate-forme, l'image d'écran de GUI peut être différent ⁷. Toutes les fonctions du menu mènent aux zones de dialogue, c'est-à-dire interactif.

Le menu complet (tree) pour le package Sim.DiffProcGUI (version 2.0) est montré cidessous. toutes les fonctions de menu mènent aux zones de dialogue c'est-à-dire interactif.

File |- Open file

|- Change working directory

|- Import data from...

|- Save

|- Save graphic

|- Quit

| ||- Quit Sim.DiffProcGUI

||- Quit R

Edit |- Eval code Ctrl+R

7. Nous emploient la version R sous Windows, cependant, le langage R est disponible sur d'autres plates-formes, des ordinateurs de Macintosh et des systèmes d'Unix/Linux. L'utilisation de R et le package Sim.DiffProcGUI sur ces autres systèmes est très semblable à leur utilisation sous Windows.

|- Clear Ctrl+T

|- Undo Ctrl+Z

|- Redo Ctrl+W

|- Cut Ctrl+X

|- Copy Ctrl+C

|- Paste Ctrl+V

|- Select All Ctrl+A

|- Delete Ctrl+D

Brownian Motion| - Creating Brownian Motion

| ||- By the normal distribution

| ||- By a random walk

|- Creating flow for Brownian Motion

| ||- By the normal distribution

| ||- By a random walk

|- Creating Arithmetic Brownian Motion

| ||- Arithmetic brownian

| ||- Flow of arithmetic brownian

|- Creating Geometric Brownian Motion

| ||- Geometric brownian

| ||- Flow of geometric brownian

|- Creating Brownian Bridge

| ||- Brownian bridge

| ||- Flow of brownian bridge

|-Brownian Motion Property

| ||- Empirical covariance for brownian motion

| ||- Limite for brownian motion

| ||- Invariance by reversal of time

| ||- Invariance by scaling

|- Brownian Trajectory in 2D

| ||- By the normal distribution

| ||- By a random walk

|- Brownian Trajectory in 3D

||- By the normal distribution

||- By a random walk

Stochastic Integral|- Stratonovitch Integral - Integral(W(s) o dW(s),0,t)

| ||- Integral(alpha o dW(s),0,t)

| ||- Integral(W(s)^{^}n o dW(s),0,t)

|- Ito Integral - Integral(W(s)dW(s),0,t)

| ||- Integral(alpha*dW(s),0,t)

| ||- Integral(W(s)^{^}n*dW(s),0,t)

|- Ito Integral vs Stratonovitch Integral ||-Integral(W(s)dW(s),0,t) vs Integral(W(s) o dW(s),0,t) ||- Integral(s*dW(s),0,t) vs Integral(s o dW(s),0,t)

||- Integral(W(s)^{^}n*dW(s),0,t)vsIntegral(W(s)^{^}n o dW(s),0,t) Stochastic Models |- Attractive Model

| ||- Attractive Model for One-System SDE

| ||- Attractive Model for Two-System SDE

| |||- Two-Dimensional Attractive Model

| |||- Three-Dimensional Attractive Model

| |||- Simulation The First Passage Time

|- Bessel Process

|- Constant Elasticity of Variance Model

|- Cox-Ingersoll-Ross Model

| ||- CIR Model

| ||- The Modified CIR and hyperbolic Process

|- Chan-Karolyi-Longstaff-Sanders Model

|- Diffusion Bridge Model

|- Double-Well Potential Model

|- Exponential Martingales Process

|- Gaussian Diffusion Model

| ||- Hull-White/Vasicek (HWV)

| |||- Hull-White/Vasicek Model

| |||- Flow of Hull-White/Vasicek Model

| ||- Ornstein-Uhlenbeck Process

| |||- Ornstein-Uhlenbeck Process

| |||- Flow of Ornstein-Uhlenbeck Process

| |||- Radial Ornstein-Uhlenbeck Process

|- Hyperbolic Process

| ||- Hyperbolic Process

| ||- General Hyperbolic Diffusion

|- Inverse of Feller Square Root Model

|- Jacobi Diffusion Process

|- Pearson Diffusions Process

|- Stochastic Process

||- Stochastic Process The Gamma Distribution ||- Stochastic Process The Student Distribution ||- Random Walk

||- White Noise Gaussian

Parametric Estimation |- Parametric Estimation of Arithmetic Brownian Motion |- Parametric Estimation of Model Black-Scholes

|- Parametric Estimation of Hull-White/Vasicek Model |- Parametric Estimation of Ornstein-Uhlenbeck Model

||- Exact likelihood Inference

||- Explicit Estimators

Numerical Solution of SDE |- One-Dimensional SDE

| ||- Euler Scheme

| ||- Predictor-Corrector Method

| ||- Milstein Scheme

| |||- Milstein Scheme

| |||- Milstein Second Scheme

| ||- Strong Ito-Taylor Scheme

| ||- Heun Scheme

| ||- Runge-Kutta Scheme

|- Two-Dimensional SDE

||- Euler Scheme

||- Predictor-Corrector Method ||- Milstein Scheme

|||- Milstein Scheme

|||- Milstein Second Scheme ||- Strong Ito-Taylor Scheme

||- Heun Scheme

||- Runge-Kutta Scheme

Statistical Analysis |- Simulation M-Samples of Random Variable |- Simulation The First Passage Time FPT |- Ploting

| ||- Histograms

| ||- Kernel Density

| ||- Empirical Distribution

|- Adjustment of Distributions | ||- Beta Distribution

| |||- Estimate of The Parameters

| |||- Kolmogorov-Smirnov Tests

| ||- Chi-Squared Distribution

| |||- Estimate of The Parameters

| |||- Kolmogorov-Smirnov Tests

| ||- Exponential Distribution

| |||- Estimate of The Parameters

| |||- Kolmogorov-Smirnov Tests

| ||- Fisher Distribution

| |||- Estimate of The Parameters

| |||- Kolmogorov-Smirnov Tests

| ||- Gamma Distribution

| |||- Estimate of The Parameters

| |||- Kolmogorov-Smirnov Tests

| ||- Log Normal Distribution

| |||- Estimate of The Parameters

| |||- Kolmogorov-Smirnov Tests

| ||- Normal Distribution

| |||- Estimate of The Parameters

| |||- Kolmogorov-Smirnov Tests

| ||- Student Distribution

| |||- Estimate of The Parameters

| |||- Kolmogorov-Smirnov Tests

| ||- Weibull Distribution

| |||- Estimate of The Parameters

| |||- Kolmogorov-Smirnov Tests

|- Density Estimation

| ||- by Histograms Methods

| ||- by Kernel Methods

|- Distribution Estimation

Help "?" |- Demos

| ||- Attractive Model for One-System SDE (2D)

| ||- Attractive Model for One-System SDE (3D)

| ||- Attractive Model for Two-System SDE (2D)

| ||- Attractive Model for Two-System SDE (3D)

| ||- Brownian Motion in 2D plane

| ||- Flow of the Diffusion Processes

| ||- Numerical Simulation of One-Dimensional SDE

| ||- Numerical Simulation of Two-Dimensional SDE

| ||- Stochastic Processes (Models)

|- Teachware

|- Start Sim.DiffProc help (.HTML) |- Start Sim.DiffProc help (.PDF)

|- Start R help system

|- About Sim.DiffProc

L'interface de Sim.DiffProc inclut quelques éléments en plus des menus et des dialogues, les expositions externes qui sont montre par la figure XIV, on a :

Plot your data : c'est pour trace votre donnée en fonction du temps (time series).

Add your data : c'est pour ajoutée votre donnée dans un graphe (par exemple pour la comparé).

Plot your solution : c'est pour tracé la solution d'une ÉDS si il existe (par exemple la solution de Black-Scholes).

Add your solution : ajouté la solution dans un graphe (par exemple simulée un flux de modèle de Black-Scholes et ajouté la solution trouvé).

Show data : les donnée importé sont sauvegardée automatiquement dans cette fenêtre. Save graphic: sauvegardée les graphes sous différent format(.pdf, .png, .jpeg, ...)

Bibliographie

[1] Allen E (2007). Modeling with Itô Stochastic Differential Equations. Springer-Verlag, New York. ISBN 13 : 978-1-4020-5953-7.

[2] Boukhetala K (1994). Simulation Study of a Dispersion About an Attractive Center. In proceeding of 11th Symposium Computational Statistics., pp. 128-130, Edited by R. Dutter and W. Grossman, Wien, Austria.

[3] Boukhetala K (1995). Identification and simulation of a communication system. Maghreb Mathematical Review, 2, pp. 55-79.

[4] Boukhetala K (1996). Computer Methods and Water Resources, Modelling and Simulation of a Dispersion Pollutant with Attractive Centre, volume 3, pp. 245-252, Computational mechanics publications edition. Boston, USA.

[5] Boukhetala K (1998). Application des Processus de Diffusion, Echantillonnage Optimal. Ph.D. thesis, University of Science and Technology Houari Boumediene, BP 32 El-Alia, U.S.T.H.B, Algeria.

[6] Boukhetala K (1998). Estimation of the first passage time distribution for a simulated diffusion process. Maghreb Mathematical Review, 7, pp. 1-25.

[7] Boukhetala K (1998). Kernel density of the exit time in a simulated diffusion. The Annals of The Engineer Maghrebian, 12, pp. 587-589.

[8] Boukhetala K, Guidoum A (2011). Sim.DiffProc : A Package for Simulation of Diffusion Processes in R. Le Hal revue du centre national de la recherche scientifique (France). Preprint submitted to Journal of Statistical Software (JSS), 25 May 2011. http://hal.

archives-ouvertes.fr/hal-00629841/fr/.

[9] Boukhetala K, Guidoum A (2011). Sim.DiffProc : Simulation of Diffusion Processes. R package version 2.0. http://CRAN.R-project.org/package=Sim.DiffProc.

[10] Boukhetala K, Guidoum A (2011). Sim.DiffProcGUI : Graphical User Interface for Simulation of Diffusion Processes. R package version 2.0. http://CRAN.R-project.org/

package=Sim.DiffProcGUI.

[11] Brown R (1828). A brief account of microscopical observations made in the months of June, July and August, 1827, on the particles contained in the pollen of plants; and on the general existence of active molecules in organic and inorganic bodies. Philosophical Magazine, 4,

pp. 161-173. http://sciweb.nybg.org/science2/pdfs/dws/Brownian.pdf.

[12] Brian J F (1992). Brownian movement in Clarkia pollen: a reprise of the first observations. The Microscope, 40(4), pp. 235-241.

[13] Dalgaard P (2001). A Primer on the R-Tcl/Tk Package. R News, 1(3), pp. 27-31.

[14] Dalgaard P (2002). Changes to the R-Tcl/Tk Package. R News, 2(3), pp. 25-71.

[15] Deheuvels P (2006). Karhunen-Loève expansions of mean-centered Wiener processes. High Dimensional Probability, 51, pp. 62-76.

[16] Deheuvels P, G. Peccati G, Yor M (2006). On quadratic functionals of the brownian sheet and related processes. Stochastic Processes and their Applications, 116, pp. 493-538.

[17] Doob J L (1942). What is a stochastic process?. The American Mathematical Monthly, 49, pp. 648-653.

[18] Douglas H, Peter P (2006). Stochastic Differential Equations in Science and Engineering. World Scientific Publishing. ISBN 981-256-296-6.

[19] Flandrin P, Borgant P, Amblard P O (2003). From Stationarity to Self-similarity, and Back: Variations on the Lamperti Transformation. in Processes with Long-Range Correlations, pp. 88-117. Springer-Verlag

[20] Forman J L, Sørensen M (2007). The Pearson Diffusions : A Class of Statistically Tractable Diffusion Processes. SSRN : Social Science Research Network, http://ssrn.com/

abstract=1150110.

[21] Franck J (2009). Modèles aléatoires et physique probabiliste. Springer-Verlag, New York. ISBN 13 : 978-2-287-99307-7.

[22] Greiner A, Strittmatter W, Honerkamp J (1988). Numerical Integration of Stochastic Differential Equations. The Journal of Statistical Physics, 51(1/2).

[23] Hadeler K, de Mottoni P, Schumacher K (1980). Dynamic Models for Animal Orientation. The Journal of Mathematical Biology, 10, pp. 307-332.

[24] Heemink A (1990). Stochastic Modelling of Dispersion in Shallow Water. Stochastic Hydrology and Hydraulics, 4, pp. 161-174.

[25] Heinz S (2011). Mathematical Modeling. Stochastic Evolution, pp. 295-334, SpringerVerlag, Berlin Heidelberg. ISBN 978-3-642-20310-7.

[26] Helland S (1983). Diffusion models for the dispersal of insects near an attractive center. The Journal of Mathematical Biology, 18, pp. 103-122.

[27] Itô K (1944). Stochastic integral. Tokyo, Proc. Jap. Acad, 20, pp. 519-529.

[28] Kolmogorov A N (1936). Math. Sbornik. N.S., 1, pp. 607-610.

[29] Kloeden P, Platen E (1989). A Survey of Numerical Methods for Stochastic Differential Equations. Stochastic Hydrology and Hydraulics, 3, pp. 155-178.

[30] Lamperti J (1962). Semi-stable stochastic processes. Transactions of the American Mathematical Society, 104, pp. 62-78. http://www.jstor.org/stable/1993933.

[31] Peter E, Eckhard P (1995). Numerical Solution of Stochastic Differential Equations. Springer-Verlag, New York. ISBN 0-387-54062-8.

[32] R Development Core Team (2011). R : A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0,

http://www.R-project.org/.

[33] R Development Core Team (2011). Writing R extensions. version 2.13.1 (2011-07-08), ISBN 3-900051-11-9.

[34] Racicot F É, Théoret R (2006). Finance Computationnelle et Gestion des Risques. Presses de l'Université du Québec. ISBN 2-7605-1447-1.

[35] Risken H (2001). The Fokker-Planck Equation : Methods of Solutions and Applications. 2nd edition, Springer Series in Synergetics. ISBN 9783540615309.

[36] Rolski T, Schmidli H, Schmidt V, Teugels J (1998). Stochastic Processes for Insurance and Finance. John Wiley & Sons.

[37] Saito Y, Mitsui T (1993). Simulation of Stochastic Differential Equations. The Annals of the Institute of Statistical Mathematics, 3, pp. 419-432.

[38] Silverman B W (1986). Density estimation for statistics and data analysis. Chapman and Hall, London.

[39] Soong T T (1973). Random differential equations in science and engineering. Academic Press, New York. LC NUMBER : QA274.23 .S58.

[40] Stefano M (2008). Simulation and Inference for Stochastic Differential Equations. Springer-Verlag, New York. ISBN 978-0-387-75838-1.

[41] Welch B (2000). Practical Programming in Tcl and Tk. 3nd edition. Prentice Hall PTR Upper Saddle River, NJ, USA. ISBN 0-13-022028-0.