Mixed POT-BM Approach for Modeling Unhealthy Air Pollution Events.

This article proposes a novel data selection technique called the mixed peak-over-threshold-block-maxima (POT-BM) approach for modeling unhealthy air pollution events. The POT technique is employed to obtain a group of blocks containing data points satisfying extreme-event criteria that are greater than a particular threshold u. The selected groups are defined as POT blocks. In parallel with that, a declustering technique is used to overcome the problem of dependency behaviors that occurs among adjacent POT blocks. Finally, the BM concept is integrated to determine the maximum data points for each POT block. Results show that the extreme data points determined by the mixed POT-BM approach satisfy the independent properties of extreme events, with satisfactory fitted model precision results. Overall, this study concludes that the mixed POT-BM approach provides a balanced tradeoff between bias and variance in the statistical modeling of extreme-value events. A case study was conducted by modeling an extreme event based on unhealthy air pollution events with a threshold u > 100 in Klang, Malaysia.

1. Introduction

In the practical application of statistical modeling, the phenomena under investigation generally involve several rare events occurring with very high or very low values, relatively corresponding to numerous data points. The existence of such rare events can affect data distribution characteristics, such as in the form of high skewness and kurtosis, which correspond to long- and heavy-tail behaviors. In environmental phenomena, the occurrence of extreme events can be observed in various fields, such as extreme precipitation and flooding [1], extreme wind speeds [2], droughts [3], extreme temperatures [4], natural hazards [5], and air pollution events [6,7].

The presence of rare and extreme events leads to a heavy-tail property that makes commonly used normal/Gaussian models unsuitable for describing the data distribution [8,9]. Thus, to address this issue, two statistical extreme-value models for extreme events have been developed, namely the generalized extreme-value (GEV) and generalized Pareto distribution (GPD) models [10,11]. Both models have been found to be useful tools for providing good insights into the behaviors of extreme air pollution events [12,13,14,15,16,17,18,19,20,21]. For example, Al-Dhurafi et al. [12] employed the GPD model to represent probability distributions in an unhealthy air pollution index for three highly polluted urban areas in Peninsular Malaysia. Martins et al. [13] employed the GEV and GPD models to investigate the behavior of air pollutants in two large urban regions (Metropolitan Area of São Paulo and Metropolitan Area of Rio de Janeiro) in South America. Masseran et al. [18] used the GPD model as a tool to investigate the risk of occurrence of unhealthy air pollution events for eight urban areas in Peninsular Malaysia. Reyes et al. [19] provided an estimation of the high ozone levels trend for seven of the meteorological stations located in the Metropolitan Area of Mexico city using the quantile function determined from a fitted GEV model. Battista et al. [21] evaluated the behavior of the air pollution level in the city of Rome (Italy) using the GEV model. Overall, all of these studies found that the results obtained with extreme-value models provided valuable information to use as a basis of managing the risk of air pollution events.

Conceptually, the GEV model is derived from the block-maxima (BM) approach, which generally uses a data selection method based on annual maximum extreme events, and this implies that only one extreme event per year can be used for GEV modeling [22,23,24]. However, this process can lead to a loss of the information contained in other large-sample values [12,25]. In reality, hazardous phenomena that occur within a year may not only be those that exert a maximum effect but also those that are the second or third most serious, which may lead to extreme events in different time periods within the year [26]. The BM approach fails to represent actual phenomena involving the occurrence of multiple extreme events within a year [27]. Therefore, a single maximum point per year definitely reduces the quality of the data representing events of interest, which could significantly contribute to the outcome bias of the extreme value analysis [28]. Although the BM size can be adjusted to include additional extreme data points, no clear rules exist to determine the optimal block size for practical applications.

GPD modeling based on the peak-over-threshold (POT) approach is conducted to determine extreme events as data points greater than a particular threshold level [29,30]. This approach can enable flexible selection of extreme data points, allowing the inclusion of a wide range of extreme events in the analysis [31,32,33]. A high number of selected data points in the extreme-value modeling process provides better precision in terms of the inference and parameter estimation of the statistical model, as well as its quantile estimates [34,35]. However, extreme data points determined with the POT approach tend to lead to the problem of dependency behaviors, which implies that there is a bias in the statistical modeling process. Although the threshold can be adjusted to overcome the dependency behavior problem, the process of determining an optimal threshold can provide only a subjective solution. This subjective threshold selection may not be able to provide an interpretable meaning for the problem under investigation [12,18]. For a problem related to an unhealthy air pollution event (air pollution index (API) > 100), a fixed threshold u = 100 can provide a meaningful analysis [36].

In summary, the BM and POT approaches have their own advantages and limitations. Therefore, this article proposes an alternative approach by combining the BM and POT approaches to provide improved tradeoff results in terms of the bias and variance of the fitted model, as seen in particular for a practical application in air pollution modeling.

2. Study Area and API Data

Klang, which is located in Peninsular Malaysia at a latitude of 101°26′44.023 E and a longitude of 3°2′41.701 N, is one of Malaysia’s largest cities, with a land area of approximately 573 km2. The population density in Klang is among the highest of any city in Malaysia, 1034 people per square kilometer [37]. Klang is also the center of important industrial and economic interests of the country. Moreover, Klang has been recognized as the 13th busiest trans-shipment port and 16th busiest container port in the world [38,39]. However, the rapid development of urban commercial and industrial areas in Klang in recent decades has elevated its risk of atmospheric pollution. Therefore, given the importance of the city’s industrial activities, monitoring and evaluating the behavior of extreme pollution events in Klang is crucial. Figure 1 presents a map of Peninsular Malaysia with the location of Klang [40].

Figure 1

(a) Map of Peninsular Malaysia (Klang is identified by the red dot); (b) map of Klang.

The API in Malaysia is derived from the five main pollutants: suspended particulate matter less than ten microns in size (PM10), nitrogen dioxide (NO2), carbon monoxide (CO), ozone (O3), and sulfur dioxide (SO2). Afroz et al. [41] reported that there are three major sources that influence the concentration levels for the CO, SO2, NO2, O3, and PM10 pollutant variables in Malaysia. They categorized thee sources as mobile sources, stationary sources, and open burning sources, which in turn can be decomposed into various factors, including industrial and construction activities, transportation exhaust emission, soil dust, open burning activity, haze events, and so on [42,43,44,45,46]. The observed data for CO, SO2, NO2, and O3 are quantified in terms of parts per million (ppm) units for the mass of a contaminant. The observed PM10 data are quantified in terms of micrograms per cubic meter (μg/m3). Thus, to measure the API indices, these pollutant variables need to be standardized to derive individual indices. Based on the observed data, the calculation of the sub-API indices can be undertaken using the mathematical formulas provided by the Department of Environment Malaysia [47]. The sub-API value for the CO pollutant can be computed using the following equation:

The sub-API value for the O3 pollutant can be computed using the following equation:

The sub-API value for the NO2 pollutant can be computed using the following equation:

The sub-API value for the SO2 pollutant can be computed using the following equation:

The sub-API value for the PM10 pollutant can be computed using the following equation:

Based on these individual indices, the API value at a particular time can be determined based on the highest value among these sub-indices [48,49]. Figure 2 shows a schematic illustration of the process of determining the API value. Table 1 describes the air quality statuses corresponding to API values at particular time [47].

Figure 2

The process of determining the API value.

Table 1

Air quality statuses corresponding to API values.

Pollution Index	Status	Health Effect
0–50	Good	Low pollution with no ill effects on health
51–100	Moderate	Moderate pollution that poses no ill effects on health
101–200	Unhealthy	Worsens the health conditions of high-risk individuals with heart and lung complications
201–300	Very unhealthy	Worsens the health conditions and reduces the tolerance to physical exercise of individuals with heart and lung complications; affects public health
>300	Hazardous	Hazardous to high-risk individuals and public health in general

The higher the API value at a particular time, the greater the threat of the occurrence of extreme pollution events. Such scenarios disrupt the economic activities of the country, lead to a high risk health problems among populations, and negatively affect environmental ecosystems. Thus, it is important to investigate the behavior of air pollution events in order to gain an understanding of these issues. This study used hourly API data from Klang for the period from 1 January 1997 to 31 August 2020 as a case study. Details of the raw API data and how to request access are available on the Department of Environment Malaysia website [50]. This study was interested in the occurrence of extreme air pollution events determined by unhealthy API indices.

3. Extreme-Value Modeling

3.1. BM Approach

In the BM approach, the extreme behavior of an air pollution event can be evaluated based on a maximum recorded API value determined from a particular period, which is defined as a block. If we let be a random variable representing API data following a particular density function F, then the probability behavior of the maximum API value from a particular block can be written as . The density function of the random variable Y is determined from the following equation:

The density of can be used as an accurate approximation model for probability distribution of the variable despite the independent and identical conditions of the original random variable not being satisfied. This flexibility makes conducting extreme-value analysis based on the BM approach for various types of problems, including those involving complex phenomena, plausible [51].

In a practical application involving a dataset, the actual density function F representing the distribution of the phenomenon under investigation is unknown. This concept implies that the distribution of F must first be estimated from the observed data before it can be used in Equation (6). However, small discrepancies in the F determination will lead to substantial discrepancies in . In parallel, subsequent analyses using the density function will contribute to large errors, consequently leading to incorrect results [52]. To overcome this problem, the density function F can be assumed to be unknown, whereas the determination of the density function can be approximated to a particular limiting distribution form as . Mathematically, this limiting distribution G(y) is valid if a sequence of constants and exists, such that: where G(y) is also a nondegenerate function [53]. According to the literature, the limiting distribution in Equation (7) provides the form of the GEV distribution as the following equation: where the location, scale, and shape parameters are denoted by the notations , , and , respectively. These parameters must be estimated to provide a practical distribution function that can be used to represent the data distribution. However, owing to the problem of small sample size deriving from the BM approach, estimation using the likelihood-based approach, such as maximum likelihood estimation, will not generate satisfactory results. As described by Hosking et al. [54], the L-moment method is an effective alternative approach to overcome this issue. The L-moment estimators of the GEV model are given as follows: where . In conjunction with Equations (9)–(11), the L-moment estimators of the terms and must be estimated from the probability weighted moments determined from the following equation: where , , and . The unbiased estimator of is determined as follows: where is the order statistic determined from the sample data (for additional details on the L-moment estimators of the GEV model, refer to [55,56,57]). The quantile function of the GEV model can be obtained by inverting its cumulative density function (CDF) [58], which is given as follow:

The quantile function knowledge in Equation (14) is crucial for determining the return period related to an extreme pollution event. The return period, referring to the probability of the time (T) block extreme event being exceeded, is in every time block period. Thus, for any particular return period , the critical value of the maximum API level can be determined using the following equation:

In parallel, the maximum API level corresponds to a particular period T, which can be computed as follows:

3.2. POT Approach

This approach uses a particular threshold value to isolate data points considered extreme from the rest of the dataset before determining a statistical model to describe an extreme phenomenon. To relate this approach to the analysis of an unhealthy air pollution event, based on Table 1, an API value exceeding the threshold of 100 is determined as the unhealthy air pollution event. Thus, a fixed threshold value indicated by u = 100 has a significant meaning in air pollution studies. Let represent a random variable in the API data following an unknown density function F. Mathematically, an unhealthy air pollution event phenomenon can be represented by a conditional event corresponding to a threshold greater than u. The conditional exceedance distribution function is determined as follows:

This conditional exceedance distribution function has a left endpoint of ; that is, equal to 0. However, if we let represent a random variable in the data above the threshold u, then the empirical cumulative density function is given by the following equation:

Equations (17) and (18) provide the same information [51]. To determine the parametric form of this density function, the limiting distribution of the normalized values exceeding the threshold must be derived. According to the literature, as the threshold approaches the endpoint of the variable, the limiting distribution of the cumulative density in Equation (18) will follow a GPD [11,59], which can be given as follows:

Meanwhile, in terms of the random variable , we obtain: where . The threshold, shape, and scale parameters are represented by , , and , respectively [60,61]. The tail behavior influenced by the existence of extreme data can be described by the shape parameter . If , then the right upper-tail distribution of the data indicates the properties of a medium-sized tail, which implies that the GPD can be approximated to an exponential distribution. On the one hand, if , then it indicates the properties of a short-tailed distribution, which leads the GPD approximation to a Pareto type II model. On the other hand, if , then it indicates the existence of long-tailed behavior that can approximate the GPD to an ordinary Pareto distribution. Moreover, the mean of the data exceeding the threshold u can be described by the scale parameter [25].

To determine the parameter estimates of the GPD model, a method based on its likelihood function can be used, as large data points can commonly be obtained for the POT approach. For the values of , where k is the total number of data points exceeding the threshold u, the log-likelihood function of the GPD is given as follows:

However, no analytical solution can be obtained from Equation (21). Thus, a numerical procedure must be utilized to obtain a final solution for the parameter estimation.

Next, similar to the GEV model, the concept of the return period can be employed in the GPD model. A return period is interpreted as the average period between extreme events exceeding a threshold u within a particular period of time [27], which is derived based on knowledge of its distribution function. The return period of the GPD model is given as the following equation [53]:

If we let , where k is the number of data points exceeding the threshold u, then the return period presented in Equation (22) can be simplified as: which implies that the API level exceeded once every m series of observations on average can be determined using the following equation:

Specifically, for an unhealthy air pollution event indicated by threshold u = 100, Equation (24) can be simplified as

In addition, by manipulating the return period formula, information on the expected return level corresponding to a particular return period can be obtained [62] as follows: where is the quantile function of the GPD model.

3.3. Mixed POT-BM Approach

The relationship between the GEV model (BM approach) and GPD model (POT approach) can be described using asymptotic model characterization [63]. As described above, in a sequence of independent random variables with a common density function F, the maximum value of data in a particular block is determined as . In a large n, the asymptotic density function of Y is , where is the GEV model presented in Equation (8). In parallel, for a particular adequately large threshold u, the density function of the random variable conditional on can be approximated to the density function of the GPD model presented in Equation (20) [53]. In a mixed POT-BM approach, for a series of random variables , the POT data are determined first by grouping them into a block containing data points satisfying the criteria of . However, some extreme events may occur in sequence over time. For example, an extreme air pollution event will likely be followed by another extreme air pollution event for several consecutive hours or days. Thus, this scenario violates the independence assumption for the selected POT data points. To overcome this problem, a declustering technique can be used to filter dependent consecutive extreme values exceeding the threshold u. Declustering is performed by selecting only the POT data blocks with a minimum separation r from one another. According to Karim et al. [64], r = 240 h is sufficient as the minimum number of separation hours between two extreme events.

After the declustering of the POT blocks has been determined, the BM concept is applied to select the maximum value of the data points in each POT block. This step implies that the extreme data points selected by the mixed POT-BM approach can satisfy or at least satisfactorily approximate the properties of an independent variable. Thus, instead of making a pragmatic choice, such as a block size length of one year—as is common practice in the BM approach, leading to a small number of selected data points and to high variance in the statistical modeling—all data points exceeding the threshold u in the POT approach can be selected; however, this choice violates the asymptotic properties of the extreme-value model, thereby implying a high bias. The mixed POT-BM approach can be a reasonable alternative for providing a balanced tradeoff between bias and variance generated by the data selection in the extreme-value statistical modeling.

In a similar vein, based on the relationship between the GEV model (BM approach) and GPD model (POT approach) described above, the argument regarding this matter leads to an intuitive justification that, for a particular adequately large threshold u, several independent maximum block data governed by the GEV model can be obtained, which implies that the mixed POT-BM approach can achieve satisfactory behaviors similar to the BM approach. Thus, the density function corresponding to the mixed POT-BM approach can be approximated to the GEV model presented in Equation (8). However, the POT-BM approach can be a highly reliable method, as it can determine more extreme data points than the common BM approach. Large independent extreme data points determined using this approach increase GEV modeling precision, as a low estimation variance can be obtained. Moreover, the model bias problem can be avoided by integrating the declustering technique in the POT-BM approach.

By taking advantage of the large number of extreme data points generated by the POT-BM approach, the maximum likelihood can be used instead of the L-moments method for parameter estimation. According to Coles [53], the likelihood-based method is equipped with a convenient set of “off-the-shelf” inference properties and can quantify uncertainties during estimation. The log-likelihood function of the GEV model for can be written as where and [57,65]. However, the log-likelihood function of the GEV model presented in Equation (27) cannot provide a simple analytical solution for parameters , , and . Thus, parameters that can maximize the log-likelihood function can be determined via a numerical optimization technique, which is available in many programming language software, such as R [66].

4. Results and Discussion

Before a detailed analysis is conducted, deriving basic information on the descriptive statistics data would be worthwhile. Figure 3 illustrates the time series plot for the observed hourly API data in the Klang area for the period from 1 January 1997 to 31 August 2020, and Table 2 presents several descriptive measures of the data.

Figure 3

Time series plot corresponding to unhealthy air pollution event threshold.

Table 2

Descriptive statistics of API data in the Klang area.

Location	Mean	Standard Deviation	Min. Value	Max. Value	Skewness	Kurtosis
Klang	55.222	20.970	0	543	4.537	65.133

The information shows that most of the time, in healthy conditions, the API values in Klang fluctuate (below a threshold of 100) around the mean, which is 55.222. However, variance is relatively high (20.970). This high variability derives from unhealthy pollution events that occur repeatedly, as shown by the spike in the data points in the time series plot. The maximum API value recorded in the 24 year data is approximately 541, which indicates a hazardous air quality status (as described in Table 1). These anomalies of API values closely relate to the occurrence of haze events in Malaysia [67]. As reported by the Department of Environment Malaysia [68], haze episodes have recurrently occurred in the years 2005, 2012, 2014, and 2015. In particular, with regard to the API value of approximately 541, this scenario stems from a severe haze event that occurred on 11 August 2005. During that period, the government of Malaysia declared a state of emergency in the Klang area. This result implies that the unhealthy and extreme air pollution event risks in Klang should not be taken lightly. Meanwhile, the skewness and kurtosis measures indicate that the data are not normally distributed. This information can lead to the reasonable application of extreme-value modeling for the data.

Moreover, as an API greater than 100 indicates the occurrence of unhealthy air pollution events, a fixed threshold u = 100 was essential in this study as the basis for the extreme-value modeling and analysis. For the POT approach, the selected threshold can affect GPD modeling precision. Thus, as this study used a fixed threshold u = 100 to assess the suitability of a POT approach, a mean residual life (MRL) plot could be utilized. The MRL plot should indicate an approximate linear behavior in the threshold u [53]. Figure 4 shows the MRL plot for the API data. In the threshold u = 100, approximate linear behavior is satisfied. Thus, u = 100 is a valid threshold for extreme-value modeling using the GPD model. Table 3 provides the results of the parameter estimation for each fitted extreme-value model. The estimated parameters of the GEV fitted model based on the BM approach have the largest standard error compared with those of the other approaches, which implies less precise results. The parameters of the GPD fitted model based on the POT approach have the smallest standard error, which implies high precision. Moreover, the estimated parameters for the GEV approximation based on the POT-BM approach can provide a considerably reduced standard error compared with those of the GEV model based on the BM approach, and it is slightly higher than the standard error obtained by the GPD model based on the POT approach. This result indicates that the application of the GEV model with the BM approach was a less precise method to deal with the data of extreme air pollution events, while both the POT method and the GEV with POT-BM approach were found to produce better results with small standard errors. However, as described above, another important criterion that needs to be evaluated for these methods is their dependency behaviors in order to avoid the bias problem in extreme-value modeling.

Figure 4

MRL plot corresponding to threshold u = 100.

Table 3

Results of parameter estimation for each fitted extreme-value model.

Model	Parameter Estimated	Standard Error
GEV based on BM	Location = 143.805	13.293
	Scale = 46.404	17.071
	Shape = 0.415	0.127
GPD based on POT	Shape = 0.2933	0.017
	Scale = 23.7401	0.511
GEV approximation based on POT-BM	Location = 110.863	1.234
	Scale = 12.911	1.376
	Shape = 0.788	0.118

Based on the results of the parameter estimation, Figure 5 presents a comparison of the fitted extreme-value models that clearly shows the most precise model in terms of the parameter estimation standard error. According to Figure 5, the fitted density of the GEV model based on the BM approach provides relatively rough results, indicating that the less precise model could represent the occurrence of extreme phenomena in the dataset. Meanwhile, GPD modeling based on POT and GEV approximation based on the POT-BM approach could provide improved model fitting with high precision. The comparison of PP plots in Figure 6 accords with the results in Figure 5 and Table 3.

Figure 5

Comparison of fitted models.

Figure 6

PP plot for each fitted model.

Apart from the precision of the fitted models in terms of their standard error, representing variance in the model estimation, another important criterion is the independence assumption for the selected extreme data. Violation of the independence assumption leads to the bias problem in the fitted models. Figure 7 shows the autocorrelation function of the data based on the different extreme-value approaches. The original observed hourly API data (top left) have a very strong dependency on their previous time lag, which clearly violates the independence assumption. Thus, the use of the BM approach (top right) could eliminate the problem of dependency behaviors, whereas the use of the POT approach (bottom left) could slightly reduce dependency behaviors, which remained noticeably significant. Moreover, the use of the mixed POT-BM approach (bottom left) made it possible to solve the problem of dependency behaviors.

Figure 7

Autocorrelation function of data based on different extreme-value approaches.

The highlighted result was obtained from the results presented in Table 3 and Figure 5, Figure 6 and Figure 7, which indicate that extreme-value modeling based on the BM approach is unbiased but less precise, with large estimation variance. However, though the POT approach tends to demonstrate high precision in extreme-value modeling, it suffers from bias due to the violation of the independence assumption. Meanwhile, the mixed POT-BM approach can provide an improved tradeoff between bias and variance. In addition, the mixed POT-BM approach can satisfy the unbiased properties of extreme-value modeling with much lower variance compared to the BM approach. In a similar vein, Figure 8 presents a comparison of the return level plot of each fitted extreme-value model. Based on Figure 8, the return level estimate determined by the GPD model based on the POT approach has a small confidence interval, but its estimation is clearly biased compared with the empirical data. For the return level estimate determined by the GEV model based on the BM approach, its precision with the empirical data is satisfactory. However, the return level plot derived from the GEV model with the POT-BM approach was found to produce a better result. The plot shows that the estimated return level curves within the range of the observed data are accurate with a small confidence interval, particularly in short return periods of time, such as 5 or 10 years. However, for a long return period, the confidence interval of the return level estimation is larger. This implies that the accuracy of the return level estimation decreases. Thus, to provide a better assessment, we suggest that the GEV (POT-BM) model needs to be re-fitted with a current set of data to obtain the most up-to-date short term (5 or 10 years of the future) evaluation of extreme pollution event over time.

Figure 8

Comparison of return level plot of fitted models.

This finding implies there is considerable uncertainty with regard to the practical use of the results. Fortunately, the mixed POT-BM approach can provide improved precision with a small confidence interval in the return level estimation. Short return periods less than 10 years can provide satisfactory estimation results for the return level estimates of extreme API events. For moderate return periods between 10 to 50 years, the precision of the estimates decreases as the range of the confidence interval increases. However, information on long-term return periods over than 50 years is not recommended for use as a reference for air pollution management planning and decisions in Klang. Overall, this study concluded that the mixed POT-BM approach is a satisfactory alternative solution for determining extreme data points for the statistical modeling of extreme events, such as unhealthy air pollution.

5. Conclusions

This article proposes the use of a novel alternative approach for determining extreme event data called the mixed POT-BM method for extreme-value modeling. A case study was conducted using the API data from Klang, Malaysia, for the period from 1 January 1997 to 31 August 2020. The extreme events of interest corresponded to unhealthy air pollution events with a threshold u > 100. Since the threshold u > 100 was used as the basis for determining the POT blocks, this approach led to the problem of dependency behaviors among the selected extreme data points, which implied a bias in the statistical modeling process. Thus, to overcome the problem of dependency behaviors present in adjacent POT blocks, a declustering technique was employed. The declustering technique was used by filtering consecutive dependent extreme data points with 240 h minimum separation from one another. It was found that the application of the declustering technique in POT blocks was able to solve the problem of dependency behaviors for selected extreme data. In addition, the BM concept was integrated to determine the maximum data points in each POT block. The extreme data points determined by the mixed POT-BM approach satisfied the independent properties of independent extreme events, with satisfactory fitted model precision results. The estimated parameters for the GEV approximation based on the POT-BM approach were found to produce a considerably lower standard error compared to the original GEV modeling with the BM approach. Apart from that, the return level plot derived from GEV model with the POT-BM approach was also found to produce a good result. The return level estimation based on short return periods of time (5 or 10 years) was found to be more accurate compared to other extreme-value models. Overall, this study concludes that the results obtained by the POT-BM method can provide an improved balance in the tradeoff between bias and variance in extreme-value modeling. Thus, the POT-BM approach is a satisfactory practical application alternative for extreme-value event analysis and modeling.