Residual-based maximum MCUSUM control chart for joint monitoring the mean and variability of multivariate autocorrelated processes

ABSTRACT Maximum multivariate cumulative sum (Max-MCUSUM) is one of the single control charts proposed for joint monitoring the mean and variability of independent observation. Since many applications yield time series data, it is important to develop Max-MCUSUM control chart for monitoring multivariate autocorrelated processes. In this paper, we propose a Max-MCUSUM control chart based on the residual of multioutput least square support vector regression (MLS-SVR). The optimal parameters of MLS-SVR model are calculated using historical in-control data and the control limit of the proposed chart is estimated using the bootstrap approach. The average run lengths of MLS-SVR-based Max-MCUSUM chart verify that the proposed chart is more sensitive to detect mean vector shift than to detect a covariance matrix shift. The illustrative examples of the proposed control chart are also provided for both simulation and real data.


Introduction
Joint monitoring mean and variability is a bivariate problem and cannot be solved as two separate problems (Gan, 1997). If the presence of special causes may lead to the changes in both mean and variability at the same time then joint monitoring of mean and variability is more reasonable (Gan, Ting, & Chang, 2004). The shift in variability can influence the control limits of the mean chart so that the mean and variability are more effective to be monitored jointly. Thaga and Sivasamy (2015) classified control charts for joint monitoring into simultaneous control chart which was consisted of two statistics plotted on one chart and single control chart which only had single statistic as a representation of mean and variability plotted on one chart. A single control chart is categorized into a single chart with traditional control limit and a single chart with twodimensional control region (McCracken & Chakraborti, 2013).
Several single control charts are introduced to monitor multivariate quality characteristics comprising Max-MEWMA chart (Xie, 1999), Max-MCUSUM chart (Cheng & Thaga, 2005b), bivariate Max-chart (Khoo, 2004), and multivariate Max-chart (Thaga & Gabaitiri, 2006). Those control charts utilize single maximum statistic from two independent transformed random variables representing the mean and variability of a process. Cheng and Thaga (2005b) proposed Max-MCUSUM chart in which the reference value was calculated from actual observation, and its upper control limit (UCL) was estimated using the Markov Chain approach. Recently, Khusna, Mashuri, Ahsan, Suhartono, and Prastyo (2018) developed a bootstrap-based Max-MCUSUM control chart with two advantages. First, the researcher can adjust the level of tightness by predetermining the reference value in Phase I monitoring process. Second, the observations do not need to meet any distribution assumption because the UCL is computed using the bootstrap method.
Many applications show that the data collected over time violate the assumption of independent observation. As a result, many researchers had developed control charts for joint monitoring auto-correlated processes. Lu and Reynolds (1999) acquainted some simultaneous control charts such as simultaneous Shewhart chart and combined Shewhart-EWMA chart based on the residual of the first-order autoregressive process. Cheng and Thaga (2005a) proposed residual-based Max-CUSUM control chart and proved that for small shifts in both mean and variability of first-order autoregressive process, the proposed chart outperforms the combined Shewhart-EWMA chart. In addition, Thaga and Yadavalli (2007) verified that the residual-based Max-EWMA chart has better performance than residual-based Max-CUSUM chart for moderate to large shifts in the process mean and variability. The Max-chart for auto-correlated processes (Thaga, Kgosi, & Gabaitiri, 2007) is more effective than the simultaneous Shewhart chart (Lu & Reynolds, 1999) in detecting shifts of mean and variability for both low and high autocorrelation levels.
In a real application, some auto-correlated quality characteristics cannot be monitored separately. A chemical process is one of the examples of the manufacturing process in which observations are highly correlated because they are measured in every unit of time order. A chemical process often has multivariate quality characteristics with a small shift of the process. Cheng and Thaga (2005b) proved that Max-MCUSUM chart outperforms Max-MEWMA chart to detect small changes in the mean vector and covariance matrix of independent observation. Therefore, this study aims to propose the residual-based Max-MCUSUM chart for joint monitoring small shift of multivariate auto-correlated processes. This study prefers to utilize multi-output least square support vector regression (MLS-SVR) model (Xu, An, Qiao, Zhu, & Li, 2013) in order to estimate the auto-correlated observation because of two reasons. First, MLS-SVR model is useful to cover the potential nonlinear relationship as well as to cope with the underlying cross-relatedness among variables. Second, the convex linear equation system in MLS-SVR model is easy to be solved and good at the computational time.
The remaining parts of this paper are organized as follows. Section 2 provides a brief explanation about MLS-SVR model. The statistic, the bootstrap control limit, and the charting procedures of MLS-SVR based Max-MCUSUM chart are explained in Section 3. Section 4.1 presents the simulation studies to determine the inputs and hyperparameters selection of MLS-SVR model. Section 4.2 displays the simulation studies to calculate the UCL of the proposed chart using the bootstrap approach. The performance evaluation, as well as the comparison of the proposed chart, is described in Section 5. Section 6 illustrates the applications of the proposed chart to monitor both simulation and real data. Finally, some conclusions and future research opportunities are given in Section 7.

Multioutput least square SVR (MLS-SVR) model
Given an observable output matrix Y ¼ ½y ij 2 R nÂm , for i ¼ 1; 2; . . . ; n, j ¼ 1; 2; . . . ; m, with n as the number of observations and m as the number of MLS-SVR outputs. The specific independent and identically distributed samples are given as fðx 1 ; y 1 Þ; ðx 2 ; y 2 Þ; . . . ; ðx n ; y n Þg, where x i 2 R d , y i 2 R m , and d explains the dimension of the input vector. Let φ : R d ! R h is a mapping function to higher dimensional Hilbert space with h dimension. The parameter w j 2 R h ðj 2 N m Þ is assumed to associate with φðxÞ and it can be written as w j ¼ w 0 þ v j , for mean vector w 0 2 R h and vector v j 2 R h ðj 2 N m Þ. The mean vector w 0 shows commonality information whereas vector v j describes specialty information. The small value of v j indicates that the output is similar to each other. The mean vector w 0 , matrix V ¼ v 1 ; v 2 ; . . . ; v m ð Þ2R hÂm , and vector b ¼ b 1 ; b 2 ; . . . ; b m ð Þ2R m are simultaneously estimated by minimizing the optimization problem as follows (Xu et al., 2013): (1) 2 R hÂm describes MLS-SVR parameter, and γ 0 ; γ 00 2 R þ are regularized parameter. The optimization problem can be formed into Lagrange function as follows: where matrix A ¼ ðα 1 ; α 2 ; . . . ; α m Þ T 2 R nÂm contains Lagrange multiplier.
The linear equations in (3) are resulted from Karush-Kuhn-Tucker (KKT) condition for optimality.  (3) and by following the idea of LS-SVR.
Since the linear equation system (4) is not easy to be solved directly, it is reformulated into following linear equation system:, Since G is a positive definite matrix, the solutions of the linear equation system (5) can be calculated as follows: (1) Solve # and η from M# ¼ N and Mη ¼ y.
Suppose that matrixα ¼α 0 1 À Á T ;α 0 2 À Á T ; . . . ;α 0 m À Á T T and vectorb are the solutions of linear equation system (5) then the decision function of MLS-SVR can be written as follows: The proper hyper-parameters of MLS-SVR model are identified using the grid search method (Hsu, Chang, & Lin, 2016). The optimal pair of hyper-parameters (γ 0 ; γ; σ) is selected based on the minimum mean square error (MSE) criterion from all possible combinations of kernel function parameter σ as well as regularized parameters γ 0 and γ 00 . Since there are m outputs, the MSE criterion is obtained from the average of MSE over each output. This research only utilizes Gaussian kernel function because its superiority had been mentioned in Xu et al. (2013). An evolutionary algorithm (Härdle, Prastyo, & Hafner, 2014) may become useful for future work to optimize the parameters of SVR.

The Proposed control chart
Let y 1 ; y 2 ; . . . ; y m are the observable output vectors that satisfy multivariate autocorrelated data, where vector y j ¼ ðy 1j ; y 2j ; . . . ; y nj Þ T for j ¼ 1; 2; . . . ; m number of MLS-SVR outputs. Each output vector is assumed to have significant partial autocorrelation function (PACF) until lag-p 1 ; p 2 ; . . . ; p m so that the inputs of MLS-SVR model are selected as follows: X ¼ y 1;ðiÀ1Þ ; . . . ; y 1;ðiÀp 1 Þ ; . . . ; y j;ðiÀ1Þ ; . . . ; y j;ðiÀp j Þ ; . . . ; y m;ðiÀ1Þ ; . . . ; y m;ðiÀp m Þ : The decision function of MLS-SVR model with optimal parametersf ðxÞ is calculated using equation (6) so that the residuals are calculated as e j ¼ y j Àf ðx j Þ. Suppose that the residual vectors e 1 ; e 2 ; . . . ; e j . . . ; e m are resulted from MLS-SVR model with y 1 ; y 2 ; . . . ; y j ; . . . ; y m output vectors. Under in-control condition and the proper selection of MLS-SVR inputs, the hyper-parameters are tuned such that those residual vectors are independent and follow a multivariate normal distribution N m ðμ e ; P e Þ. Let e 1j ; e 2j ; . . . ; e ij ; . . . ; e nj are the residuals of MLS-SVR model for j-th output with n observations. Thus, e i ¼ ðe i1 ; e i2 ; . . . ; e ij ; . . . ; e im Þ is the m Â 1 vector of residual for iÀth observation. Vectors of residual e i are assumed to have either good mean vector μ eðgÞ , when the residuals are resulted from the in-control process, or bad mean vector μ eðbÞ , when the residuals are obtained from the out-of-control process. The vectors of residual e i are assumed to have known in-control covariance matrix P eðgÞ . The MLS-SVR-based Max-MCUSUM statistics are developed from the Max-MCUSUM statistic for independent observation (Khusna et al., 2018). The statistics of MLS-SVR-based Max-MCUSUM are obtained by forming a random variable Z i as follows: eðgÞ ðe i Àμ eðgÞ Þ ½ðμ eðbÞ À μ eðgÞ Þ T P À1 eðgÞ ðμ eðbÞ À μ eðgÞ Þ 1=2 : Since the random variable Z i is a linear combination from a normal distribution, it also follows a normal distribution. If the residual vector e i is resulted from the in-control process, which has a mean vector μ eðgÞ , then the random variable Z i follows the standard normal distribution. Otherwise, random variable Z i follows the normal distribution with mean λ. The non-centrality parameter λ is calculated as: λ ¼ ½ðμ eðbÞ À μ eðgÞ Þ T X À1 eðgÞ ðμ eðbÞ À μ eðgÞ Þ 1=2 : If the residual vector e i is resulted from a shifted process, then it may have a shifted covariance matrix P eðbÞ . The random variable W i illustrating the shift of a covariance matrix can be obtained as follows:, where ΦðzÞ ¼ PðZ zÞ; for Z,Nð0; 1Þ. Function Φ À1 ð:Þ defines the inverse of cumulative distribution function of standard normal distribution and Hðx; mÞ ¼ PðX x m j Þ; for X,χ 2 m . The random variable Z i in equation is transformed into two MCUSUM statistics to monitor mean vector formulated as: As for monitoring the covariance matrix, the random variable W i in equation (10) is transformed into two MCUSUM statistics calculated as: where C ðrÞ0 ¼ 0 and S ðrÞ0 ¼ 0 are starting point. The reference value k > 0 is firstly determined in Phase I monitoring process.
Since the multivariate quality control technique is monitoring the significance of the shift magnitude, the MCUSUM statistics in equations (11-12) can be transformed as follows: Random variable Z i and W i are independent and follow the normal distribution. Therefore, two statistics in equation (13) can be formed into the single statistic of MLS-SVR-based Max-MCUSUM control chart written as: Since M ðrÞi ! 0, the MLS-SVR-based Max-MCUSUM statistic is only compared with the UCL. The OC signal occurs when M i statistic is greater than the UCL. The statistics of MLS-SVR-based Max-MCUSUM control chart in equation (14) is developed for monitoring the mean vector and covariance matrix of auto-correlated processes. However, it is not easy to determine the distribution of M ðrÞi statistics as in equation (14). Thus, the UCL of the proposed chart is estimated using the bootstrap method, one of the well-known resampling methods to estimate the parameter of unknown distribution (Efron, 1979;Efron & Tibshirani, 1994). The bootstrap procedures to estimate the UCL of MLS-SVR-based Max-MCUSUM chart are briefly explained in Figure 1. Moreover, the specified procedures to build the proposed control chart are presented at Algorithm 1.

Algorithm 1. The Charting Procedures of the Proposed Control Chart
Step I Specify the significance level α and reference value k.
Step II Build normal profile using these following steps: 1. Let y 1 ; y 2 ; . . . ; y i ; . . . ; y n are n random vectors taken from m quality characteristics of historical in-control data. Those random vectors satisfy multivariate autocorrelated processes, with y i ¼ ðy i1 ; y i2 ; . . . ; y ij ; . . . ; y im Þ is a m Â 1 vector. Save those random vectors as the outputs of MLS-SVR model. 2. For each quality characteristic, calculate the significance of PACF then define the inputs of MLS-SVR model as in equation (7). 3. For all possible combinations of kernel function parameter σ 2 2 À15 ; 2 À13 ; . . . ; 2 3 f gand regularized parameters γ 0 2 2 À5 ; 2 À3 ; . . . ; 2 15 f g and γ 00 2 2 À10 ; 2 À8 ; . . . ; 2 10 f g , follow these steps: A. Solve the MLS-SVR parameters α and b from linear equation system (5). B. Calculate the decision function of MLS-SVR model as in equation (6). C. Compute the residuals of MLS-SVR model, the difference between output and its decision function. D. Calculate the average of MSE over each quality characteristics. E. Save the optimal parameters as well as the optimal residuals producing a minimum average of MSE criterion. 4. Calculate the mean vector and covariance matrix of the optimal residuals obtained from step II.3.E, μ eðbÞ and P eðbÞ , from sample mean vector and sample covariance matrix, respectively. 5. For i ¼ 1; 2; . . . ; n samples, follow these steps: A. Calculate statistics Z i using equation (8). B. Calculate statistics W i using equation (10). C. Calculate statistics C þ ðrÞi and C À ðrÞi using equation (11). D. Calculate statistics S þ ðrÞi and S À ðrÞi using equation (12). E. Calculate statistics C ðrÞi and S ðrÞi using equation (13). F. Calculate statistics M ðrÞi using equation (14) 6. Estimate the UCL with the following procedures: A. For , ¼ 1; 2; . . . ; N 0 replications, follow these steps: a. Resample statistics M ðrÞi ; i ¼ 1; 2; . . . ; n; for B times, where B is a large number. b. Sort B bootstrap samples obtained in step II.6.A.a from small to large value. c. Compute ð100ð1 À αÞÞÀth percentile of bootstrap: M , Calculate the UCL, the average of ð100ð1 À αÞÞÀth percentile over N 0 replica- Step III Monitor new observations using these following steps: 1. Save the random vectors y nþ1 ; y nþ2 ; . . . as the outputs of MLS-SVR model. 2. Define the inputs of MLS-SVR model using the significance of PACF obtained in step II.2. 3. Train the inputs and outputs of MLS-SVR model using optimal parameters obtained in step II.3.E. 4. Repeat steps II.3.B and II.3.C with the new observations. 5. For i ¼ n þ 1; n þ 2; . . . samples, repeat steps II.5. 6. Plot statistics M ðrÞi along with its UCL obtained from step II.6.B against the number of samples. 7. Compare statistics M ðrÞi with the UCL. If any out-of-control (OC) signal is found, identify the source of shifts as follows:

Simulation studies
This section is aimed to present two kinds of simulation studies. First, the simulation studies to obtain optimal inputs and hyper-parameters of MLS-SVR model if the historical in-control data are assumed to follow a vector autoregressive (VAR), vector moving average (VMA), and VARMA model. Second, the simulation studies to calculate the UCL of the proposed chart if the optimal residual is assumed to follow a standard multivariate normal distribution.

The inputs and hyper-parameters selection of MLS-SVR model
The simulation studies are carried out in order to obtain the proper inputs and the optimal hyper-parameters of MLS-SVR model. These two conditions are addressed to obtain residuals which satisfy multivariate normal distribution and white noise assumption. These simulation studies use linear multivariate time series data such as VAR (1), VMA (1), and VARMA (1,1) which are generated as follows: VARð1Þ : This research sets up for μ 1 ¼ 5 and μ 2 ¼ 10. Table 1 summarizes the proper inputs of MLS-SVR model for some multivariate linear time series data selected using equation (7). Since each series of VAR (1) has PACF which cuts off after lag 1, the inputs of MLS-SVR model for the different parameter Φ 1 are selected as y 1ðiÀ1Þ ; y 2ðiÀ1Þ . These inputs produce MLS-SVR residuals that follow a multivariate normal distribution and satisfy white noise assumption. However, the inputs of MLS-SVR model for VMA and VARMA data are not easy to be selected due to the dies down pattern of their PACF plots (see Appendix 1). The different number of observation and the different matrix of parameter Θ q yield different significant lag of PACF such that the input selection is not straightforward. The VMA (1) data with different parameter Θ 1 as well as the VARMA (1,1) data with different parameter Φ 1 and Θ 1 need different number of MLS-SVR inputs to obtain the normality and white noise condition for their residuals. The smaller number of observation needs the fewer MLS-SVR input because of the wider significance limit of PACF plot as exhibited in Appendix 1. Furthermore, the VMA (1) as well as VARMA (1,1) data with parameter Θ 1 ¼ 0:80 0:05 0:05 0:90 ! and n = 1000 observations are displayed with three kinds of input selection in order to show that by setting the fewer input, they can produce the residuals which violate the assumptions. The input selections presented in Table 1 with n = 1000 observations use the optimal combination of hyper-parameters γ 0 , γ, and σ as summarized in Table 2. The optimal pair of hyper-parameters is obtained using a grid search method resulting minimum   ! mean of MSE criterion over 10-fold cross-validation. The residuals of MLS-SVR model for all datasets are independent and follow a multivariate normal distribution with mean vector μ e and covariance matrix P e . The mean of MSE for all datasets is close to 1because the residuals of linear multivariate time series data are generated from identity covariance matrix.

The UCL of proposed control chart
The simulation studies are conducted in order to estimate the UCL of the proposed chart using the bootstrap method. Under the in-control condition, the residuals of MLS-SVR model resulted from the proper inputs and the optimal hyper-parameters are independent and follow a multivariate normal distribution (see Table 1-2). Therefore, it is possible to conduct the simulation studies with the optimal residuals follow a standard multivariate normal distribution. In this simulation study, n ¼ 1; 000 samples are generated from standard multivariate normal distribution then statistics M ðrÞi are calculated using equation (14). Table 3 provides the UCL of proposed chart calculated using steps II.6 of Algorithm 1 for significance level α = 0.0027, B ¼ 10; 000 bootstrap samples and N 0 = 500 replications. The larger reference value k, the tighter control limit of MLS-SVR based Max-MCUSUM chart.

The performance evaluation
The performance of the proposed control chart is evaluated using ARL, the average number of sample until finding the first out-of-control observation. For various combinations of multivariate linear time series data and reference value k, the simulation studies based on Algorithm 2 are carried out to obtain the ARL of the proposed control chart. The in-control ARL (ARL 0 ) is calculated for VAR (1), VMA (1), or VARMA (1,1) data in which residuals follow a multivariate normal distribution with the mean vector μ aðgÞ ¼ 0 and covariance matrix P aðgÞ ¼ I. Otherwise, the out-of-control ARL (ARL 1 ) is calculated either by adding the mean vector shift from μ aðgÞ to μ aðbÞ ¼ μ aðgÞ þ δ or/ and by adding the covariance matrix shift from P aðgÞ to P aðbÞ ¼ P aðgÞ þ v1 mÂm . Since the magnitude of mean vector shift δ is assumed to be the same for each variable, the residuals of MLS-SVR model are in-control if δ ¼ 0 and v ¼ 0.

Algorithm 2. The Simulation of ARL for MLS-SVR-based Max-MCUSUM Control Chart
Step 1.Specify the significance level α, reference value k, number of observations n, and the number of quality characteristics m.
Step 2.Set the UCL of MLS-SVR based Max-MCUSUM chart.
Step 3.Define the parameter of multivariate linear time series model, Φ p or/and Θ q , in which residuals follow a multivariate normal distribution N m ðμ aðgÞ ¼ 0; P aðgÞ ¼ IÞ.
Step 4. Generate n samples of multivariate linear time series data using parameters setting from step 3.
Step 5. Train time series data generated in step 4 using MLS-SVR algorithm according to equation (6), then save the optimal parametersαandb along with the optimal hyper-parameters γ 0 ; γ; and σ.
Step 6. Let μ aðbÞ ¼ μ aðgÞ and P aðbÞ ¼ P aðgÞ then follow these steps: A. For N 0 replications, do these steps: a. Generate n samples of multivariate linear time series data using the parameter Φ p or/and Θ q defined in step 3, which residuals follow N m μ aðbÞ ; P aðbÞ . b. Train time series data generated in step 6.A.a using MLS-SVR parameters α andb along with the hyper-parameters γ 0 ; γ; and σ yielded from step 5. c. Calculate the statistics M ðrÞi , i ¼ 1; 2; . . . ; n; using equation (14). d. Calculate the Run Length (RL), the number of sample until finding the first M i statistic that greater than the UCL.
B. Calculate the ARL 0 , the average of RL over N 0 replications.
In this research, the ARLs of the proposed control chart are calculated for α= 0.0027, n = 1,000 samples, and N 0 = 1,000 replications. The multivariate linear time series data are generated from VAR (1), VMA (1), as well as VARMA (1,1) with the parameters shown in Table 4. It is important to be highlighted that VAR (1) model with scenario 1 has higher main diagonal of parameter Φ 1 than that with scenario 2. For m = 2 variables, the inputs of MLS-SVR model are selected as in Table 1 and the hyperparameters are tuned as in Table 2. Furthermore, the UCLs of the proposed control chart are obtained from Table 3. The value of Δδ for each variable and the value of Δv are equal to 0.25.
The ARLs of the proposed control chart for VAR (1) and VMA (1) data are presented in Tables 5-6 respectively. If the covariance matrix remains in-control and each component of the mean vector is shifted from 0 to 0.25, the ARL 1 of the proposed control chart for VAR (1) data diminishes from 373 to about 136 while the ARL 1 for VMA (1) data reduces from 374 to about 105. If each component of the covariance matrix has shifted from 0 to 0.25 and the mean vector is in-control, the ARL 1 of the proposed control chart for VAR (1) data reduces from 373 to about 208. However, the ARL 1 for VMA (1) data diminishes from 374 to about 220 to detect the same changes. Hence, the proposed chart has better performance to detect the mean vector shift of VMA (1) data than that of VAR (1) data. On the contrary, the proposed chart is more sensitive to detect the covariance matrix shift of VAR (1) data than that of VMA (1) data. Table 7 presents the ARLs of the proposed control chart using for VARMA (1,1) data. Generally, both ARL 0 and ARL 1 of the proposed chart for VARMA (1,1) data have a similar pattern to the ARLs for both VAR (1) and VMA (1) data. The stable ARL 0 at about 370 indicates that under the in-control process, the proposed control chart will detect the first out-of-control signal at about 370 samples. Moreover, the ARL 1 is   decreasing as the increasing of either mean vector shift, covariance matrix shift, or the combination of both. This evidence points out that the larger the mean vector and/or the covariance matrix shifts, the faster the ability of the proposed control chart to detect those actual shifts. The ARLs of the proposed chart for several combinations of multivariate linear time series data, reference values k, and the number of quality characteristics m are presented in Appendix 2. Figure 2 displays the ARL's plots of MLS-SVR-based Max-MCUSUM control chart for VAR (1) data. The plotted ARLs for mean vector shift are presented with the covariance matrix remains in-control and vice versa. Based on the ARLs comparison in Figure 2, some characteristics of the proposed control chart can be highlighted as follows: (a) For a certain parameter of VAR (1) data and a certain number of quality characteristics, the proposed control chart with the smaller the reference value is more sensitive to detect mean vector shift. However, the proposed control chart is more sensitive to detect the covariance matrix shift if the reference value is higher.  (b) For equal reference value and equal parameter of VAR (1) data, the higher the number of quality characteristics, the more sensitive of the proposed chart to detect the shift of mean vector. (c) For the same value of m and the same value of k, the proposed chart for VAR (1) data with scenario 2 is slightly more sensitive than that with scenario 1.
The plots of the ARLs for VMA (1) and VARMA (1,1) data are exhibited in Figure 3. Generally, the effects of k and m to the ARLs of both VMA (1) and VARMA (1,1) data are not significantly different with those to the ARLs of VAR (1) data. If each component of the covariance matrix has shifted to 0.25, then the ARLs for both VMA (1) and VARMA (1,1) data tend to be diminished at about 200. If each component of the mean vector has increased from 0 to 0.25, then the ARLs for both VMA (1) and VARMA (1,1) data are decreasing to about 100. Therefore, the proposed control chart is more sensitive to detect mean vector shift than to detect a covariance matrix shift. The existing multivariate control charts for monitoring auto-correlated data are designed to monitor the mean and variability of a process separately. The researchers cannot find any single control chart for monitoring multivariate auto-correlated processes. Therefore, it is impossible to compare the proposed chart with any existing chart. In this research, the proposed control chart is compared with Max-MCUSUM based on the residuals of the VAR model. The comparison is carried out for bivariate time series data following VAR (1) model scenario 1 with parameter Φ 1 displayed in Table 4. The ARLs of VAR-based Max-MCUSUM chart are calculated by replacing the MLS-SVR algorithm with the Box-Jenkins procedures (Box & Jenkins, 1976;Tiao & Tsay, 1989).   Table 8. If the mean vector remains in-control (δ ¼ 0) and each component of the covariance matrix has shifted to certain value of v then VAR based Max-MCUSUM chart has smaller ARL 1 than the proposed chart. However, the proposed control chart has smaller ARL 1 than VAR based Max-MCUSUM chart if each component of mean vector has shifted to  about 0.25 until 1.50 and each component of covariance matrix has shifted to some values of v. Hence, the proposed control chart tends to have better performance than VAR based Max-MCUSUM chart if the mean vector is shifted equal to or smaller than 1.50.

Applications
This section provides some illustrative examples of the proposed control chart to monitor both simulation and real data. The simulation data generated using VAR (1) model illustrate the shifts in the mean vector, covariance matrix, or combination of both shifts. In the real case, the hourly water quality data are utilized to give an example about the actual implementation of the proposed chart to monitor multivariate autocorrelated processes.

The simulation data
Four simulation datasets are generated from the VAR (1) model using equation . The parameters of VAR (1) residual for all datasets are summarized in Table 9. Each dataset contains 70 samples of training data and 30 samples of testing data. The time series plots of dataset 2 (see Figure 5(b)) and dataset 3 (see Figure 5(c)) for the last 30 observations show the shifts in the mean vector and covariance matrix, respectively. Furthermore, the time series plots of dataset 4 displayed in Figure 5(d) exhibit the changes in both mean vector and covariance matrix in the last 30 samples.
Since the simulation data are generated from the VAR (1) model with m = 3 variables, y 1;ðiÀ1Þ ; y 2;ðiÀ1Þ ; y 3;ðiÀ1Þ are selected as the inputs of MLS-SVR model. For each dataset, the first 70 observations are trained using MLS-SVR algorithm in order to obtain the optimal parametersαandb along with the optimal hyper-parameters γ 0 ; γ; and σ presented in Table  9. The last 30 observations of each dataset are fitted using those optimal parameters and hyper-parameters. Figure 6 presents monitoring results of dataset 1 and dataset 2 using MLS-SVR-based Max-MCUSUM control chart for various reference values k. There is no out-ofcontrol signal issued by the proposed control chart (see Figure 6(a-b)) because dataset 1 is generated without adding a shift either in mean vector or covariance matrix. Since the vertical axes of Figure 6(a-b) have the same scale, it is easy to understand that by setting the smaller reference value k yields the higher statistic and the higher UCL of the proposed control chart.
The proposed control chart for dataset 2 with k = 0.5 (see Figure 6(c)) finds the first OC signal at 74th sample. The first OC signal is marked with m þ because based on statistics C ðrÞ74 and S ðrÞ74 , only statistic C ðrÞ74 that is larger than the UCL. The proposed control chart with k = 1.25 (see Figure 6(d)) also detects the first OC signal marked with m þ at 73rd sample. This points out that the OC signals are caused by the mean vector shift. Moreover, the MLS-SVRbased Max-MCUSUM control charts for dataset 3 using k = 0.5 (see Figure 7(a) and k = 0.75 (see Figure 7(b)) detect the first OC signal at 73rd sample marked with v þ . This result indicates that the first OC signal at 73rd sample is only caused by a shift in the covariance matrix. The proposed control charts for dataset 4 presented in Figure 7(c-d) also detect the first OC signal at 73rd sample. Those OC signals are marked with m þ v þ indicating the changes in both mean vector and covariance matrix. The summary of the monitoring results is also presented in Table 9.

The real water quality data
This section provides the applicability of the proposed control chart to monitor the real water quality data in Surabaya. Final water turbidity (Y 1 ), water turbidity after affixation process (Y 2 ), and chlorine residual (Y 3 ) are three important quality characteristics in the water manufacturing process. Those quality characteristics satisfy the multivariate autocorrelated processes because they are monitored hourly. This study uses the hourly water quality data started from 17 th of February 2016 until 3 rd of March 2016. The first 10 days of hourly water quality data are employed as training data (monitoring Phase I), while the rests are utilized as testing data (monitoring Phase II). There are 240 observations of Phase I and 144 observations of Phase II. The water turbidity data are measured in Nephelometric Turbidity Units (NTU) while chlorine residual data are measured in part per million (ppm) unit. The time series plots of water quality data displayed in Figure 8 show a decreasing pattern started from 29 th of February 2016, especially in the final water turbidity data.
In order to convince the applicability of the proposed chart, the hourly water quality data are monitored using three types of control chart, which are MCUSUM (Crosier, 1988) chart, MLS-SVR-based MCUSUM (Crosier, 1988) chart, and the proposed chart. The MCUSUM control chart displayed in Figure 9(a) shows that most of the observations are detected as the out-of-control signals. The MCUSUM control chart presents such a bad process in the water manufacturing process and it will confuse the engineer. However, this condition actually explains that the conventional control chart is not appropriate to monitor the autocorrelated processes such as the hourly water quality data. As an alternative solution, the residual-based control charts are utilized. As explained in the charting procedures in Algorithm 1, the inputs of MLS-SVR model are selected using historical in-control data. In this case, the PACF is calculated using training data. The PACF of Y 1 is significant until lag three whereas the PACFs of Y 2 and Y 3 are significant until lag two. Therefore, the inputs of MLS-SVR model are selected as y 1ðiÀ1Þ ; y 1ðiÀ2Þ ; y 1ðiÀ3Þ ; y 2ðiÀ1Þ ; y 2ðiÀ2Þ ; y 3ðiÀ1Þ ; y 3ðiÀ2Þ . The training data are trained using MLS-SVR algorithm such that the optimal combination of hyperparameters is resulted as γ 0 ¼ 2 À5 , γ ¼ 2 À4 , and σ ¼ 2 3 . The residuals of MLS-SVR model are then monitored using both MCUSUM chart and Max-MCUSUM chart. The residual-based MCUSUM chart and residual-based Max-MCUSUM chart are displayed in Figure 9(b-c), respectively. Both residual-based control charts detect the out-ofcontrol signals started from 29 th of February 2016 indicating an improvement in the water manufacturing process should be carried out. This result can ensure the engineer about the time to take the corrective actions.

Conclusions and future research
This study proposes an MLS-SVR-based Max-MCUSUM control chart for simultaneously monitoring the mean and variability of multivariate autocorrelated processes. The optimal parameters of MLS-SVR model are calculated using historical in-control observations then the UCL of the proposed chart is estimated using the bootstrap approach. Simulation studies show that the larger the reference value, the smaller the UCL of the proposed chart. The proposed control chart is more sensitive to detect the mean vector shift if the reference value is smaller. However, the larger the reference value, the more sensitive of the proposed chart to detect the covariance matrix shift. Furthermore, the proposed control chart has a better performance than VAR-based Max-MCUSUM chart for small mean vector shift. The proposed control chart is sensitive to monitor the shifted processes in both simulation and real water quality data. This study can be extended by considering not only negative shifts but also unequal shifts in each component of mean vector and covariance matrix.
The ARLs of the proposed chart