# This code computes the actual ARL0 values of the EWMA chart
# when the process observations are correlated. It is used in
# Example 5.3 and Table 5.2.

# Inputs
# la --- lambda value (i.e., the weighting parameter)
# rho --- the critical value chosen such that the pre-specified
#         ARL0 value is reached.
# phi --- the coefficient in the AR(1) model

arl0 = function(la,rho,phi){

set.seed(100)

rep=10000 # number of replicated simulations
N=2000    # 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

crit=rep(0,N)        # rho*sqrt(la/(1-la))*sigma
for (i in 1:N){
   crit[i]=rho*sqrt(la/(2-la))
}

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

x=rep(0,N)           # the observation sequence
en=rep(0,N)          # the EWMA charting statistic

x[1] = mu0+rnorm(1) # autocorrelated observations from AR(1) model
for(i in 2:N){
   x[i] = mu0+phi*(x[i-1]-mu0)+rnorm(1)
}

en[1]=la*x[1]+(1-la)*mu0

if(abs(en[1]-mu0)>crit[1])
{
   AARL0=AARL0+1
   num=num+1
   next
}

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

en[i]=la*x[i]+(1-la)*en[i-1]

if(abs(en[i]-mu0)>crit[i])
{
   AARL0=AARL0+i
   num=num+1
   break
}

}                    # end of the j-th replicated simulation

}                    # end of all "rep" simulations

AARL0=AARL0/num

}


#la=0.05
#rho=2.216   # nominal ARL0=200 in this case

#la=0.1
#rho=2.454   # nominal ARL0=200 in this case

la=0.2
rho=2.635   # nominal ARL0=200 in this case

phi=c(-0.5,-0.25,-0.1,0,0.1,0.25,0.5)
ARL0=rep(0,length(phi))

for(i in 1:length(phi)){
   print(i)
   ARL0[i]=arl0(la,rho,phi[i])   
   print(ARL0[i])
}