##### R-code for making Figure 4.10. It also generates the data 
##### file "example49.dat" used in Example 4.9

set.seed(10)

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

N = length(x1)

k1=0.5
h1=0.881

k2=1.848
h2=2.533

postscript("fig410.ps",width=5,height=7.5,horizontal=F)

par(mfcol=c(3,2), mar=c(4,4,2,2))

x1bar = rep(0,(N/5))
x1sn2 = rep(0,(N/5))

for(i in 1:(N/5)){
   x1bar[i] = mean(x1[((i-1)*5+1):((i-1)*5+5)])
   x1sn2[i] = sd(x1[((i-1)*5+1):((i-1)*5+5)])^2
}

cn1 = rep(0,(N/5))
cn2 = rep(0,(N/5))

cn1[1] = max(0,x1bar[1]-k1)
cn2[1] = max(0,x1sn2[1]-k2)

for(i in 2:(N/5)){
   cn1[i] <- max(0,cn1[i-1]+x1bar[i]-k1)
   cn2[i] <- max(0,cn2[i-1]+x1sn2[i]-k2)
}

ii <- seq(1,(N/5))

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

plot(ii,cn1,type="o",lty=1,pch=16,xlab="n",
     ylab=expression(C[n]^{"+"}),mgp=c(2,1,0),xlim=c(0,(N/5)), 
     ylim=c(0,6),cex=0.8)
lines(ii,rep(h1,(N/5)),lty=2,cex=0.8)
title(xlab="(b)",cex=0.9)

plot(ii,cn2,type="o",lty=1,pch=16,xlab="n",
     ylab=expression(C[n]^{"+"}),mgp=c(2,1,0),xlim=c(0,(N/5)), 
     ylim=c(0,16.5),cex=0.8)
lines(ii,rep(h2,(N/5)),lty=2,cex=0.8)
title(xlab="(c)",cex=0.9)

x2bar = rep(0,(N/5))
x2sn2 = rep(0,(N/5))

for(i in 1:(N/5)){
   x2bar[i] = mean(x2[((i-1)*5+1):((i-1)*5+5)])
   x2sn2[i] = sd(x2[((i-1)*5+1):((i-1)*5+5)])^2
}

cn1 = rep(0,(N/5))
cn2 = rep(0,(N/5))

cn1[1] = max(0,x2bar[1]-k1)
cn2[1] = max(0,x2sn2[1]-k2)

for(i in 2:(N/5)){
   cn1[i] <- max(0,cn1[i-1]+x2bar[i]-k1)
   cn2[i] <- max(0,cn2[i-1]+x2sn2[i]-k2)
}

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

plot(ii,cn1,type="o",lty=1,pch=16,xlab="n",
     ylab=expression(C[n]^{"+"}),mgp=c(2,1,0),xlim=c(0,(N/5)), 
     ylim=c(0,6),cex=0.8)
lines(ii,rep(h1,(N/5)),lty=2,cex=0.8)
title(xlab="(e)",cex=0.9)

plot(ii,cn2,type="o",lty=1,pch=16,xlab="n",
     ylab=expression(C[n]^{"+"}),mgp=c(2,1,0),xlim=c(0,(N/5)), 
     ylim=c(0,16.5),cex=0.8)
lines(ii,rep(h2,(N/5)),lty=2,cex=0.8)
title(xlab="(f)",cex=0.9)

graphics.off()