# 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)
}

# Approximation de la VaR à horizon h eleve pour des rendements HMM 
# sig<-c(1,5); eps<-sim; h<-50; niv.alpha<-0.05

VaR.HMM.h.grand<-function(P,sig,V,h,niv.alpha){
d<-length(sig)
n<-length(V)
pstatio<-piP(P)
omega.bar<-sqrt(sum(sig^2*pstatio))
q<-sqrt(h)*omega.bar*qnorm(niv.alpha)
VaR<-V*(1-exp(q))
Ltt.h<--(V[(1+h):n]-V[1:(n-h)])
list(VaR=VaR,Ltt.h=Ltt.h)
}

# VaR à horizon 1 pour des rendements HMM 

VaR.HMM.h.un<-function(P,sig,V,eps,niv.alpha){
h<-1
d<-length(sig)
n<-length(V)
vecphi<-rep(0,d)
pstatio<-piP(P)
pitm1<-pstatio
pit<-pstatio

f<-function(x,p1=p1){
pnorm(x/sig[1])*p1+
pnorm(x/sig[2])*(1-p1)-niv.alpha}

xmin<-sig[1]*qnorm(niv.alpha)
xmax<-sig[2]*qnorm(niv.alpha)
VaR<-rep(0,n)

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<-t(P)%*%pit
  p1<-pitp1[1]
  q<-uniroot(f,c(xmin,xmax),p1=p1)$root
  VaR[t]<-V[t]*(1-exp(q))
                 }




Ltt.h<--(V[(1+h):n]-V[1:(n-h)])
list(VaR=VaR,Ltt.h=Ltt.h)
}


# Calcul la VaR à horizon h eleve pour des rendements  GARCH(1,1) 
# 
VaR.GARCH11.h.grand<-function(V,gamma,h,niv.alpha){
d<-length(sig)
n<-length(V)
sigma.bar<-sqrt(gamma)
q<-sqrt(h)*sigma.bar*qnorm(niv.alpha)
VaR<-V*(1-exp(q))
list(VaR=VaR)
}

# VaR à horizon 1 pour des rendements GARCH

VaR.GARCH11.h.un<-function(V,gamma,alpha,kappa,eps, niv.alpha){
n<-length(V)
ht<-rep(0,(n+1))
q<-rep(0,n)
ht[1]<-gamma
for (t in 2:(n+1)) ht[t]<-ht[t-1]+kappa*(gamma-ht[t-1])+alpha*(eps[t-1]**2-ht[t-1])
q<-sqrt(ht[2:(n+1)])*qnorm(niv.alpha)                     
VaR<-V*(1-exp(q))
list(VaR=VaR)
}


# Approximation de la VaR à horizon h eleve pour des rendements HMM avec le
# filtre de Hamilton 
# sig<-c(1,5); eps<-rep(0.1,10); h<-50; niv.alpha<-0.05

VaR.HMM.Hamilton<-function(P,sig,V,eps,h,niv.alpha){
d<-length(sig)
n<-length(eps)
vecphi<-rep(0,d)
pstatio<-piP(P)
pitm1<-pstatio
pit<-pstatio
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
}
p1<-pitp1[1]
f<-function(x){
pnorm(x/sig[1],lower.tail=FALSE)*p1+
pnorm(x/sig[2],lower.tail=FALSE)*(1-p1)-niv.alpha}

xmin<-sig[1]*qnorm(1-niv.alpha)
xmax<-sig[2]*qnorm(1-niv.alpha)


uniroot(f,c(xmin,xmax))$root
}


# Calcul la VaR à horizon h avec modele GARCH(1,1) 
# avec bootstrap dans les résidus 
# 
VaR.GARCH11<-function(eta,gamma,alpha,kappa,eps,h,niv.alpha,N=1000){
n1<-length(eta)
n<-length(eps)
htm1<-gamma
for (t in 2:n)  {
ht<-htm1+kappa*(gamma-htm1)+alpha*(eps[t-1]**2-htm1)
htm1<-ht
                     }
scenar<-rep(0,N)
sigi<-rep(0,h)
epssim<-rep(0,h)
sigi[1]<-ht
for(iscen in 1:N){
etaboot<-sample(eta, h, replace=TRUE)
#etasim<-rnorm(h)
epssim[1]<-sqrt(sigi[1])*etaboot[1]
if(h>1){
for(i in 2:h) {
sigi[i]<-sigi[i-1]+kappa*(gamma-sigi[i-1])+alpha*(epssim[i-1]**2-sigi[i-1])
epssim[i]<-sqrt(sigi[i])*etaboot[i]
                }
         }
scenar[iscen]<-epssim[h]
                  }
quantile(scenar,1-niv.alpha)
}

# 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)
eps[1:n]<-sig[chaine[1:n]]*eta[1:n]+m[chaine[1:n]]
eps
}

####
# estime un GARCH(1,1) par QMLE
#### 
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)

n <- length(eps0)
sigma2<-as.numeric(n)
sigma2[1]<-eps0[1]^2
gamma<-res$par[1]
alpha<-res$par[2]
kappa<-res$par[3]
for (t in 2:n) {
sigma2[t]<-abs(sigma2[t-1]+kappa*(gamma-sigma2[t-1])+alpha*(eps0[t-1]**2-sigma2[t-1]))
                }    
eta<-eps0[10:n]/sqrt(sigma2[10:n])
list(theta=res$par,ln=res$objective,eta=eta)
}


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

n <- length(eps0)
sigma2<-rep(0,n)
sigma2[1]<-eps0[1]^2
gamma<-var(eps0)
alpha<-res$par[1]
kappa<-res$par[2]
for (t in 2:n) {
sigma2[t]<-abs(sigma2[t-1]+kappa*(gamma-sigma2[t-1])+alpha*(eps0[t-1]**2-sigma2[t-1]))
                }    
eta<-eps0[10:n]/sqrt(sigma2[10:n])
list(theta=c(var(eps0),res$par),ln=res$objective,eta=eta)
}

####
# 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
#niv.alpha<-0.05
backtest<-rep(0,3)
VaR.moy<-rep(0,3)
P<-matrix(c(0.90,0.10,0.10,0.90),nrow=2,ncol=2)
sig<-c(1,5)/200
#sig<-c(1,5)
m<-c(0,0)
sim<-simHMM(n, P, m, sig)
std<-sd(sim)
sim<-sim/std
alphainit<-0.10
betainit<-0.75
q<-length(alphainit)
p<-length(betainit)
kappainit<-1-alphainit-betainit
restargetGARCH11<-estimgarch11.target(alphainit,kappainit,sim[1:n])
restqmlGARCH11<-estimgarch11(restargetGARCH11$theta[1],
restargetGARCH11$theta[2],restargetGARCH11$theta[3],sim[1:n])
restargetGARCH11$theta[1]<-restargetGARCH11$theta[1]*std**2
restqmlGARCH11$theta[1]<-restqmlGARCH11$theta[1]*std**2
restargetGARCH11$theta
restqmlGARCH11$theta

set.seed(8)

niv.alpha<-0.1
n2<- 50000
sim<-simHMM(n2, P, m, sig)
plot(ts(sim))
#sim<-rnorm(n2,sd=sqrt(mean(sig**2)))

V<-sim
V0<- 1000
V[1]<-exp(sim[1])*V0
for(t in (2:n2))V[t]<-exp(sim[t])*V[t-1]
plot(ts(V))

op <- par(mfrow = c(1, 2),cex.main=0.8) # 2 x 2 pictures on one plot


h<-1

Varshortterm<-VaR.HMM.h.un(P,sig,V,sim,niv.alpha)
Varshortterm.VaR<-Varshortterm$VaR 
Varshortterm.Ltt.h<-Varshortterm$Ltt.h
ntot<-length(Varshortterm.Ltt.h)

Varshortterm.QML<-VaR.GARCH11.h.un(V,restqmlGARCH11$theta[1],
                         restqmlGARCH11$theta[2],restqmlGARCH11$theta[3],sim,niv.alpha)

Varshortterm.QML.VaR<-Varshortterm.QML$VaR 

Varshortterm.VTE<-VaR.GARCH11.h.un(V,restargetGARCH11$theta[1],
                         restargetGARCH11$theta[2],restargetGARCH11$theta[3],sim,niv.alpha)
Varshortterm.VTE.VaR<-Varshortterm.VTE$VaR 

backtest[1]<-length(which(Varshortterm.VaR[1:ntot]<Varshortterm.Ltt.h[1:ntot]))/ntot
VaR.moy[1]<-mean(Varshortterm.VaR[1:ntot]/V[1:ntot])

backtest[2]<-length(which(Varshortterm.QML.VaR[1:ntot]<Varshortterm.Ltt.h[1:ntot]))/ntot
VaR.moy[2]<-mean(Varshortterm.QML.VaR[1:ntot]/V[1:ntot])

backtest[3]<-length(which(Varshortterm.VTE.VaR[1:ntot]<Varshortterm.Ltt.h[1:ntot]))/ntot
VaR.moy[3]<-mean(Varshortterm.VTE.VaR[1:ntot]/V[1:ntot])

VaR.moy.h.un<-VaR.moy
backtest.h.un<-backtest
min<-min(Varshortterm.Ltt.h)
max<-max(Varshortterm.Ltt.h,Varshortterm.QML.VaR)
plot(ts(Varshortterm.Ltt.h),ylim=c(min,max),ylab='P&L and VaR',
main='VaR at horizon h=1',xlab=NULL)
ntot<-length(Varshortterm.Ltt.h)
lines(Varshortterm.VaR[1:ntot],lwd=2)
lines(Varshortterm.QML.VaR[1:ntot],lty=3,lwd=2,col='red')
lines(Varshortterm.VTE.VaR[1:ntot],,lty=2,lwd=2,col='green')

h<-10

Varlongterm<-VaR.HMM.h.grand(P,sig,V,h,niv.alpha)
Varlongterm.VaR<-Varlongterm$VaR 
Varlongterm.Ltt.h<-Varlongterm$Ltt.h
ntot<-length(Varlongterm.Ltt.h)

Varlongterm.QML<-VaR.GARCH11.h.grand(V,
                         restqmlGARCH11$theta[1],h,niv.alpha)
Varlongterm.QML.VaR<-Varlongterm.QML$VaR 

Varlongterm.VTE<-VaR.GARCH11.h.grand(V,
                         restargetGARCH11$theta[1],h,niv.alpha)
Varlongterm.VTE.VaR<-Varlongterm.VTE$VaR 

backtest[1]<-length(which(Varlongterm.VaR[1:ntot]<Varlongterm.Ltt.h[1:ntot]))/ntot
VaR.moy[1]<-mean(Varlongterm.VaR[1:ntot]/V[1:ntot])

backtest[2]<-length(which(Varlongterm.QML.VaR[1:ntot]<Varlongterm.Ltt.h[1:ntot]))/ntot
VaR.moy[2]<-mean(Varlongterm.QML.VaR[1:ntot]/V[1:ntot])

backtest[3]<-length(which(Varlongterm.VTE.VaR[1:ntot]<Varlongterm.Ltt.h[1:ntot]))/ntot
VaR.moy[3]<-mean(Varlongterm.VTE.VaR[1:ntot]/V[1:ntot])

VaR.moy.h.10<-VaR.moy
backtest.h.10<-backtest
min<-min(Varlongterm.Ltt.h)
max<-max(Varlongterm.Ltt.h,Varlongterm.QML.VaR)

plot(ts(Varlongterm.Ltt.h),ylim=c(min,max),ylab='P&L and VaR',
main='VaR at horizon h=10',xlab=NULL)
ntot<-length(Varlongterm.Ltt.h)
lines(Varlongterm.VaR[1:ntot],lwd=2)
lines(Varlongterm.QML.VaR[1:ntot],lty=3,lwd=2,col='red')
lines(Varlongterm.VTE.VaR[1:ntot],lty=2,lwd=2,col='green')

round(100*VaR.moy.h.un,digit=2)
round(100*backtest.h.un,digit=2)
round(100*VaR.moy.h.10,digit=2)
round(100*backtest.h.10,digit=2) 
par(op)



###########################"""