# This code computes the ARL0 values of the upward CUSUM chart
# when the data are auto-correlated and follow the AR(1) model
# and when they are grouped into batches with batch size m. It
# is used in Table 4.4.

# Inputs
# k --- allowance
# h --- control limit (k and h are chosen such that a nominal
#       ARL0 value is reached).
# phi --- the coefficient in the AR(1) mode
#         (X(i)-mu0)=phi*(X(i-1)-mu0)+epsilon(i)
# m --- batch size

arl0 = function(k,h,phi,m){

set.seed(100)

rep=10000 # number of replicated simulations
N=10000    # the longest length of the observation sequence

AARL0=0   # AARL0 denotes the actual ARL0 value
num=0     # count of simulations with OC signals

mu0=0     # mu0 is the IC mean

for(j1 in 1:rep){     # the j1-th replicated simulation

x=rep(0,N)           # the original observation sequence
y=rep(0,N/m)           # the sequence of the batch means
cn=rep(0,N/m)          # the CUSUM charting statistic

x[1]=mu0+rnorm(1)
for(i in 2:m){ x[i]=mu0+phi*(x[i-1]-mu0)+rnorm(1) }
y[1]=mean(x[1:m])
cn[1]=max(0,y[1]*sqrt(m)-k)

if(cn[1]>h)
{
   AARL0=AARL0+1
   num=num+1
   next
}

for(i in 2:(N/m)){       # monitoring at the i-th batch

   for(j in 1:m){
       x[(i-1)*m+j]=mu0+phi*(x[(i-1)*m+j-1]-mu0)+rnorm(1)  
   }

   y[i]=mean(x[((i-1)*m+1):(i*m)])
   cn[i]=max(0,cn[i-1]+y[i]*sqrt(m)-k)

if(cn[i]>h)
{
   AARL0=AARL0+i
   num=num+1
   break
}

}                    # end of the j1-th replicated simulation

}                    # end of all "rep" simulations

AARL0=AARL0/num

}


k=0.25
h=5.597   # nominal ARL0=200 in this case

#k=0.5
#h=3.502   # nominal ARL0=200 in this case

#k=0.75
#h=2.481   # nominal ARL0=200 in this case

# phi is the lag-1 auto-correlation coefficient

m=5
phi=c(-0.5,-0.25,-0.1,0,0.1,0.25,0.5)  
ARL0=rep(0,7)

#phi=c(0.1,0.25)
#ARL0=rep(0,2)

for(i in 1:7){
   print(i)
   ARL0[i]=arl0(k,h,phi[i],m)
   print(ARL0[i])
}