# bruit blanc melange de 3 normales
bruit<- function (n, m0=0,  sig0=0.5, m1=3, sig1=1., p=0.05) 
{
eps <- rep(0,n)
for (i in 1:n) { 
u <- runif(1)
u
if ( u < p ) {eps[i] <- rnorm(1,mean=m1,sd=sig1)} 
if ( (p <= u) & (u < 2*p) ) {eps[i] <- rnorm(1,mean=-m1,sd=sig1)} 
if ( u > 2*p ) {eps[i] <- rnorm(1,mean=m0,sd=sig0)} 
}
as.ts(eps)
}

# fonction qui simule un STAR(1) logistique
star.sim <- function (n, a1=-0.95,  b1=0.9,  g0=5.0,  burnin=100) 
{
y <- rep(0,n+burnin)
eps <- bruit(n+burnin)

for (i in 2:(n+burnin)) 
{y[i] <- eps[i] + a1*y[i-1]+(b1-a1)*y[i-1]/(1+exp(-g0*y[i-1])) 
as.ts(y[(burnin+1):(n+burnin)])
}}

# simulation
n <- 2000
set.seed(7)
a1<- -0.95
b1<- 0.9
g0<- 5.
y <- star.sim(n, a1=a1,  b1=b1,  g0=g0) 

# estimation frequentiste d'un modele LSTAR(1)
library(tsDyn)
#vignette("tsDyn")
# fit a STAR model to the time series y
fittedstar<-lstar(y,m=1,thDelay=0,control=list(maxit=3000))
names(fittedstar$coefficients)<-c("a0","a1","b0","b1","gam","c")
summary(fittedstar)

# valeurs initiales
theta.init0<-{as.vector(c(fittedstar$coefficients[1:5],
                   var(residuals(fittedstar),na.rm=T)))}

library(MCMCpack)

p3<-function (gam,beta,sig2,y,n) {if ( gam<=0 ) return(-Inf)
Zt <- matrix(ncol = n, nrow = 4); Gt<-1/(1+exp(-gam*y))
Zt[1,]<- 1-Gt; Zt[2,]<- y*(1-Gt); Zt[3,]<- Gt; Zt[4,]<- y*Gt
eps<-c(rep(0,n)); eps[2:n]<-y[2:n]-beta%*%Zt[1:4,1:(n-1)]
-gam-sum(eps^2)/(2*sig2) }

gibbsMH<-function (y,nbatch=1000,theta.init=0.1,tun=1) {
n<-length(y); eps<-c(rep(0,n)); Zt <- matrix(nrow=4, ncol=n)
mat <- matrix(nrow=6, ncol=nbatch); mat[,1] <- theta.init
 for (i in 2:nbatch) { Gt<-1/(1+exp(-mat[5,i-1]*y))
  Zt[1,]<- 1-Gt; Zt[2,]<- y*(1-Gt); Zt[3,]<- Gt; Zt[4,]<- y*Gt
  ZtZ <- Zt[1:4,1:(n-1)]%*%t(Zt[1:4,1:(n-1)])
  Zty<-{c(sum(y[2:n]*Zt[1,1:(n-1)]),sum(y[2:n]*Zt[2,1:(n-1)]),
         sum(y[2:n]*Zt[3,1:(n-1)]),sum(y[2:n]*Zt[4,1:(n-1)]))}
  bhat<-solve(ZtZ,Zty);  Sigmaninv<-ZtZ/mat[6,i-1] 
  Sigma<-solve(diag(rep(1,4))+Sigmaninv)
  mu<-Sigma%*%(Sigmaninv%*%bhat)
  mat[1:4,i]<-mvrnorm(n=1,as.vector(mu), Sigma)
  eps[2:n]<-y[2:n]-bhat%*%Zt[1:4,1:(n-1)]
  mat[6,i]<- rinvgamma(1,shape=1+(n-1)/2,scale=1+sum(eps^2)/2)
  gamstar <- rnorm(1, mean=mat[5,i-1], sd=tun)
  unif<-runif(1,0,1)
  lognum<-p3(gamstar,bhat,mat[6,i],y,n)
  logden<-p3(mat[5,i-1],bhat,mat[6,i],y,n)
  if(lognum==-Inf) ratio<-0 else {if (logden==-Inf) 
           ratio<-1 else ratio<-exp(lognum-logden)}
  mat[5,i] <- if(unif <= ratio) gamstar else mat[5,i-1] 
   }
t(mat)}

mat<-gibbsMH(y,nbatch=100,theta.init=theta.init0,tun=1)

theta.init0<-mat[length(mat[,1]),]

mat<-gibbsMH(y,nbatch=1000,theta.init=theta.init0,tun=1)


chaine<-ts(mat)
dimnames(chaine)[[2]]<-c("a0","a1","b0","b1","gam","sig2")
plot(chaine,main="Markov chain simulated by Metropolis")
acf(chaine) # plus de non stationnarité flagrante

apply(mat, 2, mean)

chaine <- as.data.frame(chaine)
rm(a1,b1)
attach(chaine)
plot(ts(a1))
str(chaine)
ls() 

op <- par(mfrow = c(2, 2), ask=T) # 2 x 2 figures par page 
plot(ts(gam),main=expression(paste('Trace of ',gamma)),xlab='iteration',
ylab=expression(gamma))
{hist(gam, col='light blue', probability=T,xlab='',
   main=expression(paste('Posterior distribution of ',gamma)))}
lines(density(gam), col='red', lwd=3)
plot(ts(sig2),main=expression(paste('Trace of ',sigma^2)),xlab='iteration',
ylab=expression(sigma^2))
{hist(sig2, col='light blue', probability=T,xlab='',
   main=expression(paste('Posterior distribution of ',sigma^2)))}
lines(density(sig2), col='red', lwd=3)
par(op) # retour aux valeurs par défaut des paramètres graphiques