#####   R-code for making Figure 5.10. It also generates
#####   the data file "example56.dat" used in Example 5.6

set.seed(1000)

x1 = c(rnorm(50,0,1),rnorm(50,1,1))
x2 = c(rnorm(50,0,1),rnorm(50,0,2))
x3 = c(rnorm(50,0,1),rnorm(50,1,2))
x4 = c(rnorm(50,0,1),rnorm(50,0,0.5))

write.table(cbind(x1,x2,x3,x4),"example56.dat",row.names=F,col.names=F)


la=0.1
rho=2.595

N = length(x1)

En1 = rep(0,N)
En2 = rep(0,N)
En3 = rep(0,N)
En4 = rep(0,N)
En1[1] = la*x1[1]^2+(1-la)*1
En2[1] = la*x2[1]^2+(1-la)*1
En3[1] = la*x3[1]^2+(1-la)*1
En4[1] = la*x4[1]^2+(1-la)*1

for(i in 2:N){
   En1[i] = la*x1[i]^2+(1-la)*En1[i-1]
   En2[i] = la*x2[i]^2+(1-la)*En2[i-1]
   En3[i] = la*x3[i]^2+(1-la)*En3[i-1]
   En4[i] = la*x4[i]^2+(1-la)*En4[i-1]
}

U = 1+rho*sqrt(2*la/(2-la))

ii <- seq(1,N)

postscript("fig510.ps",width=5.5,height=8.5,horizontal=F)

par(mfrow=c(4,2), mar=c(4,4,1,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(-5,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(0,9),cex=0.8)
lines(ii,rep(U,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(-5,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(0,9),cex=0.8)
lines(ii,rep(U,N),lty=2,cex=0.8)
title(xlab="(d)",cex=0.9)

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

plot(ii,En3,type="o",lty=1,pch=16,xlab="n",
     ylab=expression(E[n]),mgp=c(2,1,0),xlim=c(0,N), 
     ylim=c(0,9),cex=0.8)
lines(ii,rep(U,N),lty=2,cex=0.8)
title(xlab="(f)",cex=0.9)

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

plot(ii,En4,type="o",lty=1,pch=16,xlab="n",
     ylab=expression(E[n]),mgp=c(2,1,0),xlim=c(0,N), 
     ylim=c(0,6),cex=0.8)
lines(ii,rep(U,N),lty=2,cex=0.8)
title(xlab="(h)",cex=0.9)

graphics.off()