##### R-code for making Figure 4.7. It reads data from
##### "soil.dat"

dat = matrix(scan("soil.dat"),ncol=6,byrow=T)

x =dat[,4]
xstd=(x-mean(x))/sd(x)

res = ar(x)$resid[2:length(x)]
y = (res-mean(res))/sd(res)

t=seq(1,length(x))
t1=seq(2,length(x))

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

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

plot(t,x,type="o",lty=1,pch=16,xlab="n",
     ylab=expression(X[n]),mgp=c(2,1,0),xlim=c(0,190), 
     ylim=c(10,15), cex=0.8)
title(xlab="(a)",cex=0.9)

plot(t1,res,type="o",lty=1,pch=16,xlab="n",
     ylab=expression(widehat(e)[n]),mgp=c(2,1,0),
     xlim=c(0,190), ylim=c(-2.1,1.5), cex=0.8)
title(xlab="(b)",cex=0.9)

k=0.25
h=5.597

cnx = rep(0,length(x))
cnres = rep(0,length(x)-1)

cnx[1] = max(0,xstd[1]-k)
cnres[1] = max(0,y[1]-k)

for(i in 2:length(x)){
   cnx[i] = max(0,cnx[i-1]+xstd[i]-k)
}

for(i in 2:(length(x)-1)){
   cnres[i] = max(0,cnres[i-1]+y[i]-k)
}

plot(t,cnx,type="o",lty=1,pch=16,xlab="n",
     ylab=expression(C[n]^{"+"}),mgp=c(2,1,0),
     xlim=c(0,length(t)), 
     ylim=c(0,20),cex=0.8)
lines(t,rep(5.597,length(t)),lty=2,cex=0.8)
title(xlab="(c)",cex=0.9)

plot(t1,cnres,type="o",lty=1,pch=16,xlab="n",
     ylab=expression(C[n,res]^{"+"}),mgp=c(2,1,0),
     xlim=c(0,length(t)), 
     ylim=c(0,20),cex=0.8)
lines(t,rep(5.597,length(t)),lty=2,cex=0.8)
title(xlab="(d)",cex=0.9)

graphics.off()