# This code computes the ARL1 values of the upward CUSUM chart
# when the data are auto-correlated and follow the AR(1) model
# and the shift occurs at the intial time point. It is used in
# Example 4.5 and Table 4.3.

# 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)
# delta --- mean shift size

arl1 = function(k,h,phi,delta){

set.seed(100)

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

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

mu0=0         # mu0 is the IC mean
mu1=mu0+delta # mu1 is the OC mean

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

x=rep(0,N)           # the observation sequence
cn=rep(0,N)          # the CUSUM charting statistic

x[1]=mu1+rnorm(1)
cn[1]=max(0,x[1]-k)

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

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

x[i]=mu1+phi*(x[i-1]-mu1)+rnorm(1)  # generate the observations from the AR(1) model
cn[i]=max(0,cn[i-1]+x[i]-k)

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

}                    # end of the j-th replicated simulation

}                    # end of all "rep" simulations

AARL1=AARL1/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

phi=c(-0.5,-0.25,-0.1,0,0.1,0.25,0.5) 
delta=0.5 
ARL1=rep(0,7)

for(i in 1:7){
   ARL1[i]=arl1(k,h,phi[i],delta)
}