########################################################################
#                                                                      #
# This code computes the ARL0 values of the regression-adjusted cusum  #
# chart, using the charting statistic MCZ.                             #
#                                                                      #
# Inputs                                                               #
# k --- allowance                                                      #
# p --- dimension of the process                                       #
# mu0 -- IC mean vector                                                #
# Sigma --- IC covariance matrix                                       #
#                                                                      #
########################################################################

### Note: In this code, we assume that the production process is
###       3-dimensional (i.e., p=3). It should be easy to modify
###       it for cases with another p value.


library(MASS) # The R library in which the command for generating multivariate
             # normal random vectors is contained.

arl0 = function(k,h,p,mu0,Sigma){

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


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

x = mvrnorm(N,mu0,Sigma)   # Generate N multivariate normal random vectors
cnp=matrix(0,N,p)          # the CUSUM charting statistic Cn+
cnm=matrix(0,N,p)          # the CUSUM charting statistic Cn-
mcz=rep(0,N)

### Transform x to z where z is the regression-adjusted residuals

z = t(diag((diag(solve(Sigma)))^(-0.5),p,p) %*% solve(Sigma) %*% (t(x)-mu0))

for(j1 in 1:p){
   cnp[1,j1]=max(0,z[1,j1]-k)
   cnm[1,j1]=min(0,z[1,j1]+k)
}

mcz[1]=max(c(cnp[1,],-cnm[1,]))

if(mcz[1]>h)
{
   AARL0=AARL0+1
   num=num+1
   next
}

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

for(j1 in 1:p){
   cnp[i,j1]=max(0,cnp[i-1,j1]+z[i,j1]-k)
   cnm[i,j1]=min(0,cnm[i-1,j1]+z[i,j1]+k)
}

mcz[i]=max(c(cnp[i,],-cnm[i,]))

if(mcz[i]>h)
{
   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=200          # You need to specify ARL0 and the allowance constant
k=0.5             # k here, and adjust the four numbers below.

p=3               # p is the dimension of the process
mu0 = rep(0,p)    # IC mean vector
Sigma = matrix(c(1,0.8,0.5,0.8,1,0.8,0.5,0.8,1),p,p) # IC covariance matrix

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

hL=5.0            # Initial value of the lower bound
hU=6.0            # Initial value of the upper bound

A0=arl0(k,hL,p,mu0,Sigma)
print(A0)
if(A0 >= ARL0){
   print("hL is too large")
}

A0=arl0(k,hU,p,mu0,Sigma)
print(A0)
if(A0 <= ARL0){
   print("hU is too small")
}


for(i in 1:100){

   print(i)
   ht = (hL+hU)/2
   A0 = arl0(k,ht,p,mu0,Sigma)

   if(abs(A0-ARL0) <= eps1 | abs(hU-ht) <= eps2){
      break
   }

   if(A0 < ARL0){
      hL=ht
   }

   if(A0 > ARL0){
      hU=ht
   }

}

#####

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