##### R-code for making Figure 4.8. It also generates the data 
##### file "example48.dat" used in Example 4.8

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),"example48.dat",row.names=F,col.names=F)

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

k=0.25

cn1 = rep(0,N)
cn2 = rep(0,N)
cn1[1] = max(0,x1[1]-k)
cn2[1] = max(0,x2[1]-k)

for(i in 2:N){
   cn1[i] <- max(0,cn1[i-1]+x1[i]-k)
   cn2[i] <- max(0,cn2[i-1]+x2[i]-k)
}

ii <- seq(1,N)

postscript("fig48.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,cn1,type="o",lty=1,pch=16,xlab="n",
     ylab=expression(C[n]^{"+"}),mgp=c(2,1,0),xlim=c(0,N), 
     ylim=c(0,13),cex=0.8)
lines(ii,rep(5.597,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,cn2,type="o",lty=1,pch=16,xlab="n",
     ylab=expression(C[n]^{"+"}),mgp=c(2,1,0),xlim=c(0,N), 
     ylim=c(0,13),cex=0.8)
lines(ii,rep(5.597,N),lty=2,cex=0.8)
title(xlab="(d)",cex=0.9)

graphics.off()