# Test l'existence d'un moment d'ordre 2 pour un GARCH(1,1)
# omega<-QMLE$coef[1]; alpha<-QMLE$coef[2]; beta<-QMLE$coef[3]; Sigma<-QMLE$Sigma; eps<-ser
garch.teststatiofaible<- function (omega,  alpha, beta, Sigma, eps) {
n<-length(eps)
lambda<-matrix(c(0,1,1),nrow=3,ncol=1)
den<-sqrt(as.numeric(t(lambda)%*%Sigma%*%lambda))
stat<-sqrt(n)*(alpha+beta-1)/den
pval<-pnorm(stat)
list(pval=pval,std=den/sqrt(n))
}

# calcul de la variances asymptique Sigma pour un GARCH(1,1)
#  parametre theta
# n<-10; eps<-rnorm(n); omega<-1; alpha<-0.1; beta<- 0.2
garch.etvar<- function (omega,  alpha, beta, eps) {
n<-length(eps)
gamma<-omega/(1-alpha-beta)
ht<-eps
ht[1]<-omega
S<-matrix(0,nrow=3,ncol=n)
D<-matrix(nrow=3,ncol=n)
D[,1]<-c(1,0,0)

for (t in 2:n) {
ht[t] <- omega+alpha*eps[t-1]**2+beta*ht[t-1]
D[,t]<-c(1,eps[t-1]**2,ht[t-1])+beta*D[,t-1]
S[,t]<-D[,t]/ht[t]}

eta<-eps/sqrt(ht)
mu4<-mean(eta^4)
cond<-(alpha+beta)^2+(mu4-1)*alpha^2

J<-S%*%t(S)/(n-1)
Jinv<-solve(J)
Sigma<-(mu4-1)*Jinv

list(Sigma=Sigma,mu4=mu4,cond=cond)
}

# simule un gararch(1,1)
garch.sim<- function (n, omega,  alpha, beta, valinit=100) {
gamma<-omega/(1-alpha-beta)
eta <- rnorm(n+valinit)
ht<-rep(0,n+valinit)
eps<-ht
for (t in 2:(n+valinit)) { 
ht[t] <- omega+alpha*eps[t-1]**2+beta*ht[t-1]
eps[t]<-sqrt(ht[t])*eta[t]}
return(eps[(valinit+1):(n+valinit)])}

# estime un GARCH(1,1) avec nlminb omega0<-1 alpha0<-0.05 beta<-0.7 eps0<-sim

objf3 <- function(x, eps){     
         omega <- x[1]
         alpha <- x[2]
         beta <- x[3]
         if(alpha+beta>0.9999)return(Inf)
         n <- length(eps)
         sigma2<-as.numeric(n)
         sigma2[1]<-omega
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 }

estimgarch11<- function(omega0,alpha0,beta0,eps0)
     {
valinit<-c(omega0,alpha0,beta0)
res <- nlminb(valinit,objf3, lower=c(0.0000001,0,0),upper=c(Inf,0.9999,0.9999), eps=eps0)
omega<-res$par[1]
alpha<-res$par[2]
beta<-res$par[3]
varasymp<-garch.etvar(omega, alpha, beta, eps0)
std<-sqrt(c(varasymp$Sigma[1,1],varasymp$Sigma[2,2],varasymp$Sigma[3,3])/n)
list(coef=c(omega,alpha,beta),std=std,cond=varasymp$cond, 
minimum=res$objective,Sigma=varasymp$Sigma)
}


#####################################################
alphainit<-0.05
betainit<-0.85
chouia<-0.01

indices<-c("cac.csv","dax.csv","dja.csv","dowjones.csv","djt.csv","dju.csv","ftse.csv",
"nasdaq.csv","nikkei.csv","smi.csv","sp500.csv")
name.indices<-c("CAC","DAX","DJA","DJI","DJT","DJU","FTSE","Nasdaq","Nikkei","SMI","SP500")
nbind<-length(indices)
tab.pval<-matrix('&',nrow=2*nbind,ncol=9) 
for(ind in (1:nbind)){
### lecture des donnees 
data <- read.table(indices[ind],header=TRUE,sep=",")
close<-rev(data$Close) 
date<-rev(data$Date) 
ntot<-length(close)
rend<-rep(0,ntot); rend[2:ntot]<-log(close[2:ntot]/close[1:(ntot-1)])*100 
ser<-rend[2:ntot]-mean(rend[2:ntot])
#ser<-ser[-which.min(ser)]
n<-length(ser)
#mean(ser) plot(ts(ser)) acf(ser) date[which.min(ser)] ser[which.min(ser)]

#close[which.min(ser)-1] close[(which.min(ser)-2):(which.min(ser)+2)]

omegainit<-var(ser)*(1-alphainit-betainit)

QMLE<-estimgarch11(omegainit,alphainit,betainit,ser)
test<-garch.teststatiofaible(QMLE$coef[1], QMLE$coef[2], QMLE$coef[3], QMLE$Sigma, ser) 
pval<-test$pval
std<-test$std

tab.pval[2*(ind-1)+1,1]<-as.character(name.indices[ind])
tab.pval[2*(ind-1)+2,1]<-""
tab.pval[2*(ind-1)+1,3]<-""
tab.pval[2*(ind-1)+2,3]<-""
tab.pval[2*(ind-1)+1,5]<-round(QMLE$coef[2]+QMLE$coef[3],3)
tab.pval[2*(ind-1)+2,5]<-round(std,3)
tab.pval[2*(ind-1)+1,7]<-""
tab.pval[2*(ind-1)+2,7]<-""
tab.pval[2*(ind-1)+1,9]<-""
tab.pval[2*(ind-1)+2,9]<-round(pval,4)


}

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