########################################################################
#                                                                      #
# This code searches the control limit h of the CUSUM based on T^2.    #
#                                                                      #
# 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
Sn=rep(0,N)               # the CUSUM charting statistic Cn+

Sn[1]=max(0,t(x[1,]-mu0)%*%solve(Sigma)%*%(x[1,]-mu0)-k)

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

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

Sn[i]=max(0,Sn[i-1]+t(x[i,]-mu0)%*%solve(Sigma)%*%(x[i,]-mu0)-k)

if(Sn[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=4.0             # 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=10.0            # Initial value of the lower bound
hU=15.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:20){

   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