VARMA-EGARCH Model for Air-Quality Analyses and Application in Southern Taiwan

: This study adopted the Exponential Generalized Autoregressive Conditional Heteroscedasticity (EGARCH) model to analyze seven air pollutants (or the seven variables in this study) from ten air quality monitoring stations in the Kaohsiung–Pingtung Air Pollutant Control Area located in southern Taiwan. Before the veriﬁcation analysis of the EGARCH model is conducted, the air quality data collected at the ten air quality monitoring stations in the Kaohsiung–Pingtung area are classiﬁed into three major factors using the factor analyses in multiple statistical analyses. The factors with the most signiﬁcance are then selected as the targets for conducting investigations; they are termed “photochemical pollution factors”, or factors related to pollution caused by air pollutants, including particulate matter with particles below 10 microns (PM 10 ), ozone (O 3 ) and nitrogen dioxide (NO 2 ). Then, we applied the Vector Autoregressive Moving Average-EGARCH (VARMA-EGARCH) model under the condition where the standardized residual existed in order to study the relationships among three air pollutants and how their concentration changed in the time series. By simulating the optimal model, namely VARMA (1,1)-EGARCH (1,1), we found that when O 3 was the dependent variable, the concentration of O 3 was not a ﬀ ected by the concentration of PM 10 and NO 2 in the same term. In terms of the impact response analysis on the predictive power of the three air pollutants in the time series, we found that the asymmetry e ﬀ ect of NO 2 was the most signiﬁcant, meaning that NO 2 inﬂuenced the GARCH e ﬀ ect the least when the change of seasons caused the NO 2 concentration to ﬂuctuate; it also suggested that the concentration of NO 2 produced in this area and the degree of change are lower than those of the other two air pollutants. This research is the ﬁrst of its kind in the world to adopt a VARMA-EGARCH model to explore the interplay among various air pollutants and reactions triggered by it over time. The results of this study can be referenced by authorities for planning air quality total quantity control, applying and examining various air quality models, simulating the allowable increase in air quality limits, and evaluating the beneﬁt of air quality improvement.


Introduction
According to air quality monitoring reports from the Environmental Protection Administration, major pollutants causing the Pollutant Standard Index (PSI) [1] to exceed the air quality standard are particulate matter (PM 10 ) and ozone (O 3 ). Taiwan's air quality data are released to the public in the form of pollution standards index (PSI) values, following the procedure identical to that of the United States Environmental Protection Agency. The nation's PSI was formulated and promulgated in 1994 by the Environmental Protection Administration (EPA) based on the daily concentrations The ten air quality monitoring stations are located on the Pingtung Plain, scattered across areas from 30 to 100 m above sea level. The mean annual temperature of the area is from 24 to 25 • C; the average temperature of the warmest months (July to August) exceeds 30 • C, while that of the coldest months (January to February) remains above 18 • C. The average annual rainfall is between 1500 and 2000 mm. We chose the Kaohsiung-Pingtung Air Pollutant Control Area because this region includes Kaohsiung City, which suffers from the most severe air pollution in Taiwan. The city, home to approximately 2.77 million people, has eight industrial areas and a large steel plant. We chose this Air Pollutant Control Area as our research subject because the air quality is usually poor due to its industrial activities in addition to the large number of cars and motorcycles in the city.
We collected 420 complete air pollutant monitoring statistics between 1 January 2019 and 30 May 2020 and analyzed seven pollutants (or seven variables): SO 2 , NO 2 , CO, PM 10 , O 3 , THC, and CH 4 (please refer to Figure 1 for the locations of the ten air quality monitoring stations chosen for this study). These statistics were acquired by analyzing data collected over 24 h by the EPA's autonomous monitoring stations in the area. The objective in studying the statistics of up to one and a half years is to reflect the temporality effect on air quality; thus, we were able to observe the EGARCH model during different periods with different observation frequencies and elucidate whether variables in different time series present different information. The tool adopted for corroboration was EVIEWS 10.0

Data Selection and Compilation
Prior to the analysis of the EGARCH model, we analyzed three factors using factor analyses in multiple statistical analyses that influences air quality in the Kaohsiung-Pingtung area: the photochemical pollution factor, the fuel pollution factor, and the organic pollution factor, employing factor analysis for multivariate statistics. Next, we selected the photochemical pollution factor, which is the most influential factor determining air quality in the Kaohsiung-Pingtung area, to study the VARMA-EGARCH model.
is ARMA (p, q) if it is covariance stationary and can be represented as where ϕ p 0, θ q 0, and ε t are the innovations with N(0, σ 2 ε ) and σ 2 ε > 0. The parameters p and q are called the autoregressive [AR(p)] and the moving average [MA(q)] orders, respectively. When a time series does not appear covariance stationary, the differencing procedure may be applied to make it stationary. Then, the ARMA (p, q) model can be applied to the stationary differenced time series and model so constructed is called ARIMA (p, d, q) model where d denotes the order of differencing [11,12]. The parameters ϕ and θ have been estimated using maximum likelihood method in the present study.
An inspection of autocorrelation function (ACF) and partial autocorrelation function (PACF) helps in identifying the orders AR(p) and MA(q). In addition, more objectively defined criterions such as Akaike Information Criterion (AIC), Hannan-Quinn Information Criterion (HQIC), Bayesian Information Criterion (BIC) and Final Prediction Error (FPE) can also be used to identify the correct orders p and q [12,13].

ARMA-EGARCH Modeling
Conditional mean formula: > 0 represents good news on the market while 0 reflects bad news on the market; where, γ i = 0 indicates that the degree of air pollution has a symmetrical effect in reaction to news impact. If γ i < 0, it suggests that the degree and fluctuation of air pollution incurred by air pollution shock are more significant than the situation when air pollution shock does not result in air pollution; in other words, a leverage effect can be observed when γ i < 0. To determine whether the three GARCH models are applicable to the time series analysis, we first needed to confirm whether the time series has an ARCH effect; this was tested in this study using the Lagrange Multiplier (LM) test proposed by Engle [14]. Since it is necessary for the residual of the conditional mean formula in the GARCH model to follow the white noise process, all possible orders (p, q) of the GARCH model must be compared with each other via trial and error. To determine whether or not residuals reached white noise, we later adopted the modified Q-statistics proposed by Ljung and Box [15] to analyze the residual of each series. However, the p and q orders of the model selected via trial and error might be subject to over-fitting. To solve this issue, we used the Akaike Information Criteria (AIC) [16] and Schwartz Bayesian Criterion (SBC) [17] to analyze the model under the principle of parsimonious parameterization.

Setting of the Model
This study used the following statistical principles and methods to conduct its simulation and predictive modeling for the selected photochemical pollution factors. The purpose was to explain the essence and meanings of the fat tail test, the Ljung-Box sequence test, and the ARCH test.

Fat Tail Test
The main argument in this section is related to the last argument of the previous section-in the light of a non-explained empirical fact, it is preferable to try to describe it first with the available tools before increasing the model complexity. There are two main approaches to explain fat tails through Gaussian-based models. First, the fat tail test appears because the volatility of asset returns is dynamic. Therefore, looking at the unconditional distribution as described by histogram plots or kernel density plots, we can see heavy tails but they are mostly due to the time varying volatility. Second, fat tails appear because asset returns depend, in a non-linear fashion, on other factors which are distributed according to a Gaussian law. Results of examining the skewness, kurtosis, and Jarque-Bera normal distribution can be used for determining whether the distribution of modeling errors has fat tails.

Ljung-Box Sequence Test
It was necessary to test whether the residual items in the regression model have sequence correlation before estimating the ARCH and (E)GARCH models. If the residual items have sequence correlation, the squared residual items will be examined to see if it has an ARCH effect. As such, it is very important to check if the residual items have sequence correlation before estimating the ARCH and (E)GARCH models.

Examination of the ARCH Effectiveness
If the standard deviation of a time series is stationary, we can say that the variance of the time series is homogeneous, or heterogeneous, or vice versa. Before applying data to a heterogeneous variance model, we should examine data to see if they contain heterogeneous variance. This study adopted the Lagrange Multiplier (LM) proposed by Engle (1982) and Ljung and Box Q statistics [14] to test whether the variance of the time series is heterogeneous.

Impact Response Analyses
The impulse response function (IRF) is used to study the impact of the structural disturbances on the variables in a VAR model over time; in other words, other variables' dynamic response pattern to an exogenous impulse when the said impulse impacts a variable. This study adopted the Cholesky method in dealing with orthogonalization, in order to analyze the intertemporal dynamic effects of indices across nations in the model, and the degree of dynamic interplay among those indices. The formulae are as follows: In Formula (3), i represents the variable responding to an impulse. Error ε t−1 can be interpreted as unexpected impulses in t − i. To examine how an impulse of a one-term delay influences the current term variables, one can refer to Formula (4).
In the above formula, the coefficient factor on row j, line k, formed a consecutive function over time, called the "impulse response function". This study adopted the AIC proposed by Akaike [16] to select suitable lagged variables, which were represented as follows: where, m indicates the number of variables in the model, T represents the number of samples of the three air pollutants, and SSR shows the sum of squares of the error terms.

Application of the Results of Factor Analysis
As mentioned above, this study divided air quality in the Kaohsiung-Pingtung area into three major factors according to the results of the factor analysis for multivariate statistics. It then chose the most important factor, the photochemical pollution factor, as its subject. This factor included three air pollutants: PM 10 , O 3 , and NO 2 (listed by the degree of their factor loadings). Moreover, since the data series had seasonality and cyclicity, the seven pollutant variables were standardized during factor analysis. The formula for standardization is as follows: where, Y v,t represents the original series, and µ t , σ t indicate the average and standard deviation in period 1~w.

Analysis of the Basic Properties of the Three Air Pollutants
The three air pollutants selected for the photochemical pollution factor of this study were PM 10 , O 3 , and NO 2 . Table 1 presents the analysis results of the three air pollutants' basic properties between 1 January 2019 and 30 May 2020, including: average, standard deviation, skewness, kurtosis, and statistics obtained from the Jarque-Bera normality test. Firstly, in terms of kurtosis, the kurtosis coefficients of the three air pollutants were larger than the coefficient of normal distribution (which was 2), suggesting that each variant had the characteristic of a seasonal time series (in other words, pollutant concentrations vary according to seasons). In terms of skewness, the three pollutants' skewness leaned to the right (meaning that their skewness coefficients were all positive). Among them, NO 2 had the highest value of 4.3 when compared with PM 10 and O 3 , which had values of 2.3 and 0.8, respectively (the smaller the value, the more stable the concentration). This result indicates that the concentration of NO 2 normally remains low. However, many statistics showed sudden spikes in the past, which corresponded to the results of Kuo and Ho [18]. This corroborated with the view that NO 2 had the strongest asymmetric effect, meaning that the concentration and degree of fluctuation of NO 2 formed in Kaohsiung were lower than those of the other two pollutants. Since NO 2 is least affected by the change of seasons, it is more difficult to predict its concentration fluctuations. Air pollutants that are non-compliant with the air quality standard in Kaohsiung and Pingtung are mainly PM 10 and O 3 ; hence, when the concentrations of PM 10 and O 3 increase, air pollution becomes more severe, especially during the winter. This result is reflected by the skewness of the two in Table 1, which is not high. Regarding NO 2 , an increase in its concentration can only be observed in Kaohsiung and Pingtung during certain periods in winter and early spring; since the usual contribution of NO 2 to air pollution was less than PM 10 and O 3 , it had a higher skewness coefficient. Kaohsiung City, home to approximately 2.77 million people, has eight industrial areas and a large steel plant. Air quality is usually poor due to the industrial activities in those areas and the large number of cars and motorcycles in the region. In addition, the Kaohsiung-Pingtung area is susceptible to pollutants from external sources brought by the north-eastern seasonal wind, and pollutants traveling from northern and central Taiwan. To make matters worse, pollutants are trapped in the area by an inversion caused by the seasonal wind as it descends from the Central Mountain Range. In addition, emissions from factories and cars concentrated in Greater Kaohsiung make poor air quality more likely in the region during winter. Base on the above reasons, the air pollution in this area has been quite serious over the years.  Table 1 also shows that the statistics obtained from the Jarque-Bera normality test are higher than the critical value (degree of freedom is 2, X 2 0.05,2 = 5.99) at a 5% significance level, reflecting the hypothesis that all variants refuse normal distribution. This means that all the three air pollutants had two fat tails, which proves that seasonality had a significant effect on them.

Examination of ARCH Effectiveness
The LM (Lagrance Multiplier) statistics [19] can be applied to examine whether the ARCH effect exists in a number sequence. The LM statistics is TR 2 with T being the number of samples, and R 2 being the determination coefficient value obtained using the ordinary least squares (OLS) regression; T × R 2 obeys the chi-square distribution with P degrees of freedom. When the model LM statistics is obvious, the ARCH effect exists in the number sequence. The statistics of three air pollutants listed in Table 2 indicates that the conditional variance of the three air pollutants shows a strong ARCH effect (that is, all T × R 2 values were significant at a 5% significant level). Hence, the ARCH effectiveness is appropriate for explaining these three air pollutants.

Ljung-Box Sequential Examination
Whether the residual in the regression model has serial correlation must be examined before the ARCH and EGARCH models are estimated. If the residual has serial correlation, the square of residual appears to have the ARCH and EGARCH effects. In this study, such examination was carried out using the Ljung-Box examination method; the results are listed in Table 3. All the examination statistics for L-B-Q (K) are smaller that the critical value so that null hypothesis, which is not conforming to alternative hypothesis, cannot be rejected. Hence, residuals of the various number sequences do not have serial correlation; this confirms to the phenomenon of white noise so that the model disposition is appropriate for these 3 air pollutants. To run the Ljung Box test by hand, we must calculate the statistic Q. For a time series γ of length n: where: γ j = the accumulated sample autocorrelations, m = the time lag. We reject the null hypothesis and say that the model shows lack of fit if where X 2 1−α,h is the 1 − α quantile of the chi-square distribution with h degrees of freedom.

Choosing the Best EGARCH Model
In order to subsequently simulate the best VARMA-EGARCH model, we tested different combinations with vector models VARMA and EGARCH. We chose the best one from the multiple VARMA (p, q)-EGARCH (p, q) models we created for the simulated analysis. Table 4 presents the analysis results. Among the models, VARMA (1,1)-EGARCH (1,1), the one with the smallest AIC and SC values, was chosen as the best model because it could best capture the fluctuations of air quality in different seasons.

Simulation Results of the Photochemical Pollution Factor VARMA (1,1)-EGARCH (1,1) Model
According to the results of an extra factor analysis for multivariate statistics, factors under the photochemical pollution factor can be listed, with their factor loadings from big to small, as P 10 > O 3 > NO 2 . Photochemical reactions in the atmosphere are mainly triggered by radiation given off by the sun. Once pollutants (or precursors) absorb photons and enter an electronically excited state, they react with other pollutants and form O 3 . In light of this phenomenon, this study designated O 3 as the dependent variable and PM 10 and NO 2 as independent variables and applied them to the VARMA (1,1)-EGARCH (1,1) model, in order to study how they change in the time series. Table 5 presents the relativity and trend of concentration of PM 10 , O 3 , and NO 2 in the time series under simulation with the VARMA (1,1)-EGARCH (1,1) model. This vector model can prevent bias when estimating conditional variance because it reflects the structural changes of the three air pollutants in different seasons.
The simulation results in Table 5 show that the concentration of O 3 (the t-statistic of b 0 was −0.75, meaning that it did not reach the significance level since it was <1.96) in the term could not be estimated according to the concentration of PM 10 in the same term. In contrast, the concentration of PM 10 with a lag time of one and two terms influenced the concentration of O 3 in the current term (the t-statistics of b 1 and b 2 were 2.56 and 2.04, respectively, meaning that they reached the significance level since they were >1.96). In terms of NO 2 , we found that the concentration of NO 2 in the current term could not be used to estimate the concentration of O 3 in the same term (the t-statistic of c 0 was 0.32, meaning that it did not reach the significance level since it was <1.96). However, the concentration of NO 2 with a lag time of one term influenced the concentration of O 3 in the current term (the t-statistic of c 1 was 3.79, meaning that it reached the significance level since it was >1.96). Moreover, the concentration of NO 2 with a lag time of two terms still influenced the concentration of O 3 in the current term (the t-statistic of C 2 was 2.92, meaning that it reached the significance level since it was >1.96). However, the effect of NO 2 with a lag time of two terms on the concentration of O 3 was not as remarkable as that of NO 2 with a lag time of one term. From the above analyses of the three air pollutants, we found that the concentration of O 3 was not affected by that of PM 10 and NO 2 in the same term, but after a lag time of one term, the concentration of O 3 was influenced by PM 10 and NO 2 with a lag time of one to two terms. The results help to explain why PM 10 and NO 2 are the most important precursors for photochemical reactions in the atmosphere. Upon being released into the air from various sources of emission, PM 10 and NO 2 do not immediately start photochemical reactions with the sun; rather, they only start reacting with the sun after a one-term lag and form the final product, O 3 . The photochemical reactions continue into the second term because PM 10 and NO 2 last longer in the atmosphere, thus being able to produce O 3 in the second term. However, from the simulation results in Table 5, it is found that the t-statistics of b 2 and c 2 (2.04 and 2.92, respectively) with a lag time of two terms reached the significance level, but their significance was less remarkable than that of b 1 and c 1 (2.56 and 3.79, respectively) with a lag time of one term. From the results, we can conclude that when PM 10 and NO 2 form in the atmosphere, we cannot use their concentrations to estimate the concentration of O 3 in the current term. Rather, we must wait until they have been involved in photochemical reactions for a while, as that is when O 3 starts to form. In other words, only PM 10 and NO 2 with one-or two-term lag can influence the formation of O 3 in the current term. Furthermore, one-term lagged PM 10 and NO 2 have stronger impacts on the concentration of O 3 than those with a lag time of two terms. Regarding O 3 , its concentration in the current term was significantly impacted by O 3 with a lag of one term (the t-statistic of a 1 was 7.03, meaning that it reached the significance level since it was >1.96), but was not significantly impacted by O 3 with a lag of two terms (the t-statistic of a 2 was 1.14, meaning that it did not reach the significance level since it was <1.96). These results suggested that although the concentration of O 3 in the current term is affected by that of O 3 with a one-term lag, the concentration in the current term is less impacted by O 3 with a two-term lag, as photochemical reactions start to decrease when the day proceeds into the evening. In this study, Kaohsiung-Pingtung Air Pollutant Control Area is located in the tropics. The annual average temperature is greater than 20 • C, and the photochemical reaction is the most obvious from noon to afternoon and in summer. In addition, the degree of air pollution generated in this area is also the most serious from noon to afternoon. As a result, O 3 and PM 10 concentration are usually higher during this time.
In regard to the analysis of γ i,t in the EGARCH model, the t-statistic of γ 1 in the VARMA (1,1)-EGARCH (1,1) model was a significant at 3.13, with a value of −0.082 (<0). This result suggests that when the concentrations of PM 10 and O 3 change drastically in winter, air quality in Kaohsiung and Pingtung will become noncompliant with the air quality standard. The Kaohsiung-Pingtung area is susceptible to air pollution because it is a populous industrial hub located in southern Taiwan. It has a large number of cars and scooters, as well as facilities that release pollutants regularly (namely chimneys). During winter, air quality is particularly prone to being noncompliant with the air quality standard set by the Environmental Protection Administration due to high PM 10 concentration. On the other hand, air pollutants quickly disperse in the air in summer; thus, there are fewer days with poor air quality than in the winter. This phenomenon reflects the leverage effect (when γ 1 < 0) in the EGARCH model [20]. Taiwan has done a magnificent job in containing COVID-19 and has received critical acclaim globally. The pandemic barely has any impact on businesses and industries in the nation. According to statistics up to the end of August 2020, there were only 493 confirmed cases, and only seven people died from COVID-19, a tiny fraction of the nation's 23 million population. The monitoring stations mentioned in this study were still operating during the pandemic and were not compromised by COVID-19 at all. That is, during this period, the monitoring results did not show any abnormal phenomena. Note: (1) O 3 = a o + a 1 O 3(t−1) + a 2 O 3(t−2) + b 0 PM 10(t) + b 1 PM 10(t−1) + b 2 PM 10(t−2) + c 0 NO 2(t) + c 1 NO 2(t−1) + c 2 NO 2(t−2) + d 1 ε t−1 h t = α 0 + α 1 ε 2 t−1 + α 2 ε 2 t−2 + β 1 h t−1 with ε t~N (0,1). A large ε 2 t−1 or h t−1 gives rise to a large σ 2 t . This means that a large ε 2 t−1 tends to be followed by another large ε 2 t , generating, again, the well-known behavior of volatility clustering in time series.
(3) If γ i = 0, it indicates that air quality has a symmetrical effect on the reaction to news shock. If γ i < 0, it suggests that air quality has a leverage effect on the reaction to news shock. (4) A t-statistic value ≥ 1.96 means that the parameter reaches the significance level. In contrast, a t-statistic value < 1.96 suggests that the parameter does not reach the significance level.

The Predictive Power of the VARMA (1,1)-EGARCH (1,1) Model on the Three Air Pollutants
Based on an impact responses analysis, we applied the VARMA (1,1)-EGARCH (1,1) model to estimate the three air pollutants of the photochemical pollution factor, in order to test the predictive power and degree of the model under the influence of leverage effect ( Table 5 shows that γ 1 < 0). Figures 2-4 display the estimated concentrations of the three air pollutants with the first 400 data from the sequence. The estimation results indicated that the correlation coefficients (r) of the predictive power of PM 10 , O 3 , and NO 2 were 0.832, 0.814, and 0.773, respectively, fully reflecting the analysis results mentioned previously (Kuo and Ho, 2018). In other words, the concentration of NO 2 did not change as drastically as the other two pollutants because NO 2 was less affected by the change of seasons, suggesting that it had the strongest asymmetric effect, or NO 2 has the least influence on the GARCH effect when its concentration fluctuations are triggered by the change of seasons, which makes it more difficult to predict its concentration fluctuations. Moreover, the EGARCH model could effectively reduce bias resulting from conditional heteroskedasticity when estimating correlation coefficients. The model enabled the presentation of diverse forms when it comes to dynamic variance modeling, and could accurately estimate conditional variances among the three air pollutants.

The Predictive Power of the VARMA (1,1)-EGARCH (1,1) Model on the Three Air Pollutants
Based on an impact responses analysis, we applied the VARMA (1,1)-EGARCH (1,1) model to estimate the three air pollutants of the photochemical pollution factor, in order to test the predictive power and degree of the model under the influence of leverage effect ( Table 5 shows that < 0). Figures 2-4 display the estimated concentrations of the three air pollutants with the first 400 data from the sequence. The estimation results indicated that the correlation coefficients (r) of the predictive power of PM10, O3, and NO2 were 0.832, 0.814, and 0.773, respectively, fully reflecting the analysis results mentioned previously (Kuo and Ho, 2018). In other words, the concentration of NO2 did not change as drastically as the other two pollutants because NO2 was less affected by the change of seasons, suggesting that it had the strongest asymmetric effect, or NO2 has the least influence on the GARCH effect when its concentration fluctuations are triggered by the change of seasons, which makes it more difficult to predict its concentration fluctuations. Moreover, the EGARCH model could effectively reduce bias resulting from conditional heteroskedasticity when estimating correlation coefficients. The model enabled the presentation of diverse forms when it comes to dynamic variance modeling, and could accurately estimate conditional variances among the three air pollutants.

Conclusions
This study adopted the VARMA-EGARCH model to explore the changes and influences of variation of air pollutants in the time series. This approach is unprecedented because there are barely any published studies on how multiple air pollutants change and influence each other in an Air Pollutant Control Area with the EGARCH model. Moreover, since the parameters of the traditional multivariate GARCH model are too complex, while the multivariate GARCH model cannot fully represent multivariate analysis results, we applied dummy variables from the EGARCH model to the variant formula to enable it to accurately estimate the variances of the three air pollutants and capture the seasonal fluctuations of air quality in the VARMA-EGARCH model. This study also applied the VARMA (1,1)-EGARCH (1,1) model to explore the relativity and trend of concentration of PM10, O3, and NO2 in the time series. When PM10 and NO2 form in the atmosphere, their concentrations cannot be used to estimate the concentration of O3 in the current

Conclusions
This study adopted the VARMA-EGARCH model to explore the changes and influences of variation of air pollutants in the time series. This approach is unprecedented because there are barely any published studies on how multiple air pollutants change and influence each other in an Air Pollutant Control Area with the EGARCH model. Moreover, since the parameters of the traditional multivariate GARCH model are too complex, while the multivariate GARCH model cannot fully represent multivariate analysis results, we applied dummy variables from the EGARCH model to the variant formula to enable it to accurately estimate the variances of the three air pollutants and capture the seasonal fluctuations of air quality in the VARMA-EGARCH model. This study also applied the VARMA (1,1)-EGARCH (1,1) model to explore the relativity and trend of concentration of PM 10 , O 3 , and NO 2 in the time series. When PM 10 and NO 2 form in the atmosphere, their concentrations cannot be used to estimate the concentration of O 3 in the current term. Rather, we must wait until they have been involved in photochemical reactions for a while, as that is when O 3 starts to form. In other words, only PM 10 and NO 2 with one-or two-term lag can influence the formation of O 3 in the current term; one-term lagged PM 10 and NO 2 have stronger impacts on the concentration of O 3 than those with a lag time of two terms. Finally, the estimation of the three air pollutants using the VARMA (1,1)-EGARCH (1,1) model indicated that the concentration of NO 2 did not change as drastically as the other two pollutants because NO 2 was less affected by the change of seasons. This suggested that since NO 2 had the strongest asymmetric effect, it was thus more difficult to predict its concentration fluctuations. We adopted the Pollution Standards Index for evaluating the degree of air pollution. On the other hand, the Air Pollution Index (API) is based on the Ambient Air Quality Standard GB 3095-1996 and only assesses SO 2 , NO 2 , and PM 10 . It was abolished and supplanted by the Air Quality Index (AQI), which is based on GB 3095-2012 and evaluates SO 2 , NO 2 , PM 10 , PM 2.5 , O 3 , and CO; the AQI is mainly adopted for studying the degree of influence of air quality on human health. We did not study the variation of PM 2.5 concentration because PM 2.5 was included in PM 10 ; the subject of this study was not air pollutants' impact on human health, but rather the formation mechanism of various air pollutants and the interplay among these pollutants over time.
In this study, the EGARCH model made it possible to track the degree of change of air pollutants in each time series in real time. It also took the heterogeneity of each air pollutant variant into consideration, an aspect ignored in past studies. The research findings can serve as a reference for the government when it comes to the application and certification of air quality models, simulation of models for a maximum allowable limit of increments, evaluation of air quality improvement schemes, and planning of strategies. In addition, the results of this study can also be referenced by authorities for planning air quality total quantity control, applying and examining various air quality models, simulating the allowable increase in air quality limits, and evaluating the benefit of air quality improvement. Furthermore, we will refer to the API and AQI methods if we are going to delve deeper into studies related to air pollution in the future and could be useful in the future extension of their work.