# simule un arch(1)
arch.sim<- function (n, omega,  alpha, valinit=100) {
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
eps[t]<-sqrt(ht[t])*eta[t]}
return(eps[(valinit+1):(n+valinit)])}

# estime un ARCH(1) avec nlminb omega0<-1 alpha0<-0.05 eps0<-sim

objf3 <- function(x, eps){     
         omega <- x[1]
         alpha <- x[2]
         n <- length(eps)
         sigma2<-as.numeric(n)
         sigma2[2:n]<-omega+alpha*eps[1:(n-1)]**2
         qml <- mean(log(sigma2[2:n])+eps[2:n]**2/sigma2[2:n])
         qml }

estimarch1.v2<- function(omega0,alpha0,eps0)
     {
valinit<-c(omega0,alpha0)
res <- nlminb(valinit,objf3, lower=c(0.0000001,0),upper=c(Inf,Inf), eps=eps0)
omega<-res$par[1]
alpha<-res$par[2]
list(coef=c(omega,alpha),minimum=res$objective)
}

# estime un ARCH(1) par targeting avec nlminb
objf4 <- function(x, eps)
     {
         alpha <- x[1]
         omega <- var(eps)*(1-alpha)
         sigma2<-as.numeric(n)
         sigma2[2:n]<-omega+alpha*eps[1:(n-1)]**2
         qml <- mean(log(sigma2[2:n])+eps[2:n]**2/sigma2[2:n])
       }
targetestimarch1.v2<- function(alpha0,eps0)
     {
valinit<-c(alpha0)
res <- nlminb(valinit,objf4, lower=0, upper=0.9999,  eps=eps0)
alpha<-res$par[1]
list(coef=c(var(eps0)*(1-alpha),alpha),minimum=res$objective)
}

# niter simulation de tailles n
file.remove("arch.dat")
file.remove("stat.arch1.dat")
n<- 500
niter<- 1000
para<-c(0.1,0.3,0.5,0.55,0.7,0.9)
npara<-length(para)
omegasim<-1
for(ipar in 1:npara){
alphasim<-para[ipar]
file.remove("arch.dat")
for (i in 1:niter) { 
sim<-arch.sim(n,omegasim,alphasim)
alpha.init<-0
TVEarch1<-targetestimarch1.v2(alpha.init,sim)
omega.init<-TVEarch1$coef[1]
alpha.init<-TVEarch1$coef[2]
QMLEarch1<-estimarch1.v2(omega.init,alpha.init,sim)

est<-c(TVEarch1$coef,QMLEarch1$coef)
write(est,file="arch.dat",ncolumns=4,append=TRUE)
     }

# representation des estimations 
toutest <- as.matrix(read.table("arch.dat"))

quant<-as.vector(quantile(toutest[,3]))
stat<-c(mean(toutest[,3]-omegasim),sd(toutest[,3]-omegasim),quant)
write(stat,
file="stat.arch1.dat",ncolumns=7,append=TRUE)

quant<-as.vector(quantile(toutest[,1]))
stat<-c(mean(toutest[,1]-omegasim),sd(toutest[,1]-omegasim),quant)
write(stat,
file="stat.arch1.dat",ncolumns=7,append=TRUE)


quant<-as.vector(quantile(toutest[,4]))
stat<-c(mean(toutest[,4]-alphasim),sd(toutest[,4]-alphasim),quant)
write(stat,
file="stat.arch1.dat",ncolumns=7,append=TRUE)

quant<-as.vector(quantile(toutest[,2]))
stat<-c(mean(toutest[,2]-alphasim),sd(toutest[,2]-alphasim),quant)
write(stat,
file="stat.arch1.dat",ncolumns=7,append=TRUE)


}
toutstat <- as.matrix(read.table("stat.arch1.dat"))
round(toutstat, digit=3)