# bruit blanc gaussien
bruit<- function (n, m=0,  sig=1.)
{
eps <-  rnorm(n,mean=m,sd=sig) 
as.ts(eps)}
# fonction h(g,y)
h <- function (g,y){1-exp(-g*y^2)}
# fonction qui simule un modele EXPAR d'ordre 1
expar.sim <- function (n, a0=0,  b0=0, g0=0) 
{
y <- as.numeric(n+1)
eps <- bruit(n+1)

for (i in 2:(n+1)) { y[i] <- eps[i]+(a0+b0*h(g0,y[i-1]))*y[i-1] }
as.ts(y[2:(n+1)])
}
# fonction V.n(g)
V.n <- function (g,n0,y0){mean(y0[1:(n0-1)]^2*(h(g,y0[1:(n0-1)])+1))}
# fonction D.n(g)
D.n <- function (g,n0,y0){
a<-mean(y0[1:(n0-1)]^2)
b<-mean(  y0[1:(n0-1)]^2*(h(g,y0[1:(n0-1)]))^2)
c<-mean(y0[1:(n0-1)]^2*(h(g,y0[1:(n0-1)])))
a*b-c^2}
# processus gaussien de covariance K(g,g')
z.proc <- function (g,eps,n0,y0) { 
sum(eps[2:n0]*y0[1:(n0-1)]*(h(g,y0[1:(n0-1)])+1))/sqrt(n0-1) }

# processus W(g)
w.proc <- function (g,eps,n0,y0) 
{
a<-V.n(0,n0,y0)*z.proc(g,eps,n0,y0)
b<-V.n(g,n0,y0)*z.proc(0,eps,n0,y0)
c<-D.n(g,n0,y0)*mean(y0[1:(n0-1)]^2)
if(abs(c) != 0) ((a-b)^2)/c else 0}

# parametres du modele EXPAR ajustes sur la serie y quand g est fixe
paraEXPAR <- function(g, y)
     {
         n <- length(y)
         dep<-y[2:n]
         exp1<-y[1:(n-1)]
         exp2<-h(g,y[1:(n-1)])*y[1:(n-1)]
         lm(dep~exp1+exp2-1)$coefficients
              }

# moyenne des carres des residus du modele EXPAR ajustes sur la serie y quand g est fixe
mseEXPAR <- function(g, y)
     {
         n <- length(y)
         dep<-y[2:n]
         exp1<-y[1:(n-1)]
         exp2<-h(g,y[1:(n-1)])*y[1:(n-1)]
         var(lm(dep~exp1+exp2-1)$residuals)
              }

# Niter simulation de tailles n g<-0.2
n<- 100; a0<-0;  b0<-0; g0<-0
niter<- 5000
nmax<- 5*n      
#nmax<- n      
set.seed(7)
file.remove("paraEXPAR.dat","statEXPAR.dat")
     low<-0.
     upp<-100
for (i in 1:niter) { 
     y <- expar.sim(n,a0=a0,b0=b0,g0=g0)

###################### W ######################

     eps<-rnorm(nmax)
#len=1000
#abs<- seq(low,upp,len=len)
#ord<- as.numeric(len)
#for(jj in (1:len))ord[jj]<-w.proc(abs[jj],eps,nmax,eps)
#plot(c(min(abs),max(abs)),c(min(ord),max(ord)),type='n')
#lines(abs,ord)

     opt<- optimize(w.proc, interval=c(low,upp),
         n0=nmax,y0=eps,maximum=TRUE,eps=eps)
     Sup<-opt$objective
     sup1at<-opt$maximum
     continue<-1
     passage<-0
while (continue==1) {
     if(min(sup1at,upp-sup1at)>(upp-low)/20) inter<-sup1at else inter<-runif(1)*(upp-low)

     opt11<- optimize(w.proc, interval=c(low,(low+inter)/2),
         n0=nmax,y0=eps,maximum=TRUE,eps=eps)
     sup11<-opt11$objective
     opt12<- optimize(w.proc, interval=c((low+inter)/2,inter),
         n0=nmax,y0=eps,maximum=TRUE,eps=eps)
     sup12<-opt12$objective
     sup1<-max(sup11,sup12)
     if(sup11>sup12)sup1at<-opt11$maximum else sup1at<-opt12$maximum 

     opt21<- optimize(w.proc, interval=c(inter,(inter+upp)/2),
         n0=nmax,y0=eps,maximum=TRUE,eps=eps)
     sup21<-opt21$objective
     opt22<- optimize(w.proc, interval=c((inter+upp)/2,upp),
         n0=nmax,y0=eps,maximum=TRUE,eps=eps)
     sup22<-opt22$objective
     sup2<-max(sup21,sup22)
     if(sup2>sup1)if(sup21>sup22)sup1at<-opt21$maximum else sup1at<-opt22$maximum 

     passage<-passage+1
     if(passage>9)continue<-0
     if(max(sup1,sup2)<=Sup)continue<-0 else Sup<-max(sup1,sup2)
                     }
#len<-200
#abs<- seq(low,upp,len=len)
#ord<- as.numeric(len)
#for(jj in (1:len))ord[jj]<-w.proc(abs[jj],eps,nmax,eps)
#     Sup<-max(Sup1,ord)

###################### W_n(0.5)

     sig2min.0.5<-mseEXPAR(0.5,y)
     sig2contraint<-var(lm(y[2:n]~y[1:(n-1)]-1)$residuals)     
     W1<-n*(sig2contraint-sig2min.0.5)/sig2min.0.5

###################### W_n

     opt<- optimize(mseEXPAR, interval=c(low,upp),y=y)
     inf<-opt$objective
     inf1at<-opt$minimum
     continue<-1
     passage<-0
while (continue==1) {
     if(min(inf1at,upp-inf1at)>(upp-low)/20) inter<-inf1at else inter<-runif(1)*(upp-low)

     opt11<- optimize(mseEXPAR, interval=c(low,(low+inter)/2),y=y)
     inf11<- opt11$objective
     opt12<- optimize(mseEXPAR, interval=c((low+inter)/2,inter),y=y)
     inf12<- opt12$objective
     inf1<- min(inf11,inf12)
     if(inf11<inf12)inf1at<-opt11$minimum else inf1at<-opt12$minimum 

     opt21<- optimize(mseEXPAR, interval=c(inter,(inter+upp)/2),y=y)
     inf21<- opt21$objective
     opt22<- optimize(mseEXPAR, interval=c((inter+upp)/2,upp),y=y)
     inf22<-opt22$objective
     inf2<-min(inf21,inf22)
     if(inf2<inf1)if(inf21<inf22)inf1at<-opt21$minimum else inf1at<-opt22$minimum 

     passage<-passage+1
     if(passage>9)continue<-0
     if(min(inf1,inf2)>inf)continue<-0 else inf<-min(inf1,inf2)
                     }
     sig2min<-inf
     if(sig2min>sig2min.0.5)sig2min<-{
       min(optimize(mseEXPAR, interval=c(low,0.5),y=y)$objective,
           optimize(mseEXPAR, interval=c(0.5,upp),y=y)$objective)}
#     param<-as.numeric(c(paraEXPAR(res$minimum,y),res$minimum))
#     write(param,file="paraEXPAR.dat",ncolumns=5,append=TRUE)
     W<-n*(sig2contraint-sig2min)/sig2min
#     LM<-n*(sig2contraint-sig2min)/sig2contraint
#     LR<-n*log(sig2contraint/sig2min)
#     write(c(sig2min,sig2min.0.5,W,W1),file="paraEXPAR.dat",ncolumns=4,append=TRUE)
#     LM1<-n*(sig2contraint-sig2min.0.5)/sig2contraint
#     LR1<-n*log(sig2contraint/sig2min.0.5)
#     write(c(W,LM,LR,W1,LM1,LR1,Sup),file="statEXPAR.dat",ncolumns=7,append=TRUE)
     write(c(W,W1,Sup),file="statEXPAR.dat",ncolumns=7,append=TRUE)
     }
lesstat <- read.table("statEXPAR.dat")
colnames(lesstat)<-c("W","W1","Sup")

plotsup<-max(lesstat$W1,lesstat$W)
plotsup<-10
plot(c(0,plotsup),c(dchisq(0.2,df=1),dchisq(plotsup,df=1)),
type="n",xlab="",ylab="",
main=expression(paste('Distribution of ',bold(W),' and of ',
bold(W)[n],' and ',bold(W)[n](0.5),' under ',H[0],' for n=100')),
cex.main=0.9,
lwd=1,lty=1)

abs<- seq(0,plotsup,len=100)
ord<- abs
for ( i in 1:length(abs)) ord[i]<-dchisq(abs[i],df=1)
lines(abs,ord, lwd=1,lty=1)

#plot(function (x) dchisq(x,df=1),0,plotsup,xlab="",ylab="",
#main=expression(paste('Distribution of ',bold(W),' and of ',bold(W)[n],' and ',bold(W)[n](0.5),' under ',H[0])),cex.main=0.9,
#lwd=1,lty=1)

dst<-density(lesstat$W,adjust=1)
dst$bw

densite.W<- function (z) {
dst$y[which(abs(dst$x-z)==min(abs(dst$x-z)))]}

densite.W.mod<- function (z) {
densite.W(z)+densite.W(-z)}

abs<- seq(0,plotsup,len=100)
ord<- abs
for ( i in 1:length(abs)) ord[i]<-densite.W.mod(abs[i])
lines(abs,ord, col='blue', lwd=2,lty=2)

dst<-density(lesstat$W1)
densite.W1<- function (z) {
dst$y[which(abs(dst$x-z)==min(abs(dst$x-z)))]}
densite.W1.mod<- function (z) {
densite.W1(z)+densite.W1(-z)}
abs<- seq(0,plotsup,len=100)
ord<- abs
for ( i in 1:length(abs)) ord[i]<-densite.W1.mod(abs[i])
lines(abs,ord, col='blue', lwd=2,lty=3)

dst<-density(lesstat$Sup,adjust=1)
densite.Sup<- function (z) {
dst$y[which(abs(dst$x-z)==min(abs(dst$x-z)))]}
densite.Sup.mod<- function (z) {
densite.Sup(z)+densite.Sup(-z)}
abs<- seq(0,plotsup,len=100)
ord<- abs
for ( i in 1:length(abs)) ord[i]<-densite.Sup.mod(abs[i])
lines(abs,ord, col='blue', lwd=2,lty=1)

legend(plotsup/2,0.6,c(expression(chi[1]^2),expression(bold(W)[n](0.5)),expression(bold(W)[n]),expression(bold(W))),
col=c("black","black","blue","blue"),lwd=c(1,2,2,2),lty=c(1,3,2,1))


#lesstat1 <- read.table("statEXPAR.dat")
#lesstat2 <- read.table("essai3.dat")
#lesstat3 <- read.table("essai6.dat")
#lesstat1-lesstat2
#lesstat1-lesstat3