# fonction qui simule un modele SETAR tres simple
setar.sim <- function (n, a01=0, a11=0.9, a02=-2, a12=-0.7, r0=-2, burnin=100) 
{
eps <- rnorm(n+burnin)
y <- eps
for (i in 2:(n+burnin)) { 
if(y[i-1]<=r0) y[i] <- y[i]+a01+a11*y[i-1] else y[i] <- y[i]+a02+a12*y[i-1]}
y[(burnin+1):(n+burnin)]
}
set.seed(7)
n<- 500
sim<-setar.sim(n)
length(which(sim < -2))
plot(ts(setar.sim(n)))

# parametres du modele SETAR ajustes sur la serie y quand r est fixe
paraSETAR <- function(r, y)
     {
         n <- length(y)
         dep<-y[2:n]
         exp<-y[1:(n-1)]
         regim1<-ifelse(exp< r,1,0)
         regim2<-1-regim1
         exp1<-exp*regim1
         exp2<-exp*regim2
         lm(dep~regim1+exp1+regim2+exp2-1)$coefficients
              }

# moyenne des carres des residus du modele SETAR ajustes sur la serie y quand r est fixe
mseSETAR <- function(r, y)
     {
         n <- length(y)
         dep <- y[2:n]
         exp <- y[1:(n-1)]
         regim1 <- ifelse(exp< r,1,0)
         regim2 <- 1-regim1
         exp1 <- exp*regim1
         exp2 <- exp*regim2
         var(lm(dep~regim1+exp1+regim2+exp2-1)$residuals)
              }

# Sous H0 niter simulation de tailles n
a01<- 0; a11<- 0.7; a02<- 0; a12<- 0.7; r0<- -2; n<- 100
niter<- 1000
set.seed(7)
file.remove("paraSETARH0.dat","testH0.dat")
for (i in 1:niter) { 
     y <- setar.sim(n,a01=a01,a11=a11,a02=a02,a12=a12,r0=r0)
     low<-as.numeric(quantile(y,0.15))
     upp<-as.numeric(quantile(y,0.85))
     res<- optimize(mseSETAR, c(low,upp),y=y)
     sig2min<-res$objective
     sig2contraint<-var(lm(y[2:n]~y[1:(n-1)])$residuals)
     param<-as.numeric(c(paraSETAR(res$minimum,y),res$minimum))
#     obj<-lapply(1:100,function (i) mseSETAR(low+i*(upp-low)/100,y))
#     plot(as.numeric(obj))
     write(param,file="paraSETARH0.dat",ncolumns=5,append=TRUE)
     W<-n*(sig2contraint-sig2min)/sig2min
     LM<-n*(sig2contraint-sig2min)/sig2contraint
     LR<-n*log(sig2contraint/sig2min)
     sig2min.1r<-mseSETAR(-1,y)
     W1<-n*(sig2contraint-sig2min.1r)/sig2min.1r
     LM1<-n*(sig2contraint-sig2min.1r)/sig2contraint
     LR1<-n*log(sig2contraint/sig2min.1r)
     write(c(W,LM,LR,W1,LM1,LR1),file="testH0.dat",ncolumns=6,append=TRUE)
     }
# Sous H1 niter simulation de tailles n
a01<- 0; a11<-0.9; a02<- -2; a12<- 0.2; r0<- -2
set.seed(7)
file.remove("paraSETARH1.dat","testH1.dat")
for (i in 1:niter) { 
     y <- setar.sim(n,a01=a01,a11=a11,a02=a02,a12=a12,r0=r0)
     low<-as.numeric(quantile(y,0.15))
     upp<-as.numeric(quantile(y,0.85))
     res<- optimize(mseSETAR, c(low,upp),y=y)
     sig2min<-res$objective
     sig2contraint<-var(lm(y[2:n]~y[1:(n-1)])$residuals)
     param<-as.numeric(c(paraSETAR(res$minimum,y),res$minimum))
#     obj<-lapply(1:100,function (i) mseSETAR(low+i*(upp-low)/100,y))
#     plot(as.numeric(obj))
     write(param,file="paraSETARH1.dat",ncolumns=5,append=TRUE)
     W<-n*(sig2contraint-sig2min)/sig2min
     LM<-n*(sig2contraint-sig2min)/sig2contraint
     LR<-n*log(sig2contraint/sig2min)
     sig2min.1r<-mseSETAR(-1,y)
     W1<-n*(sig2contraint-sig2min.1r)/sig2min.1r
     LM1<-n*(sig2contraint-sig2min.1r)/sig2contraint
     LR1<-n*log(sig2contraint/sig2min.1r)
     write(c(W,LM,LR,W1,LM1,LR1),file="testH1.dat",ncolumns=6,append=TRUE)
     }

plot(function (x) dchisq(x,df=2),0,max(lesstat$W,lesstat$W0),xlab="",ylab="",
main=expression(paste('Distribution of ',W[n],' and ',W[n](-1),' under ',H[0])),cex.main=0.9,
lwd=1,lty=1)

dst<-density(lesstat$W)
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)}
max(lesstat$W)
abs<- seq(0,max(lesstat$W),len=500)
ord<- abs
for ( i in 1:length(abs)) ord[i]<-densite.W.mod(abs[i])
lines(abs,ord, col='red', 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)}
max(lesstat$W1)
abs<- seq(0,max(lesstat$W1),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)

legend(max(lesstat$W)/2,0.3,c(expression(chi[2]),expression(W[n](1)),expression(W[n])),
col=c("black","blue","red"),lwd=c(1,2,2),lty=c(1,3,2))


op <- par(mfrow = c(1, 2),  cex=0.7) # 1 x 2 pictures on one plot

lesstat <- read.table("testH0.dat")
colnames(lesstat)<-c("W","LM","LR","W1","LM1","LR1")
boxplot(list(lesstat$W,lesstat$LM,lesstat$LR,lesstat$W1,
lesstat$LM1,lesstat$LR1),
names=c(expression(W[n]),expression(LM[n]),expression(LR[n]),
expression(W[n](-1)),expression(LM[n](-1)),expression(LR[n](-1))))
title(expression(paste('Test statistics under ',H[0])))

lesstat <- read.table("testH1.dat")
colnames(lesstat)<-c("W","LM","LR","W1","LM1","LR1")
boxplot(list(lesstat$W,lesstat$LM,lesstat$LR,lesstat$W1,
lesstat$LM1,lesstat$LR1),
names=c(expression(W[n]),expression(LM[n]),expression(LR[n]),
expression(W[n](-1)),expression(LM[n](-1)),expression(LR[n](-1))))
title(expression(paste('Test statistics under ',H[1])))


     par(op)