#####   R-code for making Figure 5.8. It also generates
#####   the data file "example55.dat" used in Example 5.5

set.seed(10)

x1 = c(rnorm(50,0,1),rnorm(50,0,2))
x2 = c(rnorm(50,0,1),rnorm(50,0,0.5))
write.table(cbind(x1,x2),"example55.dat",row.names=F,col.names=F)

dat = matrix(scan("example55.dat"),ncol=2,byrow=T)
x1 = dat[,1]
x2 = dat[,2]
N = length(x1)

la=0.2
rho=2.635

En1 = rep(0,N)
En2 = rep(0,N)
En1[1] = la*x1[1]+(1-la)*0
En2[1] = la*x2[1]+(1-la)*0

for(i in 2:N){
   En1[i] <- la*x1[i]+(1-la)*En1[i-1]
   En2[i] <- la*x2[i]+(1-la)*En2[i-1]
}

ii <- seq(1,N)

crit = rho*sqrt(la/(2-la))*1

postscript("fig58.ps",width=6.5,height=6.5,horizontal=F)

par(mfrow=c(2,2), mar=c(4,4,2,2))

plot(ii,x1,type="o",lty=1,pch=16,xlab="n",
     ylab=expression(X[n]),mgp=c(2,1,0),xlim=c(0,N), 
     ylim=c(-4.5,4.5), cex=0.8)
title(xlab="(a)",cex=0.9)

plot(ii,En1,type="o",lty=1,pch=16,xlab="n",
     ylab=expression(E[n]),mgp=c(2,1,0),xlim=c(0,N), 
     ylim=c(-1.5,1.5),cex=0.8)
lines(ii,rep(crit,N),lty=2,cex=0.8)
lines(ii,rep(-crit,N),lty=2,cex=0.8)
title(xlab="(b)",cex=0.9)

plot(ii,x2,type="o",lty=1,pch=16,xlab="n",
     ylab=expression(X[n]),mgp=c(2,1,0),xlim=c(0,N), 
     ylim=c(-4.5,4.5), cex=0.8)
title(xlab="(c)",cex=0.9)

plot(ii,En2,type="o",lty=1,pch=16,xlab="n",
     ylab=expression(E[n]),mgp=c(2,1,0),xlim=c(0,N), 
     ylim=c(-1.5,1.5),cex=0.8)
lines(ii,rep(crit,N),lty=2,cex=0.8)
lines(ii,rep(-crit,N),lty=2,cex=0.8)
title(xlab="(d)",cex=0.9)

graphics.off()