# 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)]
}

# 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)
              }

# computation of the critical value of the supremum test 
# using niter simulations of length  n
n<- 300
niter<- 5000
LMsup<-rep(0,niter)
set.seed(5)
for (i in 1:niter) { 
     y <- rnorm(n)
     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)
#     Wsup[i]<-n*(sig2contraint-sig2min)/sig2min
#     LR<-n*log(sig2contraint/sig2min)
     LMsup[i]<-n*(sig2contraint-sig2min)/sig2contraint
     }
critval.LMsup<-quantile(LMsup,0.95)

# Computation of the test statistics under the null
# 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("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
     res.h0<-as.vector(lm(y[2:n]~y[1:(n-1)])$residuals)
     sig2contraint<-var(res.h0)
# supremum test statistics
     W<-n*(sig2contraint-sig2min)/sig2min
     LM<-n*(sig2contraint-sig2min)/sig2contraint
     LR<-n*log(sig2contraint/sig2min)
# pointwise test statistics for r=0
     sig2min.0r<-mseSETAR(0,y)
     W0<-n*(sig2contraint-sig2min.0r)/sig2min.0r
     LM0<-n*(sig2contraint-sig2min.0r)/sig2contraint
     LR0<-n*log(sig2contraint/sig2min.0r)
# pointwise test statistics  for r=-1
     sig2min.1r<-mseSETAR(-1,y)
     W1<-n*(sig2contraint-sig2min.1r)/sig2min.1r
     LM1<-n*(sig2contraint-sig2min.1r)/sig2contraint
     LR1<-n*log(sig2contraint/sig2min.1r)
# pointwise test statistics  for r=-2
     sig2min.2r<-mseSETAR(-2,y)
     W2<-n*(sig2contraint-sig2min.2r)/sig2min.2r
     LM2<-n*(sig2contraint-sig2min.2r)/sig2contraint
     LR2<-n*log(sig2contraint/sig2min.2r)
# pointwise test statistics  for r=1
     sig2min.3r<-mseSETAR(1,y)
     W3<-n*(sig2contraint-sig2min.3r)/sig2min.3r
     LM3<-n*(sig2contraint-sig2min.3r)/sig2contraint
     LR3<-n*log(sig2contraint/sig2min.3r)
#1-LST
auxi.v1<-y[1:(n-1)]
auxi.v2<-auxi.v1^2
regauxi<- lm(res.h0~auxi.v1+auxi.v2)
sig2.1.LST<-var(regauxi$residuals)
W.1.LST<-n*(sig2contraint-sig2.1.LST)/sig2.1.LST
LM.1.LST<-n*(sig2contraint-sig2.1.LST)/sig2contraint
LR.1.LST<-n*log(sig2contraint/sig2.1.LST)
#3-LST
auxi.v1<-y[1:(n-1)]
auxi.v2<-auxi.v1^2
auxi.v3<-auxi.v1^3
auxi.v4<-auxi.v1^4
regauxi<- lm(res.h0~auxi.v1+auxi.v2+auxi.v3+auxi.v4)
sig2.3.LST<-var(regauxi$residuals)
W.3.LST<-n*(sig2contraint-sig2.3.LST)/sig2.3.LST
LM.3.LST<-n*(sig2contraint-sig2.3.LST)/sig2contraint
LR.3.LST<-n*log(sig2contraint/sig2.3.LST)


 write(c(W,LM,LR,W0,LM0,LR0,W1,LM1,LR1,W2,LM2,LR2,
         W3,LM3,LR3,W.1.LST,LM.1.LST,LR.1.LST,W.3.LST,LM.3.LST,LR.3.LST),
       file="testH0.dat",ncolumns=21,append=TRUE)
     }

# Computation of the test statistics under the alternative
a01<- 0; a11<-0.9; a02<- -2; a12<- 0.2; r0<- -2; 
set.seed(7)
file.remove("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)
# supremum test statistics
     W<-n*(sig2contraint-sig2min)/sig2min
     LM<-n*(sig2contraint-sig2min)/sig2contraint
     LR<-n*log(sig2contraint/sig2min)
# pointwise test statistics for r=0
     sig2min.0r<-mseSETAR(0,y)
     W0<-n*(sig2contraint-sig2min.0r)/sig2min.0r
     LM0<-n*(sig2contraint-sig2min.0r)/sig2contraint
     LR0<-n*log(sig2contraint/sig2min.0r)
# pointwise test statistics  for r=-1
     sig2min.1r<-mseSETAR(-1,y)
     W1<-n*(sig2contraint-sig2min.1r)/sig2min.1r
     LM1<-n*(sig2contraint-sig2min.1r)/sig2contraint
     LR1<-n*log(sig2contraint/sig2min.1r)
# pointwise test statistics  for r=-2
     sig2min.2r<-mseSETAR(-2,y)
     W2<-n*(sig2contraint-sig2min.2r)/sig2min.2r
     LM2<-n*(sig2contraint-sig2min.2r)/sig2contraint
     LR2<-n*log(sig2contraint/sig2min.2r)
# pointwise test statistics  for r=1
     sig2min.3r<-mseSETAR(1,y)
     W3<-n*(sig2contraint-sig2min.3r)/sig2min.3r
     LM3<-n*(sig2contraint-sig2min.3r)/sig2contraint
     LR3<-n*log(sig2contraint/sig2min.3r)
#1-LST
auxi.v1<-y[1:(n-1)]
auxi.v2<-auxi.v1^2
regauxi<- lm(res.h0~auxi.v1+auxi.v2)
sig2.1.LST<-var(regauxi$residuals)
W.1.LST<-n*(sig2contraint-sig2.1.LST)/sig2.1.LST
LM.1.LST<-n*(sig2contraint-sig2.1.LST)/sig2contraint
LR.1.LST<-n*log(sig2contraint/sig2.1.LST)
#3-LST
auxi.v1<-y[1:(n-1)]
auxi.v2<-auxi.v1^2
auxi.v3<-auxi.v1^3
auxi.v4<-auxi.v1^4
regauxi<- lm(res.h0~auxi.v1+auxi.v2+auxi.v3+auxi.v4)
sig2.3.LST<-var(regauxi$residuals)
W.3.LST<-n*(sig2contraint-sig2.3.LST)/sig2.3.LST
LM.3.LST<-n*(sig2contraint-sig2.3.LST)/sig2contraint
LR.3.LST<-n*log(sig2contraint/sig2.3.LST)


 write(c(W,LM,LR,W0,LM0,LR0,W1,LM1,LR1,W2,LM2,LR2,
         W3,LM3,LR3,W.1.LST,LM.1.LST,LR.1.LST,W.3.LST,LM.3.LST,LR.3.LST),
       file="testH1.dat",ncolumns=21,append=TRUE)
     }

# relative frequency of rejection of the different tests under H0
lesstat <- read.table("testH0.dat")
critchi1<-qchisq(0.95,df=1)
critchi2<-qchisq(0.95,df=2)
critchi3<-qchisq(0.95,df=3)
colnames(lesstat)<-c("W","LM","LR","W0","LM0","LR0",
                  "W1","LM1","LR1","W2","LM2","LR2","W3","LM3","LR3",
"W1LST","LM1LST","LR1LST","W3LST","LM3LST","LR3LST")
100*c(length(which(lesstat$W > critval.LMsup))/niter,
length(which(lesstat$W3 > critchi2))/niter,
length(which(lesstat$W0 > critchi2))/niter,
length(which(lesstat$W1 > critchi2))/niter,
length(which(lesstat$W2 > critchi2))/niter,
length(which(lesstat$W1LST > critchi1))/niter,
length(which(lesstat$W3LST > critchi3))/niter,
length(which(lesstat$LM > critval.LMsup))/niter,
length(which(lesstat$LM3 > critchi2))/niter,
length(which(lesstat$LM0 > critchi2))/niter,
length(which(lesstat$LM1 > critchi2))/niter,
length(which(lesstat$LM2 > critchi2))/niter,
length(which(lesstat$LM1LST > critchi1))/niter,
length(which(lesstat$LM3LST > critchi3))/niter,
length(which(lesstat$LR > critval.LMsup))/niter,
length(which(lesstat$LR3 > critchi2))/niter,
length(which(lesstat$LR0 > critchi2))/niter,
length(which(lesstat$LR1 > critchi2))/niter,
length(which(lesstat$LR2 > critchi2))/niter,
length(which(lesstat$LR1LST > critchi1))/niter,
length(which(lesstat$LR3LST > critchi3))/niter)

# relative frequency of rejection of the different tests under H1
lesstat <- read.table("testH1.dat")
colnames(lesstat)<-c("W","LM","LR","W0","LM0","LR0",
                  "W1","LM1","LR1","W2","LM2","LR2","W3","LM3","LR3",
"W1LST","LM1LST","LR1LST","W3LST","LM3LST","LR3LST")
100*c(length(which(lesstat$W > critval.LMsup))/niter,
length(which(lesstat$W3 > critchi2))/niter,
length(which(lesstat$W0 > critchi2))/niter,
length(which(lesstat$W1 > critchi2))/niter,
length(which(lesstat$W2 > critchi2))/niter,
length(which(lesstat$W1LST > critchi1))/niter,
length(which(lesstat$W3LST > critchi3))/niter,
length(which(lesstat$LM > critval.LMsup))/niter,
length(which(lesstat$LM3 > critchi2))/niter,
length(which(lesstat$LM0 > critchi2))/niter,
length(which(lesstat$LM1 > critchi2))/niter,
length(which(lesstat$LM2 > critchi2))/niter,
length(which(lesstat$LM1LST > critchi1))/niter,
length(which(lesstat$LM3LST > critchi3))/niter,
length(which(lesstat$LR > critval.LMsup))/niter,
length(which(lesstat$LR3 > critchi2))/niter,
length(which(lesstat$LR0 > critchi2))/niter,
length(which(lesstat$LR1 > critchi2))/niter,
length(which(lesstat$LR2 > critchi2))/niter,
length(which(lesstat$LR1LST > critchi1))/niter,
length(which(lesstat$LR3LST > critchi3))/niter)