##### R-code for making Figure 4.9. It reads data from
##### "example48.dat"

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)

k1=1.848
h1=7.416

k2=0.462
h2=-2.446

cn1 = rep(0,N)
cn2 = rep(0,N)
cn1[1] = max(0,x1[1]^2-k1)
cn2[1] = min(0,x2[1]^2-k2)

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

ii <- seq(1,N)

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

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

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,100),cex=0.8)
lines(ii,rep(h1,N),lty=2,cex=0.8)
title(xlab="(a)",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(-12,0),cex=0.8)
lines(ii,rep(h2,N),lty=2,cex=0.8)
title(xlab="(b)",cex=0.9)

graphics.off()