# bruit blanc melange de 3 normales bruit<- function (n, m0=0, sig0=0.5, m1=3, sig1=1., p=0.05) { eps <- rep(0,n) for (i in 1:n) { u <- runif(1) u if ( u < p ) {eps[i] <- rnorm(1,mean=m1,sd=sig1)} if ( (p <= u) & (u < 2*p) ) {eps[i] <- rnorm(1,mean=-m1,sd=sig1)} if ( u > 2*p ) {eps[i] <- rnorm(1,mean=m0,sd=sig0)} } as.ts(eps) } # fonction qui simule un modele EXPAR d'ordre 1 expar.sim <- function (n, a0=0, b0=0, g0=0, m0=0, sig0=0.1, m1=3, sig1=0.1, p=0.05) { y <- rep(0,n) eps <- bruit(n) for (i in 2:n) { y[i] <- eps[i]+(a0+b0*exp(-g0*y[i-1]^2))*y[i-1] } as.ts(y) } # objective fonction pour estimer un modele EXPAR d'ordre 1 sur la serie y objf <- function(x, y) { a <- x[1] b <- x[2] g <- abs(x[3]) mse <- 0 grd <- c(0,0,0) n <- length(y) for (i in 2:n) { mse <- mse +((y[i]-(a+b*exp(-abs(g)*y[i-1]^2))*y[i-1])^2)/n grd <- grd + (2./n)*(y[i]-(a+b*exp(-abs(g)*y[i-1]^2))*y[i-1])* c(-y[i-1],-y[i-1]*exp(-abs(g)*y[i-1]^2),b*(y[i-1]^3)*exp(-abs(g)*y[i-1]^2)) g<- abs(g) } attr(mse, "gradient") <- grd mse } # niter simulation de tailles n n<- 500 niter<- 100 set.seed(7) for (i in 1:niter) { y <- expar.sim(n, a0=-0.8, b0=2., g0=0.1) ipass <- 0 icode <- 0 while( ipass <5 & icode==0) { ipass <- ipass+1 res <- nlm(f=objf, p=c(-0.8+0.1*(runif(1)-0.5),2.+0.1*(runif(1)-0.5),0.1+0.05*(runif(1)-0.5)), y=y) if(res$code <3) icode=res$code} param <- res$estimate param[3]<-abs(param[3]) if ( i == 1) write(c(param,ipass,res$code),file="expar.dat",ncolumns=5) else write(c(param,ipass,res$code),file="expar.dat",ncolumns=5,append=TRUE) } # les estimations coefsexpar <- read.table("expar.dat") colnames(coefsexpar)<-c("a","b","g","ipass","code") op <- par(mfrow = c(3, 2)) # 3 x 2 pictures on one plot # tracer la regression set.seed(7) y <- expar.sim(n=500, a0= -0.8, b0=2., g0=0.1) plot(ts(y[c(1:100)]),ylab=expression(y[t]),xlab="Time t") title(expression(paste('(a) Simulation (',Y[t],') of an EXPAR'))) abs<- seq(-6,6,len=100) a0 <- -0.8 b0 <- 2.0 g0 <- 0.1 n <- 500 ord<-(a0+b0*exp(-g0*abs^2))*abs plot(abs,ord,type='l', ylab=expression(paste('E(',Y[t],'|',Y[t-1],'=y)')), xlab='y') x<-y[c(2:(n-1))] z<-y[c(3:n)] points(x,z,pch=20) title(expression(paste(' (b) Points (',Y[t-1],',',Y[t],') and regression function'))) coefsexpar$a <- -0.8-coefsexpar$a coefsexpar$b <- 2.-coefsexpar$b coefsexpar$g <- 0.1-coefsexpar$g boxplot(list(coefsexpar$a,coefsexpar$b,coefsexpar$g), names=c(expression(a-hat(a)),expression(b-hat(b)),expression(gamma-hat(gamma))),xlab = expression(paste('estimation errors for n=500'))) title(expression(paste(' (c) Parameter estimations'))) qqnorm(coefsexpar$a,main="", xlab='Normal quantiles', ylab=expression(a-hat(a))) qqline(coefsexpar$a) title(expression(paste(' (d) Normal Q-Q Plot of ',hat(a)))) qqnorm(coefsexpar$b,main="", xlab="", ylab=expression(b-hat(b))) qqline(coefsexpar$b) title(expression(paste(' (e) Normal Q-Q Plot of ',hat(b)))) qqnorm(coefsexpar$g,main="",xlab="", ylab=expression(gamma-hat(gamma))) qqline(coefsexpar$g) title(expression(paste(' (f) Normal Q-Q Plot of ',hat(gamma)))) ## At end of plotting, reset to previous settings: par(op)