# This code computes the actual ARL0 values of the EWMA chart
# when the actual process distribution is chi-square. It is
# used in Figure 5.5.

# 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,df){

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)[1-(1-la)^(2n)])*sigma
for (i in 1:N){
   crit[i]=rho*sqrt(la/(2-la))
#   crit[i]=rho*sqrt(la/(2-la)*(1-(1-la)^(2*i)))
}

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 = (rchisq(N,df)-df)/sqrt(2*df)      # standardized observations

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.01
#rho=1.500    # nominal ARL0=200 in this case

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

df=seq(1,10)
ARL0=rep(0,length(df))

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