source(file='coef.R') # simulation d'une EDS à sauts # deux fonctions: ajustement à (presque) tous les sauts. # ou simulation en une fois de tous les sauts sur l'intervalle Delta # si la simulation de chaque saut est impossible. # Fonction permettant de simuler une EDS à sauts. La dernière fonction # modifie les données pour # - les trier. # - ne garder que les valeurs entre -alpha et alpha # - calculer (X2-X1)^2/Delta # - calculer g(X1) # - à modifier si on veut estimer autre chose/changer l'intervalle d'estimation. #n: nombre de simulations. #Delta (multiple de delta) #delta: précision du schéma d'Euler. #ndrift,nsigma,nxi,nsauts: paramètres du modèle (drift, sigma, xi, type sauts) #lambda: intensité des sauts (ici toujours égale à 1) #beta: paramètre des sauts (cas nsauts=stable) #alpha: intervalle d'estimation #renvoie: #X #Y: approximation #T: calcul exact de g(X). euler<-function(n,Delta,delta,ndrift,nsigma,nxi,nsauts,lambda,beta,alpha) { ## Précision du modèle ###################################################################### ## Simulation de tous les sauts avec leur intensité. if (nsauts=='bin'|nsauts=='uni'|nsauts=='gau'|nsauts=='exp'|nsauts=='stable') { #simulation des temps de sauts if (nsauts=='stable') { e<- delta; ee<-e^(-beta); l<- (2-beta)/beta*(ee-1) # si nsauts='stable': beaucoup de petits sauts à simuler (on ne simule #que les sauts plus grands qu'un certain seuil), #leur nombre dépend de Delta (il faut augmenter la précision quand #Delta est petit.) lambda n'intervient pas. } else l<-lambda; # sinon, l'intensité du PP est donnée par lambda. N<-rpois(1,(n+2)*Delta*l);# nombre de sauts P<-sort(runif(N,0,(n+2)*Delta)); # temps des sauts # simulation du Brownien et des sauts en une fois. t<-seq(0,(n+3)*Delta,delta) # temps où on veut observer la simulation. # on simule à l'avance le brownien et les sauts NN<-length(t)+length(P) W<-rnorm(NN,0,1) # simulation MB (variance 1). if (nsauts=='stable') { Ws<-rnorm(NN,0,e^(1-beta/2))} # Si nsauts=stable, les # petits sauts sont remplacés par un MB de petite variance. S<-c(0,rsauts(length(P),nsauts,lambda,ee,beta))# simulation des sauts. # ee et beta n'interviennent que pour nsauts=stable, lambda pour tout le reste. # il faut maintenant bien ordonner le vecteur t. On note quand on a un saut # dans le vecteur I. TT<-sort(c(t,P),index.return=TRUE) # vecteur où il faut simuler le processus I<-(TT$ix>length(t)); # indique si on a un saut T<-TT$x; X<-seq(0,NN-1); k=1; # compteur de sauts #Simulation proprement dite. for (i in 2:NN) { if (I[i]==1) k<-k+1; # cas avec saut d<-T[i]-T[i-1]; if (nsauts=='stable') {X[i]<-X[i-1]+d*derive(X[i-1],ndrift)+sqrt(d)*sigma(X[i-1],nsigma)*W[i]+ sigma(X[i-1],nxi)*(S[k]*I[i]+sqrt(d)*Ws[i])} else #cas sans saut {X[i]<-X[i-1]+d*derive(V[i-1],ndrift)+sqrt(d)*sigma(X[i-1],nsigma)*W[i]+ sigma(X[i-1],nxi)*S[k]*I[i]} # à chaque instant( saut ou sim demandée), on recalcule Xi: Xi-1+dérive+MB+sauts (+petits sauts(cas stable)) } X<-X[I==0];# on ne garde que les points où on voulait la simulation. } ###################################################################### #deuxième fonction: simulation des sauts en une fois. if (nsauts=='gamma'|nsauts=='nig') { # simulation des sauts. if (nsauts=='gamma') {C<- rgamma((n+3)*Delta/delta,shape=delta,scale=1 ); S<-rnorm((n+3)*Delta/delta,0,sqrt(C))} if (nsauts=='nig') { beta<-alpha*delta; Y<-rnorm((n+3)*Delta/delta,0,1); N<-rnorm((n+3)*Delta/delta,0,1)^2; NN<-runif((n+3)*Delta/delta,0,1); W<-delta+delta^2*N/(2*beta^2)-delta/(2*beta^2)*sqrt(4*delta*beta^2*N+delta^2*N^2); Z<- W*(NN<=(delta/(delta+W)))+delta^2/W*(NN>(delta/(delta+W))); S<-sqrt(Z)*Y; } #simulation du MB. W<-rnorm((n+3)*Delta/delta,0,sqrt(delta)); # Simulation de l'EDS à sauts. NN=(n+3)*Delta/delta; X<-seq(0,NN); for (i in 2:NN) { X[i]=X[i-1]+derive(X[i-1],ndrift)*delta+W[i]*sigma(X[i-1],nsigma) +S[i]*sigma(X[i-1],nxi) } } #on ne garde que les valeurs pour Delta. J=seq(1,length(X),Delta/delta); X=X[J]; return(X); } trans<-function(X,alpha) { X1<-X[1:(n+2)]; X2<-X[2:(n+3)]; Y<-(X2-X1)^2/Delta; #Estimation de g Z<-Y*(abs(X2-X1)<(Delta^(1/2))); #estimation de sigma^2 T<-(X2-X1)^4/Delta; #Estimation de xi^4 V<-sigma(X1,nxi); # on ne conserve que les valeurs sur [-alpha,alpha]. J<- (abs(X1)