##### R-code for Example 7.2. It writes the related data to 
##### the file "example72.dat." It also makes "fig75.ps"

#### Generate 30 random vectors from N_3(mu,Sigma)

library(MASS)   ### R package for generating multivariate normal random vectors.

set.seed(100)

mu1 = c(0,0,0)  ### IC mean vector
mu2 = c(1,0,0)  ### OC mean vector

Sigma = matrix(c(1,0,0,0,1,0,0,0,1),3,3)  ### IC covariance matrix

x = mvrnorm(20,mu2,Sigma)

n = length(x[,1])

T0n2 = rep(0,n)

for(i in 1:n){
   T0n2[i] = t(x[i,]-mu1) %*% solve(Sigma) %*% (x[i,]-mu1)
}

write.table(round(cbind(x,T0n2),digits=3),"example72.dat",row.names=F)


postscript("fig75.ps",width=6,height=8,horizontal=F)

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

ii = seq(1,n)

### Plots of X1

plot(ii,x[,1],type="o",lty=1,pch=16,xlab="n",
     ylab=expression(X[n1]),mgp=c(2,1,0),xlim=c(0,n), 
     ylim=c(-3,4), cex=0.8)
title(xlab="(a)",cex=0.9)

k=0.5

cn1 = rep(0,n)
cn2 = rep(0,n)
cn1[1] = max(0,x[,1]-k)
cn2[1] = min(0,x[,1]+k)

for(i in 2:n){
   cn1[i] <- max(0,cn1[i-1]+x[i,1]-k)
   cn2[i] <- min(0,cn2[i-1]+x[i,1]+k)
}

plot(ii,cn1,type="o",lty=1,pch=16,xlab="n",ylab="",mgp=c(2,1,0),
     xlim=c(0,n),ylim=c(-5,15),cex=0.8)
lines(ii,cn2,type="o",lty=3,pch=1,cex=0.8)
lines(ii,rep(3.892,n),lty=2,cex=0.8)
lines(ii,rep(-3.892,n),lty=2,cex=0.8)
title(xlab="(b)",cex=0.9)

### Plots of X2

plot(ii,x[,2],type="o",lty=1,pch=16,xlab="n",
     ylab=expression(X[n2]),mgp=c(2,1,0),xlim=c(0,n), 
     ylim=c(-3,4), cex=0.8)
title(xlab="(c)",cex=0.9)

k=0.5

cn1 = rep(0,n)
cn2 = rep(0,n)
cn1[1] = max(0,x[,2]-k)
cn2[1] = min(0,x[,2]+k)

for(i in 2:n){
   cn1[i] <- max(0,cn1[i-1]+x[i,2]-k)
   cn2[i] <- min(0,cn2[i-1]+x[i,2]+k)
}

plot(ii,cn1,type="o",lty=1,pch=16,xlab="n",ylab="",mgp=c(2,1,0),
     xlim=c(0,n),ylim=c(-5,15),cex=0.8)
lines(ii,cn2,type="o",lty=3,pch=1,cex=0.8)
lines(ii,rep(3.892,n),lty=2,cex=0.8)
lines(ii,rep(-3.892,n),lty=2,cex=0.8)
title(xlab="(d)",cex=0.9)

### Plots of X3

plot(ii,x[,3],type="o",lty=1,pch=16,xlab="n",
     ylab=expression(X[n3]),mgp=c(2,1,0),xlim=c(0,n), 
     ylim=c(-3,4), cex=0.8)
title(xlab="(e)",cex=0.9)

k=0.5

cn1 = rep(0,n)
cn2 = rep(0,n)
cn1[1] = max(0,x[,3]-k)
cn2[1] = min(0,x[,3]+k)

for(i in 2:n){
   cn1[i] <- max(0,cn1[i-1]+x[i,3]-k)
   cn2[i] <- min(0,cn2[i-1]+x[i,3]+k)
}

plot(ii,cn1,type="o",lty=1,pch=16,xlab="n",ylab="",mgp=c(2,1,0),
     xlim=c(0,n),ylim=c(-5,15),cex=0.8)
lines(ii,cn2,type="o",lty=3,pch=1,cex=0.8)
lines(ii,rep(3.892,n),lty=2,cex=0.8)
lines(ii,rep(-3.892,n),lty=2,cex=0.8)
title(xlab="(f)",cex=0.9)

graphics.off()