# simulation d'un GARCH(p,q) 
simARGARCH <- function(n,omega,alpha,beta,ninit=100){
eta<-rnorm(n+ninit)
ht<-as.numeric(n+ninit)
epsilon<-as.numeric(n+ninit)
p<-length(beta)
q<-length(alpha)
epsilon[1:max(p,q)]<-sqrt(omega)*eta[1:max(p,q)]
ht[1:max(p,q)]<-epsilon[1:max(p,q)]^2
for (t in (max(p,q)+1):(n+ninit)) {
 ht[t]<-omega+sum(alpha[1:q]*(epsilon[(t-1):(t-q)]^2))+sum(beta[1:p]*ht[(t-1):(t-p)])
 epsilon[t]<-sqrt(ht[t])*eta[t]}
 epsilon[(ninit+1):(n+ninit)]
} 
# estime un GARCH(1,1) par QMLE 

objf.garch11<- function(x, eps){
         omega <- x[1]
         alpha <- x[2]
         beta <- x[3]
         n <- length(eps)
         sigma2<-as.numeric(n)
         sigma2[1]<-eps[1]^2
for (t in 2:n) sigma2[t]<-omega+alpha*eps[t-1]^2+beta*sigma2[t-1]
         qml <- mean(log(sigma2[2:n])+eps[2:n]**2/sigma2[2:n])
         qml }

garch11<- function(omega0,alpha0,beta0,eps0)
     {
valinit<-c(omega0,alpha0,beta0)
res <- nlminb(valinit,objf.garch11, lower=c(0.00000001,0,0),upper=c(Inf,1,1), eps=eps0)
res$par
}

# estime un GARCH(1,1) par QMLE 

objf.qml <- function(x, eps,q){
         omega <- x[1]
         alpha <- x[2]
         beta <- x[3]
         n <- length(eps)
         sigma2<-as.numeric(n)
         sigma2[1:q]<-eps[1:q]^2
for (t in (q+1):n) sigma2[t]<-omega+alpha*eps[t-1]^2+beta*sigma2[t-1]
         qml <- mean(log(sigma2[(q+1):n])+eps[(q+1):n]**2/sigma2[(q+1):n])
         qml }

estimgarch11<- function(omega0,alpha0,beta0,eps0,q)
     {
valinit<-c(omega0,alpha0,beta0)
res <- nlminb(valinit,objf.qml, lower=c(0.00000001,0,0),upper=c(Inf,1,1), eps=eps0,q=q)
list(theta=res$par,ln=res$objective)
}

#########################################
#          GARCH(1,q)
########################################

# estime J et kappa pour un GARCH(1,q)
MatJ <- function(omega,alpha,beta,eps){
         q<-length(alpha)
n <- length(eps)
sigma2<-as.numeric(n)
dersigma2<-matrix(nrow=(q+2),ncol=n)
score<-as.vector(rep(0,(q+2)))
Jmat<-matrix(0,nrow=(q+2),ncol=(q+2))
sigma2[1:q]<-eps[1:q]^2
dersigma2[1:(q+2),1:q]<-0
for (t in (q+1):n){
sigma2[t]<-omega+sum(alpha[1:q]*(eps[(t-1):(t-q)]^2))+beta*sigma2[t-1]
dersigma2[1:(q+2),t]<- c(1,eps[(t-1):(t-q)]^2,sigma2[t-1])+beta*dersigma2[1:(q+2),(t-1)]
score<-score+(1-eps[t]^2/sigma2[t])*dersigma2[1:(q+2),t]/sigma2[t]
Jmat<-Jmat+(dersigma2[1:(q+2),t]/sigma2[t])%*%t(dersigma2[1:(q+2),t]/sigma2[t])
} 
score<-score/(n-q)
Jmat<-Jmat/(n-q)
kappa<-mean((eps[(q+1):n]^4)/(sigma2[(q+1):n]^2))
list(score=score,Jmat=Jmat,kappa=kappa)
}

# estime un GARCH(1,q)
objfgarch1q.qml <- function(x, eps){
         q<-length(x)-2
         omega <- x[1]
         alpha <- x[2:(q+1)]
         beta <- x[q+2]
         n <- length(eps)
         sigma2<-as.numeric(n)
         sigma2[1:q]<-eps[1:q]^2
for (t in (q+1):n) sigma2[t]<-omega+sum(alpha[1:q]*(eps[(t-1):(t-q)]^2))+beta*sigma2[t-1]
         qml <- mean(log(sigma2[(q+1):n])+eps[(q+1):n]**2/sigma2[(q+1):n])
         qml }

estimgarch1q<- function(omega0,alpha0,beta0,eps0)
     {
valinit<-c(omega0,alpha0,beta0)
q<-length(alpha0)
res <- nlminb(valinit,objfgarch1q.qml, lower=c(0.00000001,rep(0,q),0),upper=c(Inf,rep(1,q),1), eps=eps0)
list(theta=res$par,ln=res$objective)
}

# solution lambda  de la minimisation 
# sous contrainte (dans theta_0 les elements de 3 à q+1 sont nuls)
objf.lambda <- function(x, Z, J){
obj<-t(x-Z)%*%J%*%(x-Z)
         obj }

lambda<- function(Z,J)
     {
q<-length(Z)-2
valinit<-c(Z[1:2],pmax(Z[3:(q+1)],0),Z[q+2]) 

res <- nlminb(valinit,objf.lambda, lower=c(-Inf,-Inf,rep(0,q-1),-Inf),upper=rep(Inf,q+2), Z=Z, J=J)
res$par
}

# pvalue estimée du test de Wald de l'hypothese de GARCH(1,1) contre GARCH(1,q)

simul.W <- function(i,Jmat,mu,Jsim,Omega) {
Z<-mvrnorm(mu=mu,Sigma=Jsim)
lamb<-lambda(Z,Jmat)
t(lamb)%*%Omega%*%lamb
}

#kappa<-res$kappa Jmat<-res$Jmat
est.pval <- function(Wn,LRn,kappa,Jmat,Jmatinv,Omega,N) {
q<-length(Jmat[1,])-2
mu<-rep(0,(q+2))
Jsim<-(kappa-1)*Jmatinv
vector.W<-sapply(1:N,simul.W,Jmat=Jmat,mu=mu,Jsim=Jsim,Omega=Omega)
W<-length(which(vector.W>Wn))/N
LR<-length(which(vector.W>LRn))/N
list(W=W,LR=LR)
}


# statistique de Wald pour tester un GARCH(1,1) contre un GARCH(1,q)
Wald.GARCH11GARCH1q <- function(omega0,alpha0,beta0,eps0,N=5000){
 q<-length(alpha0)
 n<-length(eps0)
 garch1q<-estimgarch1q(omega0,alpha0,beta0,eps0)
 theta<-garch1q$theta
 garch11<-estimgarch11(theta[1],theta[2],theta[q+2],eps0,q)
 thetagarch11<-garch11$theta
 LRn<-n*(garch11$ln-garch1q$ln)
 if(LRn<0){
 omegainit<-thetagarch11[1]
 alphainit<-c(thetagarch11[2],rep(0,q-1))
 betainit<-thetagarch11[3]
 garch1q<-estimgarch1q(omegainit,alphainit,betainit,eps0)
 theta<-garch1q$theta
 LRn<-n*(garch11$ln-garch1q$ln)
 }
 theta2<-as.vector(theta[3:(q+1)])
 K<-cbind(rep(0,(q-1)),rep(0,(q-1)),diag(rep(1,(q-1))),rep(0,(q-1)))
 res<-MatJ(theta[1],theta[2:(q+1)],theta[q+2],eps0)
 omegacont<-thetagarch11[1]
 alphacont<-c(thetagarch11[2],rep(0,q-1))
 betacont<-thetagarch11[3]
 cont<-MatJ(omegacont,alphacont,betacont,eps0)
# LRn<-2*LRn/(res$kappa-1)
 LRn<-2*LRn/(cont$kappa-1)
 Jcontinv<-try(solve(cont$Jmat),silent=TRUE)
 if(is.matrix(Jcontinv)){
  Rn<-as.numeric(n*t(cont$score)%*%Jcontinv%*%cont$score/(cont$kappa-1))}
 Jmatinv<-try(solve(res$Jmat),silent=TRUE)
 if(is.matrix(Jmatinv))dum<-K%*%Jmatinv%*%t(K)
 duminv<-try(solve(dum),silent=TRUE)
 if(is.matrix(duminv))Omega<-t(K)%*%duminv%*%K/(res$kappa-1)
 Wn<-0
 if(is.matrix(duminv))Wn<-as.numeric(n*t(theta2)%*%duminv%*%theta2/(res$kappa-1))
 if(is.matrix(Jmatinv))pval<-est.pval(Wn,LRn,res$kappa,res$Jmat,Jmatinv,Omega,N)
 pvalR<-1-pchisq(Rn,df=q-1)     
 list(theta=theta,Wn=Wn,LRn=LRn,Rn=Rn,pvalW=pval$W,pvalLR=pval$LR,pvalR=pvalR)
 }

#########################################
#          GARCH(p,1)
########################################

# estime J et kappa pour un GARCH(p,1)
# omega<-theta[1];alpha<-theta[2];beta<-theta[3:(p+2)]
# eps<-eps0
MatJ2 <- function(omega,alpha,beta,eps){
p<-length(beta)
n <- length(eps)
sigma2<-as.numeric(n)
dersigma2<-matrix(nrow=(p+2),ncol=n)
score<-as.vector(rep(0,(p+2)))
Jmat<-matrix(0,nrow=(p+2),ncol=(p+2))
sigma2[1:p]<-eps[1:p]^2
dersigma2[1:(p+2),1:p]<-0
for (t in (p+1):n){
sigma2[t]<-omega+alpha*eps[(t-1)]^2+sum(beta[1:p]*sigma2[(t-1):(t-p)])
dersigma2[1:(p+2),t]<-c(1,eps[(t-1)]^2,sigma2[(t-1):(t-p)])
for(j in 1:p){dersigma2[1:(p+2),t]<-dersigma2[1:(p+2),t]+
beta[j]*dersigma2[1:(p+2),(t-j)]}
score<-score+(1-eps[t]^2/sigma2[t])*dersigma2[1:(p+2),t]/sigma2[t]
Jmat<-Jmat+(dersigma2[1:(p+2),t]/sigma2[t])%*%t(dersigma2[1:(p+2),t]/sigma2[t])
} 
score<-score/(n-p)
Jmat<-Jmat/(n-p)
kappa<-mean((eps[(p+1):n]^4)/(sigma2[(p+1):n]^2))
list(score=score,Jmat=Jmat,kappa=kappa)
}

# estime un GARCH(p,1)
objfgarchp1.qml <- function(x, eps){
         p<-length(x)-2
         omega <- x[1]
         alpha <- x[2]
         beta <- x[3:(p+2)]
         n <- length(eps)
         sigma2<-as.numeric(n)
         sigma2[1:p]<-eps[1:p]^2
for (t in (p+1):n) {sigma2[t]<-omega+alpha*eps[t-1]^2+
                     sum(beta[1:p]*sigma2[(t-1):(t-p)])}
         qml <- mean(log(sigma2[(p+1):n])+eps[(p+1):n]**2/sigma2[(p+1):n])
         qml }

estimgarchp1<- function(omega0,alpha0,beta0,eps0)
     {
valinit<-c(omega0,alpha0,beta0)
p<-length(beta0)
res <- nlminb(valinit,objfgarchp1.qml, 
  lower=c(0.00000001,rep(0,p+1)),upper=c(Inf,rep(1,p+1)), eps=eps0)
list(theta=res$par,ln=res$objective)
}


# solution lambda  de la minimisation sous contrainte (dans theta_0 les elements de 4 à p+2 sont nuls)

lambda2<- function(Z,J)
     {
p<-length(Z)-2
valinit<-c(Z[1:3],pmax(Z[4:(p+2)],0)) 
res <- nlminb(valinit,objf.lambda, lower=c(-Inf,-Inf,-Inf,rep(0,(p-1))),upper=rep(Inf,(p+2)), Z=Z, J=J)
res$par
}
# pvalue estimée du test de Wald de l'hypothese de GARCH(1,1) contre GARCH(p,1)

simul2.W <- function(i,Jmat,mu,Jsim,Omega) {
Z<-mvrnorm(mu=mu,Sigma=Jsim)
lamb<-lambda2(Z,Jmat)
t(lamb)%*%Omega%*%lamb
}
est.pval2 <- function(Wn,LRn,kappa,Jmat,Jmatinv,Omega,N) {
p<-length(Jmat[1,])-2
mu<-rep(0,(p+2))
Jsim<-(kappa-1)*Jmatinv
vector.W<-sapply(1:N,simul2.W,Jmat=Jmat,mu=mu,Jsim=Jsim,Omega=Omega)
W<-length(which(vector.W>Wn))/N
LR<-length(which(vector.W>LRn))/N
list(W=W,LR=LR)
}

# statistiques de Wald, score et LR pour tester un 
#  GARCH(1,1) contre un GARCH(p,1)
# omega0<-omegainit alpha0<-alphainit
# beta0<-betinit eps0<-ser

tests.GARCH11GARCHp1 <- function(omega0,alpha0,beta0,eps0,N=5000){
 p<-length(beta0)
 n<-length(eps0)
 garchp1<-estimgarchp1(omega0,alpha0,beta0,eps0)
 theta<-garchp1$theta
 garch11<-estimgarch11(theta[1],theta[2],theta[3],eps0,p)
 thetagarch11<-garch11$theta
 LRn<-n*(garch11$ln-garchp1$ln)
 if(LRn<0){
 omegainit<-thetagarch11[1]
 alphainit<-thetagarch11[2]
 betainit<-c(thetagarch11[3],rep(0,(p-1)))
 garchp1<-estimgarchp1(omegainit,alphainit,betainit,eps0)
 theta<-garchp1$theta
 LRn<-n*(garch11$ln-garchp1$ln)
 }
 theta2<-as.vector(theta[4:(p+2)])
 K<-cbind(rep(0,(p-1)),rep(0,(p-1)),rep(0,(p-1)),diag(rep(1,(p-1))))
 res<-MatJ2(theta[1],theta[2],theta[3:(p+2)],eps0)
 omegacont<-thetagarch11[1]
 alphacont<-thetagarch11[2]
 betacont<-c(thetagarch11[3],rep(0,p-1))
 cont<-MatJ2(omegacont,alphacont,betacont,eps0)
# LRn<-2*LRn/(res$kappa-1)
 LRn<-2*LRn/(cont$kappa-1)
 Jcontinv<-try(solve(cont$Jmat),silent=TRUE)
 if(is.matrix(Jcontinv)){
  Rn<-as.numeric(n*t(cont$score)%*%Jcontinv%*%cont$score/(cont$kappa-1))}
 Jmatinv<-try(solve(res$Jmat),silent=TRUE)
 if(is.matrix(Jmatinv))dum<-K%*%Jmatinv%*%t(K)
 duminv<-try(solve(dum),silent=TRUE)
 if(is.matrix(duminv))Omega<-t(K)%*%duminv%*%K/(res$kappa-1)
 Wn<-0
 if(is.matrix(duminv))Wn<-as.numeric(n*t(theta2)%*%duminv%*%theta2/(res$kappa-1))
 if(is.matrix(Jmatinv))pval<-est.pval2(Wn,LRn,res$kappa,res$Jmat,Jmatinv,Omega,N)
 pvalR<-1-pchisq(Rn,df=p-1)     
 list(theta=theta,Wn=Wn,LRn=LRn,Rn=Rn,pvalW=pval$W,pvalLR=pval$LR,pvalR=pvalR)
 }


##########################
library(tseries)
library(MASS)
alphainit<-0.05
betainit<-0.85
chouia<-0.01

tab.pval<-matrix('&',nrow=11,ncol=37) 
# lecture des donnees du cac
tab.pval[1,1]<-as.character('CAC')
cacdata <- read.table("cac.csv",header=TRUE,sep=",")
cac<-rev(cacdata$Close) 
datecac<-rev(cacdata$Date) 
ntot<-length(cac)
rend<-rep(0,ntot); rend[2:ntot]<-log(cac[2:ntot]/cac[1:(ntot-1)])*100 

ser<-rend[2:ntot]-mean(rend[2:ntot])
#ser<-rend[2:ntot]
n<-length(ser)
omegainit<-var(ser)/(1-alphainit-betainit)

res<-Wald.GARCH11GARCH1q(omegainit,cbind(alphainit,chouia),betainit,ser)
tab.pval[1,3]<-round(res$pvalW,3)
tab.pval[1,5]<-round(res$pvalR,3)
tab.pval[1,7]<-round(res$pvalLR,3)
tab.pval[1,9]<-round((1-pchisq(res$Wn,df=1))/2,3)
tab.pval[1,11]<-round(res$pvalR,3)
tab.pval[1,13]<-round((1-pchisq(res$LRn,df=1))/2,3)

for(q in 3:4){
res<-Wald.GARCH11GARCH1q(res$theta[1],c(res$theta[2:q],chouia),res$theta[q+1],ser)
tab.pval[1,15+6*(q-3)]<-round(res$pvalW,3)
tab.pval[1,17+6*(q-3)]<-round(res$pvalR,3)
tab.pval[1,19+6*(q-3)]<-round(res$pvalLR,3)}

res<- tests.GARCH11GARCHp1(omegainit,alphainit,c(betainit,chouia),ser)
tab.pval[1,27]<-round(res$pvalW,3)
tab.pval[1,29]<-round(res$pvalR,3)
tab.pval[1,31]<-round(res$pvalLR,3)
tab.pval[1,33]<-round((1-pchisq(res$Wn,df=1))/2,3)
tab.pval[1,35]<-round(res$pvalR,3)
tab.pval[1,37]<-round((1-pchisq(res$LRn,df=1))/2,3)

#print(tab.pval[1,],digits=2,quote=FALSE)


############################################################
##################### SP 500 ###############################

# lecture des donnees du sp 500
tab.pval[2,1]<-as.character('SP 500')
spdata <- read.table("sp500.csv",header=TRUE,sep=",")
sp<-rev(spdata$Close) 
datesp<-rev(spdata$Date) 
ntot<-length(sp)
rend<-rep(0,ntot); rend[2:ntot]<-log(sp[2:ntot]/sp[1:(ntot-1)])*100 

ser<-rend[2:ntot]-mean(rend[2:ntot])
#plot(ts(ser))
n<-length(ser)
omegainit<-var(ser)/(1-alphainit-betainit)

res<-Wald.GARCH11GARCH1q(omegainit,cbind(alphainit,chouia),betainit,ser)
tab.pval[2,3]<-round(res$pvalW,3)
tab.pval[2,5]<-round(res$pvalR,3)
tab.pval[2,7]<-round(res$pvalLR,3)
tab.pval[2,9]<-round((1-pchisq(res$Wn,df=1))/2,3)
tab.pval[2,11]<-round(res$pvalR,3)
tab.pval[2,13]<-round((1-pchisq(res$LRn,df=1))/2,3)

for(q in 3:4){
res<-Wald.GARCH11GARCH1q(res$theta[1],c(res$theta[2:q],chouia),res$theta[q+1],ser)
tab.pval[2,15+6*(q-3)]<-round(res$pvalW,3)
tab.pval[2,17+6*(q-3)]<-round(res$pvalR,3)
tab.pval[2,19+6*(q-3)]<-round(res$pvalLR,3)}

res<- tests.GARCH11GARCHp1(omegainit,alphainit,c(betainit,chouia),ser)
tab.pval[2,27]<-round(res$pvalW,3)
tab.pval[2,29]<-round(res$pvalR,3)
tab.pval[2,31]<-round(res$pvalLR,3)
tab.pval[2,33]<-round((1-pchisq(res$Wn,df=1))/2,3)
tab.pval[2,35]<-round(res$pvalR,3)
tab.pval[2,37]<-round((1-pchisq(res$LRn,df=1))/2,3)

#print(tab.pval[2,],digits=2,quote=FALSE)
#print(tab.pval,digits=2,quote=FALSE)



####################################################################
################### Dax (Germany)#################################

# lecture des donnees du DAX
tab.pval[3,1]<-as.character('DAX')
daxdata <- read.table("dax.csv",header=TRUE,sep=",")
dax<-rev(daxdata$Close) 
datedax<-rev(daxdata$Date) 
ntot<-length(dax)
rend<-rep(0,ntot); rend[2:ntot]<-log(dax[2:ntot]/dax[1:(ntot-1)])*100 

ser<-rend[2:ntot]-mean(rend[2:ntot])
n<-length(ser)
omegainit<-var(ser)/(1-alphainit-betainit)

res<-Wald.GARCH11GARCH1q(omegainit,cbind(alphainit,chouia),betainit,ser)
tab.pval[3,3]<-round(res$pvalW,3)
tab.pval[3,5]<-round(res$pvalR,3)
tab.pval[3,7]<-round(res$pvalLR,3)
tab.pval[3,9]<-round((1-pchisq(res$Wn,df=1))/2,3)
tab.pval[3,11]<-round(res$pvalR,3)
tab.pval[3,13]<-round((1-pchisq(res$LRn,df=1))/2,3)

for(q in 3:4){
res<-Wald.GARCH11GARCH1q(res$theta[1],c(res$theta[2:q],chouia),res$theta[q+1],ser)
tab.pval[3,15+6*(q-3)]<-round(res$pvalW,3)
tab.pval[3,17+6*(q-3)]<-round(res$pvalR,3)
tab.pval[3,19+6*(q-3)]<-round(res$pvalLR,3)}

res<- tests.GARCH11GARCHp1(omegainit,alphainit,c(betainit,chouia),ser)
tab.pval[3,27]<-round(res$pvalW,3)
tab.pval[3,29]<-round(res$pvalR,3)
tab.pval[3,31]<-round(res$pvalLR,3)
tab.pval[3,33]<-round((1-pchisq(res$Wn,df=1))/2,3)
tab.pval[3,35]<-round(res$pvalR,3)
tab.pval[3,37]<-round((1-pchisq(res$LRn,df=1))/2,3)

#print(tab.pval[3,],digits=2,quote=FALSE)

####################################################################
################### FTSE (London)#################################

# lecture des donnees du FTSE
tab.pval[4,1]<-as.character('FTSE')
ftsedata <- read.table("ftse.csv",header=TRUE,sep=",")
ftse<-rev(ftsedata$Close) 
dateftse<-rev(ftsedata$Date) 
ntot<-length(ftse)
rend<-rep(0,ntot); rend[2:ntot]<-log(ftse[2:ntot]/ftse[1:(ntot-1)])*100 

ser<-rend[2:ntot]-mean(rend[2:ntot])
#plot(ts(ser))
#acf(ser^2)
#pacf(ser^2)
n<-length(ser)
omegainit<-var(ser)/(1-alphainit-betainit)

res<-Wald.GARCH11GARCH1q(omegainit,cbind(alphainit,chouia),betainit,ser)
tab.pval[4,3]<-round(res$pvalW,3)
tab.pval[4,5]<-round(res$pvalR,3)
tab.pval[4,7]<-round(res$pvalLR,3)
tab.pval[4,9]<-round((1-pchisq(res$Wn,df=1))/2,3)
tab.pval[4,11]<-round(res$pvalR,3)
tab.pval[4,13]<-round((1-pchisq(res$LRn,df=1))/2,3)

for(q in 3:4){
res<-Wald.GARCH11GARCH1q(res$theta[1],c(res$theta[2:q],chouia),res$theta[q+1],ser)
tab.pval[4,15+6*(q-3)]<-round(res$pvalW,3)
tab.pval[4,17+6*(q-3)]<-round(res$pvalR,3)
tab.pval[4,19+6*(q-3)]<-round(res$pvalLR,3)}

res<- tests.GARCH11GARCHp1(omegainit,alphainit,c(betainit,chouia),ser)
tab.pval[4,27]<-round(res$pvalW,3)
tab.pval[4,29]<-round(res$pvalR,3)
tab.pval[4,31]<-round(res$pvalLR,3)
tab.pval[4,33]<-round((1-pchisq(res$Wn,df=1))/2,3)
tab.pval[4,35]<-round(res$pvalR,3)
tab.pval[4,37]<-round((1-pchisq(res$LRn,df=1))/2,3)

#print(tab.pval[4,],digits=2,quote=FALSE)

####################################################################
################### SMI (Swiss)#################################

# lecture des donnees du SMI
tab.pval[5,1]<-as.character('SMI')
smidata <- read.table("smi.csv",header=TRUE,sep=",")
smi<-rev(smidata$Close) 
datesmi<-rev(smidata$Date) 
ntot<-length(smi)
rend<-rep(0,ntot); rend[2:ntot]<-log(smi[2:ntot]/smi[1:(ntot-1)])*100 

ser<-rend[2:ntot]-mean(rend[2:ntot])
n<-length(ser)
omegainit<-var(ser)/(1-alphainit-betainit)

res<-Wald.GARCH11GARCH1q(omegainit,cbind(alphainit,chouia),betainit,ser)
tab.pval[5,3]<-round(res$pvalW,3)
tab.pval[5,5]<-round(res$pvalR,3)
tab.pval[5,7]<-round(res$pvalLR,3)
tab.pval[5,9]<-round((1-pchisq(res$Wn,df=1))/2,3)
tab.pval[5,11]<-round(res$pvalR,3)
tab.pval[5,13]<-round((1-pchisq(res$LRn,df=1))/2,3)

for(q in 3:4){
res<-Wald.GARCH11GARCH1q(res$theta[1],c(res$theta[2:q],chouia),res$theta[q+1],ser)
tab.pval[5,15+6*(q-3)]<-round(res$pvalW,3)
tab.pval[5,17+6*(q-3)]<-round(res$pvalR,3)
tab.pval[5,19+6*(q-3)]<-round(res$pvalLR,3)}

res<- tests.GARCH11GARCHp1(omegainit,alphainit,c(betainit,chouia),ser)
tab.pval[5,27]<-round(res$pvalW,3)
tab.pval[5,29]<-round(res$pvalR,3)
tab.pval[5,31]<-round(res$pvalLR,3)
tab.pval[5,33]<-round((1-pchisq(res$Wn,df=1))/2,3)
tab.pval[5,35]<-round(res$pvalR,3)
tab.pval[5,37]<-round((1-pchisq(res$LRn,df=1))/2,3)

#print(tab.pval[5,],digits=2,quote=FALSE)


####################################################################
################### NIKKEI (Japan)#################################

# lecture des donnees du NIKKEI
tab.pval[6,1]<-as.character('NIKKEI')
nikkeidata <- read.table("nikkei.csv",header=TRUE,sep=",")
nikkei<-rev(nikkeidata$Close) 
datenikkei<-rev(nikkeidata$Date) 
ntot<-length(nikkei)
rend<-rep(0,ntot); rend[2:ntot]<-log(nikkei[2:ntot]/nikkei[1:(ntot-1)])*100 

ser<-rend[2:ntot]-mean(rend[2:ntot])
n<-length(ser)
omegainit<-var(ser)/(1-alphainit-betainit)

res<-Wald.GARCH11GARCH1q(omegainit,cbind(alphainit,chouia),betainit,ser)
tab.pval[6,3]<-round(res$pvalW,3)
tab.pval[6,5]<-round(res$pvalR,3)
tab.pval[6,7]<-round(res$pvalLR,3)
tab.pval[6,9]<-round((1-pchisq(res$Wn,df=1))/2,3)
tab.pval[6,11]<-round(res$pvalR,3)
tab.pval[6,13]<-round((1-pchisq(res$LRn,df=1))/2,3)

for(q in 3:4){
res<-Wald.GARCH11GARCH1q(res$theta[1],c(res$theta[2:q],chouia),res$theta[q+1],ser)
tab.pval[6,15+6*(q-3)]<-round(res$pvalW,3)
tab.pval[6,17+6*(q-3)]<-round(res$pvalR,3)
tab.pval[6,19+6*(q-3)]<-round(res$pvalLR,3)}

res<- tests.GARCH11GARCHp1(omegainit,alphainit,c(betainit,chouia),ser)
tab.pval[6,27]<-round(res$pvalW,3)
tab.pval[6,29]<-round(res$pvalR,3)
tab.pval[6,31]<-round(res$pvalLR,3)
tab.pval[6,33]<-round((1-pchisq(res$Wn,df=1))/2,3)
tab.pval[6,35]<-round(res$pvalR,3)
tab.pval[6,37]<-round((1-pchisq(res$LRn,df=1))/2,3)

#print(tab.pval[6,],digits=2,quote=FALSE)

####################################################################
################### Dow Jones  (DJI)    ############################

# lecture des donnees du Dow Jones
tab.pval[7,1]<-as.character('DJI')
dowjonesdata <- read.table("dowjones.csv",header=TRUE,sep=",")
dowjones<-rev(dowjonesdata$Close) 
datedowjones<-rev(dowjonesdata$Date) 
ntot<-length(dowjones)
rend<-rep(0,ntot); rend[2:ntot]<-log(dowjones[2:ntot]/dowjones[1:(ntot-1)])*100 

ser<-rend[2:ntot]-mean(rend[2:ntot])
n<-length(ser)
omegainit<-var(ser)/(1-alphainit-betainit)

res<-Wald.GARCH11GARCH1q(omegainit,cbind(alphainit,chouia),betainit,ser)
tab.pval[7,3]<-round(res$pvalW,3)
tab.pval[7,5]<-round(res$pvalR,3)
tab.pval[7,7]<-round(res$pvalLR,3)
tab.pval[7,9]<-round((1-pchisq(res$Wn,df=1))/2,3)
tab.pval[7,11]<-round(res$pvalR,3)
tab.pval[7,13]<-round((1-pchisq(res$LRn,df=1))/2,3)

for(q in 3:4){
res<-Wald.GARCH11GARCH1q(res$theta[1],c(res$theta[2:q],chouia),res$theta[q+1],ser)
tab.pval[7,15+6*(q-3)]<-round(res$pvalW,3)
tab.pval[7,17+6*(q-3)]<-round(res$pvalR,3)
tab.pval[7,19+6*(q-3)]<-round(res$pvalLR,3)}

res<- tests.GARCH11GARCHp1(omegainit,alphainit,c(betainit,chouia),ser)
tab.pval[7,27]<-round(res$pvalW,3)
tab.pval[7,29]<-round(res$pvalR,3)
tab.pval[7,31]<-round(res$pvalLR,3)
tab.pval[7,33]<-round((1-pchisq(res$Wn,df=1))/2,3)
tab.pval[7,35]<-round(res$pvalR,3)
tab.pval[7,37]<-round((1-pchisq(res$LRn,df=1))/2,3)

#print(tab.pval[7,],digits=2,quote=FALSE)




####################################################################
################### Nasdaq      #################################

# lecture des donnees du Nasdaq
tab.pval[8,1]<-as.character('Nasdaq')
nasdaqdata <- read.table("nasdaq.csv",header=TRUE,sep=",")
nasdaq<-rev(nasdaqdata$Close) 
datenasdaq<-rev(nasdaqdata$Date) 
ntot<-length(nasdaq)
rend<-rep(0,ntot); rend[2:ntot]<-log(nasdaq[2:ntot]/nasdaq[1:(ntot-1)])*100 

ser<-rend[2:ntot]-mean(rend[2:ntot])
n<-length(ser)
omegainit<-var(ser)/(1-alphainit-betainit)

res<-Wald.GARCH11GARCH1q(omegainit,cbind(alphainit,chouia),betainit,ser)
tab.pval[8,3]<-round(res$pvalW,3)
tab.pval[8,5]<-round(res$pvalR,3)
tab.pval[8,7]<-round(res$pvalLR,3)
tab.pval[8,9]<-round((1-pchisq(res$Wn,df=1))/2,3)
tab.pval[8,11]<-round(res$pvalR,3)
tab.pval[8,13]<-round((1-pchisq(res$LRn,df=1))/2,3)

for(q in 3:4){
res<-Wald.GARCH11GARCH1q(res$theta[1],c(res$theta[2:q],chouia),res$theta[q+1],ser)
tab.pval[8,15+6*(q-3)]<-round(res$pvalW,3)
tab.pval[8,17+6*(q-3)]<-round(res$pvalR,3)
tab.pval[8,19+6*(q-3)]<-round(res$pvalLR,3)}

res<- tests.GARCH11GARCHp1(omegainit,alphainit,c(betainit,chouia),ser)
tab.pval[8,27]<-round(res$pvalW,3)
tab.pval[8,29]<-round(res$pvalR,3)
tab.pval[8,31]<-round(res$pvalLR,3)
tab.pval[8,33]<-round((1-pchisq(res$Wn,df=1))/2,3)
tab.pval[8,35]<-round(res$pvalR,3)
tab.pval[8,37]<-round((1-pchisq(res$LRn,df=1))/2,3)

#print(tab.pval[8,],digits=2,quote=FALSE)
#print(tab.pval,digits=2,quote=FALSE)



####################################################################
################### DJA (Dow Jones Composite Index)     ############

# lecture des donnees du Dow Jones Composite Index
tab.pval[9,1]<-as.character('DJA')
djadata <- read.table("dja.csv",header=TRUE,sep=",")
dja<-rev(djadata$Close) 
datedja<-rev(djadata$Date) 
ntot<-length(dja)
rend<-rep(0,ntot); rend[2:ntot]<-log(dja[2:ntot]/dja[1:(ntot-1)])*100 

ser<-rend[2:ntot]-mean(rend[2:ntot])
n<-length(ser)
omegainit<-var(ser)/(1-alphainit-betainit)

res<-Wald.GARCH11GARCH1q(omegainit,cbind(alphainit,chouia),betainit,ser)
tab.pval[9,3]<-round(res$pvalW,3)
tab.pval[9,5]<-round(res$pvalR,3)
tab.pval[9,7]<-round(res$pvalLR,3)
tab.pval[9,9]<-round((1-pchisq(res$Wn,df=1))/2,3)
tab.pval[9,11]<-round(res$pvalR,3)
tab.pval[9,13]<-round((1-pchisq(res$LRn,df=1))/2,3)

for(q in 3:4){
res<-Wald.GARCH11GARCH1q(res$theta[1],c(res$theta[2:q],chouia),res$theta[q+1],ser)
tab.pval[9,15+6*(q-3)]<-round(res$pvalW,3)
tab.pval[9,17+6*(q-3)]<-round(res$pvalR,3)
tab.pval[9,19+6*(q-3)]<-round(res$pvalLR,3)}

res<- tests.GARCH11GARCHp1(omegainit,alphainit,c(betainit,chouia),ser)
tab.pval[9,27]<-round(res$pvalW,3)
tab.pval[9,29]<-round(res$pvalR,3)
tab.pval[9,31]<-round(res$pvalLR,3)
tab.pval[9,33]<-round((1-pchisq(res$Wn,df=1))/2,3)
tab.pval[9,35]<-round(res$pvalR,3)
tab.pval[9,37]<-round((1-pchisq(res$LRn,df=1))/2,3)

#print(tab.pval[9,],digits=2,quote=FALSE)



####################################################################
################### DJT (Dow Jones Transportation Average)  ########

# lecture des donnees du Dow Jones Transportation Average
tab.pval[10,1]<-as.character('DJT')
djtdata <- read.table("djt.csv",header=TRUE,sep=",")
djt<-rev(djtdata$Close) 
datedjt<-rev(djtdata$Date) 
ntot<-length(djt)
rend<-rep(0,ntot); rend[2:ntot]<-log(djt[2:ntot]/djt[1:(ntot-1)])*100 

ser<-rend[2:ntot]-mean(rend[2:ntot])
n<-length(ser)
omegainit<-var(ser)/(1-alphainit-betainit)

res<-Wald.GARCH11GARCH1q(omegainit,cbind(alphainit,chouia),betainit,ser)
tab.pval[10,3]<-round(res$pvalW,3)
tab.pval[10,5]<-round(res$pvalR,3)
tab.pval[10,7]<-round(res$pvalLR,3)
tab.pval[10,9]<-round((1-pchisq(res$Wn,df=1))/2,3)
tab.pval[10,11]<-round(res$pvalR,3)
tab.pval[10,13]<-round((1-pchisq(res$LRn,df=1))/2,3)

for(q in 3:4){
res<-Wald.GARCH11GARCH1q(res$theta[1],c(res$theta[2:q],chouia),res$theta[q+1],ser)
tab.pval[10,15+6*(q-3)]<-round(res$pvalW,3)
tab.pval[10,17+6*(q-3)]<-round(res$pvalR,3)
tab.pval[10,19+6*(q-3)]<-round(res$pvalLR,3)}

res<- tests.GARCH11GARCHp1(omegainit,alphainit,c(betainit,chouia),ser)
tab.pval[10,27]<-round(res$pvalW,3)
tab.pval[10,29]<-round(res$pvalR,3)
tab.pval[10,31]<-round(res$pvalLR,3)
tab.pval[10,33]<-round((1-pchisq(res$Wn,df=1))/2,3)
tab.pval[10,35]<-round(res$pvalR,3)
tab.pval[10,37]<-round((1-pchisq(res$LRn,df=1))/2,3)

#print(tab.pval[10,],digits=2,quote=FALSE)


####################################################################
################### DJU (Dow Jones Utilities)  ########

# lecture des donnees du Dow Jones Utilities
tab.pval[11,1]<-as.character('DJU')
djudata <- read.table("dju.csv",header=TRUE,sep=",")
dju<-rev(djudata$Close) 
datedju<-rev(djudata$Date) 
ntot<-length(dju)
rend<-rep(0,ntot); rend[2:ntot]<-log(dju[2:ntot]/dju[1:(ntot-1)])*100 

ser<-rend[2:ntot]-mean(rend[2:ntot])
n<-length(ser)
omegainit<-var(ser)/(1-alphainit-betainit)

res<-Wald.GARCH11GARCH1q(omegainit,cbind(alphainit,chouia),betainit,ser)
tab.pval[11,3]<-round(res$pvalW,3)
tab.pval[11,5]<-round(res$pvalR,3)
tab.pval[11,7]<-round(res$pvalLR,3)
tab.pval[11,9]<-round((1-pchisq(res$Wn,df=1))/2,3)
tab.pval[11,11]<-round(res$pvalR,3)
tab.pval[11,13]<-round((1-pchisq(res$LRn,df=1))/2,3)

for(q in 3:4){
res<-Wald.GARCH11GARCH1q(res$theta[1],c(res$theta[2:q],chouia),res$theta[q+1],ser)
tab.pval[11,15+6*(q-3)]<-round(res$pvalW,3)
tab.pval[11,17+6*(q-3)]<-round(res$pvalR,3)
tab.pval[11,19+6*(q-3)]<-round(res$pvalLR,3)}

res<- tests.GARCH11GARCHp1(omegainit,alphainit,c(betainit,chouia),ser)
tab.pval[11,27]<-round(res$pvalW,3)
tab.pval[11,29]<-round(res$pvalR,3)
tab.pval[11,31]<-round(res$pvalLR,3)
tab.pval[11,33]<-round((1-pchisq(res$Wn,df=1))/2,3)
tab.pval[11,35]<-round(res$pvalR,3)
tab.pval[11,37]<-round((1-pchisq(res$LRn,df=1))/2,3)

#print(tab.pval[11,],digits=2,quote=FALSE)

print(as.table(tab.pval[,c(1:2,9:26,33:37)]),digits=2,quote=FALSE)