##### R-code for Example 7.3. It writes the related data to 
##### the file "example73.dat." It also makes "fig76.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.8,0.5,0.8,1,0.8,0.5,0.8,1),3,3)  ### IC covariance matrix

x = rbind(mvrnorm(10,mu1,Sigma),mvrnorm(20,mu2,Sigma))

n = length(x[,1])
p = 3
k = 0.5                    # Specify the allowance constant

cnp=matrix(0,n,p)          # the CUSUM charting statistic Cn+
cnm=matrix(0,n,p)          # the CUSUM charting statistic Cn-
mcz=rep(0,n)               # Charting statistic MCZ
zno=rep(0,n)               # Charting statistic ZNO

### Transform x to z where z is the regression-adjusted residuals

z = t(diag((diag(solve(Sigma)))^(-0.5),p,p) %*% solve(Sigma) %*% (t(x)-mu1))

for(j1 in 1:p){
   cnp[1,j1]=max(0,z[1,j1]-k)
   cnm[1,j1]=min(0,z[1,j1]+k)
}

mcz[1] = max(c(cnp[1,],-cnm[1,]))
zno[1] = sum((cnp[1,]+cnm[1,])^2)

for(i in 2:n){       # monitoring at the i-th time point

for(j1 in 1:p){
   cnp[i,j1]=max(0,cnp[i-1,j1]+z[i,j1]-k)
   cnm[i,j1]=min(0,cnm[i-1,j1]+z[i,j1]+k)
}

mcz[i] = max(c(cnp[i,],-cnm[i,]))
zno[i] = sum((cnp[i,]+cnm[i,])^2)

}

write.table(round(cbind(x,mcz,zno),digits=3),"example73.dat",row.names=F)


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

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

ii = seq(1,n)

### Plots of X1, X2, X3

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(-2.5,3.5), cex=0.8)
title(xlab="(a)",cex=0.9)

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(-2.5,3.5), cex=0.8)
title(xlab="(b)",cex=0.9)

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(-2.5,3.5), cex=0.8)
title(xlab="(c)",cex=0.9)

### Plots of Z1, Z2, Z3

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

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

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

### Plots of three individual CUSUMs

plot(ii,cnp[,1],type="o",lty=1,pch=16,xlab="n",ylab=expression(C[list(n,1)]),
     mgp=c(2,1,0),xlim=c(0,n),ylim=c(-15,25),cex=0.8)
lines(ii,cnm[,1],type="o",lty=3,pch=1,cex=0.8)
lines(ii,rep(5.014,n),lty=2,cex=0.8)
lines(ii,rep(-5.014,n),lty=2,cex=0.8)
title(xlab="(g)",cex=0.9)

plot(ii,cnp[,2],type="o",lty=1,pch=16,xlab="n",ylab=expression(C[list(n,2)]),
     mgp=c(2,1,0),xlim=c(0,n),ylim=c(-15,25),cex=0.8)
lines(ii,cnm[,2],type="o",lty=3,pch=1,cex=0.8)
lines(ii,rep(5.014,n),lty=2,cex=0.8)
lines(ii,rep(-5.014,n),lty=2,cex=0.8)
title(xlab="(h)",cex=0.9)

plot(ii,cnp[,3],type="o",lty=1,pch=16,xlab="n",ylab=expression(C[list(n,3)]),
     mgp=c(2,1,0),xlim=c(0,n),ylim=c(-15,25),cex=0.8)
lines(ii,cnm[,3],type="o",lty=3,pch=1,cex=0.8)
lines(ii,rep(5.014,n),lty=2,cex=0.8)
lines(ii,rep(-5.014,n),lty=2,cex=0.8)
title(xlab="(i)",cex=0.9)

### Plots of MCZ and ZNO

plot(ii,mcz,type="o",lty=1,pch=16,xlab="n",ylab=expression(C[n]),mgp=c(2,1,0),
     xlim=c(0,n),ylim=c(0,25),cex=0.8)
lines(ii,rep(5.014,n),lty=2,cex=0.8)
title(xlab="(j)",cex=0.9)

plot(ii,zno,type="o",lty=1,pch=16,xlab="n",ylab=expression(widetilde(C)[n]),
     mgp=c(2,1,0),xlim=c(0,n),ylim=c(0,800),cex=0.8)
lines(ii,rep(42.031,n),lty=2,cex=0.8)
title(xlab="(k)",cex=0.9)

graphics.off()