Air Quality Time Series Based GARCH Model Analyses of Air Quality Information for a Total Quantity Control District

Air quality data collected at 8 monitoring stations located in the central Taiwan Air Quality Total Quantity Control District were analyzed using multivariate statistical factor analyses. Based on the results thus obtained, 2 major factors, i.e. photochemical pollution factor and fuel factors, were selected for the purpose of evaluating their variations and the pattern of mutual influences for the various air pollution species with respect to time series. The evaluation was conducted using a vector time series coordinated with the ARCH (Autoregressive Conditional Heteroscedacity) and GARCH (Generalized Autoregressive Conditional Heteroscedacity) models in addition to being combined with dynamic impact response analyses using a multiple time series model. The results reveal that the current O3 value is affected by the PM10 values of both a one time lag and a two times lag, as well as the NO2 value of one time lag. When the current SO2 is produced, its concentration can be used to estimate the current CO concentration, and the one time lag SO2 concentration also influences the CO concentration. Additionally, results of impact response analyses show that current CO concentration responds to variations in current SO2; this indicates that the existence of SO2 due to incomplete combustion at the pollution source is immediately reflected by the current production of CO without lagging. In this paper, the vector time series is coupled with the (G)ARCH model to convert simple data series into valuable information so that raw data are better and more completely presented for the purpose of revealing future variation trends. Additionally, the results can be referenced by authorities for planning air quality total quantity control, applying and examining various air quality models, simulating the allowable increase of air quality limits, and evaluating the benefit of air quality improvement.


INTRODUCTION
The development of (time series) theory, especially for financial series, was pioneered by Box and Jenkins (1976), who proposed the ARIMA (autoregressive integrated moving average) model for performing series (time series), especially financial series analyses.The classic regression analysis model assumes that variables of residual values are constants and that the expected values for every period remain the same.However, the results of many studies have indicated the time-varying nature of sequence order variables.Hence, Engle (1982) proposed the ARCH (autoregressive conditional heteroscedasticity) model that has been further modified by Carson et al. (2008) to the System-GARCH Model.Following Engle's ground breaking idea, many alternatives have been proposed to model conditional variances, forming an immense ARCH family; for example, the survey of Bollerslev et al. (1992), and Li et al. (2002).Of these models, the most popular is undoubtedly the generalized autogressive conditional heteroskedasticity (GARCH) model of Bollerslev (1986).Some multivariate extensions of these models have been proposed: for example, Ling and Deng (1993), Engle and Kroner (1995), Wong and Li (1997) and Li et al. (2001).In most of these multivariate extensions, the primary purpose has been to investigate the structure of the model, as in Engle and Kroner (1995), and to report empirical findings.
At present, three of the most popular models to capture the time-varying volatility of financial time series are the Generalised Autoregressive Conditional Heteroscedasticity (GARCH) model of Engle (1982) and Bollerslev (1986), the GJR (Glosten, Jagannathan and Runkle) model of Glosten et al. (1992), and the Exponential GARCH (EGARCH) model of Nelson (1991).Multivariate extensions of GARCH models are also available in the literature, such as the Constant Conditional Correlation (CCC) GARCH model of Bollerslev (1990), the Vector Autoregressive Moving Average GARCH (VARMA-GARCH) model of Ling and McAleer (2003), and the VARMA Asymmetric GARCH (VARMA-AGARCH) model of Hoti et al. (2002).
Various studies have also used neural-networks for forecasting air quality (Tsai et al., 2009).Kumar and DeRidder (2010) took heteroskedasticity of the O 3 time series explicitly into account to show how to forecast O 3 with improved confidence intervals.Moreover, their method is capable of making more accurate probability forecasts of ozone episodes in urban areas.Some studies that compare forecasting performances of ARIMA and neural-networks have been carried out by Tang et al. (1991), Shabri (2001), Ho et al. (2002), and Choon and Chuin (2008) but without a single conclusive result.
Numerous papers have been published on the pollution potential and distribution characteristics of various air pollutants.Zhang et al. (2010) has conducted successive measurements and samplings of Asian dust particles in the spring in Beijing since 2001.In this study, the chemical element composition of aerosol particles was characterized through the ground-based samples of size-segregated aerosols collected in Beijing during dust events in the spring of 2004.Huboyo et al. (2011) investigated the characteristics of fine particles (PM 2.5 ) associated with cooking, particularly temporal variations in the mass and number concentrations in a kitchen and an adjoining room.Chou et al. (2006) have studied the relationships among ambient levels of O 3 , NO and NO 2 to improve the understanding of the chemical coupling occurring among them.Cheng and Lin (2010) showed that the low PM 2.5 -to-PM 10 ratio at the concourse in the Taipei main station was likely the result of coarse PM being re-suspended in the station concourse due to passenger movement.Additionally, Han et al. (2011) used O 3 , NO and NO 2 concentrations to forecast the O 3 concentration.Furthermore, the difference in concentrations between the weekdays and weekends, and O 3 concentrations under different meteorological conditions were also discussed.
In recent years, the increasing number of industrial plants, automobiles, and motorcycles in Taiwan has made it difficult to improve the air quality of regions with concentrated sources of pollution regardless of regulations implementing more stringent air quality.Hence, promoting an air quality total quantity control is a strategy for further improving the air quality in Taiwan.The Environmental Protection Administration (Taiwan) divides Taiwan into seven air quality administration districts: Northern Taiwan District, Hsinchu-Miaoli District, Central Taiwan District, Yunlin-Chiayi-Tainan District, Kaohsiung-Pingtung District, Hualien-Taitung District, and Yilan District.Based on actual situations and requirements, the air quality total quantity limitations will be phased in for each district.The first plan to reduce the air quality total quantity will be implemented in two districts that have the worst air quality total quantity, i.e. the Central Taiwan District and the Kaohsiung-Pintung Districts.Therefore, the Central Taiwan District, which includes metropolitan Taichung, Changhua County and Nantou County, were selected for the purpose of conducting this study using the data collected by 8 air quality monitoring stations established by the Environmental Protection Administration (Taiwan) in these districts in order to carry out a multivariate statistical factor analyses.Based on the results, the air quality situations were categorized into three factors, and the two most important factors, i.e. the photochemical pollution factor and the fuel factor, were targeted for further analyses.Because dynamic variations exist among the various air pollution species, the second moment information needs to be completely grasped so that the tendency of time-dependent variations can be analyzed using the ARMA-GARCH model.Hence, patterns and results of the mutual influence among the various air pollution species can be investigated.The results will instantaneously reflect the response and correlation of various air pollution parameters and can be referenced by the government for review and certification of various air quality models in order for agencies to simulate allowable increase limits and evaluate the benefits of air quality improvement.

Selection of Air Quality Monitoring Stations
The Taichung Thermal Power Plant located in this study region was completed in 1989; it is the largest CO 2 emitter in Taiwan.With the development of the Changhua Coastal Industrial Park in this region, the air quality of the whole central region has been made worse.According to monitored air quality data collected by the Environmental Protection Administration (Taiwan), the air quality in this region is categorized in the third air quality protection class with PM 10 (particulate matter with particle size below 10 microns) and O 3 (ozone) seriously exceeding the air quality standards.1 shows the geographical locations of these stations in the central area of Taiwan.

Screening and Manipulation of Data
Prior to performing GARCH simulation studies, the air quality data collected by the 8 air quality monitoring stations were first analyzed using multivariate statistical factor analyses to look for common factors and to investigate the mutual correlations among the 7 air pollution species.In addition to the 5 major air pollutants Standard Index (PSI) as published by Environmental Protection Administration (Taiwan), including sulfur dioxide (SO 2 ), nitrogen dioxide (NO 2 ), carbon monoxide (CO), PM 10 and ozone (O 3 ), two more indices, i.e. total hydrocarbon compounds (THC) and non-methane hydrocarbon compounds (NMHC) were added in this study.The air quality index (AQI) proposed by US EPA uses the pollution potential of these 5 major air pollutants to classify the degree of their influence on human health.In literature, many authors have investigated the distribution of these 5 major air pollutants in the atmosphere (Bhaskar and Mehta, 2010;Hussein et al., 2011).Using all 7 of these indices will ensure the completeness and variety of the air quality data used in the subsequent analyses.Results of factor analyses show that three factors have eigenvalues greater than "1"; they are the organic pollution factor (i.e.NMHC, THC), the photochemical pollution factor NO 2 , PM 10 , O 3 ) and the fuel factor (i.e.SO 2 , CO).Among the above 3 factors, NO 2 , PM 10 and O 3 , which represent the photochemical pollution factor, are subjected to multivariate statistical factor analyses.The loading degrees are 0.866 for NO 2 , 0.751 for PM 10 , and 0.691 for O 3 , whereas the loading degrees for all other air-borne pollutants are not high.In the case of SO 2 and CO, which are represented in the fuel factor, the factor loading degrees are 0.872 and 0.754, respectively, and all other factors in the fuel factor are not high either.If some air pollutants in one factor have relatively high factor loading degrees, these pollutants have relatively high dependence among one another.Only the photochemical pollution factor and fuel factor were used to carry out the GARCH model simulation because the air pollutants included in these two factors cover the 5 air pollutants in the PSI published by Environmental Protection Administration (Taiwan).It can be surmised that using these two factors will fully represent the air quality pollution situation in the Central Taiwan District.The two air pollutants included in the organic pollution factor do not contribute significantly to air pollution problems in this district; they are not stipulated by the Taiwan Environmental Protection Administration as currently evaluated air pollutant standard indices either.Hence, further discussions on these two air pollutants are not included in this paper.
The series of data used in this study consist of 610 sets collected daily by the Environmental Protection Administration (Taiwan) from January 1, 2010 to September 30, 2011, and published in http://www.epa.gov.tw.During this period, some data were incompletely collected because of un-expected instrument down time for repair and maintenance, and all of these incomplete data sets are deleted, resulting in 610 complete sets of data.All statistical analyses were carried out with E-Views for Windows, version 6.0.

ARIMA Modeling (Shumway and Stoffer, 2006)
A time series {x t ; t = 0, ±1, ±2, ..…} is ARMA(p, q) if it is covariance stationary and can be represented as:   > 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 to be covariance stationary, the differencing procedure may be applied to make it stationary.Subsequently, the ARMA(p,q) model can be applied to the stationary differenced time series; the model so constructed is called the ARIMA(p, d, q,) model, where d denotes the order of differencing (Shumway and Stoffer, 2006;Brockwell and Davis, 2002).Parameters  and θ are estimated using the maximum likelihood method (Brockwell and Davis, 2002) in this study.
An inspection of both the autocorrelation function (ACF) and the partial autocorrelation function (PACF) assists in identifying the orders AR(p) and MA(q).In addition, more objectively defined criterions, such as the Akaike information criterion (AIC), the Hannone-Quinn Information Criterion (HIC), the Bayesian Information Criterion (BIC) and the Final Prediction Error (FPE) can also be used to identify the correct of p and q (Brockwell and Davis, 2002;Kumar et al., 2009).

GARCH Modeling (Kumar and De Ridder, 2010)
If values for ε t denote a real valued discrete-time stochastic process, ε t are the innovations of the ARMA process in Eq. ( 1).Engle (1982) defined these parameters as an autoregressive conditional heteroskedastic process in which all ε t are expressed as: where z t is an identically independent distributed process with a zero mean and unit variance.Although values for ε t are serially uncorrelated by definition, their conditional variances equal 2 t  , which might be autocorrelated and, therefore, may change over time.
The variance equation of the GARCH (p, q) can be expressed as (Bollerslev, 1986;Aradhyula and Holt, 1988;Shumway and Stoffer, 2006;Brockwell and Davis, 2002): where α(B) and β(B) are the appropriate polynomials of the lag operator B; D θ (0, 1) is the probability density function of the innovations or residuals with zero mean and unit variance; p ≥ 0; q ≥ 0 (6) If p equals 0, the process reduces to an ARCH(q) process.
Also, if both p and q equal 0, the conditional variance is constant.As in ARMA, the innovation ε t simply reduces to white noise.Serletis and Shahmoradi (2006) proposed an extended version of the VARMA (vector autoregressive moving average)-GARCH model to simulate natural gas price changes (g t ), and electricity price changes (e t ).The model is expressed as follows:

VARMA-GARCH Modeling
where Ω t-1 denotes the available information set in period t-1, and , for j = g, e.
Parameters h t-j and z t-k are introduced to take the anticipated and unanticipated volatilities into account.
is the error correction term from the long run cointegrating regression.

Impact response Analyses Functions
An impulse response function measures the time profile of the effect of shocks at a given point in time on the (expected) future values of variables in a dynamic system.The augmented vector autoregressive model is used for carrying out the impact response analyses: where X t = (x it , x 2t , …, x mt )' is an m × 1 vector of jointly determined dependent variables; w t is a q × 1 vector of deterministic and/or exogenous variables; and {Φ i = 1, 2, …, p}, and Ψ are m × m and m × q coefficient matrices.The following standard assumptions are made (Lütkepohl, 1991): Assumption 2: All the roots of Under Assumption 3, X t would be covariance-stationary, and Eq. ( 10) can be rewritten as the infinite moving average representation: where the m × m coeffcient matrices Ai can be obtained using the following recursive relations: with A 0 = I m and A i = 0 for I < 0, and

Essence of Model Development
The model used to simulate and predict the two selected factors was calibrated using statistical principles and methods in order to expound the significance of the fat tail test, the Ljung-Box series examination, and the ARCH examination.

Fat Tail Test
The investigation of time series empirical distribution often leads to characteristics of the fat tail test.Hence, the assumption of a normal distribution for the air quality data is not the optimal choice.An examination of the results of the skewness, kurtosis, and the Jarque-Bera normal distribution test can be used to determine whether the distribution of modeling errors has fat tails.

Ljung-Box Series Examination
The residual needs to be examined to determine whether it shows a sequencial correlation before evaluating the ARCH and GARCH models.The square of a residual that shows sequential correlations will have the ARCH effect.

Examination of the ARCH Effectiveness
Before conducting simulations using the time series combination of the ARCH and GARCH models, the process of model calibration must be carried out beforehand in order to confirm that the residual series is not related to the first order series, known as white noise, to assure an appropriate model.Next, the residual square examination is used to determine whether the model has the (G)ARCH effect.In this paper, the Q statistics proposed by Ljung-Box are used to examine whether the residual has high order autocorrelation.Only a model that has ARCH effectiveness can be used to carry out iterative non-linear calculations for estimating model parameters.

Evaluating the Photochemical Pollution Factor Simulation Model
Photochemical pollution factors include NO 2 , O 3 , and PM 10.

Analyses of Basic Characteristics
Table 1 lists the basic characteristics of these 3 air pollutants, including average, standard deviation, skewness, kurtosis, and the Jarque-Bera normal distribution examination statistics.All data show the phenomenon of "skewed on right" because they have positive skewness; NO 2 that has the highest skewness (4.86) indicates that several NO 2 data in this series experience the phenomenon of sudden increase.In contrast, PM 10 has a relatively low skewness of only 0.86 because the central air quality total quantity control district is located near two major air pollution sources, i.e. the Taichung Thermal Power Plants and the Changhua Coastal Industrial Park, causing higher PM 10 concentration and the resulting air pollution problem.This district has a relative higher PM 10 concentration during winter, so the skewness is not quite significant.The O 3 skewness of 3.89 also indicates that O 3 is a major species that causes air pollution; however, the number of days for relatively higher O 3 is relatively lower than that for relatively higher PM 10 as indicated by the relatively higher PM 10 skewness as opposed to that of O 3 .As for NO 2 , the original data indicate that it is higher only during specific periods in winter and early spring; its contribution to air quality problem is relatively insignificant, as shown by its high skewness.All three pollution factors have kurtosis values greater than the kurtosis value of 3 for a normal distribution; hence, the data series have the characteristics of seasonal series.Additionally, the statistical results from examining the Jarque-Bera normal distribution show that the 5% significance level is greater than the critical value (degree of freedom = 2 and 2 0.05,2


 5.99).This observation rejects the hypothesis of normal distribution.All three air pollution factors show the characteristics of double fat-tails, indicating that the data series is actually influenced by seasons.

Ljung-Box Series Examination
Results of using the Ljung-Box method to perform series examinations are listed in Table 2.All examination statistical values for L-B-Q(K) are smaller than the critical value so that the critical value so that null hypothesis, which does not conform to the alternative hypothesis, and therefore, cannot be rejected.This indicates that series residuals conform to white noise because they do not show serial correlations, indicating that the model disposition is quite adequate.

Examination of the ARCH Effect
Whether the ARCH effect exists in a series can be examined using the LM (Lagrance Multiplier) statistical quantity, which is expressed as LM, in which T is the number of samples, and R 2 is the determination coefficient of the results obtained by using the OLS (Ordinary Least Squares) regression.T × R 2 obeys the chi-square distribution that has a P degree of freedom.When the LM statistical value is higher than the 5% significance level, the series shows the ARCH effect.The statistical results of the three selected air pollutants shown in Table 3 reveal that all T × R 2 values are greater than the 5% significance level and that the conditional variances of all three parameters have a strong ARCH effect.Hence, using the ARCH effect to explain the photochemical pollution factor is quite appropriate.

Selecting ARCH and GARCH Models
The vector model EACF paired with various combinations of the ARCH and GARCH models has been tested in order to select the most appropriate VARMA(p,d,q)-GARCH(p,q) model for carrying out the simulation analyses.Table 4 lists the results; the VARMA(2,0,1)-GARCH(2,1) combination is selected because it has the lowest AIC and SC values.Table 5 shows the evaluation of parameters used in the VARCH(2,0,1)-GARCH(2,1) model.

Simulation Results
Table 5 shows the resulting equations and other relevant information obtained by carrying out the VARCH(2,0,1)-GARCH(2,1) model simulation.It indicates that when the current PM 10 is produced, the current O 3 concentration cannot be estimated based on the PM 10 concentration because the b 0 t-statistic of 1.29 is less than 1.96, indicating a lack of significance.However, the one time lag and two times Lag PM 10 concentrations do appear to influence the formation of O 3 concentrations.The t-statistic values of b 1 for the one time lag PM 10 (2.56), and b 2 for the two times lag PM 10 (1.99) are greater than 1.96; both indicating significance.Although some environmental engineering textbooks suggest that the production of PM 10 is not directly related to the production of O 3 , Chou (2010) has pointed out that the secondary aerosols of photochemical reactions cause high levels of atmospheric PM 2.5 and PM 10 in Taiwan.Additionally, atmospheric PM 10 includes primary and secondary aerosols, and the secondary aerosols may be classified as either secondary inorganic aerosols or secondary organic aerosols.Chen and Lee (1999) and Chang and Lee (2007) have observed that secondary organic aerosols and photochemical reactions of VOCs are closed related to the formation of atmospheric O 3 .Based on these statements, the production of PM 10 should be closely related to the photochemical reactions characteristic of atmospheric pollutants.In this research, the model simulation results also reveal that the one lag time atmospheric PM 10 concentration is actually related to the prediction of O 3 concentration.Thus, the relationship between atmospheric PM 10 concentration and O 3 production has been confirmed by measured and simulated data.As for NO 2 , its current concentration cannot be used to predict the current O 3 concentration (c 0 t-statistic of 1.00 is less than 1.96, indicating a lack of significance).However, the one time lag NO 2 concentration influences the formation of O 3 (c 1 t-statistic of 2.05 is greater than 1.96, indicating significance), whereas the two times lag NO 2 concentration becomes insignificant (c 2 t-statisticof 1.45 is less than 1.96, indicating a lack of significance).The above analyses indicate that the current O 3 concentration is affected by the one time lag and two times lag PM 10 concentrations as well as the one time lag NO 2 concentration.A possible explanation for this is that the emission of PM 10 and NO 2 to the atmosphere after they are generated by polluting activities will not lead to immediate photochemical reactions; they are subject to a photochemical reaction in order to produce O 3 only during the one time lag and two times lag.PM 10 may stay in the atmosphere for a long time (Chou, 2010); however, through photochemical reaction mechanisms, it will cause the production of O 3 even after the two times lag.Hence, the proceeding of photochemical reactions behind the time when the pollutants are emitted into atmosphere by one time lag or even two times lag.The two times lag PM 10 concentration is capable of affecting the formation of the current O 3 , and the one time lag NO 2 concentration will affect the formation of current O 3 .O 3 concentration is influenced by its own one time lag concentration (the a 1 t-statistic of 37.5 is greater than 1.96, indicating significance), but it is not quite influenced by its own two times lag concentration (the a 2 t-statistic of -0.87 is less than 1.96, indicating a lack of significance).Hence, the current O 3 concentration is influenced by its own one time lag concentration but not its own two times lag concentration.

Analyses of Impact Responses
The AIC (Akaike Information Criterion) proposed by Akaike (Posada and Buckley, 2004) is used in this paper for selecting appropriate lagged variables: where m is the number of variables in a model; T is the number of samples; SSR is the sum of error squares.
The AIC of lagged variables are first tested in order to select the one with the smallest AIC value as the most appropriate lagged variable to be used as the basis for subsequent analyses.The results of analyzing lagged variables listed in Table 6 show that the lagged variables of the eighth period have the smallest AIC, so these lagged variables are selected for carrying the impact response analyses in this study.Fig. 2 shows the responses to one-unit variations of the two photochemical air pollution parameters.Mark "1" on the X-axis represents "current", and mark "2" represent "lag one time".The Y-axis represents the degree of influence by the production or concentration change of an air pollutant on the "current" or "lag one time" concentration of another air pollutant.Fig. 2(A) shows responses of PM 10 and NO 2 to variations in O 3 .When O 3 varies, the current PM 10 and NO 2 show responses; this indicates that the involvement of PM 10 and NO 2 in photochemical reactions will lead to the production of O 3 , so the O 3 concentration is influenced by PM 10 and NO 2 .Therefore, the variation of O 3 concentration can be

Y-axis
generally determined based on concentrations of PM 10 and NO 2 .As for O 3 , it is affected by its own one time lag concentration; hence the current O 3 shows response when the current variation occurs; the influence of the two times lag becomes less and less over time.
Responses of O 3 and NO 2 to variations in PM 10 are shown in Fig. 2(B).The current O 3 does not show any response to the initial variation of current PM 10 until after the one time lag; this is similar to the simulated results.As discussed in the above simulation results section, secondary aerosols produced by photochemical reactions cause high PM 2.5 and PM 10 levels in the atmosphere of Taiwan.Additionally, the one lag time PM 10 is sufficient to predict and influence the atmospheric ozone concentration.Afterward, the O 3 concentration is affected by PM 10 .Fig. 2(B) also show significant responses of current PM 10 to its one lag time; the current PM 10 responds immediately to the production of the current PM 10 .
Fig. 2(C) showing the responses of O 3 and PM 10 to variations in NO 2 indicates that the current O 3 or PM 10 does not have responses to the variations of current NO 2 .When NO 2 is produced by polluting activities, the current O 3 concentration cannot be estimated based on NO 2 concentration until after a one time lag.The figure also reveals that NO 2 has significant response to its own lag time; the current NO 2 responds immediately to variations in its own lag period.

Investigating the Simulation of Fuel Factor
The fuel factors include two air pollutants, i.e.SO 2 and CO.

Analyses of Basic Characteristics
The basic characteristics of statistical parameters such as average, standard deviation, skewness and the Jarque-Bera normal distribution of these two air pollutants are listed in Table 7.All skewness values are positive so that they skew on the right, indicating that several data in the series show the phenomenon of sudden increase.The presence of high PM 10 and O 3 in the atmosphere is the major cause of deteriorating air quality, whereas high concentrations of SO 2 and CO cause very few days of deteriorating air quality problems.Hence, unlike the skewness of 0.86 for PM 10 , the listed skewness of 3.89 for O 3 is somewhat low.In this research, the skewness values are found to be 5.71 for SO 2 , and 17.94 for CO.The kurtosis values are 36.82for SO 2 and 568.98 for CO; both values are greater the maximum coefficient of 3 for a normal distribution.This indicates that these two air pollutants have characteristics of time-dependent series.As for the statistics from the Jarque-Bera normal distribution examination, the 5% significance level for both air pollutants are greater than the critical value, with a degree of freedom of 2 ( 2 0.05,2   5.99).Hence, the hypothesis that these two series reject the normal distribution indicates the characteristics of double fat-tail.In other words, concentrations of both SO 2 and CO are significantly influenced by seasons in addition to being time-dependent so that different concentrations are observed.

Examination of the Ljung-Box Series
The results of the Ljung-Box Examination shown in Table 8 reveal that all L-B-Q(K) examination results are smaller than the critical values.Hence, the null hypothesis cannot be rejected, indicating that the residual for either series does not have sequential correlations.This observation conforms to the phenomenon of white noise so that the model disposition is quite appropriate.

Examination of ARCH Effect
Table 9 shows the results of examining the ARCH effect; both the conditional variance values for SO 2 and CO show the ARCH effect as seen by all T × R 2 values being less than the 5% significance level that indicates significance.Hence, the ARCH effect is appropriate for explaining the fuel factor.

Selecting the ARCH and GARCH Models
Table 10 lists the testing results of using the EACF vector model combined with various combinations of the ARCH and GARCH models for evaluating the fuel factors in order to select the most appropriate combination of models.Both the  AIC and SC values for the VARMA(1,0,1)-GARCH(1,1) combination are the smallest; hence, the VARMA(1,0,1)-GARCH(1,1) combination is considered the most appropriate.Values of the various estimated parameters for the most appropriate VARMA(1,0,1)-GARCH(1,1) model are listed in Table 11.

Simulation Results
Table 11 shows the resulting equations and relevant results of carrying out the VARCH(1,0,1)-GARCH(1,1) model simulation.It indicates that when the current SO 2 is produced, the SO 2 concentration can be used for estimating the current CO concentration because the b 0 statistic of 11.88 is greater than 1.96, indicating significance.Additionally, the one time lag SO 2 concentration also influences the current CO concentration as evidenced by the b 1 statistic of 5.47 exceeding 1.96, to indicate significance.These observations can be explained based on the incomplete combustion of fuel that leads to emission of SO 2 from industries and vehicles.The incomplete combustion, which is much faster than photo-oxidation reactions, causes the current CO.Additionally, the one time lag SO 2 concentration also influences the current CO because SO 2 in the atmosphere is not involved in photochemical reactions; thus, it does not affect CO formation after its release into the atmosphere.Additionally, CO in the atmosphere is not easily dispersed so that its accumulated concentration is rather high.The CO itself is also affected by its own one time lag because the a 1 statistic of 25.68 is greater than 1.96, which indicates significance.However, the current CO concentration is not significantly influenced by the two times lag value (a 2 statistic of -1.38 being less than -1.96 in a positive number indicates a lack of significance.

Analyses of Impact Responses
The analyses results for the lagged variables are listed in Table 12.The fourth term lagged variables having the lowest AIC value are selected for carrying out the impact response analyses for the fuel factor.Fig. 3 displays the impact responses for these two air pollutants to one unit variation.
Fig. 3(A) shows the impact response for SO 2 to variations of CO; the current SO 2 shows immediate response to variations in current CO.This indicates that the production of CO is mainly caused by incomplete combustion of SO 2containing fuel that leads to an increase of atmospheric CO concentration.This conclusion is similar to the simulated  results shown in previous paragraphs that indicated that when the current SO 2 is produced, the SO 2 concentration can be used to estimate the current CO concentration.Fig. 3(A) also reveals that CO is significantly affected by the current CO and the one time lag CO concentrations.Fig. 3(B) shows the impact response for CO to variations in SO 2 ; the current CO shows a response when variations in SO 2 occur.This indicates that the incomplete combustion of SO 2 will lead to the production of current CO with lagging The figure also shows SO 2 obviously responds to the current and one time lag concentrations so that when the current variation occurs, the current concentration responds with a significant influence even on one time lag or two time lag concentrations.It can be seen in Table 5 that the b 1 statistics for the one time lag SO 2 of 2.56 is greater than 1.96.

CONCLUSIONS
The factor analyses carried out using the multivariate statistical analysis in this research leads to the selection of two relatively important factors, i.e. the photochemical pollution factor and the fuel factor because the air pollutants contained in these two factors cover the 5 air pollutants in the air pollution indices promulgated by the Environmental Protection Administration (Taiwan).Results obtained by using these two factors will therefore fully reflect the air quality conditions for the Central Taiwan Air Quality Total Quantity Control District.
The analysis results show that PM 10 has a low skewness of 0.86 because PM 10 is the major air pollutant in the Central Air Quality Total Quantity Control District.Hence, higher PM 10 concentration leads to more serious air pollution especially during the winter time when the atmospheric PM 10 increases.Additionally, the model simulation results reveal that the current O 3 concentration is influenced by the first time lag and the two time PM 10 concentrations as well as the one time lag NO 2 concentration.In other words, the current O 3 concentration is influenced by the previous PM 10 and

Y-axis
NO 2 concentrations.When the current SO 2 is produced, the current CO concentration can be estimated based on the SO 2 concentration, and the current CO concentration is influenced by the one time lag SO 2 concentration.Results of impact response analyses show that O 3 concentration does not respond to variations in current O 3 and NO 2 concentrations until after the one time lag.When PM 10 and NO 2 are produced, their concentrations cannot be used for predicting the immediate current O 3 concentration until the one time lag is reached.The SO 2 concentration responds to the variations in the current CO, indicating that the production of CO is caused by incomplete combustion of SO 2 -containing fuel, thus increasing CO concentration.Responses of CO concentration to changes in current SO 2 concentration also indicate that during incomplete combustion, SO 2 leads to the production of current CO without lagging.
The integrated VARMA-GARCH model used in this research is capable of evaluating the degree of instantaneous variations of each air pollutant in the air quality total quantity control district.Its implementation is expected to improve the model simulation results significantly by considering the heteroscedastic characteristics of data series that have been ignored in previous research.In this paper, the vector time series is coupled with the (G)ARCH model to convert simple data series into valuable information so that raw data are better and more completely presented for the purpose of revealing future trends in variation.
As mentioned earlier, the classic regression analysis model assumes that variables of residual values are constants and that the expected values for every period remain the same.However, the results of many studies indicate the time-varying nature of sequence order variables.These observations have been confirmed in various fields such as economics, banking and financing, but not in environmental engineering applications.The authors are currently conducting empirical research intended to verify the subject in question with experimental data, and we will publish the results in the near future.

Fig. 1 .
Fig. 1.Geographic locations of the 8 air quality monitoring stations established by the Environmental Protection Administration (Taiwan) in the study region.
Analyses of impact response for the various photochemical pollution factors.

Fig. 3 .
Fig. 3. Analyses of impact response for the various fuel factors.

Table 1 .
Basic statistical characteristics of the various photochemical pollution factors.

Table 2 .
Determination of Ljung-Box series correlation for the various photochemical pollution factors.

Table 3 .
Results of ARCH(q) effect examination for the various photochemical pollution factors.

Table 4 .
Results of VARMA-GARCH examination for the various photochemical pollution factors.

Table 5 .
Parameters obtained using the combined vector model and GARCH(2,1) models for the various photochemical pollution factors.

Table 6 .
AIC values of the lagged variables for the various photochemical pollution factors.

Table 7 .
Basic statistical characteristics of the various fuel factors.

Table 8 .
Determination of Ljung-Box series correlation for the various fuel factors.

Table 9 .
Results of ARCH(q) effect examination for the various fuel factors.

Table 10 .
Results of the VARMA-GARCH examination for the various fuel factors.

Table 11 .
Parameters obtained using the combined vector model and GARCH(2,1) models for the various fuel factors.

Table 12 .
AIC values of the lagged variables for the various fuel factors.