The performance of a modified EWMA control chart for monitoring autocorrelated PM2.5 and carbon monoxide air pollution data

PM2.5 (particulate matter less than or equal to 2.5 micron) is found in the air and comprises dust, dirt, soot, smoke, and liquid droplets. PM2.5 and carbon monoxide emissions can have a negative impact on humans and animals throughout the world. In this paper, we present the performance of a modified exponentially weighted moving average (modified EWMA) control chart to detect small changes when the observations are autocorrelated with exponential white noise through the average run length evaluated (ARLs) by explicit formulas. The accuracy of the solution was verified with a numerical integral equation method. The efficacy of the modified EWMA control chart to monitor PM2.5 and carbon monoxide air pollution data and compare its performance with the standard EWMA control chart. The results suggest that the modified EWMA control chart is far superior to the standard one.


INTRODUCTION
One of the world's greatest health and environmental problems is air pollution, which has steadily increased worldwide due to global industrialization and urbanization (Manisalidis et al., 2020). Air quality monitoring in many countries has shown that the levels of common pollutants have increased from 1980 until now. During 2010-2016, an estimated 55.3% of the global population were exposed to dangerous levels of air pollution (Shaddick et al., 2020). Therefore, a major global health risk to the world's population is air pollution, with more than 90% of people living in areas that do not meet the 2017 World Health Organization's (WHO) recommended threshold of air quality levels (Health Effects Institute, 2019). The WHO data also reveals that in 2016, 6.1 million deaths worldwide (around 11% of the total global deaths) are attributable to air pollution.
A widespread air pollutant consisting of a mixture of solid and liquid particles suspended in air is called particulate matter (PM), in which the physical and chemical characteristics depend on the location. The common chemical components of PM include nitrates, sulfates, ammonium, other inorganic, and biological compounds (Huang et al., 2014). PM is used as an indicator that is relevant to health and is classified in terms of the diameter or width of the particles, e.g., PM10 means the diameter is 10 µm or less and PM2.5 is 2.5 µm or less (Saarikoski et al., 2018). PM2.5 can be human-made and naturally occurring. Air pollution from vehicle emissions and traffic congestion produces both PM2.5 and carbon monoxide (CO) pollution, with major causes being diesel engine exhaust fumes and traffic congestion (Pui, Chen & Zuo, 2014;Tian et al., 2015). Emissions from factories and power plants are also major causes of PM2.5 pollution due to burning fossil fuels, especially coal (Hua et al., 2015). Moreover, the burning of agricultural waste, especially in South East Asia, and forest fires are other major causes of PM2.5 pollution (He, Lui & Zhou, 2020;Lee et al., 2018).
The small size of PM2.5 particles makes them particularly deleterious to health by directly cause respiratory and cardiovascular morbidity. Since the lungs are the primary organs affected by PM2.5, exposure can lead to lung injury by perturbations of the lung microbiome and its associated metabolome mechanism (Li et al., 2020), as well as causing respiratory diseases and lung cancer (Xing et al., 2016). However, PM2.5-induced neuroinflammation and metabolic turbulence may be mitigated by the anti-inflammatory and anti-oxidative activities of fisetin (Xu et al., 2020). Long-term exposure to PM pollution in the air causes extrinsic skin aging, (wrinkles and changes in pigmentation). Moreover, atmospheric pollutants also lead to skin diseases such as atopic dermatitis (Liao, Eie & Sun, 2020).
The trend in Central and Southern Asia has been a rise in PM2.5 levels between 2010and 2016(Shaddick et al., 2020Cheong et al., 2019;Johnston et al., 2019). Recently in Thailand, especially in Bangkok and the surrounding area, the population has been exposed to dangerously high PM2.5 pollution levels due to agricultural waste burning, foreign sources, and industry. The levels of these small particles have increased over the past few years (Wimolwattanapun, Hopke & Pongkiatkul, 2011;Oanh et al., 2013;Pinichka et al., 2017;Mahidol University, 2020). Self-care practices should be the first priority to manage this problem, while PM2.5 pollution levels should be kept under control and closely monitored periodically.
Statistical process control is a quality control approach for carrying out statistical methods to monitor and control process change. In particular, the control chart is a tool that is widely used to monitor processes to detect any changes in them, thereby preventing the occurrence of faults. In many ecological and environmental sources of data, the values are not independent over time and often comprise autocorrelated observations (autocorrelation or serial correlation is a measure of the correlation between current variable values and their past ones). Hence, the assumption in traditional control chart methodology that the observations taken from the process are independent and normally distributed does not hold. Thus, traditional control charts such as the Shewhart control chart have drastically reduced efficiency when applied to serially correlated observations. Therefore, many researchers have proposed the cumulative sum (CUSUM) (Page, 1954) and exponentially weighted moving average (EWMA) (Roberts, 1959) control charts as suitable alternatives for when the observations are autocorrelated. The EWMA control chart is recommended when the observations are not normally distributed or are autocorrelated, as has been determined using the average run length (ARL) statistic (Srivastava & Wu, 1993;Wardell, Moskowitz & Plante, 1994;Jiang, Tsui & Woodall, 2000;Carson & Yeh, 2008;Han et al., 2010). The usefulness of control chart techniques has been investigated for air pollution data in a time series. The CUSUM control chart has been employed for identifying changes in the mean air pollution level (Barratt et al., 2007), as well as for identifying important change-points in the time series of air pollutants measured at a busy roadside location in central London (Carslaw, Ropkins & Bell, 2006). The applicability of the CUSUM control chart for detecting changes in air pollutant concentrations in Delhi was investigated by Chelani (2011), who found that they have been significantly increasing over time. Moreover, standard CUSUM and EWMA control charts have been used to detect a change in air pollutant time series data in Kuwait (Al-Rashed, Al-Mutairi & Attar, 2019); the results reveal that both procedures provided early alarms in the detection of changes in air pollutants. Afterward, a newly structure control statistic was proposed by Patel & Divecha (2011) and was also generalized structure of control statistic namely modified EWMA control chart. Several researchers investigated the performance of modified EWMA chart by different situations of non-normal distribution including Aslam et al. (2017), Herdiani, Fandrilla &Nurtiti (2018) andNoiplab &Mayureesawan (2019). The application of modified EWMA procedures was demonstrated using real-life samples by comparing with the existing charts. The results showed that the proposed control chart is efficient in quick detection of the out-of-control process (Khan et al., 2018;Saghir, Ahmad & Aslam, 2019;Saghir et al., 2020;Aslam & Anwar, 2020;Supharakonsakun, Areepong & Sukparungsee, 2020).
According to these prior studies on control chart performance, the ARL is utilized to measure the robustness of the chart. In this paper, a modified EWMA control chart, which was newly developed from the traditional EWMA procedure, is presented for monitoring and detecting small and abrupt changes in autocorrelated data by evaluating the ARL. Its performance was studied comparatively with the standard EWMA chart for detecting changes in PM2.5 and CO gas level. Roberts (1959) first proposed the EWMA control chart that is very effectively at detecting small process changes. This chart's design parameters are the multiples of widths of the control limit and smoothing parameter. The EWMA statistic is defined as

Modified EWMA chart
where 0 < λ ≤ 1 is smoothing parameter and X t , t = 1,2,3,... isa sequence of autocorrelated observations with the starting value is the process target mean Z 0 = µ 0 . The respective upper and lower control limits for the EWMA chart are where H is the width of the control limit and σ is the process standard deviation. The term 1 − (1 − λ) 2i is close to 1 as i becomes large. Hence, the respective upper and lower control limit will approach their respective steady-state value as and Later, modified EWMA control chart was introduced by Patel & Divecha (2011). They corrected the inertia problem caused by errors in the EWMA statistic by considering past changes as well as the latest change in the process. The modified EWMA chart is useful for detecting process changes in observations that are autocorrelated or independent normally distributed. Recently, Khan, Aslam & Jun (2017) proposed a new control statistic structure that developed from the modified EWMA control chart as follow: where λ is the smoothing parameter; X t , t = 1,2,3,... is a sequence of autocorrelated observations; k is a constant; and the starting value is the process target mean Z 0 = µ 0 . It is similar to the EWMA statistic but with the last term extended. This chart generates a false alarm when the Z n value violates the specified control limit. In general, the upper and lower control limit are respectively given by and where L and σ is the width of the control limit of modified EWMA procedure and process standard deviation, respectively.

Method of evaluating ARL for moving average order q model
A time series is a series of data points ordered in time, and the goal of a time series analysis is usually to make a forecast of the future using time as an independent variable. A usual characteristic of a time series is autocorrelation, which is correlation between observations in the same dataset at different points in time. In other words, autocorrelated data portrays the similarity between observations as a function of the time lag between them. In a time series analysis, an MA(q) process is a common approach for modeling a univariate time series for which the error depends linearly on the current and numerous past values of the error term.
where µis the mean of the series, ε t is a white noise process assumed to be exponentially distributed, θ i is a process coefficient, and the starting value of ε 0 = s is given.
Let L(u) denote the ARL for an MA(q) process with exponential white noise on a modified EWMA control chart (see Appendix A1 for the proof) that there exists only one solution of the integral equation (see Appendix A2). It is obtained by deriving a Fredholm integral equation of the second kind as follows: with in-control process parameter α 0 and out-of-control process parameter α 1 > α 0 .

Numerical results
The ARLs of the explicit formulas are derived using a Fredholm integral equation of the second type, those of the numerical integral equation method are approximated using the Gauss-Legendre quadrature rule with 1,000 nodes, and the control ARL is set to ARL 0 = 500. The numerical approximation of the numerical integral equation and the exact result of the explicit formulas to measure the accuracy in the comparative study according to the relative error is defined as where L(u) is derived from the ARL using the explicit formulas andL(u) is an approximation of the ARL with the numerical integral equation. The numerical results are reported in Tables 1, 2 and 3. Computation of the ARL by using the explicit formulas and the numerical integral equation method on the modified EWMA control chart were carried out with a varied smoothing parameters (λ = 0.05, 0.10, 0.15 and 0.2); constant k = 1; in-control process parameter α 0 = 1; out-of-control process parameter α 1 = (1 + δ)α 0 , where δ is the shift size set as 0.001, 0.003, 0.005, 0.01, 0.05, 0.10, 0.50 or 1.00; and in-control process was ARL 0 = 500 (Tables 1 and 2). The parameters were set as µ = 2; coefficients θ 1 = −0.3, θ 2 = 0.5 with λ = 0.05, 0.1 and θ 1 = 0.1, θ 2 = 0.3 with λ = 0.15,0.2 for an MA(2) process; and coefficients θ 1 = 0.3, θ 2 = 0.5, θ 3 = 0.7 with λ = 0.05,0.1 and θ 1 = −0.5, θ 2 = −0.3, θ 3 = −0.5 with λ = 0.15,0.2 for an MA(3) process. The results indicate that when the smoothing parameter was increased, the value of ARL 1 was reduced.
The results in Tables 1 and 2 show that the analytically explicit expression of the ARL was in excellent agreement with the approximated ARL obtained from the numerical integral equation (NIE) method. The computational time for the numerical integral equation method was around 21 and 23 s for the MA(2) and MA(3) processes, respectively, while the explicit formulas required a computational time of less than one second for both. Table 3 reports the ARL values obtained by using the explicit formulas and numerical integral equation method. The parameters were set as µ = 2; coefficient parameters θ 1 = −0.3, θ 2 = 0.7, θ 3 = −0.5 for MA(3) process; λ = 0.1; and k was varies as 5λ, 10λ, 20λ, or 50λ. The results revealed that when the constant k was large, ARL 1 was reduced.
The performances of the standard and modified EWMA control charts were also compared. These were obtained by the explicit expression when varying λ (0.05 and 0.10)

Notes.
Where λ is a smoothing parameter, θ i is a process coefficient, b is UCL, δ is the shift size and ε is the relative error.

Notes.
Where λ is a smoothing parameter, θ i is a process coefficient, b is UCL, δ is the shift size and ε is the relative error.

Notes.
Where λ is a smoothing parameter, θ i is a process coefficient, b is UCL, δ is the shift size and ε is the relative error.

Notes.
Where λ is a smoothing parameter, b is UCL of the modified chart, and h is UCL of the EWMA chart.
(RMI ) defined as where ARL δ i denotes the ARLs of the EWMA and modified EWMA control charts obtained via the explicit formulas for each shift size and ARL smallest δ i denotes the smallest of the ARLs for each shift size.
The results in Table 4 show that when λ = 0.05, the performance of the modified EWMA control chart was better than the standard one for shift sizes of 0.001, 0.003, 0.005, 0.01,

Notes.
Where λ is a smoothing parameter, b is UCL of the modified chart, and h is UCL of the EWMA chart. 0.05, 0.10 and 0.30, whereas for shift size of 0.50, 1.00, and 2.00, the small RMI of the modified EWMA chart was 0.057529 while that the RMI of the EWMA control chart was 3.266348. When λ = 0.1, 0.15 and 0.2, the modified EWMA control chart was more powerful than the standard one for all cases of shift size with the zero RMI. The results indicate that overall, the modified EWMA control chart was better than the standard one at detecting process changes.

Application of the modified EWMA chart
PM2.5 and CO gas air pollutants are being constantly emitted, which is likely to increase over time in the winter and summer seasons. When the levels of these air pollutants are high (>50 µg/m 3 for PM2.5 (https://en.wikipedia.org/wiki/Air_quality_guideline) and >10,000 ppm for CO (https://www.airqualitynow.eu/download/CITEAIR-Comparing_ Urban_Air_Quality_across_Borders.pdf)), the quality of the ambient air is unhealthy to humans. Increasing PM2.5 concentration can lead to coughing, breathing difficulties, and eye irritation and can be deadly to humans. Table 5 contains a comparison of the ARLs of the modified and standard EWMA control charts obtained via the explicit formulas. PM2.5 and CO measurements were taken every day in January and May, respectively, 2020 by the Pollution Control Department, Thailand. There were small and abrupt changes in the PM2.5 and CO level data in the Din Daeng district of Bangkok (where there is a high volume of traffic) from measurements near a busy road (Chuersuwan et al., 2008). The PM2.5 and CO air pollution level data were tested for autocorrelation in the observations. The Box-Jenkins technique was applied to the two air pollution datasets to determine whether they fit forecast time series data models. The models with the lowest Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) values were considered as optimal. Moreover, t -test statistics proved that the two datasets were autocorrelated. The parameter values for the MA(1) and MA(2) processes were fitted and provided 51.163 for the mean and -0.723 and -0.380 for the coefficients, respectively. The PM2.5 level was found to be significant for the MA(2) process.
The efficiency of the modified EWMA procedure was also emphasized by its performance with the CO level data for the Din Daeng district, Bangkok, Thailand. Table 6 displays the ARLs of the modified and traditional EWMA control charts. The explicit formulas were used for measuring the ARLs of the CO gas level. The data were collected every day in May 2020. The analysis for an MA autocorrelated process resulted in a mean of 1.198 and −0.662, −0.479, and −0.495 for the coefficients of the MA(1), MA(2), and MA (3) processes. The results of the PM2.5 and CO air pollutant data indicate that the modified EWMA control chart was more effective than the standard one for detecting small shifts, and so confirms that it is excellent for monitoring unusual observations with undesirable values in a timely manner for all cases of exponential smoothing parameter.
The efficacy of the control charts was visualized by plotting graphs to showcase the effective results obtained from the proposed procedure in a comparative study. Figure 1 shows that the modified EWMA control chart detected upper PM2.5 shifts at the 7th to 11th and 17th to 21st observations. On the other hand, the standard EWMA control chart only detected shifts at the 10th to 26th observations, as illustrated in Fig. 2. Figure 3 exhibits that the modified EWMA chart detected upper CO level shifts at the 12th, 25th to 26th, and 30th to 31st observations. On the contrary, the original EWMA chart only detected upper CO level shifts at the 30th to 31st observations, as shown in Fig. 4.
The modified EWMA control chart detected the upper change in PM2.5 level at the 7th observation (i.e., the 7th January), which marked the beginning of extreme changes in PM2.5 emissions at the upper level. Meanwhile, the standard EWMA control chart detected the change at the 10th observation (i.e., the 10th January). Although the CO gas level emissions were low and harmless to the human body, the performance of the modified EWMA control chart for detecting the change in CO gas emissions was exemplary.

Notes.
Where λ is a smoothing parameter, b is UCL of the modified chart, and h is UCL of the EWMA chart.

DISCUSSION
The findings reveal that the modified EWMA control chart performed well for the case of smoothing parameter is greater than or equal to 0.1 due to the RMI of the modified EWMA chart being less than the RMI of the EWMA chart. When applied to real data, the modified EWMA control chart performed excellently for detecting shifts in the PM2.5 and CO pollution levels in all cases of smoothing parameter value. The smoothing parameter value of 0.1 is recommended in applications using the modified EWMA control chart. It is a good choice as it is easier to employ and performed better than the original EWMA control chart in all situations tested.

CONCLUSIONS
The exact ARL was provided by deriving explicit formulas that saved significantly on computational time. Therefore, it is an excellent alternative for evaluating the ARL as a measure of the effectiveness of the modified EWMA control chart. The technique showed good aptitude in monitoring and detecting small process shifts, as illustrated by changes in PM2.5 and CO gas levels examples comprising autocorrelated observations fitted to MA(2) and MA(3) models with exponential white noise. The empirical ARL shows that a smoothing parameter value of 0.1 to 0.2 supported the modified EWMA control chart far better than the EWMA control chart for all cases. Therefore, determination of the correct smoothing parameter of the chart should not be disregarded.

ADDITIONAL INFORMATION AND DECLARATIONS Funding
This research was funded by King Mongkut's University of Technology North Bangkok Contract no. KMUTNB-62-KNOW-018. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.