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