# CHAPTER 7: R COMMANDS
# 08/01/2007

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

# R PROGRAM TO ESTIMATE THE COEFFICIENTS OF AN AR(p) MODEL

p=5
x=matrix(scan("Eurostoxx50.txt"),ncol=5,byrow=T)[,4]
T=length(x)
W=10000

mu=mean(x)
sig2=var(x)

lambdareal=2*runif(p)-1
preal=p
pcomp=0
lambdacomp=0

llog=function(pr,pc,lr,lc,compsi=T,pepsi=0){

# LIKELIHOOD REPRESENTATION

if (compsi)
{

Psi=matrix(0,ncol=p,nrow=p+1)
Psi[1,]=1

if (pr>0)
{
Psi[2,1]=-lr[1]

if (pr>1)
{
for (i in 2:pr)
Psi[2:(i+1),i]=Psi[2:(i+1),i-1]-lr[i]*Psi[1:i,i-1]
}
}
if (pc>0)
{
if (pr>0)
{
Psi[2,pr+2]=-2*lc[1]+Psi[2,pr]
Psi[3:(pr+3),pr+2]=(lc[1]^2+lc[2]^2)*Psi[1:(pr+1),pr]-2*lc[1]*Psi[2:(pr+2),pr]+Psi[3:(pr+3),pr]
}
else
{
Psi[2,2]=-2*lc[1];
Psi[3,2]=(lc[1]^2+lc[2]^2);
}
if (pc>2)
{
for (i in seq(4,pc,2))
{
pri=pr+i
prim=pri-2
Psi[2,pri]=-2*lc[i-1]+Psi[2,prim]
Psi[3:(pri+1),pri]=(lc[i-1]^2+lc[i]^2)*Psi[1:(pri-1),prim]-2*lc[i-1]*Psi[2:pri,prim]+Psi[3:(pri+1),prim]
}
}
}
Psi=Psi[1:(p+1),p]
}
else
{
Psi=pepsi
}

# LOGLIKELIHOOD

x=x-mu
loglike=-x[1]^2

for (i in 2:p) loglike=loglike-(t(Psi[1:i])%*%x[i:1])^2

for (i in (p+1):T) loglike=loglike-(t(Psi)%*%x[i:(i-p)])^2

loglike=(loglike/sig2-T*log(sig2))/2

x=x+mu

list(ll=loglike,ps=Psi)
}

llo=llog(pr=preal,pc=pcomp,lr=lambdareal,lc=lambdacomp)

# MCMC REVERSIBLE SCHEME

psis=matrix(0,ncol=p,nrow=W)
mus=rep(0,W)
sigs=rep(0,W)
ncomp=rep(0,W)
llik=rep(0,W)
indacpt=rep(0,5)

for (m in 1:W){

if (runif(1)<.1{
  ind=sample(1:p,1)

  if (ind<=pcomp){ 
     ind=ind-(ind%%2==0) 
     ppropreal=preal
     ppropcomp=pcomp
     lambpropreal=lambdareal
     lambpropcomp=lambdacomp
     lambpropcomp[ind]=lambdacomp[ind]+.05*rnorm(1)
     lambpropcomp[ind+1]=lambdacomp[ind+1]+.05*rnorm(1)

     }else{

      ppropreal=preal
      ppropcomp=pcomp
      lambpropreal=lambdareal
      lambpropcomp=lambdacomp
      lambpropreal[ind-pcomp]=lambdareal[ind-pcomp]+.05*rnorm(1) 
      }

  lloprop=llog(pr=ppropreal,pc=ppropcomp,lr=lambpropreal,lc=lambpropcomp)

  if (log(runif(1))<lloprop$ll-llo$ll){

    llo=lloprop
    preal=ppropreal
    pcomp=ppropcomp
    lambdacomp=lambpropcomp
    lambdareal=lambpropreal
    indacpt[1]=indacpt[1]+1
    }
}

   if (preal<2){

       ppropreal=preal+2
       ppropcomp=pcomp-2

       ind=sample(1:pcomp,1)
       ind=ind-(ind%%2==0) 

       if (preal==0){
         lambpropreal=2*runif(2)-1
       }else{

	 lambpropreal=c(lambdareal,lambdacomp[ind]+0.01*rnorm(2)) #r2*runif(2)-1)
       }

       lambpropcomp=lambdacomp[((1:pcomp)!=ind)&((1:pcomp)!=(ind+1))]

       coef=(1/2) * (pi/4)
       }

   if (pcomp==0){

       ppropreal=p-2
       ppropcomp=2

       ind=sample(1:p,2)  # real roots to remove
       lambpropreal=lambdareal[((1:p)!=ind[1])&((1:p)!=ind[2])] 

       lambpropcomp=c(mean(lambdareal[ind])+.01*rnorm(1),.01*rnorm(1)) #rho*c(cos(theta),sin(theta))

       coef=(1/2) * (4/pi)
       }


   if ((preal>1)&&(pcomp>0)){
      
       if (runif(1)<.1){
       
         ppropcomp=pcomp-2
         ppropreal=preal+2

         ind=sample(1:pcomp,1)
         ind=ind-(ind%%2==0) 

         if (ppropcomp>0){
           lambpropcomp=lambdacomp[((1:pcomp)!=ind)&((1:pcomp)!=(ind+1))]
         }else{
           lambpropcomp=0
         }

         lambpropreal=c(lambdareal,lambdacomp[ind]+0.05*rnorm(2))

         coef=9*(1+(preal<2)) * (pi/4)

       }else{ 

         ppropreal=preal-2
         ppropcomp=pcomp+2


         ind=sample(1:preal,2) 	#real roots to remove
                                                                                                                                         
         if (ppropreal>0){

           lambpropreal=lambdareal[((1:preal)!=ind[1])&((1:preal)!=ind[2])]
	   }else{
          
              lambpropreal=0
           }
	
        lambpropcomp=c(lambdacomp,mean(lambdareal[ind])+.05*rnorm(1),.05*rnorm(1))

        coef=(4/pi) * (1+(ppropcomp<p-1))/9.
      }
    }

  lloprop=llog(pr=ppropreal,pc=ppropcomp,lr=lambpropreal,lc=lambpropcomp)

  if (log(runif(1))<log(coef)+lloprop$ll-llo$ll){

    llo=lloprop
    preal=ppropreal
    pcomp=ppropcomp
    lambdacomp=lambpropcomp
    lambdareal=lambpropreal
    indacpt[2]=indacpt[2]+1
    }

  psis[m,]=llo$ps[2:(p+1)]

if (runif(1)<.1){ #Only once in a while
  muold=mu
  mu=rnorm(1,mean=muold,sd=sqrt(2*sig2))
    
  lloprop=llog(pr=preal,pc=pcomp,lr=lambdareal,lc=lambdacomp,compsi=F,pepsi=llo$ps)

  if (log(runif(1))>lloprop$ll-llo$ll){
 
    indacpt[3]=1+indacpt[3]
    mu=muold
    }else{

       llo=lloprop
    }

  sig2old=sig2
  sig2=exp(rnorm(1,mean=log(sig2),sd=sqrt(2)*sig2))

  lloprop=llog(pr=preal,pc=pcomp,lr=lambdareal,lc=lambdacomp,compsi=F,pepsi=llo$ps)

  if (log(runif(1))>lloprop$ll-llo$ll){

    indacpt[4]=1+indacpt[4]
     sig2=sig2old
     }else{

        llo=lloprop
     }
  }

  psis[m,]=llo$ps[2:(p+1)]
  mus[m]=mu
  sigs[m]=sig2
  llik[m]=llo$ll
  ncomp[m]=pcomp
  }

# GRAPHICAL REPRESENTATION FOR THE AR(p) MODEL

par(mfrow=c(3,3),mar=c(4,4,2,1))
if (p>3)
{
hist(ncomp,main="",xlab="p",ylab="",col="gold4",breaks=seq(-1,p+1,2));par(new=T);plot(ncomp,axes=F,cex=.3,xlab="",ylab="")
axis(side=4)
}
if (p<4)
{
plot(ncomp,cex=.3,xlab="Iterations",ylab="Complex roots")
}
plot(mus,type="l",col="steelblue4",xlab="Iterations",ylab=expression(mu))
plot(500:W,sigs[500:W],type="l",col="steelblue4",xlab="Iterations",ylab=expression(sigma^2))
for (i in 1:min(3)) plot(500:W,psis[500:W,i],type="l",col="steelblue4",xlab="Iterations",ylab=expression(psi))
plot(llik,type="l",col="sienna4",xlab="Iterations",ylab="log-likelihood")

pst=matrix(1,ncol=(p+1),nrow=W);pst[,2:(p+1)]=psis
lame=apply(pst,1,polyroot)
plot((1/lame)[Mod(lame)>1],col="gold",cex=.3,xlab=expression(Re(lambda)),ylab=expression(Im(lambda)))
lines(seq(-1,1,.01),sqrt(1-seq(-1,1,.01)^2),col="sienna",lty=2,lwd=2)
lines(seq(-1,1,.01),-sqrt(1-seq(-1,1,.01)^2),col="sienna",lty=2,lwd=2)

pses=apply(psis,2,mean)
mimine=mean(mus)
predo=mimine-pses[1]*(x[p:(T-1)]-mimine)
for (i in 2:p) predo=predo-pses[i]*(x[(p-i+1):(T-i)]-mimine)
plot(x[(p+1):T],type="l",col="steelblue4",xlab="t",ylab="x")
lines(predo,lty=2,col="sienna4",lwd=1.8)

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

# R PROGRAM TO ESTIMATE THE COEFFICIENTS OF AN AR(q) MODEL

p=9
q=p
T=550

x=matrix(scan("Dyndata"),ncol=5,byrow=T)[1:(T+50),5] # observations
x=x[1:T]
W=10000

mu=mean(x)
varef=var(x)
sig2=varef

lambdareal=2*runif(p)-1
preal=p
pcomp=0
lambdacomp=0
eps=rnorm(p,sd=sqrt(sig2)) #Just about anythin'

llog=function(pr,pc,lr,lc,compsi=T,pepsi=0,eps)
{

# LIKELIHOOD REPRESENTATION

if (compsi)
{

Psi=matrix(0,ncol=p,nrow=p+1)
Psi[1,]=1

if (pr>0)
{
Psi[2,1]=-lr[1]
if (pr>1)
{
for (i in 2:pr) Psi[2:(i+1),i]=Psi[2:(i+1),i-1]-lr[i]*Psi[1:i,i-1]
}
}
if (pc>0)
{
if (pr>0)
{
Psi[2,pr+2]=-2*lc[1]+Psi[2,pr]
Psi[3:(pr+3),pr+2]=(lc[1]^2+lc[2]^2)*Psi[1:(pr+1),pr]-2*lc[1]*Psi[2:(pr+2),pr]+Psi[3:(pr+3),pr]
}
else
{
Psi[2,2]=-2*lc[1];
Psi[3,2]=(lc[1]^2+lc[2]^2);
}
if (pc>2)
{
for (i in seq(4,pc,2))
{
pri=pr+i
prim=pri-2
Psi[2,pri]=-2*lc[i-1]+Psi[2,prim]
Psi[3:(pri+1),pri]=(lc[i-1]^2+lc[i]^2)*Psi[1:(pri-1),prim]-2*lc[i-1]*Psi[2:pri,prim]+Psi[3:(pri+1),prim]
}
}
}
Psi=Psi[2:(p+1),p]
}
else
{
Psi=pepsi
}
x=x-mu

heps=rep(0,T+p)
heps[1:p]=eps
for (i in 1:T)
heps[p+i]=x[i]+sum(rev(Psi)*heps[i:(p+i-1)])

loglike=-((sum(heps^2)/sig2)+(T+p)*log(sig2))/2

x=x+mu

list(ll=loglike,ps=Psi)
}

llo=llog(pr=preal,pc=pcomp,lr=lambdareal,lc=lambdacomp,eps=eps)

# MCMC REVERSIBLE SCHEME

psis=matrix(0,ncol=p,nrow=W)
mus=rep(0,W)
sigs=rep(0,W)
ncomp=rep(0,W)
llik=rep(0,W)
indacpt=rep(0,10)
epsrec=matrix(0,ncol=p,nrow=W)
preds=matrix(0,ncol=p,nrow=W)

for (m in 1:W){

  ind=sample(1:p,1)

  if (ind<=pcomp){
     ind=ind-(ind%%2==0)
     ppropreal=preal
     ppropcomp=pcomp
     lambpropreal=lambdareal
     lambpropcomp=lambdacomp
     lambpropcomp[ind]=lambdacomp[ind]+.05*rnorm(1)
     lambpropcomp[ind+1]=lambdacomp[ind+1]+.05*rnorm(1)

     }else{

      ppropreal=preal
      ppropcomp=pcomp
      lambpropreal=lambdareal
      lambpropcomp=lambdacomp
      lambpropreal[ind-pcomp]=lambdareal[ind-pcomp]+.05*rnorm(1)
      }

  lloprop=llog(pr=ppropreal,pc=ppropcomp,lr=lambpropreal,lc=lambpropcomp,eps=eps)

  if (log(runif(1))<lloprop$ll-llo$ll){

    llo=lloprop
    preal=ppropreal
    pcomp=ppropcomp
    lambdacomp=lambpropcomp
    lambdareal=lambpropreal
    indacpt[1]=indacpt[1]+1
    }

   if (preal<2){

       ppropreal=preal+2
       ppropcomp=pcomp-2

       ind=sample(1:pcomp,1)
       ind=ind-(ind%%2==0) 

       if (preal==0){
         lambpropreal=2*runif(2)-1
       }else{

	 lambpropreal=c(lambdareal,lambdacomp[ind]+0.01*rnorm(2))
       }

       lambpropcomp=lambdacomp[((1:pcomp)!=ind)&((1:pcomp)!=(ind+1))]

       coef=(1/2) * (pi/4)
       }

   if (pcomp==0){

       ppropreal=p-2
       ppropcomp=2

       ind=sample(1:p,2) 
       lambpropreal=lambdareal[((1:p)!=ind[1])&((1:p)!=ind[2])] 

       #theta=2*pi*runif(1)
       #rho=sqrt(runif(1)) 
       lambpropcomp=c(mean(lambdareal[ind])+.01*rnorm(1),.01*rnorm(1)) #rho*c(cos(theta),sin(theta))

       coef=(1/2) * (4/pi)
       }

   if ((preal>1)&&(pcomp>0)){
      
       if (runif(1)<.1){
       
         ppropcomp=pcomp-2
         ppropreal=preal+2

         ind=sample(1:pcomp,1)
         ind=ind-(ind%%2==0) 

         if (ppropcomp>0){
           lambpropcomp=lambdacomp[((1:pcomp)!=ind)&((1:pcomp)!=(ind+1))]
         }else{
           lambpropcomp=0
         }

         lambpropreal=c(lambdareal,lambdacomp[ind]+0.05*rnorm(2))

         coef=9*(1+(preal<2)) * (pi/4)

       }else{

         ppropreal=preal-2
         ppropcomp=pcomp+2

         ind=sample(1:preal,2) 	#real roots to remove
                                                                                                                                         
         if (ppropreal>0){

           lambpropreal=lambdareal[((1:preal)!=ind[1])&((1:preal)!=ind[2])]
	   }else{
          
              lambpropreal=0 
           }
	
        lambpropcomp=c(lambdacomp,mean(lambdareal[ind])+.05*rnorm(1),.05*rnorm(1))

        coef=(4/pi) * (1+(ppropcomp<p-1))/9.
      }
    }

  lloprop=llog(pr=ppropreal,pc=ppropcomp,lr=lambpropreal,lc=lambpropcomp,eps=eps)

  if (log(runif(1))<log(coef)+lloprop$ll-llo$ll){

    llo=lloprop
    preal=ppropreal
    pcomp=ppropcomp
    lambdacomp=lambpropcomp
    lambdareal=lambpropreal
    indacpt[2]=indacpt[2]+1
    }

  psis[m,]=llo$ps

  muold=mu
  mu=rnorm(1,mean=muold,sd=sqrt(5*sig2))
    
  lloprop=llog(pr=preal,pc=pcomp,lr=lambdareal,lc=lambdacomp,compsi=F,pepsi=llo$ps,eps=eps)

  if (log(runif(1))>lloprop$ll-llo$ll){
 
    indacpt[3]=1+indacpt[3]
    mu=muold
    }else{

       llo=lloprop
    }

  if (runif(1)<.3){

  Psi=llo$ps
  heps=rep(0,2*p)
  heps[1:p]=eps  
  for (j in (p+1):(2*p))
      heps[j]=x[j]-sum(Psi*heps[(j-p):(j-1)])
  varheps=(2*p)*var(heps)

  sig2old=sig2
  sig2=varheps/rgamma(1,2*p)
  difdens=dgamma(sig2/varheps,2*p,log=T)-dgamma(sig2old/varheps,2*p,log=T)+log(sig2)-log(sig2old)

  lloprop=llog(pr=preal,pc=pcomp,lr=lambdareal,lc=lambdacomp,compsi=F,pepsi=llo$ps,eps=eps)

  if (log(runif(1))>lloprop$ll-llo$ll-difdens){		

    sig2=sig2old
    }else{

        indacpt[4]=3+indacpt[4]
        llo=lloprop
     }
  }

  if (runif(1)<.3){

  sig2old=sig2
  sig2=exp(rnorm(1,mean=log(sig2),sd=sqrt(.1*varef)))

  lloprop=llog(pr=preal,pc=pcomp,lr=lambdareal,lc=lambdacomp,compsi=F,pepsi=llo$ps,eps=eps)
  difdens=log(sig2)-log(sig2old)

  if (log(runif(1))>lloprop$ll-llo$ll-difdens){

     sig2=sig2old
     }else{

        llo=lloprop
        indacpt[5]=3+indacpt[5]
     }
  }

  if (runif(1)<.3){

  sig2old=sig2
  thismean=varef/(1+sum(Psi^2))
  sig2=exp(rnorm(1,mean=thismean,sd=sqrt(varef)))

  difdens=dnorm(log(sig2),mean=thismean,sd=sqrt(.5*varef),log=T)-dnorm(log(sig2old),mean=thismean,sd=sqrt(.5*varef),log=T)

  lloprop=llog(pr=preal,pc=pcomp,lr=lambdareal,lc=lambdacomp,compsi=F,pepsi=llo$ps,eps=eps)

  if (log(runif(1))>lloprop$ll-llo$ll-difdens){

     sig2=sig2old
     }else{

        llo=lloprop
        indacpt[6]=3+indacpt[6]
     }
  }

  Psi=llo$ps

  if (runif(1)<.5){

  heps=rep(0,2*p+1)
  keps=rep(0,p)


  for (i in 1:q){

    x = x-mu
    heps[1:p]=eps
    for (j in (p+1):(2*p+1))
      heps[j]=x[j]+sum(rev(Psi)*heps[(j-p):(j-1)])

    heps[i]=0
    for (j in 1:(q-i+1))
      keps[j]=x[j]+sum(rev(Psi)*heps[j:(j+p-1)])
    x = x+mu

    epsvar = 1/sum(c(1,Psi[i:q]^2))
    epsmean = sum(Psi[i:q]*keps[1:(q-i+1)])*epsvar
    epsmean = epsmean/epsvar
    epsvar = sig2*epsvar

    propeps = rnorm(1,mean=epsmean,sd=sqrt(epsvar))
    epspr=eps
    epspr[i]=propeps
    lloprop=llog(pr=preal,pc=pcomp,lr=lambdareal,lc=lambdacomp,compsi=F,pepsi=Psi,eps=epspr)
    propsal1=dnorm(propeps,mean=epsmean,sd=sqrt(epsvar),log=T)

    x = x-mu
    heps[i]=propeps
    for (j in (p+1):(2*p+1))
      heps[j]=x[j]+sum(rev(Psi)*heps[(j-p):(j-1)])

    heps[i]=0
    for (j in 1:(q-i+1))
      keps[j]=x[j]+sum(rev(Psi)*heps[j:(j+p-1)])
    x = x+mu

    epsvar = 1/sum(c(1,Psi[i:q]^2))
    epsmean = sum(Psi[i:q]*keps[1:(q-i+1)])
    epsmean = epsmean*epsvar
    epsvar = sig2*epsvar 
    propsal0=dnorm(eps[i],mean=epsmean,sd=sqrt(epsvar),log=T)

    if (log(runif(1))<lloprop$ll-llo$ll-propsal1+propsal0){

      indacpt[7]=(2/p)+indacpt[7]
      eps[i]=propeps;
      llo=lloprop
      }
  }
  }

  if (runif(1)<.5){

  heps=rep(0,2*p+1)
  keps=rep(0,p)

  for (i in 1:q){ 

    propeps = rnorm(1,mean=eps[i],sd=0.1*sqrt(sig2))

    heps[1:p]=eps
    heps[i]=propeps
    keps=heps[1:p]

    allop=llog(pr=preal,pc=pcomp,lr=lambdareal,lc=lambdacomp,compsi=F,pepsi=Psi,eps=keps)

    if (log(runif(1))<allop$ll-llo$ll){ 

      indacpt[8]=(2/p)+indacpt[8]
      eps[i]=propeps;
      llo=allop
      }
  }
  }

  psis[m,]=llo$ps
  mus[m]=mu
  sigs[m]=sig2
  llik[m]=llo$ll
  ncomp[m]=pcomp
  epsrec[m,]=eps

  heps=rep(0,T+2*p)
  heps[1:p]=eps
  hatx=heps
  hatx[1:T]=x-mu

  for (j in (p+1):(p+T))
      heps[j]=x[j-p]+sum(rev(Psi)*heps[(j-p):(j-1)])
  for (j in (p+T+1):(p+T+p)){
      hatx[j-p]=-sum(rev(Psi)*heps[(j-p):(j-1)])
      }

  preds[m,]=hatx[(T+1):(T+p)]+mu
  }

indacpt/W

# GRAPHICAL REPRESENTATION

sube=function(x,n)
{
y=x
if (is.matrix(x))
{
t=dim(x)
if (t[1]>1000) y=y[seq(1,t[1],length=1000),]
if (t[2]>1000) y=y[,seq(1,t[2],length=1000)]
}
else y=y[seq(1,length(x),length=1000)]
y
}


start=1000

par(mfrow=c(3,3),mar=c(4,4,2,1))
if (p>3)
{
hist(ncomp,main="",xlab="p",ylab="",col="gold4",breaks=seq(-1,p+1,2));par(new=T);
plot(sube(1:W),sube(ncomp),axes=F,cex=.3,xlab="",ylab="")
axis(side=4)
}
if (p<4)
{
plot(sube(1:W),sube(ncomp),cex=.3,xlab="Iterations",ylab="Complex roots")
}
plot(sube(start:W),sube(mus[start:W]),type="l",col="steelblue4",xlab="Iterations",ylab=expression(mu))
plot(sube(start:W),sube(sigs[start:W]),type="l",col="steelblue4",xlab="Iterations",ylab=expression(sigma^2))
for (i in 1:min(3)) plot(sube(start:W),sube(psis[start:W,i]),type="l",col="steelblue4",xlab="Iterations",ylab=expression(psi))
plot(sube(start:W),sube(llik[start:W]),type="l",col="sienna4",xlab="Iterations",ylab="log-likelihood")

pst=matrix(1,ncol=(p+1),nrow=W);pst[,2:(p+1)]=-psis
lame=apply(pst,1,polyroot)
plot(sube((1/lame)[Mod(lame)>1]),col="gold",cex=.3,xlab=expression(Re(lambda)),ylab=expression(Im(lambda)))
lines(seq(-1,1,.01),sqrt(1-seq(-1,1,.01)^2),col="sienna",lty=2,lwd=2)
lines(seq(-1,1,.01),-sqrt(1-seq(-1,1,.01)^2),col="sienna",lty=2,lwd=2)

plot(sube(start:W),sube(epsrec[start:W,1]),ylim=range(epsrec[start:W,]),type="l",ylab=expression(epsilon),col="steelblue")
for (i in 2:p) lines(sube(start:W),sube(epsrec[start:W,i]),col="steelblue")