Multivariate control charts based on Bayesian state space models

Multiv ariate Con trol Charts based on Ba y esian State Space Mo dels K. T rian tafyllop oulos ∗ 31 Jan uary 2006 Abstract This paper develops a new multiv ariate control c har t ing metho d for vector auto cor- related and s erially corre la ted proce sses. The main idea is to pro p o se a Bay esia n mul- tiv aria te lo cal level mo del, which is a g eneraliza tio n of the Shewha rt-Deming mo del fo r auto corr elated pro cesses, in o rder to pr ovide the pr edictive error distr ibution of the pro- cess and then to apply a univ ariate mo dified EWMA control chart to the logar ithm of the Bay es’ factors of the predictive error densit y versus th e target error density . The resulting chart is pr op osed as ca pable to deal with b oth the non-norma lity and the auto cor r ela- tion structure of the log Bay es’ factors . The new control charting scheme is general in application and it ha s the adv antage to control sim ultaneous ly not only the pro cess mean vector and the disp ersio n cov a riance matrix, but also the ent ire tar get distribution o f the pro cess. Two examples of London metal exchange data and o f pro ductio n time ser ies da ta illustrate the capabilities of the new co nt ro l chart. Some key wor ds: time series, SPC, m ultiv a riate con tro l c hart, sta te space model, EWMA. 1 In tro duc tion In the last decades m ultiv ariate Statistical Pro cess Control (SPC) h as receiv ed consider ab le atten tion, since in practice m any pro cesses are observe d in a v ector form (Montg omery 1 ). Univ ariate con trol charts ha v e b een extensiv ely discuss ed in the literature (Montg omery 1 , Bo x and Luce ˜ nno 2 , d el Castilo 3 ) and man y efforts hav e b een devo ted to upgrad in g the con- trol charts for: (a) cases of correlated un iv ariate pro cesses; and (b) cases of multiv ariate uncorrelated pro cesses. Multiv ariate control charting has b een discussed in many studies, e.g. T racy et al. 4 , Liu 5 , Kourti and MacGregor 6 , Mason et al. 7 , V argas 8 , Y e et al. 9 and Pan 10 among man y others. Review pap er s on multiv ariate control charts include Lo wry and Mon tgomery 11 , Su lliv an and W o o dall 12 , Mon tgomery and W o o dall 13 , Bersimis et al. 14 and Y eh et al. 15 . Most of the current researc h has b een fo cused on the Hotelling’s T 2 con trol chart and the multiv ari- ate EWMA con trol c hart for control ling the pro cess mean. Y eh et al. 16 , Su rtihadi et al. 17 , Cheng and Thaga 18 and Costa and Rahim 19 prop ose and study multiv ariate EWMA and CUSUM con trol c harts to con trol the disp ers ion of a multiv ariate pro cess. As stated b efore univ ariate cont rol charts for auto correlated pro cesses ha ve b een d iscussed in the literature ∗ Department of Probabili ty and Statistics, Hicks Building, Univ ersity of Sheffield, Sheffield S3 7RH, UK, Email: k.triantafyllopo ulos@sheffield.ac.uk , T el: +44 114 222 3741, F ax: + 44 114 222 3759. 1 (Mon tgomery 1 , Bo x and L u ce ˜ nno 2 ), ho we ver, f or multiv ariate pro cesses the general fo cus h as b een placed to uncorrelated pro cesses. Dy er et al. 20 , Jiang 21 , K algonda and Kulk arn i 22 and No orossana and V aghefi 23 consider multiv ariate con trol c harting for auto correlated pro cesses based on autoregressiv e-mo ving-a verag e (ARMA) time series mo dels and the T 2 and mul- tiv ariate CUS UM con trol c harts are illustrated. P an and Jarrett 24 build a multiv ariate T 2 con trol c hart for the forecast errors of the pro cess. Th ey consider a state-space appr oac h for mo delling the und erlying pr o cess and they p oint out that th e p r oblem of monitoring m ulti- v ariate pro cesses is a pr oblem of multiv ariate time series forecasting as we ll as a problem of con trol c harting. S ome forms of Ba y esian con trol charts, kn o wn also as adaptiv e or dynamic con trol c harts, are discussed in T agaras 25 , T agaras and Niko laidis 26 , de Magalh˜ aes et al. 27 and in references therein. Adaptiv e cont rol c harts offer the flexibilit y and versatilit y to dy- namically change the samp ling size and the sampling interv al of a S h ewhart con trol c hart, but they are disadv an taged in that the complexit y is increased and usu ally th e mo deller has to resort to Monte Carlo simulatio n. Our aim in this pap er is to constru ct a multiv ariate con trol c hart for auto corr elated pro cesses in suc h a wa y that the sc heme will b e capable to monitor the pro cess mean v ector only , the pr o cess disp er s ion cov ariance matrix only , or b oth the pr o cess mean vect or and the pro cess d isp ersion co v ariance matrix. W e p rop ose a new control chart based on the theory of sequentia l Ba yes’ factors (W est and Harrison 28 ). First w e fit a lo cal leve l mo del to the multiv ariate p ro cess and then we app ly a univ ariate mo dified EWMA con trol chart to the loga rithm of the Ba yes’ factor to monitor the disp ersion of the predictiv e distribution of the data f rom the target distribution. Our mo del mak es use of a generalizati on of the Shewhart-Deming m o del f or multiv ariate auto correlated pro cesses (Deming 29 , del Castilo 3 , T riant afyllop oulos et al. 30 ). Section 2 give s the necessary time series b ac kground. The pr op osed con trol chart is dis- cussed in detail in S ection 3. In Sections 4 and 5 tw o examples, consisting of data from th e London metal exc hange and from a pro d u ction of a plastic mould, illustr ate the metho dol- ogy and giv e light to the design and implementat ion of the n ew control c hart. Concluding commen ts are giv en in Section 6 and the app end ix d etails a pro of of an argument in S ection 3. 2 Bac kground The conv ent ional con tr ol c harts are based on the Shewh art-Deming mo del, e.g. for a p × 1 pro cess v ector y t this mo del sets y t = µ + ǫ t , ǫ t ∼ N p (0 , Σ) , (1) where µ is the p ro cess mean v ector and Σ is th e pro cess disp ersion co v ariance matrix, known also as the measuremen t co v ariance matrix. Here N p (0 , Σ) indicates the p -dimensional n ormal distribution with mean vect or zero and co v ariance matrix Σ. The measurement drift sequence { ǫ t } is assumed uncorrelated and this mak es the generating p ro cess { y t } an uncorrelated sequence to o. In this pap er we extend the ab ov e mo d el by consid ering equ ation (1), but no w µ is replaced b y a time-dep endent µ t , whic h follo w s a multi v ariate random wa lk m o del, kno wn also as local lev el m o del (Durbin and K o opman 31 ). Discoun t W eigh ted Regression (D WR), which originated in the path-br eaking w ork of Bro wn 32 , is a metho d for forecasti ng auto correlated time series. Considering univ ariate time 2 series Ameen and Harrison 33 dev elop ed further D WR for more complex time series. The reviews of Ameen 34 and Goo dwin 35 suggest that D WR can mo del efficien tly time s er ies in a wide range of situations. T r ian tafyllop oulos and Pik oulas 36 dev elop ed a multiv ariate v ersion of DW R and these authors fo cused on the estima tion of the measuremen t co v ariance matrix. In th is pap er we consider the D WR metho d of T rian tafyllop oulos 37 for m ultiv ariate lo cal lev el mo dels defined b y y t = µ t + ǫ t and µ t = µ t − 1 + ω t , (2) where ǫ t ∼ N p (0 , Σ) an d ω t ∼ N p (0 , Ω t Σ). The scalar Ω t is sp ecified with the aid of a discoun t factor δ and the sequences { ǫ t } and { ω t } are mutually and individu ally uncorrelated, e.g. E ( ǫ i ǫ ′ j ) = E ( ω k ω ′ ℓ ) = E ( ǫ r ω ′ s ) = 0, for all i 6 = j , k 6 = ℓ and for all r , s . Here E ( · ) denotes exp ectatio n and ǫ ′ j denotes the r o w ve ctor of ǫ j . Th e mo del definition is complete b y sp ecifying a p rior d istribution p ( µ 0 | Σ), which is usually the p -dimen s ional n orm al distribu tion, e.g. µ 0 | Σ ∼ N p ( m 0 , P 0 Σ), for some kno wn prior mean v ector m 0 and a p ositiv e scalar P 0 > 0. It is fu rther assu med that µ 0 is u ncorrelated of all ω t . F or some p ositiv e in teger N > 0, let y t = ( y 1 , y 2 , . . . , y t ) b e the information set comprising data up to and in cluding time t , for t = 1 , 2 , . . . , N . With the pr ior µ 0 | Σ ∼ N p ( m 0 , P 0 Σ), the p osterior densit y of µ t | Σ , y t is µ t | Σ , y t ∼ N p ( m t , P t Σ), where m t and P t are up dated b y m t = m t − 1 + P t − 1 δ + P t − 1 e t = δ m t − 1 + P t − 1 y t δ + P t − 1 and P t = 1 δ + P t − 1 , (3) with e t = y t − E ( y t | y t − 1 ) = y t − m t − 1 b eing the one-step forecast err or v ector at time t − 1. Define the residual error v ector r t = E ( ǫ t | y t ) = y t − m t . F or eac h time t the estimator S t of Σ is achiev ed by least squares estimatio n as S t = 1 t t X i =1 r i e ′ i = 1 t t X i =1 δ e i e ′ i δ + P i − 1 , (4) after observing that r t = y t − m t = y t − m t − 1 − P t − 1 e t δ + P t − 1 = e t − P t − 1 e t δ + P t − 1 = δ e t δ + P t − 1 . Details of the d er iv ations of m t , P t and S t app ear in T rian tafyllop oulos an d Pik oulas 36 and T riant afyllop oulos 37 . F rom the ab o ve it follo ws th at the one-step forecast density is y t +1 | Σ = S t , y t ∼ N p  m t , ( δ + P t ) S t δ  and the corresp ondin g one-step forecast error d ensit y is e t +1 | Σ = S t , y t ∼ N p  0 , ( δ + P t ) S t δ  , (5) where e t +1 = y t +1 − E ( y t +1 | y t ) = y t +1 − m t . The adequacy of the mo del is ev aluated via the mean of squ ared stand ard one-step forecast error ve ctor (MSSE), the mean of absolute p er centag e one-step forecast er r or ve ctor (MAPE) 3 and the mean of absolute one-step forecast err or v ector (MAE). Th ese statistics are discussed in Chatfield 38 and for data y 1 , y 2 , . . . , y N they are defined b y M S S E = 1 N N X t =1  ( e ∗ 1 t ) 2 ( e ∗ 2 t ) 2 · · · ( e ∗ pt ) 2  ′ , e ∗ t =  ( δ + P t − 1 ) S t − 1 δ  − 1 / 2 e t , M AP E = 1 N N X t =1  | e 1 t | y 1 t | e 2 t | y 2 t · · · | e pt | y pt  ′ , M AE = 1 N N X t =1 [ | e 1 t | | e 2 t | · · · | e pt | ] ′ , where e ∗ t is the stand ard one-step forecast er r or, y t = [ y 1 t y 2 t · · · y pt ] ′ , e t = [ e 1 t e 2 t · · · e pt ] ′ and { δ − 1 ( δ + P t − 1 ) S t − 1 } − 1 / 2 denotes the in verse of th e symmetric square r o ot of the matrix δ − 1 ( δ + P t − 1 ) S t − 1 based on the sp ectral decomp osition of symmetric matrices (Gupta and Nagar 39 ; pages 6-7). If the mo del fit is go o d th e MSSE should b e close to the v ector [1 1 · · · 1] ′ , while MAPE and MAE should b e as small as p ossible in absolute v alue. Note that the MAPE, as a p ercent age statistic, mak es sense only for a p ositiv e v alued pro cess y t , for all t . If this is not the case, then MAPE can n ot hav e a meaningful in terpretation and it sh ou ld b e excluded from the statistical analysis (Chatfield 38 ). 3 The Ba y esian Con trol Chart 3.1 The Main I dea Ba yes’ factors ha ve b een extensiv ely discussed in the statistics literature and recen tly they ha v e b een applied sequenti ally for time series, see e.g. W est and Harrison 28 (Chapter 11). Salv ador and Gargallo 40 prop ose a monitoring sc heme, based on Ba y es’ fact ors, for m ultiv ari- ate time series, b ut this approac h is not suitable for con trol c harting, b ecause it is applied in a mo del selection pr oblem. In add ition to this, most of the Ba yesia n time series monitoring (including th e w ork of Salv ad or and Gargallo 40 ) relies up on sim u lated based metho d s and in particular Mon te Carlo simulatio n. I n this pap er we fa vour non-iterativ e tec h n iques, b ecause they are f aster, more flexible and easier to apply . Once w e hav e the distribu tion (5 ) we can constru ct a target d istribution for the d isp ersion of y t from the target mean and then compare these t w o distributions. It is well kno wn (see e.g. P an and Jarrett 24 ) that the f orecast errors e i and e j ( i 6 = j ) are appro ximately uncorrelated and the appro ximation is so go o d as S t is closer to Σ. Sup p ose no w th at the target mean of { y t } is denoted b y µ and th e p ro cess disp ersion co v ariance matrix is denoted by V . This notation is consistent with th e S hewhart-Deming mo del as in equation (1), with V = Σ so that E ( y t ) = µ and V ar( y t ) = Σ, w here V ar( y t ) denotes the co v ariance matrix of y t . Is is assumed that µ is a generally u nkno wn v ector, bu t not sto chastic . In our m o del of equation (2) w e ha v e E ( y t | µ t ) = µ t and V ar( y t | µ t ) = Σ , bu t no w µ t is sto chastic and it also c hanges with time according to the random walk m o del of (2 ). W e p ostulate that, if the p ro cess is in con trol, the one step forecast mean of y t will b e close to the target mean v ector µ and the forecast co v ariance matrix of y t will b e close to the target disp ersion co v ariance matrix V . Thus we can define the target error d istr ibution b y ε t ∼ N p (0 , V ), where ε t = y t − µ is the p ro cess error, also kno w n in the pro cess adju stmen t literature (del Castillo 3 ) as d isturbance drift. Here we assume that V is p ositiv e definite matrix, although the pr op osed app roac h can b e mo difi ed when V is p ositiv e semi-definite. According to the ab o v e p ostulate, if m o del (1 ) describ es w ell the in-con trol pro cess, d ensit y (5) should b e close to the ab o v e target distribu tion. In 4 order to find out “ho w close” it is, w e form the Ba y es’ factor at time t : B F ( t ) = f e ( t ) f ε ( t ) = f e ( e t | Σ = S t − 1 , y t − 1 ) f ε ( ε t ) , t = 1 , 2 , . . . , N , where f e ( t ) and f ε ( t ) denote the p r obabilit y densit y fu n ctions of e t and ε t , resp ectiv ely . F or consistency in the ab o ve equation w e n eed to make the conv entio n y 0 = ∅ (the n ull or empt y set). Since b oth densities f e ( t ) and f ε ( t ) are normal we ha ve B F ( t ) = s δ p det ( V ) ( δ + P t − 1 ) p det ( S t − 1 ) exp  ( y t − µ ) ′ V − 1 ( y t − µ ) / 2 − δ ( y t − m t − 1 ) ′ S − 1 t − 1 ( y t − m t − 1 ) / (2 δ + 2 P t − 1 )  , (6) where det( · ) denotes the d eterminan t of a square matrix. The Ba yes’ factor B F ( t ) tak es v alues from 0 to + ∞ . W e will sa y that the pro cess y t is in con trol at time t , if f e ( t ) = f ε ( t ), or if B F ( t ) = 1; otherw ise the pr o cess will b e out of con trol, at this time p oin t. An out of con trol signal migh t b e caused b ecause of a m ean shift (e.g. when E ( y t | y t − 1 ) = m t − 1 is significan tly differen t th an µ ) or b ecause of a disp ersion shif t (e.g. V ar ( y t | Σ = S t − 1 , y t − 1 ) = ( δ + P t − 1 ) S t − 1 /δ is significan tly different than V ). 3.2 The Mo dified EWMA Con trol Chart for Correlated Data A control chart for the Ba y es’ factor B F ( t ) can conclude whether B F ( t ) is close to 1 and thus whether the pro cess is in control or not. Since B F ( t ) is p ositiv e v alued, it is more conv enient to w ork with the logarithm of the Ba y es’ f actor LB F ( t ) = log B F ( t ) = p log δ / 2 + { log det ( V ) } / 2 − p { log( δ + P t − 1 ) } / 2 − { log det ( S t − 1 ) } / 2 + ( y t − µ ) ′ V − 1 ( y t − µ ) / 2 − δ ( y t − m t − 1 ) ′ S − 1 t − 1 ( y t − m t − 1 ) / (2 δ + 2 P t − 1 ) (7) and so we can construct an appr op r iate univ ariate con trol chart for LB F ( t ). In ord er to prop ose suc h a c hart we need to deal w ith t w o issu es: (a) the v alues of LB F ( t ) will b e serially correlate d and (b) the distribution of LB F ( t ) m ight not b e normal. Considering (a), in our d evelo pm ent it is clear that, from the d efinition of th e B F ( t ), either the original data y t are i.i.d. or auto-correlate d, the resulting data B F ( t ) (or LB F ( t )) will b e correlated and hence, if the Sh ewhart or any other control c hart is to b e used successfully , they should b e mo dified appropriately to accommo date for correlated observ ations. Many authors ha v e demonstrated that the Shewhart control c harts need to b e mo d ified in ord er to cater for s er ially correlat ed observ ations (V asilop oulos an d S tam b oulis 41 ; Sc hmid 42 ). Similarly , the EWMA needs also to b e mo dified and the resulting mo d ified EWMA control c hart has b een discussed in many articles including Schmid 43 and V anBrac kle and Reynolds 44 . Acco rd ing to Harris and Ross 43 ignoring serial correlation has a stronger effect in EWMA than in the Shewhart cont rol c hart, but as we will see later the EWMA con trol chart is pr eferable to Shewhart, b ecause it is more robust to the assumption of normalit y . On e could also consider the mo difi ed CUSUM c hart for correlated obs erv ations, but we will n ot further d iscuss this in the presen t pap er. Pro ceeding with (b) one needs to c hec k the assump tion of normalit y , b efore app lying a mo dified EWMA (or S hewhart or CUSUM) con trol chart. Borror et al. 46 studied th e ARL 5 p erforman ce of the EWMA and they su ggested that the EWMA with a smo othing parameter equal to 0.05 is very effectiv e, ev en in the presence of non-normalit y of the observ ations. This result agrees with Mon togomery 1 who states for th e EWMA “It is almost a p erfectly non-parametric (d istr ibution free) pro cedu re”. Mara ve lakis et al. 47 study th e robu stness to normalit y of th e EWMA b y tabulating characte ristics of the run length d istributions (e.g. ARL) for observ ations generated by sev eral gamma distribu tions. These results conclude that, for relativ ely lo w v alues of the damping parameter of the EWMA and for shifts in th e mean the EWMA con trol c hart can b e used, ev en in th e absence of normalit y . Mo reov er, if the pro cess is in-con trol follo w ing a symmetrical, but n ot n ormal, distribu tion, then the EWMA can b e applied s u ccessfully . T o the follo w ing w e look at the emp irical d istr ibution of LB F ( t ) wh en the p ro cess is in con trol and when it is out of con trol. W e generate 1000 v ectors from a b iv ariate normal distribution N 2 ( µ, V ) with µ =  0 0  and V =  1 2 2 5  and we generate 1000 v ectors for three out of control scenarios. In scenario 1 w e sim ulate data from N 2 ( µ d , V ) (deviations from the mean µ ); in scenario 2 w e sim ulate data from N 2 ( µ, V d ) (deviations from the cov ariance matrix V ); in scenario 3 we simula te data from N 2 ( µ d , V d ) (deviations from b oth µ and V ), where µ d =  0 . 5 0  and V d =  1 2 . 5 2 . 5 8  . Figure 1 sh ows the histograms of the LB F ( t ) for the ab o ve four s cenarios (one in con trol and three out of con trol scenarios). F rom th is figure we obs erv e that, although the d istribution of the LB F ( t ) for the in-con trol p r o cess (panel (a) in Figure 1) is not-normal, it is roughly symmetric. The distributions of the LB F ( t ) for the out of cont rol pro cesses app ear to b e sligh tly skew ed, b u t the histograms are n ot conclusive . Th e imp ortan t p oint is th e non- normalit y of the LB F ( t ) an d the symm etry of the d istribution of the in-control pro cess. This enables us to m ak e use of the mo dified EWMA con trol c hart, b ut w e note that the mo dified CUSUM con trol c hart can also b e u sed. A more formal confirm ation of th e non- normalit y of th e distribution of LB F ( t ) can b e carried out by the usin g s tandard tests of normalit y , ho wev er, here the histograms are deemed sufficien t to declare the non-normalit y of the distribution of LB F ( t ). W e u se a tw o phase control sc heme; in Ph ase I the mean µ and th e co v ariance matrix Σ are estimated and adjustmen ts are applied if necessary , while in Phase I I the EWMA cont rol c hart is app lied to detect an y changes in the mean of LB F ( t ). Thus we pr op ose the algorithm: Algorithm 1. Ther e ar e two phases: Phase I: We fit the DWR mo del (2) for a set of historic al data t = 1 , 2 , . . . , N ∗ , with N ∗ < N . We che ck the p erformanc e and ade quacy of the mo del via the M SSE, M APE and MAE over al l t = 1 , 2 , . . . , N ∗ and we p ossible apply adjustments to the DWR mo del, (e.g. adjustments in the me an level) so that we obtain optimal values m opt = m N ∗ , S opt = S N ∗ , δ = δ opt ensuring that in Phase I the mo del matches the in- c ontr ol pr o c ess. The mo difie d EWMA c ontr ol chart is applie d so that c ontr ol limits ar e ade quately define d ac c or ding to pr e- sp e cifie d ARL curves. F or this to b e designe d, a state-sp ac e mo del for the pr o c ess LB F ( t ) ne e ds to b e i dentifie d and her e simple AR and ARMA mo del ling wil l b e gener al ly ac c eptable. 6 Histogram of LBF (a) in control Frequency −6 −4 −2 0 2 4 6 0 100 300 500 Histogram of LBF (b) out of control (mean) Frequency −6 −4 −2 0 2 4 6 0 100 200 300 Histogram of LBF (c) out of control (variance) Frequency −6 −4 −2 0 2 4 6 0 100 200 300 400 Histogram of LBF (d) out of control (mean & variance) Frequency −5 0 5 10 0 100 300 500 Figure 1: Histograms of th e log Ba y es’ factor LB F ( t ) for an in-con trol pro cess (panel (a)) and out-of-con trol p ro cesses (panels (b)-(d)). Th e out of con trol scenarios considered are deviations f r om the m ean vec tor (pan el (b)), deviations from the co v ariance matrix (panel (c)) and deviations from b oth the mean ve ctor and the co v ariance matrix (panel (d)). Phase I I: We fit the DWR mo del with the mo del c omp onents fr om Phase I (e.g. δ = δ opt , m t = m opt , Σ = S opt and we apply a mo difie d E WMA c ontr ol chart at observations LB F ( t ) with the c ontr ol limits identifie d at Phase I, for t = N ∗ + 1 , N ∗ + 2 , . . . , N . In ord er to apply the mo dified EWMA control c hart w e first calculate the series z t with observ ations x t = LB F ( t ) as z t = λx t + (1 − λ ) z t − 1 , 0 < λ ≤ 1 . (8) The parameter λ is the EWMA smo othing p arameter and as it is ment ioned ab o v e, for λ = 0 . 05 or λ = 0 . 1 the control chart is robu s t to normalit y . Then, the con trol limits of th e mo dified EWMA con trol chart are µ z ± cσ z , (9) where µ z = E ( z t ), σ 2 z = lim t →∞ V ar( z t ) (asymptotic v ariance of z t ) and c > 0 is determined according to the required ARL. F or AR(1) dep end ence x t = φx t − 1 + ν t and for large t , the 7 asymptotic v ariance σ 2 z is σ 2 z = σ 2 λ { 1 + φ (1 − λ ) } (1 − φ 2 )(2 − λ ) { 1 − φ (1 − λ ) } , where ν t ∼ N (0 , σ 2 ) an d σ 2 , φ are assu med kno w n. In practice these parameters are estimated at Phase I. According to Schmid 43 the asymptotic v ariance σ 2 z p erforms b etter than the exact v ariance of z t , whic h is giv en in Schmid 43 and which pro d uces time-dep endent control limits. Most of th e literature on this topic fo cuses on deriving the v ariance σ 2 z assuming simple time series mo d els for x t , e.g . as in the abov e AR(1) or as in the ARMA(1,1) mo del considered in V anBrac kle and Reynolds 44 . Algorithm 1 can b e simplifi ed , if at Ph ase I, the quan tities P t and S t con v erge to stable v alues and these v alues are determined in Phase I f or b oth ph ases. Th is b r ings up a we ll kno wn pr ob lem, whic h has r eceiv ed considerable atten tion in the time series literature (see e.g. Durbin and Ko opman 31 ). How eve r, for the D WR and s imilar multiv ariate mo dels limiting results for P t and S t ha v e not b een y et established. The n ext theorem (wh ic h pro of is in the app end ix) states that P t and S t con v erge to stable limiting v alues. Theorem 1. In the DWR mo del (2) the estimator S t of the me asur e ment c ovarianc e matrix Σ c onver g es in pr ob ability to Σ and the non-sto chastic sc alar p ar ameter P t c onver ges to the limit P = ( √ δ 2 + 4 − δ ) / 2 , i.e. S t P − → Σ and P t − → P . F rom Theorem 1 the estimator S t is consisten t and from the pro of of this theorem (giv en in the app en d ix), S t is also unbiased estimator. Theorem 1 suggests th at P t − 1 in the calculation of LB F ( t ) of equation (7) can b e r eplaced by its limit P . F rom equation (3) and Theorem 1 , the forecast of y t , m t − 1 can b e appro ximated b y m t − 1 = m 0 + P δ + P t − 1 X i =1 e i = m 0 + √ δ 2 + 4 − δ √ δ 2 + 4 + δ t − 1 X i =1 e i , where P t − 1 of equation (7) is rep laced b y P . Figure 2 shows how fast { P t } con v erges to its limit P , for a p r ior P 0 = 1 / 10 00 and three v alues of δ . This figure p oin ts out that P t is b ound ed ab ov e b y 1, but for δ = 0 . 2, this b ound is only achiev ed after t > 13 (solid lin e in Figure 2), while f or δ = 0 . 9, this b ound is ac hiev ed for an y t > 1 (dotted line in Figure 2). This giv es an emp irical indication of the sp eed of con v ergence of { P t } , for several v alues of δ . The limit P is kno wn b efore the algorithm starts (e.g. P dep ends only on δ ) and, giv en enough data in Phase I, the limit Σ can b e app ro ximated by Σ ≈ S N ∗ , in the end of Ph ase I. T his can h a v e an additional b enefit on computational sa vings, b ut more imp ortan tly it giv es a theoretical ju s tification that th e DWR pro duces a goo d cop y of the pro cess { y t } and therefore this mo del is appr op r iate for the monitoring part at Ph ase I I of Algorithm 1. F or example, if P t and S t w ere not con verging to stable v alues, no matter ho w man y data we collect ed at Phase I, the co v ariance matrix of y t and thus its uncertain ty w ould c hange o v er time resulting in an uns table time series mo del. F alse alarms are probab le in the fr amew ork of suc h unstable mo dels, whic h should b e a vo ided. In the design and application of the con trol c hart it is imp ortant to suggest v alues of m 0 , P 0 , δ and S 0 and to stud y their sensitivit y an d infl uence to the p erformance of the pr op osed con trol c hart. Since these suggestions are related to f orecasting as in equation (5), r esults on the sensitivit y of su ch prior parameters follo w f rom T rian tafyllop oulos and Piko ulas 36 and 8 Convergence of Pt Time Pt 0 10 20 30 40 1 2 3 4 5 Figure 2: Rate of conv ergence for the sequence { P t } of Theorem 1; the solid line p lots { P t } for δ = 0 . 2, the dashed line plots { P t } for δ = 0 . 5, the dotted line plots { P t } for δ = 0 . 9 and the dashed/dotted line is the critical b oun d of 1. T riant afyllop oulos 37 . It is worth while noting that, giv en en ou gh d ata in Phase I, the v alues of m 0 , P 0 and S 0 are not critical to th e forecast p erformance, as in time series mo delling prior information is d eflated o ve r time. Th is is in dicated in T h eorem 1 from th e fact that P do es not dep en d on P 0 . The v alue of δ ca n b e critical in forecasting and a general r ecommendation is that s everal v alues of δ (in the range of (0 , 1)) are applied in Phase I and ac cording to the forecast p erformance (see Section 2) a v alue of δ is decided. One should n ote that h igh v alues of δ (e.g. δ = 0 . 9) yield smo oth forecasts with lo w forecast v ariances, bu t these f orecasts are sometimes unable to forecast abrupt changes in the d ata; lo w v alues of δ (e.g. δ = 0 . 1) yield more p r ecise forecasts in the presence of “wild d ata”, but these forecasts come with increased forecast v ariances. Our prop osal for the mo dified EWMA cont rol c hart for the LB F ( t ) pro cess is motiv ated from the fact that the obser v ations LB F ( t ) p ossess auto correlation an d non-normalit y . The approac h is mo d el-based, and so a comparison with traditionally used multiv ariate control c harts, su ch as the Hotelling’s T 2 and the M-EWMA (which are b oth data-based con tr ol c harts), is difficult and in m an y o ccasions it can not giv e j ustice. Within the mo d el-based con trol charting metho ds , it app ears that our approac h can b e compared w ith the residu al c hart (P an and Jarrett 24 ), bu t again the comparisons need to mak e sure that mo del uncer- tain t y (w h ether for example th e D WR is a goo d mo del or an alternativ e time series mo del p erforms b etter) should b e ideally remov ed b efore any comparison is attempted. F or example a m iss-sp ecification of a time series mo del migh t result to a false result in the comparison 9 of the comp eting con trol c harts. F r om our exp erience the DWR works generally well (since it is a generaliza tion of the Shewh art-Deming mo del), b ut this m ight not b e the case for ev ery multiv ariate pro cess. W e b eliev e that such a comparison should d eserv e the length and the detail of a whole pap er and thus h ere we do n ot pu rsue th is pro ject. Next we giv e t wo examples illustrating the design and application of the prop osed con trol c hart. T o the ab o ve we ha v e assu med that give n a pr o cess { y t } the int erest is in bu ilding a con trol c hart f or monitoring sim ultaneously the p ro cess mean and the disp ersion co v ariance matrix. How ev er, in some cases the interest is p laced on monitoring the disp ersion co v ariance matrix only . In this case we can mo d ify the con trol sc heme b y consider in g a mo d ified EWMA con trol c hart of the log-Ba yes’ f actors of the fir s t order difference pro cess z t = y t − y t − 1 , whic h from equation (2) has zero mean. Con trol c h arts based on { z t } w ill b e more r obust as compared to those for { y t } , sin ce the uncertain ty of monitoring the p ro cess mean of { y t } has b een remo v ed. 4 London M etal Exc hange Data London metal exc hange (LME) is the w orld’s premier n on-ferrous metals market tr ad in g currentl y alum in ium, copp er, lead and zinc, among other n on -fer r ous metals. Information on the LME and its functions can b e foun d in its web site: ht tp://www .lme.co. uk . The review of W atkins and McAleer 48 explores th e recen tly gro wing literature on the LME marke t and T riant afyllop oulos 37 discusses the correlation of sp ot and fu ture cont ract prices of aluminium based on the D WR mo del of Section 2. In this p ap er we d iscuss data of sp ot prices for th e four metals aluminium (v ariable { y 1 t } ), copp er (v ariable { y 2 t } ), lead (v ariable { y 3 t } ) and zinc (v ariable { y 4 t } ). The data are collect ed from January 2005 until Octob er 2005 f or eve ry trading da y ex- cluding wee kends and bank holida ys; Figure 3 p lots the data. W e form the observ ation vect or y t = [ y 1 t y 2 t y 3 t y 4 t ] ′ and we are interested in kno wing whether volat ilit y is apparent , for t = 151 until t = 220. In other words w e wan t to kno w whether from t to t + 1, the v ariabilit y of the observ ations y t and y t +1 has c hanged. Th is is a ma jor concern to econometricians, b ecause if there is evidence for vo latilit y , this means there is u ncertain t y in inv estments and ideally the v olatilit y should b e understo o d and explained. In order to answ er this imp ortant question w e form the first order difference of the series { y t } , defined by x t = y t − y t − 1 , for t > 1 (Figure 4). Adopting the u s ual forecasting str ategy of commo d it y forecasting, giv en data up to time t − 1, the forecast mean of y t at time t is ju s t the v alue of y t − 1 and so w e can write E ( y t | y t − 1 ) = y t − 1 . W e note th at the tru e mean of x t ma y not b e zero (un less in mo d el (2) it is µ t = µ + ω t ), but it is tru e that conditionally on y t − 1 or y t − 1 w e ha v e E ( x t | y t − 1 ) = E ( y t ) − y t − 1 = 0, since E ( y t | y t − 1 ) = y t − 1 . F r om Figure 4 we obser ve that th e series { x t } fluctuates around zero and vola tilit y can b e detected as significan t deviations from the zero target ; s uc h deviations can b e detecte d w ith the aid of a con trol chart of Section 3. First w e need to mak e sure th at th e D WR mo del fits the differenced series { x t } wel l. W e tak e t = 1 − 150 as Phase I, in wh ic h the adequacy of the D WR mo del is ev aluated. The p erforman ce statistics of S ection 2 are: M S S E = [0 . 993 1 . 486 0 . 866 1 . 323 ] ′ and M AE = [18 . 93 2 45 . 187 14 . 569 19 . 082 ] ′ , suggesting an acceptable fit. Of course the MAPE is not a v ailable, since { x t } is n ot a p ositiv e v alued pro cess (Secti on 2). W e ha v e designed a mo dified EWMA con trol c hart for the LB F ( t ) of the p ro cess { x t } according to the discu s sion of Section 3. Figure 5 sho ws four co ntrol c harts corresp ond ing to 10 1700 1800 1900 2000 3200 3600 4000 850 900 950 1000 1200 1300 1400 1500 0 50 100 150 200 London metal exchange data Time y 4 y 3 y 2 y 1 Figure 3: LME data y t = [ y 1 t y 2 t y 3 t y 4 t ] ′ , consisting of aluminium ( { y 1 t } ), copp er ( { y 2 t } ), lead ( { y 3 t } and zinc ( { y 4 t } ) sp ot prices (in US dollars p er tonne of eac h metal). four v alues of the EWMA smo othing parameter λ . Typically the con trol c hart is robu st to normalit y for small v alues of λ , but for these v alues the con trol chart is only d etecting v ery small d rifts in the mean this migh t not b e desirable. As λ in creases th e mo dified EWMA con trol c hart is losing its robustness o ver normalit y , bu t for symmetric pro cess distributions, suc h as the empirical d istr ibution of the LB F ( t ) sho wn in Figure 1, the E WMA con trol c hart might still b e used for λ = 0 . 5. The correlation of the LB F ( t ) is accoun ted b y th e autoregressiv e mo del of S ection 3 and an analysis inv olving the d ata at Phase I sho ws that an the autoregressiv e parameter φ = 0 . 1 is adequate to capture the auto correlation of LB F ( t ). According to T ables for the ARL of the mo dified EWMA con trol c hart (see e.g. Shiau and Hsu 49 ) we c ho ose the v alue of c in equation (9) so that AR L = 370 . 4, e.g. for λ = 0 . 05 and φ = 0 . 1 w e ha v e c = 2 . 469. The remainder of the con trol limits are calculated as in equation (9). Figure 5 shows that the pro cess in Ph ase I I app ears to b e in con trol, for λ = 0 . 05 and λ = 0 . 1, while f or λ = 0 . 2 and λ = 0 . 5 the cont rol c hart returns an out of con trol p oint at t = 172 (with v alues z 172 = − 1 . 852 and z 172 = − 2 . 999, resp ectiv ely). Th e mean of th e EWMA z t is sligh tly lo w er th an zero, wh ic h indicates that, for the entire pr o cess { x t } , there will b e some deviation of the pr edictiv e d ensit y f e from the target density f ε . It is up to the mo deller to d ecide whether such deviatio n fr om the target distribution is w orth of declaring the pro cess out of control . In searc h of a m ore automatic appr oac h , one can lift up the whole con trol chart so that in Phase I the mean of z t is exac tly zero. This can b e p er f ormed automatica lly , in the end of Phase I, and this will declare the pro cess in cont rol in Phase I I, 11 −80 −40 0 20 40 60 −150 −50 0 50 100 150 −60 −40 −20 0 20 40 −50 0 50 0 50 100 150 200 London metal exchange differenced process Time x 4 x 3 x 2 x 1 Figure 4: LME differenced pro cess x t = [ x 1 t x 2 t x 3 t x 4 t ] ′ , consisting of alum in ium ( { x 1 t } ), copp er ( { x 2 t } ), lead ( { x 3 t } ) and zinc ( { x 4 t } ). T h e horizont al lines, placed at zero, indicate no v olatilit y . for λ = 0 . 05 , 0 . 1, wh ile for λ = 0 . 2 , λ = 0 . 5 there is an out of con tr ol p oin t at t = 172. I n Figure 5 the v alue of λ = 0 . 5 is rather high to ensuring correct con trol limits of the mo dified EWMA c hart (see the relev ant d iscussion in p age 7); here the c hart with λ = 0 . 5 is mainly sho wn for comparison p urp oses with the charts w ith lo w er v alues of λ , but in practice w e suggest th at λ do es n ot exceed 0.2, u nless there is strong evidence to su pp ort the assu mption of n ormalit y for the distribution of LB F ( t ). It is w orth p oin ting out th at th e concen tration of consecutiv e E WMA v alues under the mean in Ph ase I I is causing warning, whic h is apparent in all c harts. The phenomenon is more apparent in the c harts for λ = 0 . 05 and λ = 0 . 1 and it can suggest the out of con trol state of the pro cess at t = 172,whic h is apparent in the c h arts with λ = 0 . 2 and λ = 0 . 5. The interpretatio n of the ou t of con trol signal at t = 172 can not b e d one ju st by lo oking at Figure 4 and more dedicated metho ds of out of con trol v ariable iden tification need to b e emplo yed, see e.g. Bersimis et al. 14 . 5 Pro d uction Time Series Data In an exp erimen t of pro d uction of a plastic mould the qualit y is cen tered on the con trol of temp erature and its v ariation. F or this pu rp ose fiv e measurements of the temp eratur e of the mould hav e b een take n, for 276 time p oin ts. The exp erimen t is fully describ ed in Pa n and Jarr ett 24 and these authors sho w that this 5-dimensional pro du ction pro cess { y t } is b oth auto correlated and serially correlat ed including b oth v ector autoregressiv e and mo ving a v er- 12 Modified EWMA chart (a) lambda = 0.05 LBF 0 50 100 150 200 −1.0 −0.5 0.0 Modified EWMA chart (b) lambda = 0.1 LBF 0 50 100 150 200 −1.5 −1.0 −0.5 0.0 Modified EWMA chart (c) lambda = 0.2 LBF 0 50 100 150 200 −2.0 −1.0 0.0 0.5 Modified EWMA chart (d) lambda = 0.5 LBF 0 50 100 150 200 −3 −2 −1 0 1 Figure 5: Mo dified EWMA cont rol chart for the log Bay es’ factor of th e LME differenced pro cess. Plots (a)-(d) show the mo dified cont rol chart for differen t v alues of the s mo othing parameter λ . In eac h plot of the p anel, the solid horizon tal lin e indicates the mean of the EWMA and the dotted h orizon tal lines indicate the con trol limits; the v ertical line separates Phase I (for t = 1 − 150) and Phase I I (for t = 151 − 210 ). age terms. These authors us e a v ector state space c harting approac h b ased on the Hotelling con trol c hart resulting on 12 out of con trol signals at Ph ase I I (time p oints from t = 181 to t = 220) and hence concluding that the pro cess falls badly out of con trol at Phase I I. W e ha v e used the d ata at Phase I (time p oin ts t = 1 − 180) in order to estimate the target mean vect or µ = [208 . 24 5 153 . 638 53 . 063 − 22 . 742 16 . 126] ′ (as the a ve rage of eac h y it : t = 1 − 180) and the disp ersion co v ariance matrix V =       0 . 168 − 0 . 001 0 . 633 − 0 . 438 0 . 015 − 0 . 001 0 . 023 − 0 . 017 0 . 006 − 0 . 002 0 . 633 − 0 . 017 25 . 621 − 15 . 658 0 . 453 − 0 . 438 0 . 006 − 15 . 65 8 14 . 181 − 0 . 596 0 . 015 − 0 . 002 0 . 453 − 0 . 596 0 . 951       (as the samp le cov ariance m atrix of eac h y t : t = 1 : 180), where y t = [ y 1 t , y 2 t , y 3 t , y 4 t , y 5 t ] ′ . Th e D WR fits w ell with M S S E = [0 . 855 0 . 950 0 . 992 1 . 161 0 . 996] ′ , whic h is close to [1 1 1 1]. T he other t w o p erf orm ance s tatistics are M AE = [1 . 378 0 . 899 4 . 450 3 . 316 0 . 945] ′ and M AP E = [0 . 007 0 . 006 0 . 089 − 0 . 0 59] ′ , where for { y 4 t } the “–” indicates that the MAPE is not av ailable, since this v ariable is not p ositive v alued (see th e relev an t discussion for MAPE in Section 2). The ab o v e p er f ormance statistics s uggest that the mo del fi t is go o d and therefore we can pro ceed with con trol c harting at Phase I I ( t = 181 − 279 ). The fir st thing to do is to fin d a su itable AR(1) mo del for the p r o cess LB F ( t ). A s uitable mo del is the AR(1): LB F ( t ) = − 4 . 624 + 0 . 062 LB F ( t − 1) + ν t . According to the d iscussion 13 Modified EWMA chart (a) lambda = 0.05 LBF 0 50 100 150 200 250 −1.0 −0.5 0.0 0.5 1.0 Modified EWMA chart (b) lambda = 0.1 LBF 0 50 100 150 200 250 −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5 Figure 6: Mo dified EWMA control c hart for the log Ba yes’ facto r of the Pro d uction p ro cess. Plots (a)-(b) sh o w tw o c harts for v alues of the smo othing parameter λ = 0 . 05 and λ = 0 . 1. F or b oth plots, the solid horizonta l line indicates the target mean 0 and the dotted horizonta l lines indicate the cont rol limits; th e solid vertica l line separates Phase I (for t = 1 − 180) and Phase I I (for t = 181 − 276). ab o ve , we r emov e the interce pt − 4 . 624 so that w e can obtain a in-con trol pro cess in Phase I. Th us we d esign the m o dified EWMA con trol c hart f or LB F ( t ) + 4 . 624. Again w e use tables for the mo d ified EWMA con trol c hart and for λ = 0 . 05 the resulting control c hart is give n in Figure 6. This figur e agrees with the residual c hart of Pa n and Jarr ett 24 , that finds the pro cess in Phase I I out of control for most of the data p oin ts. In Ph ase I chart of panel (b) of Figure 6 gives one out of con tr ol p oint , whic h is in agreemen t with P an and Jarrett 24 , but in panel (a) of Figur e 6 the control c hart detects more out of control p oint s in Ph ase I. The EWMA con trol c hart is robu st to non-normalit y f or the low v alues of λ = 0 . 05 and λ = 0 . 1, but f or λ = 0 . 05 th e chart is more sensitiv e to small shifts in the mean of LB F ( t ), resu lting to the detect ion of out of con trol p oin ts in Phase I. Any out of con trol p oin ts in Ph ase I should b e immediately inv estigated and usual SPC pro cedures of r emo ving influence of these p oints in the calc ulation of the con trol limits should b e applied (Mon tgomery 1 ). 6 Conclusions This pap er d ev elops a new multiv ariate con trol chart based on Ba ye s’ factors. This con trol c hart is sp ecifically aimed at m ultiv ariate auto correlated and serially correlate d pro cesses. The general idea is to form a target distrib ution, to construct a pr edictiv e densit y with go o d forecast ab ility and then to apply a univ ariate con trol c hart for the logarithm of th e Ba yes’ factor of the predictiv e error dens it y against the target error densit y . Although in this 14 pap er, for simplicit y , we ha v e consid ered normal distribu tions for the target and the pr edictiv e densities, in general applicatio n the p rop osed con trol c harts can b e applied considerin g other densities to o as long as they are a v ailable in analytic form. W e ha ve restricted our d iscussion to th e mo d ified EWMA con trol c hart, but other con trol c harts suc h as the mo dified CUSUM and n on-parametric con trol c harts can b e app lied. A ma jor adv antag e of our approac h as compared to other multi v ariate con trol charts is that once w e h a v e obtained the log Ba y es’ factors we can app ly an y appr opriate univ ariate con trol c hart. A difficult y app ears to b e that the resulting Ba y es’ factors pro cess is b oth auto correlated and non-normal, but w e b eliev e th e design of the p rop osed chart is a c hallenge that can attract and motiv ate further researc h in this so imp ortant area of statistica l p ro cess control. Ac kno wledgemen ts I sh ou ld lik e to thank the editor E r ik Mønness and t wo anonymous referees for m aking several v aluable suggestions, whic h considerably impro ved the pap er. App endix Pr o of of The or em 1. First w e pr o v e S t P − → Σ. It suffices to pr o v e that S t is unbiased estimator and that its co v ariance matrix conv erges to zero. F rom equations (4) and (5) we obtain E ( S t ) = 1 t t X i =1 δ E ( e i e ′ i ) δ + P i − 1 = 1 t t X i =1 δ ( δ + P i − 1 )Σ ( δ + P i − 1 ) δ = 1 t ( t Σ) = Σ and so S t is u n biased for Σ . F or the con vergence, let vec h( · ) denote the column stac king op erator of a low er p ortion of a co v ariance matrix and let k · k denote a matrix n orm defined in a suitable linea r space. F rom equation (5) w e ha ve V ar { v ec h ( S t ) } = 1 t 2 t X i =1  δ δ + P i − 1  2 V ar { v ec h ( e i e ′ i ) } . (A-1) F rom equation (5) e i follo w s a p -v ariate normal distr ibution and so by writing e i = [ e i 1 e i 2 · · · e ip ] ′ , w e ha v e that Co v( e ij , e ik ) = E ( e ij e ik ) are b ounded, since these exp ectations are expressed as momen ts of the m ultiv ariate normal distribution (T rian tafyllop oulos 50 ). Hence V ar { v ec h ( e i e ′ i ) } has finite elemen ts and so we can write k V ar { v ec h ( e i e ′ i ) } k < M , f or s ome M > 0. F or an y ǫ > 0 define t 0 = [ ǫM ] (the in tegral part of ǫM ). F r om P i − 1 > 0 we ha v e that δ / ( δ + P i − 1 ) < 1, for all i = 1 , 2 , . . . t . Then k V ar { v ec h ( S t ) }k = 1 t 2      t X i =1  δ δ + P i − 1  2 V ar { v ec h ( e i e ′ i ) }      ≤ M t 2      t X i =1  δ δ + P i − 1  2      ≤ tM t 2 = M t < ǫ, for an y t > t 0 . This sho w s that lim t →∞ V ar { v ec h ( S t ) } = 0 and so S t P − → Σ. 15 Pro ceeding no w with { P t } we sho w that { P t } is a Cauc hy s equence in the real line and hence lim t →∞ P t = P exists. T o prov e that { P t } is a Cauc hy sequence, it suffices to pro ve that lim t →∞ | P t − P t − 1 | = 0, where | · | denotes absolute v alue. First we sh o w that exists p ositiv e int eger t 0 suc h that for all t > t 0 it is P t < 1. The pro of of th is is by con tradiction. Supp ose that for all t 0 exists t > t 0 suc h that P t ≥ 1. Without loss in generalit y take t 0 = t ∗ and P t ∗ = 1. Then we see that P t ∗ +1 = 1 / ( δ + P t ∗ ) = 1 / ( δ + 1) < 1, P t ∗ +2 = 1 / ( δ + P t ∗ +1 ) = ( δ + 1) / ( δ 2 + δ + 1) < 1 and like wise P t ∗ + k < 1, for all k ≥ 1. So we can pic k t 0 = t ∗ + 1 so that w e can not fi nd any t > t 0 with P t ≥ 1, whic h contradicts th e hypothesis. Thus exists t 0 > 0 so that for all t > 0 it is P t < 1. This in turn implies that δ + P t − 1 > 1 , ∀ t > t 0 . (A-2) F rom the definition of P t of equation (3), w e obtain P t − P t − 1 = 1 δ + P t − 1 − 1 δ + P t − 2 = − P t − 1 − P t − 2 ( δ + P t − 2 )( δ + P t − 2 ) = · · · = ( − 1) t − 1 ( P 1 − P 0 ) Q t − 1 i =1 ( δ + P t − i )( δ + P t − i − 1 ) . No w pic k t 0 as in (A-2) and define M = min { δ + P t − 1 , ( δ + P t − 2 ) 2 , . . . , ( δ + P t 0 +1 ) 2 } so that M > 1. Then | P t − P t − 1 | = | 1 − δ P 0 − P 2 0 | Q t 0 i =0 ( δ + P i ) 2 Q t − t 0 − 2 i =1 ( δ + P t − 1 )( δ + P t − i − 1 ) < | 1 − δ P 0 − P 2 0 | Q t 0 i =0 ( δ + P i ) 2 M t − t 0 − 1 → 0 , since lim t →∞ M t − t 0 − 1 = + ∞ . This pro v es that lim t →∞ | P t − P t − 1 | = 0 and so { P t } is a Cauc hy sequence. Thus lim t →∞ P t = P exists and from equ ation (3) we ha ve P = 1 / ( δ + P ), for whic h w e derive P = ( √ δ 2 + 4 − δ ) / 2, after rejecting the negat ive root P = ( − √ δ 2 + 4 − δ ) / 2. References [1] Mon tgomery DC. Intr o duction to Statistic al Quality Contr ol. Wiley , 4th edition: New Y ork, 2000. [2] Bo x GEP , Luce ˜ no A. Statistic al Contr ol b y M onitoring and F e e db ack A djustment. Wiley: New Y ork, 1997. [3] Del Castillo E. Statistic al Pr o c ess A djustment for Qu ality Contr ol. Wiley: New Y ork, 2002. [4] T racy ND, Y oun g, JC, Mason RL. Mu ltiv ariate control charts for ind ividual observ ations. Journal of Quality T e chnolo gy 1992 24 : 88-95. [5] Liu R Y. Con trol c harts for multiv ariate p ro cesses. Journal of the Americ an Statistic al Asso ciation 1995 90 : 1380 -1387. [6] Kourti T, MacGregor JF. Mu ltiv ariate SPC metho ds for p ro cess and pr o duct monitoring. Journal of Quality T e chnolo gy 1996 28 : 409-42 8. [7] Mason RL, Chou YM, Y oung JC . Applying Hotelling’s T-2 statistic to b atc h pro cess. Journal of Quality T e chnolo gy 2001 33 : 466-47 9. [8] V argas NJA. Robust estimation in multiv ariate control c harts f or individu al observ ations. Journal of Quality T e chnolo gy 2003 35 : 367-37 6. 16 [9] Y e N, Borror CM, Pa rmar D. Scalable c hi-square d istance v ersus con ve ntional statistical distance f or p r o cess monitoring with uncorrelated data v ariables. Quality and R eliability Engine ering International 2003 19 : 505-515. [10] P an X. An alternativ e app roac h to m ultiv ariate EWMA con trol chart. Journal of Applie d Statistics 200 5 32 : 695-70 5. [11] Lo wry CA, Mon tgomery DC. A r eview of m ultiv ariate con trol c harts. IIE T r ansactions 1995 27 : 800-810. [12] Sulliv an JH, W o o dall WH. A comparison of multi v ariate control c harts for individual observ ations. Journal of Q uality T e chnolo gy 1996 28 : 398-40 8. [13] Mon tgomery DC, W o o dall WH. A discussion on statisticall y-based pr o cess monitoring and con tr ol. Journal of Quality T e chnolo gy 1997 29 : 157-162 . [14] Bersimis S , P s arakis S, Panareto s J. Multiv ariate statistical pro cess con trol charts: an o v erview. Qu ality and R eliability Engine ering International 200 6 (to app ear). [15] Y eh, AB, Lin, DK-J, McGrath, RN. Multiv ariate con trol c harts for monitoring co v ariance matrix: a review. Quality T e chnolo gy and Quantitative Management 2006 (to app ear). [16] Y eh AB, Lin DKJ , Zhou H, V enk ataramani C. A multiv ariate exp onential ly w eigh ted mo ving av erage con trol c hart for monitoring pro cess v ariabilit y . J ournal of Applie d Statistics 2003 30 : 507-536. [17] Surtihadi J, Raghav ac hari M, Ru nger G. Multiv ariate con trol c harts for pro cess disp er- sion. International Journal of Pr o duction R ese ar ch 2004 42 : 2993-3009 . [18] Cheng, SW, Thaga, K . Multiv ariate max-CUSUM chart. Quality T e chnolo gy and Quan- titative M anagement 200 5 2 : 221-235. [19] Costa AFB, Rahim MA. Monitoring Pro cess Mean and V ariabilit y with One Non-cen tral Chi-square Chart. Journal of A pplie d Statistics 2004 31 : 1171-1183. [20] Dy er JN, C on er ly MD, Adams BM. A sim ulation stud y and ev aluation of multiv ariate forecast based control charts applied to ARMA pro cesses. Journal of Statistic al Computa- tion and Simulation 2003 73 : 709- 724. [21] Jiang W. Multiv ariate cont rol charts for monitoring auto correlated pr o cesses. Journal of Quality T e chnolo gy 2004 36 : 367-379. [22] Kalgonda, AA, K u lk arni, SR. Multiv ariate qualit y con trol c hart for auto correlated pro- cesses. Journal of Applie d Statistics 2004 31 : 317- 327. [23] No orossana R, V aghefi SJ M. Effect of auto correlation on p erformance of the MCUSUM con trol chart. Quality and R eliability Engine ering International 2006 DOI: 10.1002 /qre.695 [24] P an X, Jarrett J. App lying state space to SPC: m on itoring m ultiv ariate time series. Journal of Applie d Statistics 2004 31 : 397-41 8. [25] T agaras G. A sur vey of recen t dev elopmen ts in the design of adaptive con trol c h arts. Journal of Quality T e chnolo gy 1998 30 : 212-23 1. 17 [26] T agaras G, Nik olaidis Y. Comp aring the effecti veness of v arious Ba ye sian x control c harts. Op er ations R ese ar ch 2002 50 : 878-888. [27] De Magalh˜ aes MS, C osta AFB, Neto FDM. Ad aptiv e con trol charts: a Marko vian ap- proac h for pro cesses su b ject to ind ep endent disturb ances. International Journal of Pr o duc- tion Ec onomics 2006 99 : 236-246. [28] W est M, Harrison PJ. Bayesian F or e c asting and Dynamic Mo dels (2nd edn). Springer: New Y ork, 1997. [29] Deming WE. Out of the Crisis. Massac husetts Insitute of T ec h nology Cent er for Ad- v anced Engineering Study: C am bridge, MA, 1986. [30] T rianta fyllop ou los K, God olphin JD, Go dolphin EJ. Pr o cess impro v ement in the micro- electronic ind ustry b y state space mo delling. Quality and R e liability Engine ering Interna- tional 2005 21 : 465-475. [31] Durbin J , Ko opman SJ. Time Series Analysis by State Sp ac e Metho ds. Oxford Univ ersit y Press: Oxford, 2001 . [32] Bro wn R G. Smo othing, F or e c asting and Pr e diction of Di sc r ete Time Series. Englewoo d Cliffs, Pren tice Hall: Ney Jersey , 1962. [33] Ameen JR M, Harrison PJ. Discount weig hted estimation. Journal of F or e c asting 1984 3 : 285-2 96. [34] Ameen JRM. Sequen tial discoun t estimation. The Statistician 1988 37 : 227-23 7. [35] Go o dwin P . Ad j usting j udgemen tal extrap olations using Theil’s method and discoun ted w eigh ted r egression. Journal of F or e c asting 1997 16 : 37-46. [36] T rianta fyllop ou los K , Pik oulas J. Multiv ariate Ba yesia n regression applied to the problem of net work securit y . Journal of F or e c asting 2002 21 : 579-594. [37] T rianta fyllop ou los K . Multiv ariate d iscoun t weig hted r egression and lo cal level mo d- els. Computational Statistics and D ata Analysis 2006 DOI: 10.1016/j.c sd a.2005 .07.003 (in press). [38] Chatfield C. Time-Series F or e c asting . Chapman and Hall: New Y ork, 2001 . [39] Gupta AK, Nagar DK. Matrix V ariate Distributions. Chapman and Hall: New Y ork, 1999. [40] Salv ador M, Gargallo P . Automatic monitoring and in terven tion in m ultiv ariate dynamic linear mo dels. Computa tional Statistics and Data Analysis 2004 47 : 401-43 1. [41] V asilop oulos, A V, Stam b oulis, AP . Mod ification of cont rol c hart limits in the pr esence of correla tion. Journal of Quality T e chnolo gy 1978 10 : 20-30. [42] Sc hmid , W. On th e run length of a S h ewhart c hart for correlated d ata. Statistic al Pap ers 1995 36 : 111-130. 18 [43] Sc hmid , W. On EWMA charts for time series. In F r ontiers in Statistic al Quality Contr ol 1997 5 : 114-137. [44] V anBrac kle, LN, Reynolds, MR. EWMA and CUSUM control charts in the presence of auto correlation. Communic ations in Statistics: Si mulation and Computatio n 1997 26 : 979-1 008. [45] Harris, TJ, Ross, WH. Statistical pro cess control pr o cedures for correlated observ ations. Canadian Journal of Chemic al Engine ering 1991 69 : 48-57. [46] Borror, CM, Champ, CW, Ridgon, SE. Poisson EWMA con trol c harts. Journal of Qu ality T e chnolo gy 1998 30 : 352-36 1. [47] Mara v elakis, PE, P anaretos, J, Psarakis, S. An examination of the robu stness to non normalit y of the EWMA control c harts for the d isp ersion. Communic ations in Statistics: Simulation and Computation 2005 34 : 1069- 1079 [48] W atkins C , McAleer M. Econometric mo delling of n on-ferrous metal p r ices. Journal of Ec onomic Surveys 2004 18 : 651-701. [49] Shiau, J-JH, Hsu, Y-C. Robustness of the EWMA con trol chart to non-norm alit y for auto correlated pro cesses. Quality T e chnolo gy and Quantitative Management 2005 2 : 125- 146. [50] T rianta fyllop ou los, K. On th e ce ntral m omen ts of the m ultidimensional Gaussian distri- bution. Mathematic al Scientist 2003 28 : 125- 128. 19