# lecture des donnees du cac
cacdata <- read.table("cac.csv",header=TRUE,sep=",")
str(cacdata); indcac<-cacdata$Close ; datecac<-rev(cacdata$Date)
n<-length(indcac); indcac[1:n]<-indcac[n:1]
cac=rep(0,n); cac[2:n]<-100*log(indcac[2:n]/indcac[1:(n-1)])
#data.entry(as.vector(datecac),indcac,cac)
#rm(list=ls())


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

phi<-function(j,omega,alpha,epsilon){
q<-length(alpha[1,])
d<-length(omega)
volat<-sqrt(abs(omega[j]+sum(alpha[j,1:q]*epsilon[2:(q+1)]^2)))
dnorm(epsilon[1],mean=0,sd=volat)
}

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

#fonction objectif calculée à l'aide de la représentation matricielle
objf <- function(x, y,d,q){
 omega<-rep(0,d); alpha <- matrix(nrow=d,ncol=q)
 p<-matrix(nrow=d,ncol=d)
 for(j in 1:d) {
   omega[j] <- x[(j-1)*(q+1)+1]
   alpha[j,1:q] <- x[((j-1)*(q+1)+2):(j*(q+1))]
                }
 for(i in 1:d) {
  for(j in 2:d) {
    p[i,j] <- x[((q+1)*d+(i-1)*(d-1)+(j-1))]
    p[i,1] <- 1-sum(p[i,2:d]) 
               }
                }
 n <- length(y)
 F<-matrix(0,nrow=d,ncol=n)
 x<-y[(q+1):1]
 pstatio<-piP(p)
 for (j in 1:d) F[j,(q+1)]<-pstatio[j]*phi(j,omega,alpha,x)
 for (k in (q+2):n)  {
  x<-y[k:(k-q)] 
  for (j in 1:d) {
   for (l in 1:d) F[j,k]<-F[j,k]+p[l,j]*F[l,k-1]
   F[j,k]<-F[j,k]*phi(j,omega,alpha,x)
                 }
                   }
 -(log(sum(F[,n])))/n
}

# fonction objectif calculée à l'aide du filtre de Hamilton
objf2 <- function(x, y,d,q){if(min(x)<0)return(Inf)
 omega<-rep(0,d); alpha <- matrix(nrow=d,ncol=q)
 p<-matrix(nrow=d,ncol=d)
 for(j in 1:d) {
   omega[j] <- x[(j-1)*(q+1)+1]
   alpha[j,1:q] <- x[((j-1)*(q+1)+2):(j*(q+1))]
                }
 for(i in 1:d) {
  for(j in 2:d) {
    p[i,j] <- x[((q+1)*d+(i-1)*(d-1)+(j-1))]
    p[i,1] <- 1-sum(p[i,2:d]) 
               }
                }
 if(min(p)<0)return(Inf)
 n <- length(y)
 pit<-matrix(0,nrow=d,ncol=n)
 pitm1<-matrix(0,nrow=d,ncol=n+1)
 vecphi<-rep(0,d)
 vrais<-0
 pstatio<-piP(p)
 pitm1[,(q+1)]<-pstatio
 for (t in (q+1):n)  {
  x<-y[t:(t-q)] 
  for (j in 1:d) vecphi[j]<-phi(j,omega,alpha,x)
  den<-sum(pitm1[,t]*vecphi)
  if(den<=0)return(Inf)
  pit[,t]<-(pitm1[,t]*vecphi)/den 
  pitm1[,t+1]<-t(p)%*%pit[,t]
  vrais<-vrais+log(den) 
                     }
 -vrais/n
}

# fonction objectif calculée à l'aide du filtre de Hamilton
# marche aussi quand q=0
objf3 <- function(x, y,d,q){if(min(x)<0)return(Inf)
 omega<-rep(0,d); if (q>0) alpha <- matrix(nrow=d,ncol=q)
 p<-matrix(nrow=d,ncol=d)
 for(j in 1:d) {
   omega[j] <- x[(j-1)*(q+1)+1]
   if (q>0) alpha[j,1:q] <- x[((j-1)*(q+1)+2):(j*(q+1))]
                }
 for(i in 1:d) {
  for(j in 2:d) {
    p[i,j] <- x[((q+1)*d+(i-1)*(d-1)+(j-1))]
    p[i,1] <- 1-sum(p[i,2:d]) 
               }
                }
 if(min(p)<0)return(Inf)
 n <- length(y)
 pit<-matrix(0,nrow=d,ncol=n)
 pitm1<-matrix(0,nrow=d,ncol=n+1)
 vecphi<-rep(0,d)
 vrais<-0
 pstatio<-piP(p)
 pitm1[,(q+1)]<-pstatio
 for (t in (q+1):n)  {
  x<-y[t:(t-q)] 
  for (j in 1:d) {
if (q>0) vecphi[j]<-phi(j,omega,alpha,x)
else
vecphi[j]<-phi2(j,omega,x)
                 }

  den<-sum(pitm1[,t]*vecphi)
  if(den<=0)return(Inf)
  pit[,t]<-(pitm1[,t]*vecphi)/den 
  pitm1[,t+1]<-t(p)%*%pit[,t]
  vrais<-vrais+log(den) 
                     }
 -vrais/n
}

# calcul de probabilités lissées à l'aide du filtre de Hamilton
# marche aussi quand q=0
liss <- function(x, y,d,q){if(min(x)<0)return(Inf)
 omega<-rep(0,d); if (q>0) alpha <- matrix(nrow=d,ncol=q)
 p<-matrix(nrow=d,ncol=d)
 for(j in 1:d) {
   omega[j] <- x[(j-1)*(q+1)+1]
   if (q>0) alpha[j,1:q] <- x[((j-1)*(q+1)+2):(j*(q+1))]
                }
 for(i in 1:d) {
  for(j in 2:d) {
    p[i,j] <- x[((q+1)*d+(i-1)*(d-1)+(j-1))]
    p[i,1] <- 1-sum(p[i,2:d]) 
               }
                }
 if(min(p)<0)return(Inf)
 n <- length(y)
 vrais<-0
 pit.t<-matrix(0,nrow=d,ncol=n)
 pit.tm1<-matrix(0,nrow=d,ncol=n+1)
 vecphi<-rep(0,d)
 pstatio<-piP(p)
 pit.tm1[,(q+1)]<-pstatio
 for (t in (q+1):n)  {
  x<-y[t:(t-q)] 
  for (j in 1:d){
if (q>0) vecphi[j]<-phi(j,omega,alpha,x)
else
vecphi[j]<-phi2(j,omega,x)}
  den<-sum(pit.tm1[,t]*vecphi)
  if(den<=0)return(Inf)
  pit.t[,t]<-(pit.tm1[,t]*vecphi)/den 
  pit.tm1[,t+1]<-t(p)%*%pit.t[,t]
  vrais<-vrais+log(den) 
                     }
 pit.n<-matrix(0,nrow=d,ncol=n)
 pit.n[,n]=pit.t[,n]
 for (t in n:(q+2))  {
  for (i in 1:d) {pit.n[i,t-1]<-
pit.t[i,t-1]*sum(p[i,1:d]*pit.n[1:d,t]/pit.tm1[1:d,t])}
                     }
pitm1et.n<-array(0,dim=c(d,d,n))
for (t in (q+2):n)  {
 for (i in 1:d) {
  for (j in 1:d) {
pitm1et.n[i,j,t]<-p[i,j]*pit.t[i,t-1]*pit.n[j,t]/pit.tm1[j,t]
                 } } }
                
liss<-list(probaliss=pit.n,probatransliss=pitm1et.n,vrais=vrais)
}

# une iteration de l'algorithm EM pour un MS-ARCH(0)

EM <- function(omega,pi0,p,y){
#if(min(omega,pi0,p)<0)return(Inf)
#for (i in 1:d) if(sum(p[i,1:d])!=1)return(Inf)
 d<-length(omega)
 n <- length(y)
 vrais<-0
 pit.t<-matrix(0,nrow=d,ncol=n)
 pit.tm1<-matrix(0,nrow=d,ncol=n+1)
 vecphi<-rep(0,d)
 pit.tm1[,1]<-pi0
 for (t in 1:n)  {
  for (j in 1:d) vecphi[j]<-dnorm(y[t],mean=0,sd=sqrt(abs(omega[j])))
  den<-sum(pit.tm1[,t]*vecphi)
  if(den<=0)return(Inf)
  pit.t[,t]<-(pit.tm1[,t]*vecphi)/den 
  pit.tm1[,t+1]<-t(p)%*%pit.t[,t]
  vrais<-vrais+log(den) 
                     }
 pit.n<-matrix(0,nrow=d,ncol=n)
 pit.n[,n]=pit.t[,n]
 for (t in n:2)  {
  for (i in 1:d) {pit.n[i,t-1]<-
pit.t[i,t-1]*sum(p[i,1:d]*pit.n[1:d,t]/pit.tm1[1:d,t])}
                     }
pitm1et.n<-array(0,dim=c(d,d,n))
for (t in 2:n)  {
 for (i in 1:d) {
  for (j in 1:d) {
pitm1et.n[i,j,t]<-p[i,j]*pit.t[i,t-1]*pit.n[j,t]/pit.tm1[j,t]
                 } } }

omega.final<-omega
pi0.final<-pi0
p.final<-p
for (i in 1:d)  {
omega.final[i]<-sum((y[1:n]^2)*pit.n[i,1:n])/sum(pit.n[i,1:n])
pi0.final[i]<-pit.n[i,1]
for (j in 1:d) {
 p.final[i,j]<-sum(pitm1et.n[i,j,2:n])/sum(pit.n[i,1:(n-1)])
               }
                }        

liss<-{list(probaliss=pit.n,probatransliss=pitm1et.n,vrais=vrais,
omega.final=omega.final,pi0.final=pi0.final,p.final=p.final)}
}

# MS(4)-ARCH(0)
#d=2
#thetainit<-c(0.81,4.46,0.023,0.898)
#d=3
#thetainit<-c(0.42,1.27,6.77,0.015,0.013,0.982,0.010,0.080,0.0882)
d=4
thetainit<-c(0.40,1.17,2.90,12.11,0.016,0.01,0.01,0.95,0.01,0.01,0.028,0.95,0.01,0.139,0.01,0.682)
# algorithme EM
omega.init<-thetainit[1:d]
p.init<-matrix(0,nrow=d,ncol=d)
for (i in 1:d) {p.init[i,2:d]<-thetainit[(d+(i-1)*(d-1)+1):(d+i*(d-1))]
p.init[i,1]<-1-sum(p.init[i,2:d])}
pi0.init<-piP(p.init)
tol<-0.001
for(iter in 1:60){
     if ( iter == 1) 
     write(c("Itération ",iter),file="MS-CAC.res",ncolumns=2)
     else
     write(c("Itération ",iter),file="MS-CAC.res",ncolumns=2,append=TRUE)
     write("paramètres initiaux ",file="MS-CAC.res",append=TRUE)
for (i in 1:d) { write(c(omega.init[i],p.init[i,]),file="MS-CAC.res",
                  ncolumns=(d+1),append=TRUE)}
resiter<-EM(omega.init,pi0.init,p.init,cac)
     write("vraisemblance initiale ",file="MS-CAC.res",append=TRUE)
     write(resiter$vrais,file="MS-CAC.res",append=TRUE)
omega.final<-resiter$omega.final
pi0.final<-resiter$pi0.final
p.final<-resiter$p.final
vrais.final<-EM(omega.final,pi0.final,p.final,cac)$vrais
     write("nouveaux paramètres ",file="MS-CAC.res",append=TRUE)
for (i in 1:d) { write(c(omega.final[i],p.final[i,]),file="MS-CAC.res",
                  ncolumns=(d+1),append=TRUE)}
     write("nouvelle vraisemblance ",file="MS-CAC.res",append=TRUE)
     write(vrais.final,file="MS-CAC.res",append=TRUE)
crit<-{max(abs(omega.init-omega.final),abs(p.init-p.final),
       abs(resiter$vrais-vrais.final))}
     write("critère ",file="MS-CAC.res",append=TRUE)
     write(crit,file="MS-CAC.res",append=TRUE)
     write("proba invariante ",file="MS-CAC.res",append=TRUE)
 pstatio<-piP(p.final)
     write(pstatio,file="MS-CAC.res",append=TRUE)
   if (crit<tol) break
omega.init<-omega.final
pi0.init<-pi0.final
p.init<-p.final
 }

theta<-rep(0,d^2)
theta[1:d]<-omega.final
for (i in 1:d) theta[(d+(i-1)*(d-1)+1):(d+i*(d-1))]<-p.init[i,2:d]

# graphiques
n<-length(cac) 
lissage<-liss(theta, y=cac,d=d,q=0)
volatplus<-rep(0,n)
volatcol<-rep('blue',n)
color<-c('blue','yellow','green','red')
for (t in 1:n) volatcol[t]<-color[which.max(lissage$probaliss[1:d,t])]
unique(volatcol)

for (t in 1:n) volatplus[t]<-2*sqrt(theta[which.max(lissage$probaliss[1:d,t])])
volatmoins<- -volatplus
Ordrerecord.moins<-order(cac)[1:2]
Rendement.moins<-cac[Ordrerecord.moins]
Date.rendement.moins<-datecac[Ordrerecord.moins]
plot(ts(cac),type="n",axes=FALSE,ylab='CAC 40',xlab=NULL)
lines(ts(cac),col=gray(0.5))
lines(volatplus,col='blue')
lines(volatmoins,col='blue')
box()
axis(2)
lab1<-c("","19/Aug/91","11/Sep/01","")
at1<-c(-500,Ordrerecord.moins[2],Ordrerecord.moins[1],n)
axis(1,at=at1, lab=lab1, cex.axis=1) 

##############