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

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).