The quality control chart for monitoring multivariate autocorrelatedprocesses

Computational Statistics &Data Analysis 51(2007)3862–3870/locate/csdaThe quality control chart for monitoring multivariateautocorrelated processesJeffrey E.Jarrett a ,∗,Xia Pan ba University of Rhode Island,USAb University of Illinois at Springfield,USAReceived 13July 2005;received in revised form 15January 2006;accepted 17January 2006Available online 30March 2006AbstractPreviously,quality control and improvement researchers discussed multivariate control charts for independent processes and univariate control charts for autocorrelated processes separately.We combine the two topics and propose vector autoregressive (V AR)control charts for multivariate autocorrelated processes.In addition,we estimate AR(p )models instead of ARMA models for the systematic cause of variation.We discuss the procedures to construct the V AR chart.We examine the effects of parameter shifts and by example present procedures to show the feasibility of V AR control charts.We simulate the average run length to assess the performance of the chart.©2006Elsevier B.V .All rights reserved.Keywords:Quality control and improvement;V AR control charts;Multivariate autocorrelated processes1.Introduction and purposeShewhart (1931)developed statistical process control (SPC)techniques for the improvement and assurance of product and service quality.The assumptions of these methods include serial independence of successive samples taken from a process.Alwan and Roberts (1988)proposed modeling and filtering to monitor residuals by a one-step-forecast of autocorrelated process possessing systematic patterns.English and Sastri (1990)applied Kalman filter to monitoring univariate autoregressive (AR)(1)process.Wardell et al.(1994)derived the run length distribution of Shewhart chart on residuals.Montgomery and Mastrangelo (1991),Lu and Reynolds (1999)discussed some SPC charts for autocorrelated data.Runger and Willemain (1995)compared Shewhart charts based on several process transformations to eliminate autocorrelation.Snoussi et al.(2005)indicated that residuals control charts provide much better shift detection properties.Testik (2005)considered using an AR (1)process to modify control limits for residual based charts.Loredo et al.(2002)provided a method for monitoring autocorrelated processes based on Hawkins (1991,1993)regression adjustment.They examined the performance of residual based quality control charts in terms of the average run length to observation based control charts via Monte Carlo simulations.Their results indicate the superiority of residual based control in comparison with the observation based control charts when detecting a mean shift in short run∗Corresponding author.Tel.:+14018744169;fax:+14018744312.E-mail addresses:jejarrett@ (J.E.Jarrett),xpan1@ (X.Pan).0167-9473/$-see front matter ©2006Elsevier B.V .All rights reserved.doi:10.1016/j.csda.2006.01.020J.E.Jarrett,X.Pan /Computational Statistics &Data Analysis 51(2007)3862–38703863autocorrelated data.This is important because the i.i.d.assumptions of standard Shewhart control charts do not hold for large amounts of historical data.Kalgonda and Kulkarni (2003,2004)considered using the Hotelling T 2statistic to monitor several quality char-acteristics.Due to difficulties of this statistic in determining the factors responsible for out of control signals when considering changes in the process mean vector and/or process covariance matrix,they proposed a new diagnostic,“the D-technique”.Their analysis including a dummy variable multiple linear regression provided some evidence of usefulness of constructing multivariate control charts.Their second study (Kalgonda and Kulkarni,2004)extended the usefulness of multivariate control charts in detecting an out of control status and identifying variable(s)responsible for the out of control situation.Another study,Jaerkaporn et al.(2003)corroborated the usefulness of using average run length methods to the standard Shewhart control charts for individuals.Their method showed promise that average run length control charts are superior for detecting shifts in the mean than the traditional charts.Sliwa and Schmid (2005)extended the use of average run length to multivariate time series.Their study focused on the autocovariance and the cross-covariance structure of financial assets.Hence,monitoring the process control for new multivariate time series resulted from their approach.Finally,Cheng and Thaga (2005)produced a CUSUM control chart capable of detecting changes in the mean and standard deviation of autocorrelated data.Based on the average run length,they claim that their new control chart (MCAP)is useful for monitoring modern production processes where one produces high quality goods with a tiny fraction of nonconforming output units.We propose to continue the study of multivariate autocorrelated processes by proposing a vector autoregressive (V AR)control chart based on ARL for precision manufactured products that are now produced en masse in high technology oriented industries.One example of such production is optical communication products manufacturing where the output is a multivariate autocorrelated time series.One bases the production of fiber optic on SiO 2rod made from condensation of silicon and oxygen gases.The preparation of the SiO 2rod needs to control temperature,pressure and the concentrations of different components.We find similar processes in semiconductors where they prepare materials.In these processes,correlated variables are plentiful and the addition of multivariate analysis is necessary to monitor the process.In addition,these variables are cross-correlated with time leads and lags.Univariate quality control charts are effective for monitoring individually correlated variables but are not effective for services,which are cross-correlated.We propose to introduce (V AR)control charts designed to monitor multivariate processes of the type described.Shewhart control charts do not emphasize parameter estimation.This lack of concern emanates from using the rough criterion of average run length (ARL)to determine effectiveness.For autocorrelated processes,estimation is the key procedure for control chart construction.Some of the proposed control charts for autocorrelation are also not applicable if we ignore estimation.Hence,we will consider estimation issues associated with ARIMA modeling.We arrange the remaining parts of this paper as follows:In Section 2,we describe the principles and design of V AR control charts.In Section 3,we examine the effects of parameter shifts.In Section 4,we show an example V AR chart construction for a multiple variable AR (3)process.In turn,we examine the st,in Section 5,we summarize.2.V AR control chartWe propose using a V AR control chart (another title is MultivariateAutoregression (MAR)control chart).The process is required either stationary or we can achieve stationarity with the appropriate filtering process,i.e.differencing.For a multivariate autoregressive process of n variables,denote y t =(y 1t ,y 2t ,...,y nt ) as a (n ×1)vector.These cross-correlated variables contain autocorrelation of order p .The process follows a V AR modely t = c + 1L + 2L 2+···+ p L py t + e t = c + (L) y t + e t ,(1)where L is the backshift operator, c =(c 1,c 2,...,c n ) is the constant vector,and e t =(e 1t ,e 2t ,...,e nt ) is the vectorfor the error term.Each j is a (n ×n)coefficient matrix for lag j . e t contain no autocorrelation but does contain cross-sectional correlation.That is,E e t e t = is invariant about time but may not be a diagonal (n ×n)matrix.Some elements of the matrix j are zero when variables are uncorrelated.In this model,the AR coefficient matrix, = 1 2··· p,and the constant term ¯c represent systematic patterns.The error term vector e t represents nonsystematic disturbances.We monitor the process to filter the systematic element.In turn,we monitor the special-cause.To filter out and estimate the systematic element for a process that is in-control3864J.E.Jarrett,X.Pan /Computational Statistics &Data Analysis 51(2007)3862–3870requires multivariate regression,which,in turn,requires huge sample sizes for asymptotic estimation.To construct control charts by these methods,we implement a two-stage routine.For in-control cases,we estimate the systematic elements.We denote =( c )being an n ×(1+np)matrix,and x t = 1, y t −1, y t −2,..., y t −pbeing a 1×(1+np)vector.Then Eq.(1)becomesy t = x t + e t .(1)Due to the usual assumptions about e t ,we estimate the coefficient matrix n by ordinary least squares (OLS).ForGaussian AR(p )process, e t ∼i.i.d N 0, ,these OLS estimates are ˆ=T t =1y t xt Tt =1x t xt −1,(2)where T is the number of observations in Stage I.The estimated systematic model isˆ y t=ˆ x t .(3)In addition,we estimate error term as ˆ e t = y t −ˆ y t.We then develop the special-cause control chart based on e t .Maximum likelihood methods estimate the covariance of e t asˆ=(1/T )T t =1ˆ e t ˆ e tp→ .(4)If the process is in-control,the model in (1)is appropriate and the error term is uncorrelated and normal,the estimatederror e t should be also approximately normal and invariant with time.Thus,it is asymptotically N 0, .We now apply the multivariate control chart e t .In Stage I the Hotelling T 2statistic is,approximately 2(n)when T −n is large.T 2=ˆ e t −¯ˆ e ˆ −1 ˆ e t −¯ˆ e ∼ 2(n),(5)where ¯ˆ e is mean of Stage I,¯ˆ e = T t =1ˆ e t /T ,which should be close to zero.Since the estimates of the AR(p )model areasymptotically available,we use chi-square instead of Beta distribution.For Stage II observations,Y i ,(i >T ),we denote the first Stage II observation as the (T +1)st observation.We substitute the estimated error by the one-step-ahead forecast error,while we compare the first Stage II observation with the one-step-ahead forecast from Stage I.That is,from the known T observations in Stage I,we predict the first Stage II status asˆ y T +1=ˆ xT +1.(6)Moreover,the forecasting error is the difference between the real observation and the forecast,ˆ e T +1=ˆ y T +1− yT +1.(7)There are several possible schemes for Stage II chart.If we treat ˆas Wishart distribution,the chart statistic is T 2T +1=ˆ e T +1−¯ˆ e ˆ −1 ˆ e T +1−¯ˆ e ∼n(T +1)(T −1)T (T −n)F ,n,T −n .(8)However,in an autoregression environment,we only have asymptotic properties,ˆ→ .Therefore,the more appropriate way is to use the approximation:T 2T +1= ˆ e T +1−¯ˆ e ˆ −1 ˆ e T +1−¯ˆe∼ 2(n).(9)J.E.Jarrett,X.Pan /Computational Statistics &Data Analysis 51(2007)3862–38703865For the following observations,the forecasting errors are acquired recursively one-step-ahead each.Generally,forthe forecasting error of the i th observation,ˆe i =ˆ x i − y i ,the chart statistic is T 2i =ˆ e i −¯ˆ e ˆ −1 ˆ e i −¯ˆ e∼ 2(n).(10)An F distribution is asymptotically close to a chi-square distribution (divided by T ).The asymptotic distributions ofthe estimates provide,under the same confidence level,more narrow confident intervals for true parameters than the exact distribution for limited sample size T.Therefore,control charts using chi-square distributions are more sensitive than those using F distributions.Selecting statistic :One statistic is given in (10).By selecting the expression of the mean item of T 2i we can have morechart schemes.One option is directly choosing ¯ˆ e as its asymptotic value,zero.Then the item ¯ˆe in (10)will disappear:T 2i =ˆ e iˆ −1ˆ e i ∼ 2(n).(11)The corresponding Stage I statistic isT 2=ˆ e tˆ −1ˆ e t ∼ 2(n).(5)Another option is using recursive mean.This is corresponding to the chart of moving centerline in univariate case(Montgomery and Mastrangelo,1991).The mean for the (i +1)th observation’s statistic can be defined as ¯ˆ e i =ˆ e i +(i −1)ˆe i −1 /i so that the statistic is T 2i =ˆ e i −¯ˆ e i ˆ −1 ˆ e i −¯ˆ e i∼ 2(n).(12)Although the recursive mean varies every time when a new observation is added,the control limit of the chart doesnot change.This is much better than the moving center line in the univariate case.Note,all the above schemes have the same chart that we determine the control limit (UCL)by 2(n)say, 20.01(n).Only the statistic values are different for different schemes.Hence,the V AR chart is actually simple in its demonstration,and model estimation is relativelysimple.Among (10)–(12),the statistic in (11)is the most convenient.Because the asymptotic mean of ˆ¯e t is zero thus(11)is theoretically justified for large sample size.Determining the AR order :In practice,determining the autocorrelation order p may be the issue before chart setting.To find a parsimonious AR order,the null hypothesis is that the AR order is p 0as opposed to being p 1>p 0.The test statistic is a likelihood ratio:2(L 1−L 0)=T log ˆ 0 −log ˆ 1 ∼ 2 n 2(p 1−p 0) ,(13)where T is the sample size in Stage I. 0is the estimated covariance matrix of the error term given in (10)under the null hypothesis being true. 1is also from (10)but under the alternative hypothesis.The chi-square distribution has a degree of freedom n 2(p 1−p 0)because each variable has (p 1−p 0)less lags under the null hypothesis and each of the n equations include n variables and totally n 2restrictions.3.The effect of shiftThere are three classes of parameters in V AR models:the process mean ,the shift on covariance matrix of error term, ,and the shift autocorrelation structure of the model, .The most important shift of interest is the shift on the mean of the process.For a mean shift,the relationship between the mean of a stationary time series vector, is the expectation of Eq.(1).The following theorem gives the effects of mean shift on the V AR chart statistic.There are possibly two ways for a shift to build on the ARMA process.One is that the shift is added on the constant item,meaning there is really a physical3866J.E.Jarrett,X.Pan /Computational Statistics &Data Analysis 51(2007)3862–3870structure that the historical statuses are added on the constant.Another way is that the shift is added on the mean of the process,which is more an observed result.We here only consider the first situation.Theorem 1.For a process model in (1)and VAR chart in (10),a shiftin the process constant adds a normal distribution with mean of ˆ kˆ −1ˆ k and variance of 4ˆ k ˆ −1ˆ k on the in-control T 2.The effects depend on the shift magnitude and the estimates of VAR parameter matrices ,ˆ =(I n −ˆ (L)) .Proof.In the following proof,we use the estimated coefficients rather than the known coefficients.Denote the in-control residual here as ˆ e t ,and the out-of-control (due to shift)residuals as ˆ t .In addition,we use z t to denote theprocess after shift.These changes in denotation are not necessary but just for an alternative expression.The results for the shift effects are the same except that one is univariate ARMA and the other is multivariate AR,because the asymptotic estimates should be just the known parameters.The forecasting error without shift occurrence is ˆ e t = y t −ˆ y t.After the shift,the process status is z t rather than y t .A sudden shift at time just before t 0on mean for some of the variables impliesz t = yt + for t t 0.Since the shift is unknown before detection,we still base the one-step-ahead forecasting on the estimated coefficient matrix ˆ= ˆ c ˆ .Thus,for a V AR(p )process,ˆz t =ˆ c +ˆ (L) z t when t t 0+p .Then the forecasting error after the shift isˆ εt t= z t −ˆ z t = z t −ˆ c −ˆ (L) z t = y t + −ˆ c −ˆ (L) y t −ˆ (L) = y t −ˆ (L) −ˆ y t+ =ˆ e + I n −ˆ (L)=ˆ e t +ˆ.Here the ˆe t was the forecasting error at time t ,had the shift not occurred.With the shifted ˆ εt =ˆ e t −ˆ ,the Hotelling T 2statistic will be a noncentral chi-square distribution:shift T 2t=ˆ ε tˆ −1ˆ εt = ˆ e t −ˆ ˆ −1 ˆ e t −ˆ =ˆ e tˆ −1ˆ e t −2ˆ ˆ −1ˆ e t +ˆ ˆ −1ˆ =noshift T 2t −2ˆ ˆ −1ˆ e t +ˆˆ −1ˆ .The first term is the in-control Hotelling statistic noshift T 2t ∼ 2(n)with mean of n and variance of 2n .Since asymptotically ˆ e t ∼N 0, ,the second term −2ˆ ˆ −1ˆ e t is a linear combination of normal distributions so that it is normal with mean of 0and variance of 4ˆˆ −1ˆ ˆ −1ˆ =4ˆ ˆ −1ˆ .Its effect,we view that the first term increases variation but without shifting the mean.The third term is a positive definite scalar.The last two terms are normaldistribution with mean of ˆˆ −1ˆ and variance of 4ˆ ˆ −1ˆ .Therefore,we measure the total effect of the process shift ˆ on T 2by ˆˆ −1ˆ .When t 0 t <t 0+p ,let k =t −t 0.We have z t = y t +for t 0 t <k and k lags shift while p −k lags do not shift.The one-step-ahead prediction isˆz t =ˆ c +ˆ 1 z t −1+ˆ 2 z t −2+···+ˆ k z t −k +ˆ k +1 y t −k −1+···+ˆ p y t −p =ˆ c +ˆ 1L +ˆ 2L 2+···+ˆ k L k z t +ˆ k +1 y t −k −1+···+ˆ p y t −p .J.E.Jarrett,X.Pan /Computational Statistics &Data Analysis 51(2007)3862–38703867Thusˆ εt t= z t −ˆ z t = z t − ˆ c +ˆ 1L +ˆ 2L 2+···+ˆ k L k z t +ˆ k +1 y t −k −1+···+ˆ p y t −p= y t + − ˆ c +ˆ 1L +ˆ 2L 2+···+ˆ k L k yt +ˆ k +1 y t −k −1+···+ˆ p y t −p − ˆ 1L +ˆ 2L 2+···+ˆ k L k=ˆ e t +ˆ k,where ˆ k =I n −ˆ1−ˆ 2−···−ˆ k and k =0,1,2,...,(t −t 0).Therefore the effect on V AR chart statistic is,for t 0 t t 0+p ,shift T 2t =noshift T 2t−2ˆ k ˆ −1ˆ e t +ˆ kˆ −1ˆ k ,where k =t −t 0.This is still a positive shift on the mean of the in-control T 2and on the variance of the in-control T 2.From the definition of ˆ and ˆ k ,obviously ˆ =ˆ pfor t t 0+p ,because k =0for k >p ,and the estimates of k should have the same property asymptotically.The effects during t 0 t t 0+p can be either larger or smaller than theeffects during t t 0+p ,depending on the specific V AR structure ˆ(L)and ˆ . Other systematic shifts,such as the shift on the covariance matrix of error term, ,or the shift on the autocorrelationstructure of the model, ,will also violate the distributions of the chart statistic,hence out of control signals will be observed after the shifts.The larger the shift scale,the more serious the violation of the chi-square distribution.The out-of-control ARL,which determines when the V AR chart detects the shift,relates to the control limit and the scale of the shift.4.A simulated example of a V AR chartFollowing is an illustration showing the procedures and properties of a V AR chart.We specify three as both the order of autocorrelation and the number of variables.We design the covariance matrix of the error term vector as= 10.5−1.50.54.25−1.75−1.5−1.752.99.In this way,we keep the positive definite.The constant term is set to be (3 =(1.12.11.1),and the three coefficient matrices are as follows:1= 0.50.2−0.08−0.10.70.20.25−0.220.54 , 2= −0.250.23−0.1−0.340.320.120.43−0.180.15 , 3= 0.21−0.1−0.13−0.09−0.32−0.210.150.050.48.By this design,all of the characteristic values ofI n p − 1p −1− 2 p −2−··· p =0,lie in the unit circle,| i |<1,i =1,2,...,p .This setting guarantees a stationary model (p =3,n =3).In addition,the variable means are = cI n − 1L − 2L 2− 3L 3 −1=(0.498211.649915.082).With the above parameters,we sampled and experimented using Matlab.Fig.1(a)and (b)are the results for the two statistics in (5)and (5 )at Stage I,where the UCL was preset at =0.005.From the figures,all the 100observations fell within the UCL.Actually,if we repeat the same simulation,sometimes we found one or two observations fallout of UCL.This may imply the in-control ARL may be suitable.In fact,in Fig.1(c)and (d)if we try 500observations,we observed from one to five out UCL samples.Fig.1also showed that the statistics in (11)and (11 )at Stage II almost have no difference.To confirm the autocorrelation order p ,we run the likelihood ratio test in (13)for the hypotheses (p 0=2vs .p 1=3),(p 0=3vs .p 1=4),(p 0=3vs .p 1=6),and (p 0=4vs .p 1=6).At =0.05the null hypothesis was rejected for first test and cannot be rejected for the other three tests.This showed that the recommended method for selecting order p3868J.E.Jarrett,X.Pan /Computational Statistics &Data Analysis 51(2007)3862–3870STAGE I VAR CHART WITH MEAN ERRORS, alpha=0.005Observations H o t e l l i n g T -s q u a r e S t a t i s t i cSTAGE I VAR CHART WITHZERO MEAN, alpha=0.0027ObservationsH o t e l l i n g T -s q u a r e S t a t i s t i c 02468101214STAGE I VAR CHART WITH ZERO MEAN, alpha=0.005ObservationsH o t e l l i n g T -s q u a r e S t a t i s t i c0102030405060708090100024********STAGE I VAR CHART WITH MEANERRORS, alpha=0.005ObservationsH o t e l l i n g T -s q u a r e S t a t i s t i c(a)(b)(c)(d)Fig.1.(a)Stage |V AR chart with mean errors,alpha =0.005.(b)Stage |V AR chart with zero mean,alpha =0.005.(c)Stage |V AR chart with mean errors,alpha =0.005.(d)Stage |V AR chart with zero mean,alpha =0.0027.is valid for this model.It should be noted that not rejecting hypothesis (p 0=4vs .p 1=6)does not mean the selection is finished.Continuing searches occur until we reject the null hypothesis for (p 0=2vs .p 1=3).We examine the in-control ARL via experimental simulation.Observe the histogram of the run length distribution based on certain estimated in-control process parameters in Fig.2(a).The ARL is 215for =0.0027,and the respective standard deviation is 227.Hence,the UCL should be set higher than anticipated.Fig.2(b)shows the run length distribution at higher UCL, =0.0020.This ARL is 282.6and the standard deviation is 305.3.When a shift occurs on the variable means of the process,say,the constant term of the first and the third variable,c 1and c 3,are doubled,then c =(2.22.12.2),and the mean vector becomes 1 =[−0.619818.195027.5453].By comparison,this is a very large shift in mean and covariance matrix (about a one error-term sigma shift on Variable 1,half error-term-sigma shift on Variable 2,and four sigma shift in Variable 3).The result of the shift is the ARL =1.24and its standard deviation is 0.56.( =0.002).The V AR chart is thus very sensitive to shifts as shown in Fig.3at observation 80.J.E.Jarrett,X.Pan /Computational Statistics &Data Analysis 51(2007)3862–38703869100300500700900110013001500010********60Histogram of Run Length DistributionRun LengthF r e q u e n c y5015025035045055065075085095010501150125013501450050100150200250300350400Histogram of Run Length DistributionRun LengthF r e q u e n c y(a)(b)Fig.2.(a)Histogram of run length distribution.(b)Histogram of run length distribution at higher UCL.102030405060708090100010********60Run Chart for ShiftObservationsT -s q u a r e S t a t i s t i cFig.3.Run chart for shift.If a small shift occurs,for example,the constant term of the third variable increases by 50%,(corresponding to about 0.4error-term-sigma shift on Variable 2and 0.9error-term-sigma shift on Variable 3),the ARL is 104and standard deviation of ARL is 117(a =0.0020).Similarly,if the constant of the third variable decreases 50%,the effect on ARL is 113.Similar to the increase case,although the shift on the means are different (half error-term-sigma shift on Variable 1and one error-term-sigma shift on Variable 3).3870J.E.Jarrett,X.Pan/Computational Statistics&Data Analysis51(2007)3862–3870Since the simulations were based on certain estimates of the parameters instead of the true underlying parameters, the ARL values are different from the results using true parameters.For example,when true parameters were used,the mean shift effect in above give ARL=84.5with8.67standard error.ARL based on estimates may be larger or smaller than the ARL from true parameters,depending on whether overvaluing or undervaluing the estimates.As the system is dynamic,a given set of data can produce overestimation on one parameter and underestimation on the other.Since the true process is unknown,there is actually no way to assess if one overvalues or undervalues the ARL.The only prudent way is to collect as many data as possible in Phase I to acquire asymptotic estimates.5.SummaryVar control charts are a combination of univariate residual chart for autocorrelation and multivariate chart for independent processes.Based on our development and carefully explained methodology,one can model a multivariate process by AR equations describing its physical properties and estimate the AR(p)model by OLS procedures,and applies V AR control chart to monitoring autocorrelated multivariate processes.Hence,the V AR control chart is both useful as long as sample size are sufficient for OLS estimation.ReferencesAlwan,B.M.,Roberts,H.V.,1988.Time-series modeling for detecting level shifts of autocorrelated processes.J.Bus.Econ.Stat.6(1),87–96. Cheng,S.W.,Thaga,K.,2005.Max-Cusum chart for autocorrelated processes.Stat.Sin.15,527–546.English,J.R.,Sastri,T.,1990.Enhanced quality control in continuousflow put.Ind.Eng.19,258–262.Hawkins,D.M.,1991.Multivariate quality control based on regression-adjusted variables.Technometrics33(1),61–75.Hawkins,D.M.,1993.Regression adjustment for variables in multivariate quality control.J.Qual.Technol.25(3),170–182.Jaerkaporn,D.,et al.,2003.Process monitoring for correlated gamma-distributed data using generalized-linear-model-based control charts.Qual.Rel.Eng.Int.19(6),477–491.Kalgonda,A.A.,Kulkarni,S.R.,2003.Diagnosis of multivariate control chart signal based on dummy variable regression mun.Stat.: Theory and Methods32(8),1665–1684.Kalgonda,A.A.,Kulkarni,S.R.,2004.Multivariate quality control chart for autocorrelated process.J.Appl.Stat.31(3),317–327.Loredo,E.N.,Jaerkpaporn,D.,Borror,C.M.,2002.Model-based control chart for autoregressive and correlated data.Qual.Rel.Eng.Int.18(6), 489–496.Lu,C.W.,Reynolds,M.R.,1999.Control charts for monitoring the mean and variance of autocorrelated processes.J.Qual.Technol.31,259–274. Montgomery,D.C.,Mastrangelo,C.M.,1991.Some statistical process control methods for autocorrelated data.J.Qual.Technol.23(3),179–193. Sliwa,P.L.,Schmid,W.,2005.Monitoring the cross-covariances of a multivariate time series.Metrika61(1),89–115.Snoussi,A.,Ghourabi,M.E.,Limam,M.,2005.On SPC for short run autocorrelated mun.Stat.:Simulation and Computation34(1), 219–234.Testik,M.C.,2005.Model inadequacy and residuals control charts for autocorrelated processes.Qual.Rel.Eng.Int.21(2),115–130.Wardell,D.G.,Moskowitz,H.,Plante,R.D.,1994.Run-length distribution of special-cause control charts for correlated processes.Technometrics 36,1084–1105.。