########################################################################
#                                                                      #
# This code searches the control limit h_SS of the self-starting MEWMA #
# chart (7.45)-(7.46).                                                 #
#                                                                      #
# 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(lambda,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
xbar=matrix(0,N,p)        # Xbar_n
Sn2=array(0,c(N,p,p))     # S_n^2
un=matrix(0,N,p)          # u_n
En=matrix(0,N,p)          # Self-starting charting statistic E_{n,SS}
Tn2=rep(0,N)              # T_{n,SS}^2

xbar[2,]=(x[1,]+x[2,])/2
Sn2[2,,]=(x[2,]-x[1,])%*%t(x[2,]-x[1,])/2

for(i in 3:4){
   xbar[i,]=((i-1)*xbar[i-1,]+x[i,])/i
   un[i,]=(x[i,]-xbar[i-1,])*sqrt((i-1)/i)
   Sn2[i,,]=Sn2[i-1,,]*(i-2)/(i-1)+(un[i,]%*%t(un[i,]))/(i-1)
}


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

xbar[i,]=((i-1)*xbar[i-1,]+x[i,])/i
un[i,]=(x[i,]-xbar[i-1,])*sqrt((i-1)/i)
Sn2[i,,]=Sn2[i-1,,]*(i-2)/(i-1)+(un[i,]%*%t(un[i,]))/(i-1)

En[i,]=lambda*un[i,]+(1-lambda)*En[i-1,]
Tn2[i]=(t(En[i,]) %*% solve(Sn2[i-1,,]) %*% En[i,])*
       (2-lambda)/(lambda*(1-(1-lambda)^(2*i)))

if(Tn2[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 smoothing parameter
lambda=0.2        # lambda 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
Sigma = matrix(c(1,0,0,0,1,0,0,0,1),p,p)

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

hL=13.0            # Initial value of the lower bound
hU=15.0            # Initial value of the upper bound

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

A0=arl0(lambda,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(lambda,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