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

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