##### R-code for computing ARL_0 values of the two-sided
##### EWMA chart consided in Table 5.5.

# This code computes the actual ARL0 values of the EWMA chart
# when the actual process distribution is N(0,1) and the parameters
# are estimated.

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


arl0 = function(la,rho,mu1,sigma1){

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


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 = rnorm(N,mu1,sigma1) 

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

if(abs(en[1])>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])>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.1
rho=2.454   # nominal ARL0=200 in this case for the one-sided EWMA

mu0=c(-0.1,-0.1,-0.1,0,0,0,0.1,0.1,0.1)
sigma0=c(0.9,1.0,1.1,0.9,1.0,1.1,0.9,1.0,1.1)
mu1=-mu0/sigma0
sigma1=1/sigma0

ARL0=rep(0,9)

for(i in 1:9){
   print(i)
   ARL0[i]=arl0(la,rho,mu1[i],sigma1[i])   
   print(ARL0[i])
}