##### R-code for Example 6.2. It writes the related data to
##### the file "example62.dat." It also creates "fig61.ps."

set.seed(100)

x1 = rnorm(10,0,0.5)
x2 = rnorm(10,1,0.5)
x3 = rnorm(10,0,0.5)

x =c(x1,x2,x3)
n=length(x)

### Compute F(1,m)

Q0 = rep(0,n)                   ## Q0[i] actually denotes Q(0,i)
for (i in 2:n){
   Q0[i] = var(x[1:i])*(i-1)
}

F1 = Q0

### Compute F(2,m)

F2 = rep(0,n)
temp2 = matrix(0,n,n)

for (i in 2:n) {

   for (j in 1:(i-2)) {
      temp2[i,j]=F1[j]+var(x[(j+1):i])*(i-j)
   }
   temp2[i,i-1]=F1[i-1]

   F2[i] = min(temp2[i,1:(i-1)])
}

### Compute F(3,m)

F3 = rep(0,n)
temp3 = matrix(0,n,n)

for (i in 3:n) {

   for (j in 1:(i-2)) {
      temp3[i,j]=F2[j]+var(x[(j+1):i])*(i-j)
   }
   temp3[i,i-1]=F2[i-1]

   F3[i] = min(temp3[i,1:(i-1)])
}

i = 1:n

write.table(round(cbind(i,x,temp3[n,],temp2[20,]),digits=3), 
            file="example62.dat",col.names=T, row.names=F)

### Make a figure

postscript("fig61.ps",width=7.2,height=3.6,horizontal=F)

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

plot(i,x,type="o",lty=1,pch=16,xlab="i",
     ylab=expression(X[i]),mgp=c(2,1,0),xlim=c(0,length(x)), 
     ylim=c(-1,2.5), cex=0.8)
title(xlab="(a)",cex=0.9)

plot(i[1:29],temp3[n,1:29],type="l",lty=1,xlab="h",
     ylab="",mgp=c(2,1,0),xlim=c(0,n), 
     ylim=c(0,13),cex=0.8)
lines(i[1:19],temp2[20,1:19],lty=2,cex=0.8)
title(xlab="(b)",cex=0.9)

graphics.off()