#####   R-code for making Figure 5.12. It also generates
#####   the data file "example57.dat" used in Example 5.7

set.seed(1000)

x1 = c(rnorm(250,0,1),rnorm(250,1,1))
x2 = c(rnorm(250,0,1),rnorm(250,0,2))
x3 = c(rnorm(250,0,1),rnorm(250,1,2))
x4 = c(rnorm(250,0,1),rnorm(250,0,0.5))

xbar1=rep(0,100)
xbar2=rep(0,100)
xbar3=rep(0,100)
xbar4=rep(0,100)

sn1=rep(0,100)
sn2=rep(0,100)
sn3=rep(0,100)
sn4=rep(0,100)

for(i in 1:100){
   xbar1[i] = mean(x1[((i-1)*5+1):(i*5)])
   sn1[i] = sd(x1[((i-1)*5+1):(i*5)])
   xbar2[i] = mean(x2[((i-1)*5+1):(i*5)])
   sn2[i] = sd(x2[((i-1)*5+1):(i*5)])
   xbar3[i] = mean(x3[((i-1)*5+1):(i*5)])
   sn3[i] = sd(x3[((i-1)*5+1):(i*5)])
   xbar4[i] = mean(x4[((i-1)*5+1):(i*5)])
   sn4[i] = sd(x4[((i-1)*5+1):(i*5)])
}

write.table(cbind(xbar1,xbar2,xbar3,xbar4,sn1,sn2,sn3,sn4),"example57.dat",
row.names=F,col.names=F)

### Define EWMA charts for detecting mean shifts

la=0.1
rho=2.731

N = length(xbar1)

EnM1 = rep(0,N)
EnM2 = rep(0,N)
EnM3 = rep(0,N)
EnM4 = rep(0,N)
EnM1[1] = la*xbar1[1]+(1-la)*0
EnM2[1] = la*xbar2[1]+(1-la)*0
EnM3[1] = la*xbar3[1]+(1-la)*0
EnM4[1] = la*xbar4[1]+(1-la)*0

for(i in 2:N){
   EnM1[i] = la*xbar1[i]+(1-la)*EnM1[i-1]
   EnM2[i] = la*xbar2[i]+(1-la)*EnM2[i-1]
   EnM3[i] = la*xbar3[i]+(1-la)*EnM3[i-1]
   EnM4[i] = la*xbar4[i]+(1-la)*EnM4[i-1]
}

U = 0+rho*sqrt(la/(2-la))
L = 0-rho*sqrt(la/(2-la))

### Define EWMA charts for detecting variance shifts

la=0.1
rho=2.836

EnV1 = rep(0,N)
EnV2 = rep(0,N)
EnV3 = rep(0,N)
EnV4 = rep(0,N)
EnV1[1] = la*sn1[1]^2+(1-la)*1
EnV2[1] = la*sn2[1]^2+(1-la)*1
EnV3[1] = la*sn3[1]^2+(1-la)*1
EnV4[1] = la*sn4[1]^2+(1-la)*1

for(i in 2:N){
   EnV1[i] = la*sn1[i]^2+(1-la)*EnV1[i-1]
   EnV2[i] = la*sn2[i]^2+(1-la)*EnV2[i-1]
   EnV3[i] = la*sn3[i]^2+(1-la)*EnV3[i-1]
   EnV4[i] = la*sn4[i]^2+(1-la)*EnV4[i-1]
}

UV = 1+rho*sqrt(2*la/(2-la)/4)

ii <- seq(1,N)

postscript("fig512.ps",width=7,height=8.5,horizontal=F)

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

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

plot(ii,EnM1,type="o",lty=1,pch=16,xlab="n",
     ylab=expression(E[n,M]),mgp=c(2,1,0),xlim=c(0,N), 
     ylim=c(-2,2),cex=0.8)
lines(ii,rep(U,N),lty=2,cex=0.8)
lines(ii,rep(L,N),lty=2,cex=0.8)
title(xlab="(b)",cex=0.9)

plot(ii,EnV1,type="o",lty=1,pch=16,xlab="n",
     ylab=expression(E[n,V]),mgp=c(2,1,0),xlim=c(0,N), 
     ylim=c(0,5.1),cex=0.8)
lines(ii,rep(UV,N),lty=2,cex=0.8)
title(xlab="(c)",cex=0.9)

plot(ii,xbar2,type="o",lty=1,pch=16,xlab="n",
     ylab=expression(X[n]),mgp=c(2,1,0),xlim=c(0,N), 
     ylim=c(-2,4.5), cex=0.8)
title(xlab="(d)",cex=0.9)

plot(ii,EnM2,type="o",lty=1,pch=16,xlab="n",
     ylab=expression(E[n,M]),mgp=c(2,1,0),xlim=c(0,N), 
     ylim=c(-2,2),cex=0.8)
lines(ii,rep(U,N),lty=2,cex=0.8)
lines(ii,rep(L,N),lty=2,cex=0.8)
title(xlab="(e)",cex=0.9)

plot(ii,EnV2,type="o",lty=1,pch=16,xlab="n",
     ylab=expression(E[n,V]),mgp=c(2,1,0),xlim=c(0,N), 
     ylim=c(0,5.1),cex=0.8)
lines(ii,rep(UV,N),lty=2,cex=0.8)
title(xlab="(f)",cex=0.9)

plot(ii,xbar3,type="o",lty=1,pch=16,xlab="n",
     ylab=expression(X[n]),mgp=c(2,1,0),xlim=c(0,N), 
     ylim=c(-2,4.5), cex=0.8)
title(xlab="(g)",cex=0.9)

plot(ii,EnM3,type="o",lty=1,pch=16,xlab="n",
     ylab=expression(E[n,M]),mgp=c(2,1,0),xlim=c(0,N), 
     ylim=c(-2,2),cex=0.8)
lines(ii,rep(U,N),lty=2,cex=0.8)
lines(ii,rep(L,N),lty=2,cex=0.8)
title(xlab="(h)",cex=0.9)

plot(ii,EnV3,type="o",lty=1,pch=16,xlab="n",
     ylab=expression(E[n,V]),mgp=c(2,1,0),xlim=c(0,N), 
     ylim=c(0,5.1),cex=0.8)
lines(ii,rep(UV,N),lty=2,cex=0.8)
title(xlab="(i)",cex=0.9)

plot(ii,xbar4,type="o",lty=1,pch=16,xlab="n",
     ylab=expression(X[n]),mgp=c(2,1,0),xlim=c(0,N), 
     ylim=c(-2,4.5), cex=0.8)
title(xlab="(j)",cex=0.9)

plot(ii,EnM4,type="o",lty=1,pch=16,xlab="n",
     ylab=expression(E[n,M]),mgp=c(2,1,0),xlim=c(0,N), 
     ylim=c(-2,2),cex=0.8)
lines(ii,rep(U,N),lty=2,cex=0.8)
lines(ii,rep(L,N),lty=2,cex=0.8)
title(xlab="(k)",cex=0.9)

plot(ii,EnV4,type="o",lty=1,pch=16,xlab="n",
     ylab=expression(E[n,V]),mgp=c(2,1,0),xlim=c(0,N), 
     ylim=c(0,5.1),cex=0.8)
lines(ii,rep(UV,N),lty=2,cex=0.8)
title(xlab="(l)",cex=0.9)

graphics.off()