# This code compute the actual ARL0 value of a 2-sided EWMA chart for detecting
# process variance shifts. It is used in Example 5.6 and Figure 5.11

# Inputs
# la --- lambda value (i.e., the weighting parameter)
# rho1 --- the critical value of the upward chart
# rho2 --- the critical value of the downward chart

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

la=0.1
rho1=3.199
rho2=1.750

U=1+rho1*sqrt(2*la/(2-la))
L=1-rho2*sqrt(2*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)         # generate N(0,1) observations

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

if(en[1]>U | en[1]<L)
{
   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]^2+(1-la)*en[i-1]

if(en[i]>U | en[i]<L)
{
   AARL0=AARL0+i
   num=num+1
   break
}

}                    # end of the j-th replicated simulation

}                    # end of all "rep" simulations

if(num==0) AARL0=N else AARL0=AARL0/num