2.1 MCMA and other Mexican cities

^{Huerta and Sansó (2007)} illustrate the behavior of the daily extreme values of O₃ for the center of the MCMA from 1990 to 2002 through a method for both time and space. Under a Bayesian approach, they define the temporal component as a dynamic linear model (DLM) and the spatial element imposed by the evolution matrix of the DLM. As a result, the maximum levels of O₃ decreased in that period, as well as in the region.

On the other hand, ^{Gong and Ordires (2015)} predict the maximum daily concentrations of ozone in the MCMA through several models of artificial intelligence: multiple linear regression, neural networks, support vector machines, random forest and two assembly techniques. They find that prediction errors in relation to the current day are around 50%. Likewise, with maximum values of O₃ in Mexico City, ^{Rodríguez et al. (2016)} present a trend analysis for this pollutant in the five regions, between 2001 and 2014, by means of a generalized distribution model of extreme values, using Open-BUGS/Winbugs. The results show that the monthly maximum of O₃ decreased during this period.

Later, ^{Aguilar and Reyes (2018)} perform an analysis through Haar wavelet transformation for extreme values of O₃ and its precursors (NO₂ and CO) in Mexico City from 2015 to 2016. During these years, notorious changes were implemented in driving restrictions by the local government. The results showed that for multi-day events of O₃ exhibit, T periods are greater than four days, while for NO₂ and CO they are greater than two. In addition, it is shown that these air pollutants are multi-temporal and comprise a correlated group.

By means of a space-time regression, ^{Yeongkwon et al. (2018)} study the behavior of the concentrations of PM_2.5, PM₁₀, O₃, NO₂, CO and SO₂ in Mexico City. Through models that consider hourly, daily, monthly, semiannual and annual averages, they deduce satisfactory results in the forecasts, except for the series of PM₁₀, PM_2.5 and SO₂.

In other cities of Mexico, the problem of air pollution and its monitoring has also become visible; in this regard, ^{Hernández et al. (2004)} present a geospatial analysis of CO with monthly and annual averages from 1995 to 2001 in Toluca (capital of the State of Mexico). They use the software Surfer 7.2 and determine how the dispersion of the CO is favored with winds and precipitations, where the norm established is only surpassed during the winter, a situation that does not happen in the MCMA.

^{Hernández (2009)} carries out a study for the Metropolitan Zona of Guadalajara (MZG), in center-west Mexico, to predict the daily maxima of O₃ monitored by seven stations from 1997 to 2006. By means of the theory of extreme values, this study verifies that with high temperatures and low wind speeds, an increase in O₃ is obtained. In addition, when using non-stationary and atmospheric information (wind direction, relative dampness, temperature, wind speed), it is possible to predict the daily maximum of O₃ in order to have a monitoring activity that protects the population from dangerous concentrations.

^{Corona and Rojas (2009)} perform an air quality diagnosis of the city of Mexicali, Baja California (in the northwest of Mexico) through concentration data of pollutants, which were interpolated with the Surfer 8 software. It was found that air quality (measured through the number of days in which the standard is exceeded) was not satisfactory in the period from 1997 to 2005.

For the city of Hermosillo, Sonora, ^{Cruz et al. (2013)} found that concentrations of heavy metals (Pb, Cd, Ni, Cu, and Cr) in samplings of suspended particles were below the maximum permissible; however, air quality was not satisfactory. ^{Hernández et al. (2017)} made estimates of long-term trends of O₃ in the Metropolitan Zone of Monterrey (MZM) and compared it with the MCMA and MZG through the analysis of hourly, daily and monthly averages. The estimates were elaborated with the tools openair, WindRose, timeVariation and TheilSen of the R software, for the period from 1993 to 2014. It was observed that the trend of O₃ for the MZM increases in spring, summer and autumn, compared to the downward trend of the MCMA. There was not a clear trend for the MZG in the reference period.

2.2 Other cities

Air quality is currently not satisfactory in various urban areas around the world, so there is a need to conduct studies that establish and anticipate the behavior of air pollutants. In this sense, ^{Gramsh et al. (2006)} analyze contaminant trends in Chile using clusters to identify patterns; in particular, PM₁₀ and O₃ are examined, evidencing that concentration levels vary with seasonality. For the case of PM₁₀, concentrations are higher in winter and O₃ in summer. The topographic and meteorological characteristics of the evaluated areas play an important role in the recognized patterns.

Later, ^{Jaramillo et al. (2007)} refer to the importance of O₃ forecast studies in the city of Cali, Colombia, using the Box-Jenkins methodology to analyze the period from April to July 2003. From 2496 pieces of weekly data analyzed, the first 2232 were used for the estimate, while the rest were used to corroborate the results of the model. The best model had the predictive capacity of 8 h, so that preventive measures could be taken beforehand.

In North America, ^{Camalier et al. (2007)} describe and exemplify a generalized linear model (GLM) in which they relate O₃ with meteorological variables to infer its trends in 39 cities of the USA. They establish a separate model for each city and propose another model with splines for the non-linear relationships that exist among the meteorological components, O₃ and seasonal changes. From this analysis, it is evident that O₃ rises when temperature increases and humidity decreases.

^{Brantley et al. (2019)} applied regression by quantiles with natural cubic splines to estimate the trends of CO, NO, NO₂, NO_y (oxides of nitrogen), BC (black carbon) and other air pollutants, in a railyard area of Atlanta (Georgia, United States). Among several findings, they observed similar trends in NO and NO_y as well as differentiated concentrations of pollutants according to their location.

^{Kumar and de Ridder (2010)} made forecasts of O₃ maximum daily concentrations using the GARCH model in association with the FFT-ARIMA model for four urban areas of Brussels and London. They found that certain factors, such as traffic, play an important role in the concentration of the pollutant, since the structure of the GARCH model depends on the site, unlike the FFT-ARIMA model; therefore, the use of the GARCH model not only improves the intervals of prediction, but also makes more accurate forecasts of O₃ maxima.

In Bulgaria, ^{Georgieva et al. (2014)} analyze the concentration of pollutants (NO, NO₂, NO_x, PM₁₀, SO₂ and O₃) from 2011 to 2012 to assess air quality, using factor analysis (FA) and the Box-Jenkins methodology. Multicollinearity was found and three factors were established: F1 = {NO₂, NO, NO_x, PM₁₀}; F2 = {O₃}; and F3 = {SO}, explaining 90.74% of the total variance. On the other hand, SARIMA models were estimated for each of the six pollutants using a forecast horizon of 72 h, with results observed to be accurate. It was detected that in the case of PM₁₀, higher concentrations are present in winter, exceeding national and European standards.

To predict average concentrations of PM₁₀ in three urban areas of Andalusia (Spain), ^{Palomares et al. (2019)} use five methods, divided into parametric (persistence model and multiple linear regression) and nonparametric (adaptive linear neuron, multilayer and radial basis function). The data correspond to the period from 2005 to 2010 and meteorological measurements are used as exogenous variables for the conformation of the models. It was found that the forecasts with non-parametric models were better, therefore including information on weather conditions improves their predictive capacity. However, such models have a disadvantage, since they require a large amount of data.

In Asia, for Hong Kong, ^{So et al. (2007)} carried out an ANOVA analysis of PM_2.5 where long-term trends and spatial variations were evaluated, considering roads, urban environments and rural environments. Samples collected every six days for 12 months were analyzed in strategic sites. Among several results, seasonal variations were found similar: high in autumn and winter and low in summer; moreover, the behavior of PM_2.5 was not similar in the areas considered.

^{Chang and Yao (2008)} recognize an association between economic development and emissions of air pollutants in densely populated cities that account for 20% of the GDP, which identifies the immediate challenges to understand and control such problems. They review long-term trends of air pollutants for the case of three cities in China, through chemical analysis and ANOVA, finding that the increase in pollutants in cities is similar: high in autumn and winter and lower in summer.

^{Zamri et al. (2009)} obtain forecasts of CO and NO₂ of maximum monthly concentrations through the Box-Jenkins methodology. The data covers the period from 1997 to 2006 in four areas of Malaysia. They emphasize that, in general, the 2016 forecasts do not exceed the permissible limits. The prediction of PM_2.5 in Beijing is seen in ^{Aditya et al. (2018)}, where several machine learning models are compared, and logistic regression is chosen. Daily data of the following variables are considered: temperature, wind speed, dew point and pressure. Likewise, for prediction purposes, an autoregressive model (AR) is used and it is affirmed that these models are efficient to predict PM_2.5 levels.

^{Jaiswal et al. (2018)} study the annual trend of CO, NO₂, SO₂, PM_2.5 and PM₁₀ from 2013 to 2016 for Varanasi, India. They use the Mann-Kendall test and perform their forecast using an ARIMA model. The results show that PM_2.5, CO, NO₂ and SO₂ have a decreasing tendency, unlike PM₁₀, which shows an upward trend.

^{Vita et al. (2018)} develop an algorithm using the fractional Kalman filter (FKF) method to predict the concentrations of the following pollutants: O₃, NO, NO₂, SO₂, PM_2.5 and PM₁₀. The Air Pollution Model and Chemical Transport Model Model (TAPM-CTM) is used to measure the concentrations of air pollutants. They compare the results with the Kalman filter and show that the model with FKF has better precision.

The diverse methodologies that have been used to estimate trends and forecasts of different atmospheric pollutants have varied according to the intended objective, being the most important GLM, DLM, extreme values, Box-Jenkins (ARIMA), clusters, splines, Bayesian approaches and machine learning algorithms. The method applied in this study could be considered as simpler, because it has no distributional assumptions to satisfy and its implementation and maintenance are economic. In fact, the analyst with this proposal has a fast description of the pollutant trends and their dynamics. Only empirical evidence of the pollutants under study is used, thus leaving additional variables aside.

3.1 Controlled smoothing and forecasts

The approach used for the Hodrick and Prescott filter through the Kalman filter is followed as suggested by ^{Guerrero (2007)} and 2008). It encapsulates a set of data observed over time (technically a time series) composed by the sum of a non-observable trend and a random component:

where τ _t is the trend, η _t is the random component from which no specific distribution is assumed with Var(η _t ) = σ_η ² and N is the number of data. Through the controlled smoothing method, the random component is gradually reduced, inducing smoothness. Smoothed values are estimated with this method using an index that will be discussed later. The filter is simple to apply for any data series. The structure of the filter comes from the following expression:

where λ is the smoothing parameter that penalizes for smoothness. In the first term, the goodness of fit of the trend with the observed data is sought, while the second is the induction of smoothness. When the parameter approaches zero (λ → 0) the estimated trend approaches the observed data, and when the value of the parameter is ∞ (λ → ∞) the trend tends to a straight line. In a matrix notation, we have

where y = (y ₁ ,…,y _N )', τ = (τ ₁ ,…,τ _N )' and K ₂ is the matrix (N - 2) × N given by

K2= 1-2100…000001-210…0000.....…....00000…01-21

The solution to the problem is

where τ̂ represents the estimated trend given a value of λ, while the variance of τ̂ in terms of mean square error is given by Var(τ̂ ) = (I _N +λK ₂ ' K ₂ )-1 σ _η ^-2, where I _N is the N dimensional identity matrix. Using the positive square root of the diagonal elements from Var(τ̂ ), it is possible to construct intervals for estimated trends.

Given that the second derivative of M (λ) evaluated in τ = τ̂ is a symmetric and positive matrix, we see where (4) produces a minimum. It should be noted that to obtain τ̂, the matrix of dimension N × N has to be inverted, so this calculation can cause instability and lack of precision of the numerical solution when N is large; in addition to obtaining (4), it is required that there is no missing data in the series. In this context, the Hodrick and Prescott filter estimated through the Kalman filter is a solution, in which it is required to formulate a state-space model as follows:

with  xt=τtτt-1, At=2-110, ct'=10 and wt=εt0 where εt and ηt are two independent random errors of zero mean, uncorrelated and identically distributed with var(εt) = σε2 and var(ηt) = ση2. Therefore, the equation of state has the following form:

where, according to ^{Guerrero (2008}), λ = σ _η ² / σ _ε ². To equate the results of the Kalman filter with smoothing with those obtained by (4), it is assumed that σ _ε ² = 1 and σ _η ² = λ.

On the other hand, in ^{Guerrero (2007)} the forecast of the trend, τ̂_N (h), for the period h of τ _N+h with origin in N for h ≥ 1 and d = 2, is

where we assume μ = 0. Likewise, in order to obtain forecasts by intervals, all the available data is considered, and both the trend and its respective estimation interval are recalculated.

3.2 Smoothness index

The smoothness index (proposed by ^{Guerrero, 2008}), which measures the smoothness of the trend, is given by

where tr(.) represents the matrix trace. Note that when λ → 0 the smoothness S(λ; N) → 0, therefore when λ → ∞ we have that (λ; 𝑁) → 1. In other words, the index is between zero and one. The index for practicality can be expressed in percentage terms. It can be shown that the maximum value of the achievable smoothness index for a set of 𝑁 data is given by 1-2/𝑁. In particular, given the amount of data per pollutant available for this study (𝑁 = 3950), the maximum smoothness achievable is 99.95%.

Note that different indexes of smoothness, S (λ; N), generate different trends; therefore, there is no single trend or, under such optic, an optimal trend. In fact, an infinite number of trends could be estimated from a single set of specific data. With this method, estimations can be made from a trend that is exactly the set of observed data to another that is a straight line at the limit. The problem focuses on deciding what percentage of smoothness is appropriate. This technique has been used for time series in economics and demographic issues; however, there is no history of its use in environmental issues.

To decide the percentage of smoothness through the said index, ^{Guerrero et al. (2017)} suggest a certain index of smoothness based on the structure of the data. However, the sample sizes of the time series in that work very different: in this document there are 3950 observations per series compared to 48, 100 and at most 400 observations in that reference. Additionally, the applications illustrated by these authors have substantially less variability than that observed in the series analyzed here. In fact, it should be noted that the index, S(λ; N), depends on two parameters: λ and N. Therefore, the larger the sample size, the greater the index of smoothness, i.e., when N → ∞ we have that S(λ; N) → 1. With this, it is clear that with the amount of data used here, significantly large smocks are required.

The proposal to impose a pertinent smoothness index in this work consists of choosing one that generates trends that are not affected by extreme observations or seasonal effects; this is, a high index of smoothness is imposed such that in the estimated trends, it mitigates the said effects. When a change in the variability of the respective pollutant is detected, a segmented analysis can be done as proposed by ^{Guerrero and Silva (2015)}. In particular, a percentage of smoothness close to 1 will be imposed by the amount of data, which will overcome the possible generation of spurious trends. After multiple tests, a smoothness index of 99.91% (which corresponds to a smoothing parameter of λ = 1 × 10¹¹) is then used for the estimation of contaminant trends.

To choose the appropriate λ, we first selected the observed time series with the highest variance and with seasonality effects. In this case, it was the O₃ time series. Secondly, we tested different values of λ. With the chosen value it was possible to appreciate the underlying trend (see Fig. SM1 in the supplementary material). Third, we estimated the rest of the trends with a similar λ, so that the same smoothness index was imposed for all the pollutants series. This enabled us to make valid comparisons between all the estimated trends.

On the other hand, for the estimation of Var (τ̂ ) and S(λ; N) the effects of any discrepancies in terms of differentiated quality of the data used are mitigated from the statistical point of view, since the estimated trends are standardized with the use of an identical smoothness index for all pollutant series. This does not necessarily occur when using other statistical techniques, where there are automatic criteria for the selection of the smoothing parameter. It is also assumed that the quality of the data for each pollutant is the same throughout the RAMA.

3.3 Precision measurements of forecasts

The accuracy of pollutant forecasts is assessed using the following precision measures (see ^{Petris et al., 2007}):

Mean absolute deviation (MAD):

Mean square error (MSE):

Mean absolute percentage error (MAPE):

where et = y _t - ŷ _t, ŷ _t = τ̂ , which refers to the forecast error in period y _t and ŷ _t. This represents the real and predicted values of the trend in period t. For all these measures, the smaller the value, the better fit the data set to the estimated trend (see ^{Klimberg et al., 2010}).

3.4 Delimitation of the MCMA and monitoring network location

The spatial behavior of air pollutants is not homogeneous for the entire MCMA (e.g., ^{SEDEMA, 2016}, ²⁰¹⁸, as well as other references available at shorturl.at/cuAG1). For this reason, the MZCM was divided into five regions according to the localization of the monitoring units (^{Rodríguez et al., 2016}) and air pollutant reports from the Atmospheric Monitoring System (Mexico City’s atmospheric monitoring system, SIMAT) (shorturl.at/yzQ03). This regions are: northwest (NW), northeast (NE), center (C), southwest (SW) and southeast (SE). It is considered that the estimated trends are important to describe the main characteristics of the five regions studied, each of which contain 34 monitoring units that, according to information provided by SIMAT and RAMA, could be located as shown in Figure 1.

Fig. 1 Location of the monitoring network. Source: own elaboration. 

The type of pollutant monitored in each of the centers was identified. However, since they are not measured constantly, an evaluation was made of each monitoring unit and pollutant by region. Thus, only those centers that recorded information during the 10 years of the study (2008-2018) were taken into account. The catalogue of the 20 stations by regions that monitor the selected pollutants is shown in Table I.

3.5 Data

RAMA monitors the concentrations of eight pollutants (O₃, SO₂, NO₂, NO, NO_x, PM₁₀, PM_2.5 and CO) every hour. Since each pollutant can be captured by more than one monitoring station in the five geographical regions, the maximum daily concentration reached was identified in each station in order to extract its daily maximum in a specific region. This was done during the period from January 1, 2008 to October 24, 2018, obtaining a total of 3950 daily data. Except for O₃ in the northeast region and SO₂ in the northwest and northeast regions, the rest of the pollutants had a certain percentage of missing data (Table II).

With missing values, it is justified to apply the Hodrick and Prescott filter through the Kalman filter as proposed by ^{Guerrero (2008)}, both for estimating their trends and for their respective forecasts. It is worth mentioning that SAS software version 13.2 and R version 3.5.2 were used to process the data and prepare the estimates.

In the MCMA, the so-called Índice Metropolitano de Calidad del Aire (metropolitan air quality index, IMECA), relies on the Official Mexican Standards to determine permissible limits for the preservation of the population’s health. This index may suffer changes due to the updating of some environmental standards such as NOM-025-SSA1- 2014O3 (^{Secretaría de Salud, 2010}), NOM-022-SSA1-2010 (^{Secretaría de Salud, 2014a}) and NOM-020-SSA1-2014 (^{Secretaría de Salud, 2014b}). The IMECA takes O₃, SO₂, NO₂, PM₁₀, PM_2.5 and CO into account. It also considers the pollutants NO_x and NO. Table III shows the criteria to establish air quality according to the concentration corresponding to each of the six criteria pollutants.

4.1 Trends

4.2 Forecasts

It is important to emphasize that no forecasts are made of the maximum observed pollutant concentrations, but of the estimated trends. Thus, these forecasts are made in a seven-day horizon. Note that the observed data for the analysis were registered from January 1, 2008 to October 31, 2018; however, the estimation of the trend was only contemplated until October 25, to then predict the period from 26 to 31 October, 2018. With the intention of verifying the validity of the forecasts, we also estimated the trend with all the data until October 31, 2018. We also obtained the precise MAD, MSE and MAPE for each pollutant by region.

In general, the forecasts of the eight pollutants trends showed a significant proximity between both series (see Table IV), since through the MAPE it is possible to observe that the average error between the estimate of the trend and the forecast was less than 1.5%. Compared to the MAD, these were less than 0.5%, and for MSE they did not exceed 0.2%. The following paragraphs describe in detail what happens with the forecasts by region.