##### R-code for computing rho_U values in Table 5.4.

# This code searches the rho value of the upward EWMA chart for detecting
# process variance shifts such that a pre-specified ARL0 value is reached.

# Inputs
# la --- lambda value (i.e., the weighting parameter)
# rho --- the critical value chosen such that the pre-specified
#         ARL0 value is reached.
# df --- degree of freedoms of the chi-square distribution

arl0 = function(la,rho){

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

U=1+rho*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)
{
   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)
{
   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
}

}

#####

ARL0=1000          # You need to specify ARL0 and the lambda value "la"
la=0.1            # here, and adjust the four numbers below.

eps1=10^(-1)      # precision requirement of the actual ARL_0
eps2=10^(-4)      # precision requirement of searched rho

rhoL=5.5          # Initial value of the lower bound
rhoU=6.5          # Initial value of the upper bound

A0=arl0(la,rhoL)
if(A0 >= ARL0){
   print("rhoL is too large")
}

A0=arl0(la,rhoU)
if(A0 <= ARL0){
   print("rhoU is too small")
}

for(i in 1:100){

   print(i)
   rhot = (rhoL+rhoU)/2
   A0 = arl0(la,rhot)

   if(abs(A0-ARL0) <= eps1 | abs(rhoU-rhot) <= eps2){
      break
   }

   if(A0 < ARL0){
      rhoL=rhot
   }

   if(A0 > ARL0){
      rhoU=rhot
   }

}

#####

# After running this code, in R, just type
# > rhot     # obtain the searched rho value
# > A0       # find the actual ARL_0 value