Max-EWMA chart for time and magnitude monitoring using generalized exponential distribution

Abstract Joint monitoring of time and magnitude is an important issue in many industrial and non-industrial fields and several proposals exist in the literature. The aim of this article is to propose a new maximum exponentially weighted moving average (Max-EWMA) chart assuming generalized exponential distribution for time as well as magnitude. As the test-statistic of Max-EWMA chart consists of two EWMA statistics, it is easy to identify the out-of-control signal raised by the chart. The performance of the chart is evaluated using different run length characteristics including, average run length (ARL), standard deviation of run length (SDRL), and quantiles of the run length distribution. Besides extensive simulation studies, a real data set is also used to show the practicality of the proposed chart. The results show that the proposed chart is efficient in detecting different size of small to moderate shifts.


Introduction
Statistical process control (SPC) tools, especially control charts are not limited only to manufacturing processes but equally applicable in many non-manufacturing fields.Monitoring the event time and its associated magnitude is important due to its broad spectrum applications, and therefore, several types of control charts have been developed to study them collectively and individually.The interval time between events (T) can be monitored by using a control chart known as the time-between-events (TBE) control chart or the T control chart.Recently, new control charts have been proposed to monitor time and its associated magnitude (Ali and Pievatolo 2018;Ali 2020).For example, T&X chart is proposed by Wu, Jiao, and He (2009) that is used to monitor the frequency and magnitude of an event (or the ratio r).
In this article a new Max-EWMA control chart for simultaneous monitoring of magnitude and time stamp is introduced assuming both magnitude, X, and time stamp, T, follow the generalized exponential (GE) distribution.The probability density functions for magnitude and time, respectively, are given below.

À Á cÀ1 e
Àt u t , t > 0: (2) where u x , u t , d, and c denote the scale and shape parameters of the magnitude and time distributions, respectively.In the literature, a max-EWMA control chart for simultaneous monitoring of magnitude and frequency of an event is proposed by Sanusi, Teh, and Khoo (2020).The chart found efficient for the detection of downward shifts in magnitude and upward shifts in time assuming gamma distribution for magnitude and exponential distribution for time-betweenevents.The EWMA chart is defined as Z i ¼ kX i þ ð1 À kÞZ iÀ1 where X i is the observed value at the time point i.The time or sample number is plotted on the horizontal axis against the EWMA statistic Z i plotted on the vertical axis.The Max-EWMA control chart was first proposed by Chen and Cheng (1998) and its extension was given by Xie (1999).In fact, two EWMA charts were combined to construct a Max-EWMA chart in such a way that it can monitor two different statistics simultaneously.Later, for jointly monitoring of mean and variance, an EWMA control chart based upon a non-central chi-squared statistic was developed by Costa and Rahim (2004).
In the Max-EWMA chart, the maximum absolute value computed by two EWMA charts, one for TBE and the other for magnitude is plotted against the samples number or time index.This chart can easily identify the direction and source of an out-of-control signal when it is detected.
To monitor the time stamps for both phase I and II, a control chart based on gamma distribution has been proposed by Yang et al. (2016).For monitoring the TBE, a one-sided generally weighted moving average (GWMA) control chart on the basis of gamma distribution has been proposed by Chakraborty, Human, and Balakrishnan (2017).Ali, Pievatolo, and G€ ob (2016) gave an overview related to high-quality processes by using TBE charts.Some recent contributions in this direction can be seen in Zhang and Chen (2004), Khoo and Xie (2009), Pascual (2010), Akhundjanov and Pascual (2015), Shafae et al. (2015) and Wang, Bizuneh, and Abebe (2017), and references cited therein.
To monitor time and magnitude jointly, several authors introduced different charts including Wu, Jiao, and He (2009), Liu et al. (2009), and Qu et al. (2014).For simultaneous monitoring of magnitude and time stamp, a distribution free upper-sided EWMA type control chart has been developed by Wu, Castagliola, and Celano (2020).Similarly, a Max-EWMA chart was proposed by Sanusi, Teh, and Khoo (2020) to monitor simultaneously time and magnitude.For autocorrelated processes, a Max-CUSUM control chart was introduced by Petcharat (2018).Sheu and Lin (2003) developed a Max-GWMA control chart and concluded that this chart is effective in the detection of shifts in the mean as well as the variance.Control charts using first passage time distribution have also been proposed by Ali and Pievatolo (2018) and Ali (2020) to jointly monitor time and magnitude using renewal reward process assumption.
The three parameter exponentiated Weibull distribution is introduced by Mudholkar and Srivastava (1993) and its special case is known as the generalized exponential distribution (GED), which is a generalization of the ordinary exponential distribution.The density function of the GED is right skewed unimodal and in many cases, it is a better substitute of Weibull distribution and gamma distributions.
Several properties of GED have been discussed by different researchers (Gupta and Kundu 1999, 2001a, 2001b).In quality control, control charts have been developed by using the GED.For example, a control chart by using the GED on the basis of renewal process is introduced by Ali (2017).This chart shows more flexibility than the previous exponential TBE charts.Likewise, Chiang et al. (2017) used GED to introduce five control charts, including the Shewhart-type chart and four parametric bootstrap charts for the GE percentiles.Considering the practical importance and increasing popularity of the GED, we use this distribution to develop a Max-EWMA control chart for simultaneous monitoring of magnitude and TBE data.Besides a comprehensive simulation study, real data examples are especially focused to show the applications of the proposed chart.
The rest of the article is organized as follows.Section 2 presents the Max-EWMA chart assuming GED.Section 3 gives a description of in-control ARL of the proposed chart while Sec. 4 presents a detailed out-of-control performance assessment of the Max EWMA chart.Two real data examples are discussed in Sec. 5 whereas Sec.6 concludes the study.

Max-EWMA control chart for GED
Let an event magnitude, X, and TBE, T, both follow the GED, i.e., X $ GEðd, u x Þ, where d is the shape parameter and u x is the scale parameter and T $ GEðc, u t Þ where c denotes the shape parameter and u t denotes the scale parameter of TBE distribution.Assuming independence between time stamp and magnitude, the process in-control parameters are denoted by (d 0 , u x 0 ) and (c 0 , u t 0 ), whereas the out-of-control (OOC) parameters are defined as where the change in magnitude shape and scale parameters is denoted by d d and d u x while the change in time stamp shape and scale parameters is denoted by d c and d u t : It is to be noted that when d d , d u x , d c and d u t are equal to one, then the process is said to be in-control.To calculate the Max-EWMA chart, we first calculate the following two independent standardized statistics.
where EðXÞ ¼ u x ½wðd þ 1Þ À wð1Þ and VarðXÞ ¼ u 2 x ½w 0 ð1Þ À w 0 ðd þ 1Þ, represent the mean and variance of the generalized exponential distribution for the magnitude, psið:Þ and its derivatives are the digamma and polygamma functions.In particular, the digamma is defined by wðxÞ ¼ d ln ðCðxÞÞ=dx: We refer to Gupta and Kundu (1999) for the details to derive these expressions.Then, the two EWMA statistics using Z i and W j are computed as follows.
where the starting values for the EWMA statistics are taken as R 0 ¼ EðZ j Þ ¼ 0, H 0 ¼ EðW j Þ ¼ 0 and x is the smoothing parameter lies within (0,1].Next, the Max-EWMA statistic for monitoring magnitude and TBE is computed by combining R j and H j as follows.
N j ¼ maxfjR j j, jH j jg: It is clear from Eq. ( 8) that if the time or magnitude shifts from in-control state to an out-of-control state, the statistic N j become greater than the statistic N j computed from in-control data.Hence, only the upper control limit is required for the Max-EWMA control chart because N j is a non-negative real number.Hence, the upper control limit is computed by where E(N) and Var(N) denote the mean and variance of the in-control N j , respectively.It is worth mentioning that explicit expressions for mean and variance are difficult to derive and these can be calculated numerically by generating in-control data.Furthermore, the constant L is a multiplier that manages the width of the control limit.
The procedure to construct the Max-EWMA control chart for simultaneously monitoring magnitude and time by using generalized exponential distribution is summarized as follows.
(1) Generate in-control time and magnitude data and construct the UCL of the Max-EWMA chart using Eq. ( 9) by calculating the values of E(N) and Var(N) from the in-control process.The ðx, LÞ are fixed to achieve the desired in-control ARL (ARL 0 ).
(2) Calculate Z j , W j , R j , H j and N j for each individual sample by setting the initial values as R 0 ¼ H 0 ¼ 0. (3) Plot N j against the sample number j on the control chart.If N j < UCL, declare the process in-control, otherwise out-of-control.By using the identification symbols given in Table 1, label the plotted points.If N j !UCL, then examine both jR j j and jH j j against the UCL.Label the plotted point with "X" if jR j j !UCL, otherwise label "T" if jH j j !UCL: For a simultaneous shift in both magnitude and time interval between events, label the plotted point with "XT" if any of the two statistics jR j j and jH j j are greater than UCL.(4) Investigating the underlying cause(s) associated to each out-of-control signal and then resume the process again from Step 3.

In-control ARL of the Max-EWMA chart
To develop a Max-EWMA control chart for simultaneously monitoring time and magnitude using GED, the design strategy is to first locate the optimal L value by satisfying the desired incontrol average run length (ARL) for a given smoothing parameter.The ARL measures on average the out-of-control signal given by the chart whereas the stability of a control chart run length is measured by the SDRL.As the computation of ARL expression is not directly possible, the ARL value is computed using 5000 simulation runs.To obtain L with the desired in-control ARL 0 (say 370) assuming u x ¼ u t ¼ 0.005, 0.01, 0.50, 1.00, 1.50 and x ¼ 0.05, 0.10, and 0.12, Table 2 shows the in-control ARL.It is clear from the table that a large L value is obtained for large x.Similarly, small L values are obtained with small shape parameters.

Out-of-control performance of the Max-EWMA chart
To assess the performance of the Max-EWMA chart, the run length characteristics like the ARL, SDRL, and quantiles of the run length distribution are used.In particular, Monte Carlo simulations are used to determine run length characteristics of the proposed Max-EWMA chart.The values of the shape parameter are fixed to have increasing, decreasing, and constant hazard functions of the distribution (Ali 2017).To this end, different shifts like d 2 f1, 0:85, 0:80, 0:70, 0:50g with the ith percentiles P i , i ¼ 10, 25, 50, 75, 90 and in-control ARL 370 are used.In particular, fd 0 ¼ 0.5 (decreasing hazard function), u x ¼ 0.005, c 0 ¼ 0.5 (decreasing hazard function), u t ¼ 0:005g, fd 0 ¼ 1.0 (constant hazard function), u x ¼ 0.005, c 0 ¼ 1.0 u t ¼ 0:005g (constant hazard function), and fd 0 ¼ 1.5 (increasing hazard function), u x ¼ 0.005, c 0 ¼ 1.5 (increasing hazard function), u t ¼ 0:005g cases with x ¼ 0:05, 0:10, 0:12 have been considered.The small values of x are taken to detect small shifts (Montgomery 2007).The run length profile of the Max-EWMA chart computed by assuming x ¼ 0:05 are reported in Tables 3-13 and A.1-A.4 whereas the results for x ¼ 0:10, 0:12 are reported in the Supplementary Materials.Table 3 lists the ARL and other run length measures assuming x ¼ 0:05 for different shifts in the scale parameters of the time and magnitude distributions assuming shape parameters equal to 0.5.It is observed from the table that the ARL(SDRL) increased by 71.9% (76.7%) when there is 15% shift in the scale parameter of the magnitude and it is increased by 73.4% (76.3%) for 15% shift in the scale parameter of the time distribution.However, when there are simultaneous shifts of ðd u x , d u t Þ ¼ ð0:85, 0:85Þ and (0.70,0.70), compared to the ARL 0 , ARL 1 increased by 520.5% and 374.2%, respectively.The percentiles show that the run length distribution is rightly skewed where the ARL is always greater than P 50 but less than P 75 .For example, when ðd u x , d u t Þ ¼ ð0:85, 0:85Þ, the median is 1935.0, the ARL is 2303.0 and P 75 ¼ 3810.5.For very large shifts, it is noticed that the ARL 1 < ARL 0 : Next, Table 4 presents the ARL values assuming the shape parameter equal to one.It is noticed from the table that the ARL increased by 23.3% (23.5%) when there is 15% shift in the scale parameter of the magnitude and increased by 23.6% (22.9%) when there is 15% shift in the scale parameter of the time distribution.However, if there are simultaneous shifts in ðd u x , d u t Þ ¼ ð0:85, 0:85Þ, the ARL increases by 70.2% whereas ARL 1 reduces by 89.3% in the case of ðd u x , d u t Þ ¼ ð0:50, 0:70Þ: Furthermore, the ARL is biased when the shape parameter is less than or equal to one.Similarly, Table 5 lists the ARL values with 1.5 as the shape parameter.From the table it is observed that the ARL reduced by 11.9% (14.1%) when there is a 15% shift in the scale parameter of the   Next, to study shifts in the scale parameter of the magnitude distribution, Table 6 lists the ARL assuming x ¼ 0:05 for different shape parameters.For the shape parameter equals to 0.5 it is observed that the ARL increased by 72.6% (72.9%) when there is a 15% shift in the scale parameter of the magnitude distribution.Furthermore, the percentiles show that the run length distribution is right skewed, i.e., the ARL is always greater than P 50 but less than P 75 .For the shape parameter equals to one, the ARL increased by 27.9% (26.0%) when there is a 15% shift in the scale parameter of the magnitude.Contrary to shape parameters less than or equal to one, it is observed that the ARL reduced by 11.8% (14.8%) when there is a 15% shift in the scale parameter of the magnitude distribution with shape parameter equals to 1.5.Thus, the chart outperforms when the shape parameter is greater than one, i.e., increasing hazard function.Similar trends are noticed for shifts in the scale parameter of the TBE (Table 7).
Table 8 lists the ARL values assuming x ¼ 0:05 with in-control shape parameter equals to 0.5 and one can see that the ARL increased by 49.0% (53.8%) when there is a 20% shift in the shape parameter of the magnitude and increased by 82.0% (83.9%) when there is a 20% shift in the scale parameter of the time distribution.Also, if there are simultaneous shifts ðd d , d u t Þ ¼ ð0:80, 0:80Þ and (0.70,0.80),ARL 1 increased by 356.2% and 207.3%, respectively.Thus, in this case, the ARL is biased.Contrary to 0.5 as the shape parameter, Table 9 tabulates the ARL values with shape parameter equals to one.It is noticed from the table that the ARL increased by 2.2% (2.5%) when there is a 20% shift in the shape parameter of the magnitude and reduced by 8.3% (12.1%) if there is a 20% shift in the scale parameter of the time distribution.Also, if there are   10 show that the run length distribution is right skewed, where the ARL is always greater than P 50 but less than P 75 .For example, when ðd d , d u t Þ ¼ ð0:80, 0:80Þ, the median is 135.0, the ARL is 180.5 and P 75 ¼ 239.0.Compared to the shape parameter equals to 0.5 and 1.0, the results of the shape parameter equals to 1.5 are always unbiased.Tables 11-13 present the ARL values assuming x ¼ 0:05 and in-control shape parameter equals to 0.5, 1.0, and 1.5, respectively.For the shape parameter 0.5, it is observed that the ARL reduced by 61.3% (70.6%) when there is 30% shift in the shape and scale parameters of the magnitude distribution.Also, if there are simultaneous shifts of ðd d , d u x Þ ¼ ð0:50, 0:50Þ and (0.80,0.80),ARL 1 reduced by 88.9% and increased by 11.8%, respectively.with in-control shape parameter equals to 1.0.Here, the ARL reduced by 87.3% (93.4%) if there is a 30% shift in the shape and scale parameters of the magnitude distribution.Also, if there are simultaneous shifts in ðd d , d u x Þ ¼ ð0:50, 0:50Þ and (0.80,0.80),ARL 1 compared to the ARL 0 reduced by 94.2% and 69.3%, respectively.The percentiles in Table 12 show that the run length distribution is right skewed where the ARL is always greater than P 50 but less than P 75 .For example, when ðd d , d u x Þ ¼ ð0:70, 0:70Þ, the median is 41.0, the ARL is 46.8 and P 75 ¼ 57.0.Table 13 tabulates the ARL values for the shape parameter equals to 1.5.Contrary to previous values of the shape parameters, it is noticed from the table that the ARL reduced by 90.7% (95.8%) when there is a 30% shift in the shape and the scale parameters of the magnitude distribution.Also, if there are simultaneous shifts in ðd d , d u x Þ ¼ ð0:50, 0:50Þ and (0.80,0.80),ARL 1 reduced by 95.5% and 90.7%, respectively.Tables A. 1-A.3 show the ARL values assuming x ¼ 0:05 with in-control shape parameters equal to 0.5, 1.0, and 1.5, respectively.For 0.5 as the shape parameter, it is noticed that the ARL increased by 61.5% (58.1%) when there is a 30% shift in the scale parameter of the magnitude and increased by 31.3% (26.2%) if there is a 30% shift in the shape parameter of the time distribution.Also, if there are simultaneous shifts in ðd u x , d c Þ ¼ ð0:50, 0:50Þ and (0.80,0.85),ARL 1 reduced by 70.5% and increased by 293.7%, respectively.Thus, large shifts are detected more quickly than small shifts.The percentiles in Table A.1 show that the run length distribution is right skewed, that is, the ARL is always greater than P 50 but less than P 75 .For example, when ðd u x , d c Þ ¼ ð0:50, 0:50Þ, the median is 88.0, the ARL is 108.1 and P 75 ¼ 136.0.For the shape parameter equals to 1.0, Table A.2 shows the ARL reduced by 60.9% (68.6%) when there is a 30% shift in the scale parameter of the magnitude and reduced by 44.3% (49.3%) if there is a 30% shift in the shape parameter of the time distribution.Also, if there are simultaneous shifts in ðd u x , d c Þ ¼ ð0:50, 0:50Þ and (0.80,0.85),ARL 1 reduced by 90.3% and increased by 11.0%,  respectively.Next, Table A.3 lists the ARL values with in-control shape parameter equals to 1.5.It is observed that the ARL reduced by 92.7% (97.6%) when there is a 50% shift in the scale parameter of the magnitude and reduced by 87.9% (92.7%) when there is a 50% shift in the shape parameter of the time distribution.However, when there are simultaneous shifts in ðd u x , d c Þ ¼ ð0:50, 0:50Þ and (0.80,0.85),ARL 1 reduced by 93.2% and 46.3%, respectively.Thus, the Max-EWMA chart has unbiased ARL when the in-control shape parameters of time and magnitude distributions are greater than one.Also, the run length distribution is noticed right skewed for the assumed shape parameters.Table A.4 shows the ARL values for shape parameters equal to 0.5, 1.0, and 1.5, respectively.For 0.5 as the shape parameter value, it is observed that when the shift in the shape parameter of the magnitude increased by 40% from 1.00, the ARL reduced by 69.4% (72.4%).Similarly, if the shift size decreased by 20% from 1.00, the ARL increased by 53.2% (52.5%).For the shape parameter equals to one, it is observed from the table that when the shift in the shape parameter of the magnitude distribution increased by 20% the ARL reduced 51.1% (52.6%).Similarly, when the shift size decreased by 20% the ARL increased by 1.7% (0.3%).For the shape parameter equals to 1.5, it is noticed that when a shift in the shape parameter of the magnitude increased by 20% the ARL reduced by 52.2% (54.6%).Also, if the shift decreased by 40% from 1.00 the ARL reduced by 22.7% (26.3%).

A real life application
To show the application of the Max-EWMA chart, in this section a real data set is discussed.et al. (2021) reported a data set of a company for monitoring the bottleneck machine from 08/01/12 to 27/12/18 (format DD/MM/YY) listed in Table S1 of the Supplementary Materials.It includes the times between two consecutive breakdowns (T i in days) and also the consequent cost (X i in euros), that include the repair and restart cost of spare parts, the cost of work force and production disruption.As the parameters of the process are unknown, the 44 records of amplitudes, X, and time between two consecutive breakdowns, T, has been taken to fit generalized exponential, exponential and gamma distributions.As Rahali et al. (2021) stated that there is a slight positive correlation between X and T, we assume that X and T are independent to implement the proposed chart.The resulting parameters are reported in Table 14 and assuming x ¼ 0.05 and L ¼ 2.783, the UCL is 0.15098.Using the same data set, the gamma and exponential distributions are also used.Moreover, using x ¼ 0.05 and L ¼ 2.783, we obtained the UCL as 0.1539.The Max-EWMA chart is plotted and shown in Figure 1b.From the figure note that an out-of-control signal is detected by the control chart at 11th observation assuming exponential and gamma distributions whereas Figure 1a detects the first out-of-control point at 8th observation using the GE distribution for time and magnitude.Thus, the proposed chart is more efficient for bottleneck data monitoring.

Conclusion
This article proposed a new Max-EWMA chart assuming the generalized exponential distribution to jointly monitor time and magnitude.Different properties of the proposed chart have been assessed assuming different shifts.The results from the simulation study suggested that the chart performs the best when the shape parameter of the time/magnitude distribution is greater than one.Furthermore, the proposed chart has the smallest ARL when the smoothing parameter is less

Figure 1 .
Figure 1.The Max-EWMA control chart for the breakdowns data of a bottleneck machine.

Table 3 .
Run length profile of the Max-EWMA when X $ GEð0:50, 0:005Þ and T $ GEð0:50, 0:005Þ assuming x ¼ 0:05 for different shifts in u x and u t : the time distribution.Also, if there are simultaneous shifts in ðd u x , d u t Þ ¼ ð0:85, 0:85Þ and (0.50,0.70)ARL 1 reduces by 21.4% and 92.7%, respectively.The percentiles in Table5show that the run length distribution is right skewed, that is, the ARL is always greater than P 50 but less than P 75 .For example, when ðd u x , d u t Þ ¼ ð0:85, 0:85Þ, the median is 215.0, the ARL is 291.1, and P 75 ¼ 394.2.

Table 4 .
Run length profile of the Max-EWMA when X $ GEð1, 0:005Þ and T $ GEð1, 0:005Þ assuming x ¼ 0:05 for different level of shifts in u x and u t :

Table 6 .
Run length profile of the Max-EWMA assuming x ¼ 0:05 for different level of shifts in u x :

Table 7 .
Run length profile of the Max-EWMA for different level of shifts in u t : respectively, compared to the ARL 0 .For the shape parameter equals to 1.5, Table10lists the ARL values.It is observed from the table that the ARL reduced by 20.1% (26.2%) when there is 20% shift in the shape parameter of the magnitude distribution and reduced by 46.5% (54.8%) when there is a 20% shift in the scale parameter of the time distribution.Also, when there are

Table 9 .
Run length profile of the Max-EWMA when X $ GEð1, 0:005Þ and T $ GEð1, 0:005Þ assuming x ¼ 0:05 for different level of shifts in d, u t :

Table 12 .
Run length profile of the Max-EWMA when X $ GEð1, 0:005Þ and T $ GEð1, 0:005Þ assuming x ¼ 0:05 for different level of shifts in d and u x :

Table 14 .
Parameter estimates using Different Distributions for Bottleneck Data.

Table A .
2. Run length profile of the Max-EWMA when X $ GEð1, 0:005Þ and T $ GEð1, 0:005Þ assuming x ¼ 0:05 for different level of shifts in u x and c.Table A.3.Run length profile of the Max-EWMA when X $ GEð1:5, 0:005Þ and T $ GEð1:5, 0:005Þ assuming x ¼ 0:05 for different level of shifts in u x and c.