# piP = proba invariate d'une chaine avec proba de transition p
piP<-function(P){
d<-length(P[1,])
A<-rbind(t(P)-diag(rep(1,d)),1)
b<-c(rep(0,d),1)
solve(A[2:(d+1),1:d],b[2:(d+1)])
}

phi2<-function(j,omega,x){
d<-length(omega)
volat<-abs(omega[j])
dnorm(x,mean=0,sd=volat)
}

# filtre de Hamilton pour prévoir les carres à horizon h
# sig<-c(1,5); eps<-rep(0.1,10); h<-50
Hamilton<-function(P,sig,eps,h){
d<-length(sig)
n<-length(eps)
vecphi<-rep(0,d)
pstatio<-piP(P)
pitm1<-pstatio
pit<-pstatio
pitph<-matrix(nrow=d,ncol=h)
sigti.t<-rep(0,h)
for (t in 1:n)  {
for (j in 1:d) vecphi[j]<-phi2(j,sig,eps[t])
  den<-sum(pitm1*vecphi)
  if(den<=0)return(Inf)
  pit<-(pitm1*vecphi)/den 
  pitm1<-t(P)%*%pit
                     }
pitp1<-pit
for(i in 1:h) {
pitp1<-t(P)%*%pitp1
sigti.t[i]<-sum((sig**2)*pitp1)
}
sqrt(sigti.t)
}

# Modele GARCH(1,1) pour prévoir les carres à horizon h
# gamma<-3; alpha<-0.5; kappa<-0.1
# 
prevGARCH11<-function(gamma,alpha,kappa,eps,h){
n<-length(eps)
sigti.t<-rep(0,h)
htm1<-gamma
for (t in 2:n)  {
ht<-htm1+kappa*(gamma-htm1)+alpha*(eps[t]**2-htm1)
htm1<-ht
                     }
for(i in 1:h) {
sigti.t[i]<-gamma+(ht-gamma)*((1-kappa)**(i+1))/(1-kappa)
}
sqrt(sigti.t)
}

# n<- 10 P<-matrix(c(0.1,0.1,0.9,0.9),nrow=2,ncol=2) fonction qui simule une chaine de markov
MarkovChain.sim <- function (n, P) 
{
d<-length(P[1,])
chaine<-rep(0,n)

u<-runif(1)

sum<-cumsum(piP(P))
chaine[1]<-min(which(u<=sum))
for (t in 2:n){
u<-runif(1)
sum<-cumsum(P[chaine[t-1],])
chaine[t]<-min(which(u<=sum))
}
chaine
}

# simule un HMM
simHMM<- function (n, P, m, sig) 
{
d<-length(P[1,])
eta <- rnorm(n)
eps<-as.numeric(n)
chaine<-MarkovChain.sim(n, P)
for (t in 1:n) eps[t]<-sig[chaine[t]]*eta[t]+m[chaine[t]]
eps
}

###
# omega<-omegasim; alpha<-alphasim;beta<-betasim
# simulation d'un GARCH(p,q) 
simGARCH <- 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
#### x<-c(2,0.5,0.1) eps<-rnorm(100)
objf.GARCH11 <- function(x,eps){  
gamma<-x[1]; alpha<-x[2]; kappa<-x[3]

petit<-0.001

if((1-alpha-kappa) > (1-petit))return(Inf)

gamma<-max(gamma,petit)
alpha<-max(alpha,0)
kappa<-max(kappa,petit)

         n <- length(eps)
         sigma2<-as.numeric(n)
         sigma2[1]<-eps[1]^2
for (t in 2:n) {
sigma2[t]<-abs(sigma2[t-1]+kappa*(gamma-sigma2[t-1])+alpha*(eps[t-1]**2-sigma2[t-1]))
                }    
qml <- mean(log(sigma2[2:n])+((eps[2:n]**2)/sigma2[2:n]))
         qml }
#eps0<-sim; alpha0<-alphainit; beta0<-betainit  

estimgarch11<- function(gamma0,alpha0,kappa0,eps0)
     {
valinit<-c(gamma0,alpha0,kappa0)
petit<-0.00000001
res <- nlminb(valinit,objf.GARCH11, lower=c(petit,0,petit),
upper=c(Inf,rep(1,2)), eps=eps0)
list(theta=res$par,ln=res$objective)
}


# estimation d'un modele GARCH(p,q) par QML sur la serie eps
#p<-2; q<-1
#max(p,q)
# x<-c(omegainit,alphainit,betainit); eps<-sim
objf <- function(x,p,q,eps){     
         omega <- x[1]
         alpha <- x[2:(q+1)]
if(p>0)  beta <- x[(q+2):(p+q+1)]
         r<-max(p,q)
         n <- length(eps)
         sigma2<-as.numeric(n)
         sigma2[1:r]<-eps[1:r]^2
for (t in (r+1):n) {
if(p>0)
{sigma2[t]<-omega+sum(alpha[1:q]*(eps[(t-1):(t-q)]^2))+sum(beta[1:p]*sigma2[(t-1):(t-p)])}
else
{sigma2[t]<-omega+sum(alpha[1:q]*(eps[(t-1):(t-q)]^2)) }
                   }    
         qml <- mean(log(sigma2[(r+1):n])+((eps[(r+1):n]**2)/sigma2[(r+1):n]))
         qml }
     
#eps0<-sim; omega0<-omegainit; alpha0<-alphainit; beta0<-betainit  
estimgarchpq<- function(omega0,alpha0,beta0,eps0)
     {
valinit<-c(omega0,alpha0,beta0)
q<-length(alpha0)
p<-length(beta0)
res <- nlminb(valinit,objf, lower=c(0.00000001,rep(0,(p+q))),
upper=c(Inf,rep(3,(p+q))), 
eps=eps0,p=p,q=q)
list(theta=res$par,ln=res$objective)
}

####
# estime un GARCH(1,1) par targeting
#### x<-c(2,0.5,0.1) eps<-rnorm(100)
objf.targetGARCH11 <- function(x,eps){  
alpha<-x[1]; kappa<-x[2]

gamma<-var(eps)
petit<-0.00000001
if((1-alpha-kappa) > (1-petit))return(Inf)

gamma<-max(gamma,petit)
alpha<-max(alpha,0)
kappa<-max(kappa,petit)

         n <- length(eps)
         sigma2<-as.numeric(n)
         sigma2[1]<-eps[1]^2
for (t in 2:n) {
sigma2[t]<-sigma2[t-1]+kappa*(gamma-sigma2[t-1])+alpha*(eps[t-1]**2-sigma2[t-1])
                }    
qml <- mean(log(sigma2[2:n])+((eps[2:n]**2)/sigma2[2:n]))
         qml }
#eps0<-sim; alpha0<-alphainit; beta0<-betainit  
estimgarch11.target<- function(alpha0,kappa0,eps0)
     {
valinit<-c(alpha0,kappa0)
petit<-0.00000001
res <- nlminb(valinit,objf.targetGARCH11, lower=c(0,petit),
upper=c(rep(1,2)), eps=eps0)
list(theta=c(var(eps0),res$par),ln=res$objective)
}

####
# x<-valinit; eps<-eps0
# estime un GARCH(p,q) par targeting
####
objf.target <- function(x,p,q,eps){     
         alpha <- x[1:q]
if(p>0){
beta <- x[(q+1):(p+q)]
omega<-var(eps)*(1-sum(alpha)-sum(beta))
       }
if(p==0) {omega<-var(eps)*(1-sum(alpha))}


if(omega<0.00000001)omega<-0.00000001
         r<-max(p,q)
         n <- length(eps)
         sigma2<-as.numeric(n)
         sigma2[1:r]<-eps[1:r]^2
for (t in (r+1):n) {
if(p>0)
{sigma2[t]<-omega+sum(alpha[1:q]*(eps[(t-1):(t-q)]^2))+sum(beta[1:p]*sigma2[(t-1):(t-p)])}
else
{sigma2[t]<-omega+sum(alpha[1:q]*(eps[(t-1):(t-q)]^2)) }
                   }    
         qml <- mean(log(sigma2[(r+1):n])+((eps[(r+1):n]**2)/sigma2[(r+1):n]))
         qml }
#eps0<-sim; alpha0<-alphainit; beta0<-betainit  
estimgarchpq.target<- function(alpha0,beta0,eps0)
     {
valinit<-c(alpha0,beta0)
q<-length(alpha0)
p<-length(beta0)
res <- nlminb(valinit,objf.target, lower=c(rep(0,(p+q))),
upper=c(rep(1,(p+q))), eps=eps0,p=p,q=q)
omega<-var(eps0)*(1-sum(res$par))
list(theta=c(omega,res$par),ln=res$objective)
}


#library(tseries)





# niter simulation de tailles n
set.seed(7)

n<- 1000
P<-matrix(c(0.90,0.10,0.10,0.90),nrow=2,ncol=2)
sig<-c(1,5)
m<-c(0,0)
sim<-simHMM(n, P, m, sig)
#plot(ts(sim)) 
#alphainit<-0.75
#omegainit<-var(sim)*(1-alphainit)
#q<-length(alphainit)
#betainit<-c()
#p<-length(betainit)
#restqmlARCH1<-estimgarchpq(omegainit,alphainit,betainit,sim)
#restargetARCH1<-estimgarchpq.target(alphainit,betainit,sim)


alphainit<-0.10
betainit<-0.75
#omegainit<-var(sim)*(1-alphainit-betainit)
q<-length(alphainit)
p<-length(betainit)
#restqmlGARCH11<-estimgarchpq(omegainit,alphainit,betainit,sim)
#restargetGARCH11<-estimgarchpq.target(alphainit,betainit,sim)
kappainit<-1-alphainit-betainit
restqmlGARCH11<-estimgarch11(var(sim),alphainit,kappainit,sim)
restargetGARCH11<-estimgarch11.target(alphainit,kappainit,sim)



#var(sim)
#sum(piP(P)*(sig**2))
#restqmlARCH1$theta[1]/(1-restqmlARCH1$theta[2])
#restargetARCH1<-estimgarchpq.target(alphainit,betainit,sim)
#restargetARCH1$theta[1]/(1-restargetARCH1$theta[2])




varlimit<-sqrt(sum(piP(P)*(sig**2)))

h<-20
op <- par(mfrow = c(1, 2)) # 2 x 2 pictures on one plot

obs<-rep(0,100)

#gamma<-restqmlARCH1$theta[1]/(1-restqmlARCH1$theta[2])
#alpha<-restqmlARCH1$theta[2]
#kappa<-1-alpha
#prevARCH1QML<-prevGARCH11(gamma,alpha,kappa,obs,h)
gamma<-restqmlGARCH11$theta[1]
alpha<-restqmlGARCH11$theta[2]
kappa<-restqmlGARCH11$theta[3]
prevGARCH11QML<-prevGARCH11(gamma,alpha,kappa,obs,h)



gamma<-restargetGARCH11$theta[1]
alpha<-restargetGARCH11$theta[2]
kappa<-restargetGARCH11$theta[3]
prevGARCH11TVE<-prevGARCH11(gamma,alpha,kappa,obs,h)


prevexact<-Hamilton(P,sig,obs,h)
maxval<-1.96*max(prevexact,varlimit,prevGARCH11QML,prevGARCH11TVE)

plot(ts(1.96*prevexact),ylim=c(-maxval,maxval),type="n",
   ylab='Prediction interval',xlab="horizon h")
lines(ts(1.96*prevexact),lwd=2)
lines(-ts(1.96*prevexact),lwd=2)
abline(1.96*varlimit,0,col='blue')
abline(-1.96*varlimit,0,col='blue')
lines(ts(1.96*prevGARCH11QML),lty=3,lwd=2,col='red')
lines(-ts(1.96*prevGARCH11QML),lty=3,lwd=2,col='red')
lines(ts(1.96*prevGARCH11TVE),lty=2,lwd=2,col='green')
lines(-ts(1.96*prevGARCH11TVE),lty=2,lwd=2,col='green')
legend(h/2,5,c("exact","asymp","QMLE","VTE"),
col=c("black","blue","red","green"),lwd=c(1.5,1,2,2),lty=c(1,1,3,2))


obs<-rep(sig[2],100)


gamma<-restqmlGARCH11$theta[1]
alpha<-restqmlGARCH11$theta[2]
kappa<-restqmlGARCH11$theta[3]
prevGARCH11QML<-prevGARCH11(gamma,alpha,kappa,obs,h)



gamma<-restargetGARCH11$theta[1]
alpha<-restargetGARCH11$theta[2]
kappa<-restargetGARCH11$theta[3]
prevGARCH11TVE<-prevGARCH11(gamma,alpha,kappa,obs,h)


prevexact<-Hamilton(P,sig,obs,h)
maxval<-1.96*max(prevexact,varlimit,prevGARCH11QML,prevGARCH11TVE)

plot(ts(1.96*prevexact),ylim=c(-maxval,maxval),type="n",
   ylab='Prediction interval',xlab="horizon h")
lines(ts(1.96*prevexact),lwd=2)
lines(-ts(1.96*prevexact),lwd=2)
abline(1.96*varlimit,0,col='blue')
abline(-1.96*varlimit,0,col='blue')
lines(ts(1.96*prevGARCH11QML),lty=3,lwd=2,col='red')
lines(-ts(1.96*prevGARCH11QML),lty=3,lwd=2,col='red')
lines(ts(1.96*prevGARCH11TVE),lty=2,lwd=2,col='green')
lines(-ts(1.96*prevGARCH11TVE),lty=2,lwd=2,col='green')
legend(h/2,6,c("exact","asymp","QMLE","VTE"),
col=c("black","blue","red","green"),lwd=c(1.5,1,2,2),lty=c(1,1,3,2))


     par(op)