# Introduction ultivariate statistical process control is often used in chemical and process industries where autocorrelation is most prevalent. Traditional multivariate statistical process control techniques are based on the assumption that the successive observation vectors are independent. In recent years, due to automation of measurement and data collection systems, a process can be sampled at higher rates, which ultimately leads to autocorrelation. Consequently, when the autocorrelation is present in the data, it can have a serious impact on the performance of classical control charts. This point has been made by numerous authors, including Berthouex, Hunter, and Pallensen (1978), Harris and Ross (1991), Montgomery and Mastrangelo (1991). Runger (1996) has presented a realistic model that generates autocorrelation and cross correlation and provides a useful approach to characterizing process data. The interpretation of these charts: charts based on modeling residuals is not as simple as the authors suggest, and the alternative engineering feedback control methods are often more appropriate with such highly auto correlated data. This article considers the problem of monitoring the mean vector of a process in which observations can be a highly first order vector autoregressive VAR (1) (Hotelling, 1947) and MEWMA chart in a simple way and has ability to detect small as well as large shift as soon as possible as required by some industrial processes with high level of first order vector autoregressive data. MMOEWMA control statistic gives weight to the past observation vectors in slightly different way than MEWMA and each current change is considered with full weight. This corrects MEWMA statistic for suffering from inertia problem. This article discusses the procedures to construct the Multivariate Modified EWMA chart. Simulate the average run length to assess the performance of the chart. The MMOEWMA vector auto correlated control chart is defined in second section and the derivation of upper control limits is kept with Appendix 1. Further, performance of MMOEWMA monitoring scheme is illustrated along with MEWMA scheme for real multivariate chemical process data in third section. ARL properties of MMOEWMA are derived and compared with MEWMA in fourth section. The comparisons reveals that MMOEWMA scheme outperform MEWMA scheme. Computation of ARL values were carried out using Markov chain approach described in Appendix 2. # II. # Multivariate Control Charts for Monitoring the Process Mean Suppose that the p x 1 random vectors Y 1 , Y 2 , Y 3 , ? each representing the p quality characteristics to be monitored, are observed over time. These vectors may represent individual observations or sample mean vectors. To study the performance of the various multivariate control charts, it will be assumed that Y n , n = 1, 2, ?, are independent multivariate normal random vectors with mean vectors ? n , n =1, 2,? For simplicity, it is assumed that each of the random vectors has the known co-variance matrix ?. # a) The Multivariate Shewhart Control Chart Hotelling's (1947) introduced a multivariate control-chart procedure based on his Hotelling's chisquare statistics defined as 0 is a specified control limit. Because this procedure is based on only the most recent observation, it is insensitive to small and moderate shifts in the mean vector. n n n Y Y 1 ' 2 ? ? = ? . At any n th stage in process, n n n Y Y 1 ' 2 ? ? = ? > h,(1) b) The Multivariate EWMA (MEWMA) control chart Lowery et al. (1992) proposed the MEWMA chart as natural extension of EWMA chart. It is a popular chart used to monitor a process with p quality characteristics for detecting small to moderate shifts. The in-control process mean is assumed without loss of generality to be a vector of zeros, and covariance matrix ?. The MEWMA control statistic is defined as vectors, X n =?Y n + (I-?)X n-1, n=1,2,?,(2) where X 0 = 0, 1 x p vector and ? = diag(? 1 , ? 2 ,?, ? p ), 0 h 1 ,(3) where h 1 (>0) is chosen to achieve a specified in control ARL and ? xn is the covariance matrix of X n given by ? xn = {?/(2?)}?, under equality of weights of past observations for all p characteristics; ? 1 =? 2 =?= ? p =?. The UCL = 2 / 1 1 2 / 1 ) ( 2 h ? ? ? ? ? ? ? ? ? . If one or more points fall beyond h 1 , the process is assumed to be out-of-control. The magnitude of the shift is reflected in the noncentrality parameter ? 1 ' ? -1 ? 1. They conclude that an assignable causes result in a shift in the process mean from ? 0 to ? 1. c) MMOEWMA control chart for monitoring the first order vector autoregressive VAR (1) process mean The MMOEWMA chart as natural extension of Patel and Divecha (2011) proposed Modified EWMA chart. The desirable properties of a multivariate auto correlated control chart are that it is easy to implement and is effective for detecting shifts of all sizes as per technical specifications. The Multivariate Modified EWMA chart that introduce considers past observations similar to MEWMA scheme and additionally considers past as well as latest change in the process. Let Y n , n = 1,2,?, are sequence of first order auto correlated normal random vectors with mean vector ? n , and common covariance matrix ?. Further it is assumed without loss of generality that the in control process mean vector is ? 0 = (0, 0,?, 0)' = 0. The MMOEWMA chart statistic is a modification in MEWMA chart statistic. To define MMOEWMA control statistic as vector X n given by, X n = ?Y n + (I-?)X n-1 + (Y n -Y n-1 ), n?1,(4) where X 0 is the p-dimensional zero vector and ? = diag(? 1 , ? 2 ,?, ? p ), 0 h 2 , (5) where h 2 (>0) is chosen to achieve a specified in control ARL. If one or more points fall beyond h 2 , the process is assumed to be out-of-control. ? xn is the covariance matrix of X n given by ? xn = ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? 2 ) 1 (2 2 , under equality of weights of past observations for all p characteristics; ? 1 =? 2 =?= ? p =?, and past and current changes. The upper control limit of MMOEWMA chart is, 1 displays the part of measurements on three temperature column taken every minute from a chemical process that is working in control and out of control situations, abrupt changes and small shifts occur. Here number of variable p=3. Three variables are temperature columns Y 1 with mean ? 01 =92.24, Y 2 with mean ? 02 =95.56, Y 3 with mean ? 03 = 100.27. To assume covariance matrix to be UCL = 2 / 1 2 ) 1 ( 2 2 ? ? ? ? ? ? ? ? + ? ? ? ? ? ? 2 / 1 2 ) (h (discussed in Appendix 1). # III. Illustration ? ? ? ? ? ? ? ? ? ? ? = ) ( ) ,( ) , ( ) , ( ) ( ) , ( ) , ( ) , ( )( 33 2 3 1 3 3 2 22 1 2 3 1 2 1 11 Y V Y Y C Y Y C Y Y C Y V Y Y C Y Y C Y Y C Y V = ? ? ? ? ? ? ? ? ? ? ? ? ? ? 1 5 . 0 5 . 0 5 . 0 1 5 . 0 5 . 0 5 . 0 1 Upper control limit of MEWMA UCL = 2 / 1 1 2 / 1 ) ( 2 h ? ? ? ? ? ? ? ? ? =0.882, where ?=0.1, h 1 =14.7, and MMOEWMA UCL = 2 / 1 2 ) 1 ( 2 2 ? ? ? ? ? ? ? ? + ? ? ? ? ? ? 2 / 1 2 ) (h =0.839, where ?=0.1, h 2 = 5.417. UCL is used in average run length to choose appropriate value of decision interval. The MMOEWMA chart gives an out-of-control signal as soon as 2 2 n T = X n ' ? xn -1 X n > h 2, h 2 = 5.417 and the MEWMA chart gives an out-of-control signal as soon as All the ARL computations were carried out using Markov chain approach described in Appendix 2. MMOEWMA is the chart for multivariate processes having autocorrelated observations. However, assuming that MEWMA chart can be applied at least to the residual vectors, to compare the ARL values of MMOEWMA with that of MEWMA having common parameters. Table 3 to table 6 showed that the control limits for MMOEWMA are quite small than those of MEWMA as well as 2 ? chart for the same in control ARL. For two, three, four variable cases, the smaller out of control ARL imply that MMOEWMA chart performs excellent in detection of shifts, be it small, moderate or large for every of in control ARLs. MMOEWMA chart with ? 1 =? 2 =?= ? p =? =1 is multivariate Shewhart chart version for multivariate autocorrelated process. 2 1 n T = X n ' ? xn -1 X n > h 1 , h 1 = 14.78 . # Conclusion A simple multivariate control chart for monitoring small as well as large shifts in highly first order vector autoregressive VAR (1) process such as multivariate chemical process is given. It is good method to monitor first order vector autoregressive process in chemical/other industries. # Appendix 1 Derivation of the Covariance Matrix of MMOEWMA Statistic x n By repeated substitution in equation X n = ?Y n + (I -?)X n-1 + (Y n -Y n-1 ), n?1, it can be shown that X n = + 0 n Y ) (I ? ? + ? = ? ? 1 - n 0 j 1 - j - n j - n j ) Y - (Y ) (I (i) The expectation of X n gives, E(X n ) = ? , (mean of Y n ). Lemma 1: If ? 1 = ? 2 = ?= ? p = ?, then the expression for Multivariate Modified EWMA statistic X n = ? ? ? = 1 0 n j (1-?) j Y n-j + (1-?) n Y 0 + ? ? = 1 0 n j (1-?) j (Y n-j -Y n-j-1 ) Taking expectation on both side, E(X n ) = ? ? ? = 1 0 n j (1-?) j E(Y n-j )+ (1-?) n E(Y 0 )+ ? ? = 1 0 n j (1-?) j E(Y n-j -Y n-j-1 ) But ] ) 1 ( 1 [ )] 1 ( 1 [ ] ) 1 ( 1 [ ) 1 ( 0 n n j n j ? ? ? ? ? ? ? ? = ? ? ? ? = ? ? = ? E (X n ) = 0 ) 1 ( ] ) 1 ( 1 [ + ? + ? ? µ ? µ ? n n = ? Lemma 2 : If the starting value of process is, X 0 = µ 0 = Y 0 and 0< ? ? 1 is a constant. The mean is, E(X n ) = E((1-?)X n-1 + ?Y n + (Y n -Y n-1 )) = µ 0 . For p = 2 and n=1, we have, X 1 = ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ? ? ? ?Y Y Y Y Y Y Y Y ? ? ? ? = ? ? ? ? ? ? 21 11 X X say, So that, Cov(X 1 ) = ? X1 = ? ? ? ? ? ? ) ( ) , ( ) ( 21 21 11 11 X V X X Cov X V (ii) Where ) ( 11 X V and ) ( 21 X V are the variance of univariate modified EWMA statistic,and) , ( 21 11 X X Cov = ? 1 ? 2 ? 12 + (1-? 1 )(1-? 2 ) ? 12 + cov( 10 11 Y Y ? , 20 21 Y Y ? ). Then as per (ii) Cov(X n ) = ? Xn = ? ? ? ? ? ? ) ( ) , ( ) ( 2 2 1 1 n n n n X V X X Cov X V Now Y n 's are autocorrelated normal with covariance matrix ?, so that the (Y n -Y n-1 ) 's (n?1 ) have covariance matrix 2(I-Rho)? with Rho = diag(? y1 , ? y2 , ?, ? yp ). Then, taking ? 1 =? 2 =?= ? p =? and ? y1 , ? y2 , ?, ? yp ?1, as n tends to infinity. Lemma 3: The variance of univariate Modified EWMA (MOEWMA) control statistic X n is, V (X n ) = V(? ? ? = 1 0 n j (1-?) j Y n-j )+ V((1-?) n Y 0 )+ V( ? ? = 1 0 n j (1-?) j (Y n-j -Y n-j-1 )) © 2013 Global Journals Inc. (US) # Global Journal of Computer Science and Technology Volume XIII Issue III Version I 29 ( D D D D D D D D ) Year V (X n ) = (1-?) 2n V(Y 0 ) + ) , () 1 ( 2 ) ( ) 1 ( 11 0 1 0 1 2 2 2 2 ? ? ? ? = ? = + ? ? ? ? + ? j n j n n j n j j j n j Y Y Cov Y V ? ? ? ? + )] ( ), [( ) 1 ( 2 ) ( ) 1 ( 2 1 1 1 0 1 0 1 2 1 2 ? ? ? ? ? ? ? ? = ? = + ? ? ? ? ? ? + ? ? ? ? j n j n j n j n n j n j j j n j n j Y Y Y Y Cov Y Y V ? ? ( ) ( ) ? ? ? = ? ? ? ? ? + ? = ? ? ? ? ? ? + ? ? + 1 0 1 1 1 2 1 0 1 2 ) ( , ) 1 ( ) ( , ) 1 ( n j j n j n j n j n j j n j n j n j Y Y Y Cov Y Y Y Cov ? ? ? ? Since Y n 's are autocorrelated normal with variance ? 2 , the variance of (Y n -Y n-1 ) (n?1) is = ) ( n X V 2 2 ) 2 ( ) 1 ( 2 ) 2 ( ?? ? ? ? ? ? ? ? ? + ? + ) 2 ( ) 1 ( 2 2 ? ? ? ? ? ? + ) 2 ( ) 1 )( 1 ( 4 2 3 ? ? ? ? ? ? ? ? ? + ) 2 ( 1 2 2 ) 1 ( ) 2 ( ) 1 ( 2 2 2 2 2 1 ? ? ? ? ? ? ? ? ? ? ? ? + ? ? (iii) In normal autocorrelated process (a) with 3 ? nearly negative half and 1 ? , 2 ? nearly equal and opposite in sign and being monitored for small shifts, ( with autocorrelation ? nearly one (??1) the above expression (iii), reduces to = ) ( n X V 2 2 ) 2 ( ) 1 ( 2 ) 2 ( ?? ? ? ? ? ? ? ? ? + ? (iv) Let 2 ) 2 ( ) 1 ( 2 ?? ? ? ? ? ? is a small value for high value of ? and small ? , sometimes even negligibly small such that modified EWMA limits equal EWMA limits. Therefore, V(X n ) = ? ? ? ) 2 ( ? ? ? + ) 2 ( ) 1 ( 2 ? ? ? ? ? ? ? ? ? 2 . (v) Therefore, Multivariate MOEWMA covariance from equation (iv, v) becomes ) ( 1n X V = ? ? ? ) 2 ( ? ? ? + ) 2 ( ) 1 ( 2 ? ? ? ? ? ? ? ? ? = ) ( 2n X V = ) , ( 2 1 n n X X Cov . In general the best approximation of covariance matrix of the MMOEWMA p-variable vectors is given by, ? = Xn ? ? ? + ? ? ? ? ? ? ? ? ) 2 ( ) 1 ( 2 ) 2 ( (vi) The UCL of MMOEWMA control chart is UCL = ( ) 2 / 1 2 2 / 1 ) 2 ( ) 1 ( 2 ) 2 ( h ? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? (vii) The MMOEWMA chart gives an out-of-control signal as soon as 2 2 n T = X n ' ? xn -1 X n > h 2, (viii) Where h 2 (>0) is chosen to achieve a specified in control ARL. If one or more points fall beyond h 2 , the process is assumed to be out-of-control. The magnitude of the shift is reflected in the non-centrality parameter ? 1 ' ? -1 ? 1. We conclude that an assignable causes result in a shift in the process mean from ? 0 to ? 1. Appendix 2 RL Computation for MMOEWMA Scheme using Markov Chain Approach Following Runger and Prabhu (1996) and Molnau et al. (2001) the Markov chain approach of ARL for MMOEWMA has been derived. Different choices of ? (weighting factor), h 2 (decision value), and p (number of variable) are considered. # In and out of-Control Case For the in or out control case, the ARL analysis can be simplified as a one dimensional Markov chain. # Global Journal of Computer Science and Technology Volume XIII Issue III Version I In this case the two dimensional range of X n is represented by the X 1 and X 2 axes, and the states used for the Markov chain are assumed as circular rings. Because Y n has a spherical distribution, the probability of transitioning from state i to state j, denoted as p(i, j), depends only on the radii of states i and j. For i = 0,1,2,?,m and j not equal to zero. p(i,j) = P (d n in state j | d n-1 in state i) = P [ (j-0.5) g < ) ( ) 1 ( 1 1 ? ? ? + ? + n n n n Y Y X Y ? ? < (j+0.5) g | d n-1 = ig]. Given that d n-1 = ig, X n-1 is distributed as igU, and j =0,1,2,?,m p(i,j) = P [ (j-0.5) g < ) ( ) 1 ( 1 ? ? + ? + n n n Y Y igU Y ? ? < (j+0.5) g | d n-1 = ig]. Let e denote the p component unit vector e' =(1,0,0,?,0). According to Runger and Prabhu (1996) Y n and U are independent spherical random variables, without loss of generality it can assume that U is identity equal to e to obtain p(i,j) = P [ {(j-0.5) g }/ ?< ? ? For any control chart that is approximated by a Markov chain, the run length performance can be determined from the transition probability matrix. Assume that a Markov chain has s states (see Brook and Evans (1972)). The transition probability matrix contains the transition probabilities for moving from state to state. Let this s x s matrix of transition probabilities be presented as P, where the process mean vector is such that the non centrality parameter is ?. Let the s x 1 vector q designate the starting state of the Markov chain. The vector q will have a one in the component corresponding to the starting state and zeros in all of the other components. The zero state ARL of a scheme modeled as a Markov chain represented by ARL=q'(I-P) -1 1. (ix) / )] ( ) 1 [( 1 ? ? + ? + n n n Y Y ige Y < {(j+0.5) g} / ?]. Let ) , ( # Steps of ARL Computation for MMOEWMA Step-1 Choose the parameter ? (Weighting factor), h 2 (decision value), p (number of variable), and shift in mean vector d. Step-2 The upper control limit of MMOEWMA chart is, UCL = Step-6 Non centrality parameter, c = [ {(1-?) i g / ?}+d] 2 , degree of freedom is p , d is the shift in mean vector. Step-7 For j not equal to zero (j ? 0), p(i,j) = P Step-8 For the case where j =0, we have p(i,0) = P [ ? 2 (p,c) < {(0.5) 2 g 2 / ? 2 }]. Step-9 Adjust the t. p.m.(R a ) such that row sums are unity. Step-10 Compute u = [I -R] -1 1 Step-11 Compute q = R a '* I ![signals that a statistically significant shift in the mean has occurred; that is process out-of-control, where h > M](image-2.png "") 1![Figure 1 and 2 depict MEWMA and MMOEWMA statistics charting.](image-3.png "CFigure 1") 1![Figure 1 : MEWWA Control Chart (p=3) Observe the shoot up bar showing abrupt shift in figure 2 which is completely missing in figure 1.](image-4.png "Figure 1 :") 1NOObservationsDeviated ObservationsY 1Y 2Y 3Y 1 -? 01Y 2 -? 02Y 3 -? 03092.2495.56100.270.000.000.00192.8395.16100.770.59-0.400.50???????1292.8495.16100.760.60-0.400.491392.8495.16100.760.60-0.400.491492.8495.16100.760.60-0.400.491592.8495.16100.760.60-0.400.49???????a) MMOEWMA chart for monitoring Multivariate Chemical Process using table 1 and table 2 2YearD D D D D D D D )(shows that, MMOEWMA vector givesbest forecasts for the process mean vector;undoubtedly better than the MEWMA prediction barringabruptchangesituation(1201 stobservation).MMOEWMA also detects all the shifts more timely ascompared to MEWMA for chemical process data. 2temperature DataNo.MEWMA vector, ?=0.1,MEWMAMMOEWMA vector,MMOEWMAh 1 = 14.78statistics?=0.1, h 2 = 5.417statisticsX 1X 2X 3T n1 2X 1X 2X 3T n2 210.06-0.040.050.240.560.460.552.91?????????120.43-0.290.3612.480.58-0.140.515.19130.44-0.300.3713.480.58-0.160.505.47*140.47-0.310.3814.430.58-0.190.505.73*150.47-0.320.3915.30*0.58-0.210.505.98*?????????0.39-0.460.2115.33*0.36-0.480.195.46*0.38-0.460.2115.17*0.36-0.480.195.400.38-0.460.2115.02*0.36-0.480.195.350.37-0.460.2114.87*0.35-0.480.195.30?????????-0.84-0.29-0.5214.26-0.86-0.32-0.555.37-0.84-0.29-0.5214.41-0.86-0.32-0.555.42*-0.84-0.29-0.5314.56-0.87-0.31-0.555.47*-0.85-0.29-0.5314.72-0.87-0.31-0.555.53*-0.85-0.29-0.5314.87*-0.87-0.31-0.565.59*?????????0.020.67-0.1314.91*0.050.70-0.115.44*0.030.67-0.1314.750.060.70-0.105.38?????????12000.430.290.434.720.460.320.461.9112010.430.290.556.491.681.541.8129.10*12020.430.290.546.330.340.190.451.54?????????14160.31-0.180.6414.730.31-0.190.635.2414170.31-0.190.6314.79*0.30-0.190.635.25?????????14340.30-0.230.6115.35*0.28-0.250.605.4114350.30-0.230.6115.38*0.28-0.250.605.43*?????????14390.30-0.240.6115.48*0.28-0.260.595.46* 3h 1 =10.7512.3413.1014.7815.1616.94ShiftP =2 and ? =0.10p =3 and ? = 0.1p =4 and ? = 0.10.050199950210074979950.539.551.145.661.252.3681.012.113.713.515.314.516.51.57.037.697.668.408.208.992.04.975.395.435.835.796.232.53.904.234.264.544.514.813.03.273.493.523.753.743.97Table 4 : Average Run Lengths of MMOEWMA ChartsARL Values for Multivariate Modified EWMA (MMOEWMA) Charts (p=2, p=3 and p=4)h 2 =3.934.5254.785.4175.5466.221Shiftp =2 and ? =0.1p =3 and ? = 0.1p =4 and ? = 0.10.05001000500100050010000.530.540.0931.5241.5732.342.771.010.0611.6310.5412.1510.8912.551.55.696.386.046.736.307.02.03.874.294.164.584.384.802.52.873.173.123.423.313.613.02.022.242.252.462.622.63 50132Year28Volume XIII Issue III Version ID D D D ) C(Global Journal of Computer Science and Technology? = Shift 0 0.5 1.0 1.5 2.0 2.5 3.0 ? = h 2 =2 ? chart h 10.6 200. 116. 42. 15.8 6.9 3.5 2.2 Table 6 : Average Run Lengths of MMOEWMA Charts MEWMA Chart 0.1 0.2 0.4 0.6 h 1 8.66 9.65 10.29 10.53 ARL values for p=2 from Lowery et. al. (1992) 200 201 199 200 28.1 35.10 51.9 73.6 10.2 10.10 13.2 19.3 6.12 5.50 5.74 7.24 4.41 3.80 3.54 3.86 3.51 2.91 2.55 2.53 2.92 2.42 2.04 1.88 MMOEWMA Chart 0.1 0.2 0.3 0.4 0.6 3.135 3.742 4.223 4.705 5.850.8 10.58 200 95.5 28.1 10.3 4.75 2.75 1.91 0.8 7.561ShiftARL values for p=202002002002002002000.521.124.7629.929.951.674.981.08.127.888.658.6514.9323.991.54.784.114.064.065.758.922.03.292.642.442.442.924.122.52.451.871.691.691.862.363.01.901.431.321.321.411.63© 2013 Global Journals Inc. (US) [i,1]<-r%*%Xn[i,1]+(1-r)*Z1+(Xn[i,j]-X0 ) * Monitoring Sewage Treatment Plants: Some quality control aspects PMBerthouex WGHunter LPallesen Journal of Quality Technology 10 1978 * An Approach to the Probability Distribution of CUSUM Run Lengths DBrook DAEvans Biometrika 59 1972 * Statistical Process Control Procedures for Correlated Observations TJHarris WHRoss Canadian Journal of Chemical Engineering 69 1991 * Multivariate Quality Control, Techniques of Statistical Analysis HHotelling Eisenhart, Hastay, and Wallis 1947 McGraw-Hill New York * A multivariate exponentially weighted moving average control chart CALowry WHWoodall CWChamp SERigdon 1992 * Technoetrics 34 * A Program for ARL Calculation for Multivariate EWMA Charts WEMolnau GCRunger DCMontgomery KRSkinner ENLoendo Journal of Quality Technology 33 4 2001 * Some Statistical Process Control methods for Autocorrelated data (with discussion) DCMontgomery CMMastrangelo Journal of Quality Technology 26 1991 * Modified exponentially weighted moving average (EWMA) control chart for an analytical process data AKPatel JDivecha Journal of Chemical Engineering and Material Science 2 1 2011 * A Markov Chain Model for the Multivariate Exponentially Weighted Moving Averages Control Chart GCRunger SSPrabhu Journal of American Statistical Association 91 436 1996 Theory and Methods * Multivariate statistical process control for autocorrelated process GCRunger International Journal of production research 34 1996. 1996 * The R Foundation for Statistical Computing Version 2.0 RLanguage 2004. 2004-11-15 * R-Program for monitoring Modified Multivariate Exponentially Weighted Moving Average Control Chart