# calcul des variances asymptique Sigma et Sigmastar 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, tol=sqrt(.Machine$double.eps)) {
garch.etvar<- function (omega,  alpha, beta, eps, tol=0.0001) {
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-alpha-beta,-gamma,-gamma)

for (t in 2:n) {
ht[t] <- omega+alpha*eps[t-1]**2+beta*ht[t-1]
D[,t]<-c(1-alpha-beta,eps[t-1]**2-gamma,ht[t-1]-gamma)+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

inv.Sigmastar<-S%*%t(S)/(n-1)
K<-matrix(inv.Sigmastar[2:3,1],nrow=2,ncol=1)
J<-inv.Sigmastar[2:3,2:3]
if(abs(det(J))>tol)Jinv<-solve(J) else Jinv<-J

if(abs(det(inv.Sigmastar))>tol)Sigmastar<-solve(inv.Sigmastar) else Sigmastar<-inv.Sigmastar


d<-as.numeric(inv.Sigmastar[1,1]-t(K)%*%Jinv%*%K)
d<-1/d


e<-(1-alpha-beta)/(1-beta)
c<-mean(ht[1:n]^2)/e^2
C<-rbind(1,-Jinv%*%K)
Sigma<-Sigmastar+(c-d)*C%*%t(C)

L<-matrix(c(1-alpha-beta,-gamma,-gamma,0,1,0,0,0,1),nrow=3,ncol=3)
Sigma<-t(L)%*%Sigma%*%L
Sigmastar<-t(L)%*%Sigmastar%*%L
list(Sigma=Sigma,Sigma.star=Sigmastar,mu4=mu4,cond=cond)
}

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

objf3 <- function(x, eps){     
         n <- length(eps)
         omega <- x[1]
         alpha <- x[2]
         beta <- x[3]
         if((alpha+beta>0.9999)|(alpha+beta=='NaN'))return(Inf)
         sigma2<-rep(0,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)
     {
n<-length(eps0)
valinit<-c(omega0,alpha0,beta0)
#borne<-2*var(eps0)
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((varasymp$mu4-1)*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)
}

# estime un GARCH(1,1) par targeting avec nlminb
objf4 <- function(x, eps)  {
         n<-length(eps)
         alpha <- x[1]
         beta <- x[2]
         if(alpha+beta>0.999)return(Inf)
         omega <- var(eps)*(1-alpha-beta)
         sigma2<-rep(0,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])
       }
targetestimgarch11<- function(alpha0,beta0,eps0)
     {
n<-length(eps0)
valinit<-c(alpha0,beta0)
res <- nlminb(valinit,objf4, lower=c(0,0), upper=c(0.9999,0.9999),  eps=eps0)
alpha<-res$par[1]
beta<-res$par[2]
omega<-var(eps0)*(1-alpha-beta)
varasymp<-garch.etvar(omega, alpha, beta, eps0)
std<-sqrt((varasymp$mu4-1)*c(varasymp$Sigma.star[1,1],varasymp$Sigma.star[2,2],varasymp$Sigma.star[3,3])/n)
list(coef=c(omega,alpha,beta),std=std,cond=varasymp$cond, minimum=res$objective)
}

############################### VTE sur series reelles
VTE.seriesrelles<-function(indices,name.indices,alphainit,betainit){
nbind<-length(indices)
tab.pval<-matrix('&',nrow=4*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])
if(name.indices[ind]=="Nasdaq")ser<-ser[-which.min(ser)]



VTE<-targetestimgarch11(alphainit,betainit,ser)

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

tab.pval[4*(ind-1)+3,1]<-""
tab.pval[4*(ind-1)+4,1]<-""
tab.pval[4*(ind-1)+3,3]<-round(VTE$coef[1],3)
tab.pval[4*(ind-1)+4,3]<-round(VTE$std[1],3)
tab.pval[4*(ind-1)+3,5]<-round(VTE$coef[2],3)
tab.pval[4*(ind-1)+4,5]<-round(VTE$std[2],3)
tab.pval[4*(ind-1)+3,7]<-round(VTE$coef[3],3)
tab.pval[4*(ind-1)+4,7]<-round(VTE$std[3],3)
tab.pval[4*(ind-1)+3,9]<-""
tab.pval[4*(ind-1)+4,9]<-""
}
tab.pval
}

############################### QMLE sur series reelles
QMLE.seriesrelles.omegainit<-function(indices,name.indices,omegainit,
                       alphainit,betainit){
nbind<-length(indices)
tab.pval<-matrix('&',nrow=4*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])
if(name.indices[ind]=="Nasdaq")ser<-ser[-which.min(ser)]



QMLE<-estimgarch11(omegainit,alphainit,betainit,ser)

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

tab.pval[4*(ind-1)+3,1]<-""
tab.pval[4*(ind-1)+4,1]<-""
tab.pval[4*(ind-1)+3,3]<-""
tab.pval[4*(ind-1)+4,3]<-""
tab.pval[4*(ind-1)+3,5]<-""
tab.pval[4*(ind-1)+4,5]<-""
tab.pval[4*(ind-1)+3,7]<-""
tab.pval[4*(ind-1)+4,7]<-""
tab.pval[4*(ind-1)+3,9]<-""
tab.pval[4*(ind-1)+4,9]<-""
}
tab.pval
}


############################### QMLE sur series reelles 
#avec omegainit=var(ser)*(1-alphainit-betainit)

QMLE.seriesrelles<-function(indices,name.indices,alphainit,betainit){
nbind<-length(indices)
tab.pval<-matrix('&',nrow=4*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])
if(name.indices[ind]=="Nasdaq")ser<-ser[-which.min(ser)]

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



QMLE<-estimgarch11(omegainit,alphainit,betainit,ser)

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

tab.pval[4*(ind-1)+3,1]<-""
tab.pval[4*(ind-1)+4,1]<-""
tab.pval[4*(ind-1)+3,3]<-""
tab.pval[4*(ind-1)+4,3]<-""
tab.pval[4*(ind-1)+3,5]<-""
tab.pval[4*(ind-1)+4,5]<-""
tab.pval[4*(ind-1)+3,7]<-""
tab.pval[4*(ind-1)+4,7]<-""
tab.pval[4*(ind-1)+3,9]<-""
tab.pval[4*(ind-1)+4,9]<-""
}
tab.pval
}

############################### VTE + QMLE 
seriesrelles<-function(indices,name.indices,alphainit,betainit){
nbind<-length(indices)
tab.pval<-matrix('&',nrow=4*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])
if(name.indices[ind]=="Nasdaq")ser<-ser[-which.min(ser)]


VTE<-targetestimgarch11(alphainit,betainit,ser)
omega.init<-VTE$coef[1]
alpha.init<-VTE$coef[2]
beta.init<-VTE$coef[3]

QMLE<-estimgarch11(omega.init,alpha.init,beta.init,ser)

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

tab.pval[4*(ind-1)+3,1]<-""
tab.pval[4*(ind-1)+4,1]<-""
tab.pval[4*(ind-1)+3,3]<-round(VTE$coef[1],3)
tab.pval[4*(ind-1)+4,3]<-round(VTE$std[1],3)
tab.pval[4*(ind-1)+3,5]<-round(VTE$coef[2],3)
tab.pval[4*(ind-1)+4,5]<-round(VTE$std[2],3)
tab.pval[4*(ind-1)+3,7]<-round(VTE$coef[3],3)
tab.pval[4*(ind-1)+4,7]<-round(VTE$std[3],3)
tab.pval[4*(ind-1)+3,9]<-""
tab.pval[4*(ind-1)+4,9]<-""
}
tab.pval
}
#####################################################
alphainit<-0.05
betainit<-0.85

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")

time.VTE<-system.time(VTE.seriesrelles(indices,name.indices,alphainit,betainit))

time.QMLE<-system.time(QMLE.seriesrelles(indices,name.indices,
                    alphainit,betainit))

time.VTE.QMLE<-system.time(seriesrelles(indices,name.indices,alphainit,betainit))

print(time.VTE,digits=2,quote=FALSE)
print(time.QMLE,digits=2,quote=FALSE)
print(time.VTE.QMLE,digits=2,quote=FALSE)


omegainit<-1.
alphainit<-0.0
betainit<-0.0

time.VTE<-system.time(VTE.seriesrelles(indices,name.indices,alphainit,betainit))

time.QMLE<-system.time(QMLE.seriesrelles.omegainit(indices,name.indices,
                    omegainit,alphainit,betainit))

time.VTE.QMLE<-system.time(seriesrelles(indices,name.indices,alphainit,betainit))



print(time.VTE,digits=2,quote=FALSE)
print(time.QMLE,digits=2,quote=FALSE)
print(time.VTE.QMLE,digits=2,quote=FALSE)

QMLE.seriesrelles.omegainit(indices,name.indices,
                    omegainit,alphainit,betainit)