A multi-dimensional non-homogeneous Markov chain of order K to jointly study multi-pollutant exceedances

In this work we consider a multivariate non-homogeneous Markov chain of order K≥0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K \ge 0$$\end{document} to study the occurrences of exceedances of environmental thresholds. In the model, d≥1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d \ge 1$$\end{document} pollutants may be observed and, according to their respective environmental thresholds, a pollutant’s concentration measurement may be considered an exceedance or not. The parameters of the model are the order of the chain, and its initial and transition distributions. These parameters are estimated under the Bayesian point of view with the maximum a posteriori and leave-one-out cross validation methods used to estimate the order. In the case of the initial and transition probabilities, the estimation is made through samples generated using their respective posterior distributions. Once these parameters are obtained, we may estimate the probability of having no, one or more pollutants exceeding the associated environmental thresholds. This is made using the Markov property as well as a recurrence formula. Results are applied to the case where d=2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d = 2$$\end{document} which will correspond to ozone and particulate matter with diameter smaller than 10 microns (PM10\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{10}$$\end{document}) measurements obtained from the Mexico City monitoring network.


Introduction
Air pollution is a common hazard to which inhabitants of large and medium sized cities around the world are exposed. Among the many pollutants present in those cities are ozone (O 3 ) and particulate matter with a diameter smaller than 10 microns (PM 10 ). The health hazard caused by these two pollutants are well known. For instance, high levels of ozone may make a susceptible part of the population staying in that environment for a long period of time to experience upper respiratory system problems as well as eyes irritation (see, for example, Bell et al. 2004Bell et al. , 2005Cifuentes et al. 2001;Dockery et al. 1992; Galizia and Kinney 1999;O'Neill et al. 2004;Gauderman et al. 2004;Gouveia and Fletcher 2000;Loomis et al. 1996;Martins et al. 2002;WHO 2006). In the case of PM 10 , they may also irritate the respiratory system, and there is evidence that longterm exposure to this pollutant may increase the risk of cardiovascular disease and lung cancer (Gauderman et al. 2004;Mauderly 1997;Mauderly and Oberdörster 1997;Thurston 1996;Itô and Thurston 1996). Additionally, if pregnant women are exposed to these two pollutants, adverse effects may be present in the newborn. Therefore, it is very important to study their behaviours.
A first step to understanding the behaviours of these, and other pollutants in general, is to monitor their levels. Hence, environmental authorities throughout the world have assembled monitoring networks in problematic and not-so-problematic cities. Also, environmental standards have been implemented along with many other actions aimed to decreasing air pollution levels and avoid population exposure.
Mexico is among the many countries adopting measures to control air pollution levels. Since the late 1980s, environmental authorities at the federal level and, in particular, in Mexico City, have been implementing measures such as restrictions on the number of vehicles circulating in the metropolitan area, regulation of more than 300 industries, restrictions on the circulation of some types of vehicles on Saturdays, a decrease of the country's environmental thresholds to declare the occurrence of an environmental alert (NOM 2014a, b). For a timeline of the implementation of the many measures, see, for instance, Achcar et al. (2011).
Nowadays, emergency alerts in Mexico City are based on exceedances of ozone and/or PM 10 with a lately implemented protocol for particulate matter with a diameter smaller that 2.5 microns (PM 2.5 ). In the present work one interest is in estimating the probability of having one of these type of exceedances in a given day, taking into account the present and past information about these pollutants. We are also interested in estimating the probability that in a given day in the year we have either no exceedances, exceedance of one of the pollutants or exceedances of both. In order to do that, we use a bivariate non-homogeneous Markov model of order K ≥ 0. The parameters of the model will be estimated using the Bayesian point of view.
The main aim of the present work is to understand the historical behaviour of the data. As a secondary aim, we are interested in providing a way of warning the possible occurrences of exceedances of the Mexican air quality thresholds for ozone and PM 10 . Hence, if the historical data suggest a high probability of an exceedance occurring, then close attention should be paid to the present behaviour of the data, and if this behaviour suggests the continuity of the past behaviour, then some measures could be taken in advance to prevent population exposure to high levels of pollution. Note that due to the non stationarity of the data, the estimated probability of occurrence of an exceedance may be higher than it really is; however, that is not a bad thing since one of the aim is just to alert the population and be on the safe side.
Non-homogeneous Markov models have already been used to study problems in air pollution. For instance, in Paroli et al. (2005) we have a non-homogeneous Markov mixture of periodic auto-regressive models to study air pollution data from the lagoon of Venice, Italy. A Gibbs sampling algorithm estimates the parameters within the Bayesian framework. In Lagona et al. (2011) non-homogeneous hidden Markov models are used to study multi-pollutants exceedances problems. In that work an E-M type algorithm estimates the parameters involved in the model. Also, Rodrigues et al. (2015) consider a non-homogenous version of the model used in Álvarez et al. (2005) to study the sequence recording the ozone exceedances in Mexico City. The maximum a posteriori method is used to estimate the order and the transition probabilities of the chain. In Rodrigues et al. (2019), a non-homogeneous Markov chain with transition probabilities having a seasonal behaviour rules the sequence recording whether or not an ozone exceedance has occurred in a given day. A similar approach was also taken by Paroli et al. (2005).
Other environmental areas also use non-homogeneous Markov chain models. For instance, Rajagopalan et al. (1996) and Hughes et al. (1999) use non-homogeneous Markov models to study the occurrence of precipitation. In Drton et al. (2003), these models are used to study occurrence of tornado activity.
Although in the present work we also consider the non-homogeneous Markov chain approach and the Bayesian point of view to estimate the parameters present in the model, the novelty here is that two or more pollutants may be studied via a multivariate non-homogeneous Markov chain of order K ≥ 0 model. We assume that the order of the chain as well as its initial and transition probabilities are unknown and need to be estimated. The chosen estimation methods are the maximum a posteriori and the leave-one-out cross validation (see, for instance, Chang 2019) for the order of the chain and, having this information, a sample of the initial and transition probabilities generated using their posterior distributions are used to estimate them.
This work is organised as follows. In Sect. 2, we describe the mathematical model. Section 3 gives the Bayesian formulation of the model. In Sect. 4 an application to the case of Mexico City ozone and PM 10 data is presented. In Sect. 5, a discussion of the results as well as some additional comments are given. An appendix, given before the list of references, presents some additional results as well as a simulation study to corroborate de choice of the order of the chain using the methods considered here. This work also has a Supplementary Material file presenting the tables with the values of the estimated probabilities of interest in the case of the selected model.

A multivariate non-homogeneous Markov chain model
Let N > 0 and T l > 0 be given natural numbers representing, respectively, the number of years and the number of days in the lth year where pollutants' measurements were taken, l = 1, 2, . . . , N . Hence, depending on whether or not we have a leap year, either T l = 366 or T l = 365, l = 1, 2, . . . , N . Following Drton et al. (2003), we take T l = T = 366, l = 1, 2, . . . , N , with the convention that for non-leap year a measurement is set to zero when t = 366. Hence, from now on, we drop the index "l" from the notation of the observational period.

Remark
We could also leave each year with its original length. In Rodrigues et al. (2019) that was performed in the model used there. The aim was to verify whether or not significant changes in the results could occur. The findings were that it does not matter if we take T l = T = 366 for all l = 1, 2, . . . , N , or if we consider the original number of days for each year. Hence, in order to simplify the model and the computer programmes, we keep the convention that T l = T = 366, l = 1, 2, . . . , N , assigning value zero to the measurement if we have t = 366 and a non-leap year.
Assume that there are d ≥ 1 pollutants we are interested in studying. Let Z i,t be the measurement of the ith pollutant on the tth day, t = 1, 2, . . . , T , i = 1, 2, . . . , d. Denote by Z t = (Z 1,t , Z 2,t , . . . , Z d,t ) the vector recording the pollutants' measurements at time t in a given year, t = 1, 2, . . . , T . Indicate by L i > 0 the environmental threshold for the ith pollutant, i = 1, 2, . . . , d. Consider . . . , E d,t ) : t ≥ 0} the process defined as follows, Hence, E i,t indicates whether or not on the tth day of a given year the ith pollutant's measurement exceeds its corresponding environmental threshold, i = 1, 2, . . . , d, t = 1, 2, . . . , T . By definition of the vector E t , it is possible to see that it assumes values on the set This could be viewed as the binary decomposition of a number in S 2 = {0, 1, . . . , 2 d − 1} since we can associate a vector in S 1 to a number in S 2 using the transformation f : S 1 → S 2 given by f (x 1 , x 2 , . . . , x d ) = d−1 i=0 x i+1 2 i . Hence, the process E, which assume values on a d-dimensional space, may be viewed as a unidimensional process assuming values in S 2 . Let W = {W t : t ≥ 0} indicate that process. Therefore, if we denote byw an element of S 2 , let x = (x 1 , x 2 , . . . , x d ) ↔w indicate thatw ∈ S 2 corresponds to the state x ∈ S 1 via the transformation f (·).
Remarks 1. On occasions the "tilde" used to indicate an element in S 2 , will be dropped. That will be made when numbers, instead of letter, are used to indicate the state in S 2 .
2. Since the aim here is to analyse whether or not exceedances of the the thresholds occurred, the processes E and W can be used without problems. Note that we are not interested in studying the possible dependence between the pollutants' measurements and the way the concentration of one of them affects the other. That is studied elsewhere.
3. Up to now we have three sequences. The sequence Z = {Z t : t ≥ 0} which records the actual pollutants' measurements, the sequence E = {E t : t ≥ 0} which is the transformed sequence Z and records whether or not the pollutants' measurements are above the associated thresholds and the sequence W = {W t : t ≥ 0} which records the univariate values associated with E through the function f (·).
As in Álvarez et al. (2005) and Rodrigues et al. (2015), we assume that W is ruled by a non-homogeneous Markov chain of order K ≥ 0. We also assume that K is a random variable with state space S = {0, 1, 2, . . . , M}, where M ≥ 0 is a given fixed integer number. Denote this chain by X (K ) = {X (K ) t : t = 1, 2, . . . T }. Therefore, for each year the sequence W is considered a realisation of the non-homogeneous Markov chain X (K ) . Note that each element X (K ) t of X (K ) has as its coordinates values of the sequence W, i.e., X Thus, X (K ) has as its state space the set S We may also see an element S as an element of S for A a set, we use the notation |A| to indicate the cardinality of A). This may be achieved by using a transformation g(·) similar to f (·), but now defined on S (K ) 3 and with image in S (K ) 4 . This function g : and m ∈ S (K ) 4 we use the notation y ↔ m to indicate their correspondence. Therefore, each K -dimensional state y of the Markov chain X (K ) may be associated with a unidimensional value m ∈ S (K ) 4 via g(·).
From Rodrigues et al. (2015) we also have that, for K ≥ 2, if the set of observed value is (z 1 ,z 2 , . . . ,z T ) ∈ (S 2 ) T , then the transition probabilities of are different of zero, if and only if, the next state is X (K ) t+1 = (z t+1 ,z t+2 , . . . ,z t+K ), i.e., if the observation following z t ,z t+1 , . . . ,z t+K −1 isz t+K . Hence, as in Álvarez et al. (2005) and Boys and Henderson (2002), instead of considering the transition from X . . ,z t+K ), we may consider the transformed state spaces S (K ) 4 and S 2 , and the transition probabilities of X (K ) may be written as (see, for instance, Álvarez et al. 2005), where m ∈ S (K ) 4 ,j ∈ S 2 , and 1 ≤ t ≤ T − K . As in Rodrigues et al. (2015), indicate by Q 2. Unless otherwise stated, from now on, we are going to use the state space S (K ) 4 and the corresponding transition probabilities (2).
In addition to estimating the order K of the Markov chain, we will also estimate its transition probabilities P (K ) mj (t), for each t, as well as the initial probabilities Q , the transition matrix at time t.
As an additional information provided by the model, we are able to estimate P(W t = i),ĩ ∈ S 2 . These values report the probability of whether or not we have exceedances of one, more than one or none of the pollutants. They are obtained by taking advantage of the Markov property and using a recursive form. A way of obtaining these probabilities is given at the end of the next section.
Remark Even though at the end we are dealing with a univariate chain, it reflects the behaviour of a multi-dimensional problem since we are studying the behaviours of several pollutants at once. The transformed univariate chain is used in order to facilitate the analysis.

A Bayesian estimation of the parameters of the model
The order of a Markov chain could be estimated using the auto-correlation function associated to the chain. An alternative method to estimate the order and consequently the transition probabilities is to use the so-called reversible jump Markov chain Monte Carlo algorithm. That was used in Álvarez and Rodrigues (2008). However, the line followed here (using maximum a posteriori and leave-one-out) is simpler and also produces suitable results. In the case of the initial and transition probabilities we have, for instance, the maximum likelihood method (see, for instance, Fleming and Harrington 1978) and the empirical estimators (see, for instance, Aalen and Johansen 1978). In the present work, we use the Bayesian approach (see, for instance, Robert 2007 andCarlin andLouis 2000) to estimate the order and the transition probabilities. Inference is performed using the information provided by the so-called posterior distributions of the parameters. The posterior distribution of a vector of parameters θ given the observed data D, indicated by P(θ | D), is such that P(θ | D) ∝ L(D | θ ) P(θ ), where L(D | θ ) is the likelihood function of the model, and P(θ ) is the prior distribution of the vector θ. These elements are given as follows.

The likelihood function and the prior and posterior distributions
In the present case, the vector of parameter is θ = (K , Q (K ) (1), P (K ) (t), t = 1, 2, . . . , T − K ) which belongs to the following sample space 2, . . . , l; l l=1 x i = 1} is the (l − 1)-dimensional simplex. (If K = 0, then the parametric space reduces to = T |S 2 | .) Let Y = {Y t : t = 1, 2, . . . , T } indicate our observed data. Instead of using the measurements Z t as our observed values, we will use the values in S 2 associated with E t which in turn is associated with both the measurements Z t , t = 1, 2, . . . , T , and the transformed sequence W. Hence, we take Y = W.
Let (x 1 ,x 2 , . . . ,x K ) be such that the observation at time t isx 1 . For a order K ≥ 1, indicate by n (K ) mĩ (t) the number of years in which the vector (x 1 ,x 2 , . . . ,x K ) corresponding to a state m ∈ S (K ) 4 is followed by the observationĩ ∈ S 2 . Also, define n Since a Markov model is assumed, the likelihood function is given by (see, for instance, Drton et al. 2003;Fleming and Harrington 1978;Boys and Henderson 2002;Fan and Tsai 1999) When K = 0, the expression (3) simplifies to where The prior distribution of the vector of parameters is given as follows. We assume a prior independence of P (K ) (t) as functions of t. Also, since the forms of P (K ) (t) and where P(Q (K ) (1) | K ) and P P (K ) (t) | K are the prior distributions of the initial distribution Q (K ) (1) and of the transition matrix P (K ) (t) given the order of the chain, respectively, and P(K ) is the prior distribution of the order K .
Remark When we have K = 0, the vector of parameters is θ = (Q (0) m (t); t = 1, 2, . . . , T ), whose prior distribution simplifies to P(θ | Given the nature of transition matrices, we assume that rows are independent and that, given the order K of the chain, each row of the transition matrix P (K ) (t) will have as prior distribution a Dirichlet distribution with appropriate hyperparameters. Therefore, given that K = k, row (P (k) for t in the appropriate range. The initial distribution Q (K ) (1), will also have a Dirichlet prior distribution, but now with hyperparameters (α (K ) m ; m ∈ S (K ) 4 ). Therefore, If K = 0, then Q (0) (t) has as prior distribution a Dirichlet distribution with hyperparameters (α (0) m (t); m ∈ S (0) 4 = S 2 ), t = 1, 2, . . . , T . Therefore, for any given t, We assume that K has as prior distribution a truncated Poisson distribution with rate λ > 0, defined on the set S; i.e., where I A (x) = 1, if x ∈ A and is zero otherwise. Therefore, from Boys and Henderson (2002); Fan and Tsai (1999); Álvarez et al. (2005); and Rodrigues et al. (2015), the conditional posterior distribution of P (K ) (t) given K , is proportional to the product of Dirichlet distributions with hyperparameters (n (K ) Additionally, the posterior distribution of the initial distribution Q (K ) (1) given K is also proportional to a Dirichlet distribution, but now with hyperparameters (α Álvarez et al. (2005), with the appropriate adaptation for the case of K = 0. Hence, the posterior distribution of the order K is Hence, we will use (6) to estimate the value of K that maximises the posterior distribution P(K | Y). We will also use the leave-one-out cross-validation method (see, for instance, Chang 2019). Once that K is estimated, samples of the posterior transition and initial distributions associated to that K are obtained. The generated values are used to estimate the corresponding distributions via the law of large numbers.

Remark
The hyperparameters of the prior distributions will be considered known and will be specified later when the model is applied to the data.

Estimation of the probabilities of exceedances
Once the posterior distributions are obtained and samples of size M > 0 of the corresponding initial and transitions probabilities are produced, we may estimate the probabilities P(W t =ĩ),ĩ ∈ S 2 , t = 1, 2, . . . , T , using the law of the large numbers in the following way. For each element l of the sample of the initial and transition probabilities it is possible to calculate the a value of P(W t =ĩ) which we indicate by where each P (l) (W t =ĩ) is obtained as follows. (We consider separately the cases where K = 0, K = 1 and K ≥ 2.) Case 1. When X (K ) has order K = 0, then P(W t =ĩ) corresponds to the probabilities If the order of X (K ) is K = 1, then we may take advantage of the Markov property as well as the estimated initial distributions. That is made as follows. Note that, for t ≥ 2 andj ∈ S 2 , we may write P(W t =j) recursively as (Morales-Morillón 2018), where (1). In order to estimate the values of P (l) (W t =j),j ∈ S 2 , we proceed as follows.
Step 1 Generate a sample of size M > 0 of Q Step 1) and generate a sample also of size M of P (1) the lth element of corresponding samples,ĩ,j ∈ S 2 , t = 1, 2, . . . , T − 1, l = 1, 2, . . . , M . Then, for each l = 1, 2, . . . , M , use (8) to recursively obtain P (l) (W t =j). Hence, Case 3. When the selected order of the chain X (K ) is K ≥ 2, then we proceed as follows. Step Step 2 For t ≥ K + 1 and m ↔ (x 1 ,x 2 , . . . ,x K ), consider the generated sample of size M > 0 of Q Step 1) and additionally generate a sample also of size M of the transition probabilities P In order to obtain P (l) (W t =ĩ),ĩ ∈ S 2 , t = K + 1, K + 2, . . . , T , note that, using the Markov property, we may write the following recursive formula, where P(X (K ) Remark As mentioned earlier, besides using the maximum a posteriori to select the order of the chain, we also use a modification to the non-homogeneous case of the leave-one-out cross validation method (LOO) -see, for instance, Chang (2019) and references therein. Once this order is estimated, if a lower order is selected, we may also use, for instance, the well known sum of the absolute values of the differences (SAD) using the estimated and observed initial and transition probabilities to see how well the estimated probability values describe the observed. These methods are described in Section A of the Appendix.

Application to ozone and PM 10 data from the monitoring network of Mexico City
In this section we consider the model presented in Sect. 2, in the particular case of d = 2. An application will be made to measurements of two pollutants collected at Mexico City's monitoring network. This network comprises of several monitoring stations placed throughout the metropolitan area. The pollutants considered are ozone (given in parts per million -ppm) and PM 10 (given in microgram per cubic meter -μg/m 3 ).
Even though ozone in Mexico City has been measured since late 1980s, systematic PM 10 monitoring started only in 1995. Hence, we have decided to use measurements from that year on for both pollutants. Therefore, we have considered measurements obtained from 01 January 1995 to 31 December 2019. In the analysis, we will use the daily maximum measurements for those pollutants. Measurements are taken minute by minute in each station and the hourly averaged result is reported every hour. In a given station the ozone daily maximum measurement is the maximum of the 24-hour averaged results. Ozone daily maximum measurement for the metropolitan area is the maximum of all daily maxima recorded in all stations of the monitoring network. PM 10 daily measurement in a given station is the average results over the 24-hour period for that station. The PM 10 daily maximum measurement for the metropolitan area is the maximum of the daily averaged results for all stations in the network. We will consider the overall measurements for the entire metropolitan area instead of splitting it into regions as it was made in the past when only ozone measurements were taken into account. The thresholds considered are those specified by NOM (2014a, b), i.e., 0.095 ppm for ozone and 75 μg/m 3 for PM 10 .

Data description
In Fig. 1 we have the figure of the Mexico City's ozone and PM 10 measurements during the observational period considered here. In that figure we also have two horizontal lines representing the valid environmental threshold for each pollutant.
Looking at Fig. 1 there are several episodes where high peaks have occurred in the PM 10 data, especially in the beginning of the observational period. We may also observe, as we move towards the end, that even though peaks occur, they are not as high as in the beginning. In the case of ozone, we also see higher concentrations in the beginning of the observational period. However, towards the end, we may see that the values vary around the line 0.095. Also, from the plots given in Fig. 1, we may see that perhaps the series of measures implemented by Mexico and, in particular, Mexico City since the late 1980s have contributed to the decrease of ozone and PM 10 levels, including the height of their peaks, in the metropolitan area.
If we analyse the maximum, minimum and mean values throughout the years, we have that the overall mean values were 86.92 μg/m 3 and 0.128 ppm in the cases of PM 10 and ozone, respectively. The corresponding minima were 14 μg/m 3 and 0.02 ppm. In the case of maximum values, we have 379 μg/m 3 and 0.35 ppm in the case of PM 10 and ozone, respectively. In Fig. 2 we may see the behaviours of these quantities. In this figure (top two plots), we have the yearly behaviours of the mean (thick continuous lines), maximum (thick dashed lines) and minimum (thin dashed lines) of the measurements.
Looking at Fig. 2 (top two plots) we see that rapid year by year decreases occur on the mean and maximum ozone measurements from 1995 to about the year 2008. That might be a consequence of the several measures implemented by the environmental authorities in order to decrease those levels. Most of those measures were implemented aiming to decrease ozone levels since this pollutant was the one still causing many of the environmental alerts in the city. Nevertheless, some of those measures also contributed to decrease PM 10 levels since they impact industries, cars, busses, taxis and trucks emissions. Hence, the decrease was also seen in the mean of the yearly PM 10 measurements. However, it seems that some peaks of emission of PM 10 still occur since the yearly maximum values do not decrease as much as the mean's. The number of yearly exceedances of the respective thresholds varies from 124 to 312 days in the case of PM 10 and from 178 to 341 in the case of ozone. The highest number of exceedances occurred in the years 1997 and 1996 in the cases of PM 10 and ozone, respectively, and the corresponding smallest values were in the years 2014 and 2012. In total we have 6659 and 5104 ozone and PM 10 exceedances of the corresponding thresholds. These values account for approximately 72.93% and 55.9% of the observed days, respectively.
In Fig. 2 we also have the plots (bottom two plots) of the number of exceedances year by year for both pollutants. Looking at these plots we may see that in the year 2000 we also had a high number of ozone exceedances. That, however, is smaller than the number in 1996 if only by one. Looking at the two bottom plots in Fig. 2, we see that the number of ozone as well as PM 10 exceedances also decrease rapidly until around 2008 and then vary around a lower value from that year on. That is also a reflex of the decrease in the measurements (maximum and mean) throughout the years shown in Fig. 1 and top two plots in Fig. 2. Hence, we see an improvement in the air quality related to ozone and PM 10 in later years of the observational period. However,

Programming details
In the model described in Sect. 2 we have the following specifications. The whole observational period has 9131 days for the total of N = 25 years of which 6 are leap and 19 are non-leap. The dimension of our vector is d = 2 corresponding to the two pollutants we are interested in studying. Hence, Z t = (Z 1,t , Z 2,t ), t = 1, 2, . . . , T = 366, where Z 1,t and Z 2,t correspond to the ozone and PM 10 measurements at day t, respectively, t = 1, 2, . . . , T . Hence, we have The thresholds considered are L 1 = 0.095 ppm and L 2 = 75 μg/m 3 in the case of ozone and PM 10 , respectively. In order to account for the possible dependence of few past days of maximum measurements, we are using M = 8. Therefore, the state space of the order K of the Markov chain X (K ) is S = {0, 1, 2, . . . , 8} and we have λ = 1 in (5).
The hyperparameters of the Dirichlet prior distributions are assigned in a similar manner as in Morales-Morillón (2018). Therefore, the values of α As it is possible to see in those works, occurrences and selection of orders from 0 to 17 were allowed depending on the data and threshold used. Therefore, even though the value of λ is one, it allows the occurrences of order larger than one whenever the data ask for it.
2. The hyperparameters of the Dirichlet prior distributions are chosen to reflect the occurrences of the events to which they are associated. For instance, in the case where K = 1, events that have not been observed, but are part of the space of possible events, will have a small probability assigned to it, whereas those that have been observed will have probabilities that absorb the contribution of both observed data and prior distributions.
3. Another option to assign the hyperparameters of the Dirichlet prior distributions is to proceed as in Álvarez et al. (2015), Álvarez and Rodrigues (2008), Rodrigues et al. (2015) where small values of the hyperparameters were assigned to rare events and larger values to those occurring more frequently. That, however, was not followed here because we wanted to see if not so specific prior distributions would produce good results as in Morales-Morillón (2018) where the model was used to study problems in genetics. Additionally, using different hyperparameters in the present case could imply in having to consider a large number of vectors of hyperparameters since we would have to take into account every possible transition at each time t.

Results
The selection of the order K of the Markov chain is made using the maximum a posteriori as well as the leave-one-out cross validation (see Appendix). Hence, in the case of the maximum a posteriori method we select the value of K that maximise its posterior distribution. The posterior probability P(K | Y), given by (6), is proportional to L(Y | K ) λ K /K ! and may be analysed using the value of log These values are shown in Table 1 where we also have the values of LOO for each K (first two columns).
It is possible to see by looking at Table 1 that, in the case of MAP, the highest probability is achieved when K = 1 followed by K = 2. Hence, the selected order is K = 1. When we use the LOO method to select the order of the chain, the order with the smallest value is selected. In the present case, that value is given by the case K = 0 closely followed by K = 1. For purpose of comparison, we have also obtained the estimated values of K using MAP and LOO for two additional vectors of hyperparameters of the Dirichlet prior distributions. In one of them all values of α appearing in those prior distributions are equal to one and in the other they are all equal to 1/2. These values are also given in Table 1. Looking at Table 1 we see that in all cases there is a consistency in the selection of the order using MAP (the selected order is always K = 1). When we use LOO, the chosen order also is K = 1 for values of α = 1/2, 1 in agreement with the choice made using MAP. Hence, the choice of K = 1 for the order of the chain. Once the order is selected, we generate a sample of the corresponding initial and transition probabilities using their associated posterior distributions. We may also obtain the probabilities P(W t =ĩ), t = 1, 2, . . . , T ,ĩ ∈ S 2 , applying (7) and using the generated samples and the recursive formula (8) since the chosen order is K = 1.
In Fig. 3 we have the estimated (blue lines) and the observed (dashed black lines) transition probabilities P (1) i j (t) as well as the 95% credible interval (green and maroon lines). Estimation of the probabilities was made using a sample of size 1000. (The values of these probabilities are given in the Supplementary Material file.) When we look at Fig. 3 it is possible to see that when the observed transition has value zero, i.e., either the transition was not actually observed or it was observed only a number of times that was negligible when compared to the number of observed years, the estimated posterior distribution gives a non-zero value to those transitions, oscillating between small and very small values. Hence, even though some transition events may not have been observed during the timeframe considered here, their posterior transition probabilities allow the possibility that they occur with small probabilities. That is useful because if a given transition has not happened during the observational period, it does not mean it will not happen in the future. Additionally, we may see that the largest value given to non-observable events is around 0.25. In the cases where we do have observed transitions, the estimated transition probabilities represent well the observed ones. We may also see from Fig. 3 that transitions from state 0 (i.e., no exceedances of either pollutant) to state 0 are higher when we consider the time from around the end of June to the beginning of December. This period corresponds to part of summer/autumn seasons. Summer season has rainy days most of the time, hence, we may have a low probability of both pollutants exceeding their corresponding threshold. Also, we see a high probability of going from a 0 state to a 1 state (corresponding to (1) (t)

Fig. 3 Estimated (blue lines) and observed (dashed black lines) posterior transition probabilities P
(1) mĩ (t), m ∈ S 2 ,ĩ ∈ S 2 , t = 1, 2, . . . , T − 1, as well as the 95% credible interval (green and maroon lines), when we have K = 1 a (1, 0) state, i.e., an exceedance of ozone only) during the same period when high values for P (1) 00 (t) occur. Hence, if an exceedance occurs in that period is more likely to be caused by high levels of ozone than PM 10 's. That may be corroborated when we look at the transition probabilities from 0 to 2 (corresponding to the state (0, 1) indicating an exceedance of PM 10 only). We may also see that if we are in the state 0, then there also exists a large probability of staying at that state in the beginning of the year and also during the summer/autumn seasons.
Other transitions that call attention because of their relatively higher values are from 1 ↔ (1, 0) to 1 and from 3 ↔ (1, 1) to 3. The high values of the former occur in a similar period of those from 0 to 1 only with higher values. The latter occurs practically the whole year. Hence, once we have exceedances of PM 10 during spring/summer/autumn, there is a large probability of continuing to have it. Similarly, once both pollutants exceeded their corresponding thresholds, they are bound to continue to exceed them.
Another interest here is in estimating the probability of occurrence of one of the states in S 2 throughout the year. The estimated values of P(W t =j),j ∈ S 2 , are given in the Supplementary Material files. In Fig. 4 we have the plots of the estimated (blue lines) and empirical/observed (black dashed lines) probabilities P(W t =ĩ),ĩ ∈ S 2 , t = 1, 2, . . . , T , with the 95% credible interval (green and maroon lines) when we have the selected order of the chain K = 1.
Looking at Fig. 4, it is possible to see that the fit of the estimated and observed probabilities is very good. Additionally, non-zero probabilities with small values were assigned to events that had observed probability zero or approximately zero as in the case of P(W t = 2) (left bottom plots) for the period ranging from mid April to mid October, and a small number of days in the case of P(W t = 0) (left top plots). In that way, even if events have not been observed up to now, the corresponding estimated Remark In addition to the results for the selected model (K = 1), we have decided to include in the appendix an analysis when K = 2 in order to illustrate how to proceed if the selected order of the chain is K ≥ 2 (see Section C of the Appendix).

Discussion and comments
In the present work a multivariate non-homogeneous Markov chain of order K ≥ 0 was used to jointly study occurrences of pollutants' exceedances. The model was applied, in its two-dimensional version, to ozone and PM 10 data from Mexico City. Estimation of the order K and the associated initial and transition probabilities is performed using the Bayesian point of view. The approach used to study this multi-dimensional problem was to transform the vector recording of whether or not exceedances have occurred into a unidimensional indicator. These will constitute the values of each element of the observed data. On these data the state spaces S (K ) 3 and S (K ) 4 associated to the Markov chain X (K ) are defined and, also, the corresponding initial and transition probabilities.
We would like to call attention to the fact that, even though the series of measurements are not stationary, the non-homogeneous Markov chain model is applied to the transformed series whose states depend on whether or not threshold exceedances have occurred. Hence, the non-stationarity is diluted a little. That can be seen in Fig. 5.
Note that some non stationarity is still present. However, as mentioned earlier, even though this may produce higher exceedance probabilities, this increasing effect is only used in order to alert about the possibility that higher levels of pollution may occur. Estimation of the order of the chain was performed using both the maximum a posteriori (MAP) and leave-one-out cross-validation (LOO) methods. The selected order was K = 1 in the case of MAP and K = 0 closely followed by K = 1 in the case of LOO. Hence, we take K = 1 as the selected order. The associated initial and transition probabilities were obtained by averaging the generated values obtained from the respective posterior distributions.
In order to see how well MAP and LOO estimate the order of the chain, we have performed a simulation study. Using the estimated initial and transition probabilities associated with the order K = 1, we have generated 100 realisations of the nonhomogeneous chain. We have considered the cases where N = 25, 50, 75, 100 and used MAP and LOO to estimate the order of the chain. Results are given in Table 2 in Section B of the Appendix. As we see both methods agree on the selection of the order K giving higher/lower values to K = 1 followed by K = 2. This supports the selected order in the application considered in the present work.
Using the estimated initial and transition probabilities, the values of the probabilities of having exceedances of one, both or none of the pollutants were obtained. That was made taking advantage of the Markov property of the model as well as a recursive formula using the generated initial and transition probabilities. Results show a good fit of the estimated values to the observed. In order to certify that, in Fig. 6 we have the plots of |P(W t =ĩ) −P(W t =ĩ)| for t = 1, 2, . . . , T ,ĩ ∈ S 2 , where P(W t =ĩ) and P(W t =ĩ) indicate the observed (empirical) and the estimated probabilities using the  Fig. 6 Absolute values of the differences between the observed and estimated P(W t =ĩ), t = 1, 2, . . . , T , i ∈ S 2 , in the cases K = 0 (blue lines), K = 1 (black lines) and K = 2 (red lines) generated values of the initial and transition probabilities obtained from their respective posterior distribution. We present the cases K = 0 (blue lines), K = 1 (black lines) and K = 2 (red lines) since the latter is the order with second largest probability, and the first case is the one showing, in general, the small differences.
Looking at Fig. 6 we may see that, when compared to the case K = 2, a better fit is achieved when we have K = 1. Even though the values of the differences in the case K = 2 are also small, we may see that when compared to the differences when K = 1, the latter presents a better overall fit. In the case of K = 0, the differences are even smaller than those in the case of K = 1. However, it is worth calling attention to the fact that when K = 0, the estimated posterior probabilities receive more weight from the data which are what rule the empirical (observed) probabilities, hence, the best adjustment. We may also see that in the case K = 2, the differences throughout the year are not large either. Larger misfits occurred for the periods ranging from June to October in the case of P(W t = 2), and from January to June and November to the end of the year in the case of P(W t = 3). However, even in these cases the differences are small.
If we take into account the information provided by the SAD, we have that the values, for K equal to zero, one and two are 3.75, 263.58 and 3393.404, respectively. Again, we may see that the case K = 1 outperforms the case K = 2. As before, the case K = 0 has the smallest SAD, but as already argued that is due to the form the empirical and estimated values of the corresponding probabilities are obtained.
Going back to Fig. 4, where the plots of the estimated P(W t =ĩ),ĩ ∈ S 2 , t = 1, 2, . . . , T , are presented, we may see that the probability of occurrence of states 0, 1 and 2 are mostly below 0.75 with values below 0.5 in the cases of states 0 and 2. Additionally, in the first half of the year the probabilities of occurrence of state 3 are larger than 0.25 and most of the time are in the interval (0.5; 0.75) which, compared to the probabilities of occurrences of the remaining states, are very high. Hence, the model and the data corroborate the fact that during winter and spring there are larger probabilities of both ozone and PM 10 exceedances due to the dry season and clear skies. Additionally, in March/April there are large numbers of forest fires in the south of Mexico, with smoke travelling to Mexico City. Hence, it is likely to have occurrences of PM 10 exceedances during that period. We may also see that exceedances of both ozone and PM 10 , i.e., the occurrence of state 3, are fewer during the rainy season and especially during the period ranging from the beginning of July to the end of October. Very low estimated and practically zero observed probabilities are also noticed when we consider state 2 (meaning an exceedance of PM 10 only) during the period ranging from the beginning of April until mid-October. Additionally, in the case of the occurrences of state 1 (exceedances of ozone only) the probabilities are above 0.5 for the period ranging from mid-June to the beginning of September. In contrast, the probabilities of having exceedances of both pollutants (state 3) are low (around 0.25) for a similar period. Additionally, we may see that for the period ranging from 01 January to approximately the beginning of June, while the probabilities of having no exceedances, exceedances of either only ozone or PM 10 are below 0.3 for most of the time, in the case of exceedances of both pollutants these probabilities are always above 0.5. Hence, is more likely to have exceedances of both pollutants at the same time than having either only of one of them or none at all.
The results presented here may be used to obtain information about the occurrences of exceedances of none, one or more pollutants in given days of the year. That is made using the estimated P(W t =ĩ),ĩ ∈ S 2 , t = 1, 2, . . . , T . We may do that as follows. If we want to know the probability of having both ozone and PM 10 exceedances at a given time t, t = 1, 2, . . . , T , we just need to look at the values of P(W t = 3). For instance, in the case of t = 90, which would correspond to the 31st of March in a non-leap year and the 30th of March in a leap year, we have that this probability is approximately 0.777 whereas the probability of having an exceedance of ozone only is P(W 90 = 1) ≈ 0.127, and of PM 10 only is P(W 90 = 2) ≈ 0.009 (see tables in the Supplementary Material file). We also have that P(W 90 = 0) ≈ 0.086.

Remark
Note that the estimated initial and transition probabilities reflect averages behaviours. Also, note that even though there might be differences in the values of these probabilities if the day, say 30 March, falls in a weekday or in a weekend, it will fall on a weekend only at intervals of about six years, hence, the expected value of the corresponding probability will appear just 1/6 of the times with little effect on the total average. Another quantity that may be obtained using the model is, as in Rodrigues et al. (2015), the probability of having the occurrence of exceedances of either none, one or both pollutants in a given time when we have information of a few days earlier.
For instance, assume we want to know the probability that at time t = 100 (a day at the beginning of April) we have exceedances of both pollutants, i.e., an occurrence of state 3, and that in the three previous days we have the following sequence of states: 2, 1, 2 (i.e., exceedances of PM 10 only, ozone only and, again, PM 10 only). Therefore, we want to know the probability of the following sequence of states, Hence, we have, Hence, the probability of having either no exceedances or exceedances of only one of the pollutants is of order ten times lower than having exceedances of both pollutants. Recall that March/April is spring and there are many forest fires and also clear skies.
If we want to know the probability of having an ozone exceedance on the day t = 250, then we just need to use the values of P(W 250 = 1) and P(W 250 = 3) since the event { an ozone exceedance } may be written as the union of the events { an exceedance of ozone only } and { exceedances of both ozone and PM 10 }, i.e., { an ozone exceedance } = { occurrence of (1, 0) } ∪ { occurrence of (1, 1)}. Therefore, P({an ozone exceedance at t = 250}) = P(W 250 = 1) + P(W 250 = 3) ≈ 0.509 + 0.126 = 0.635.
(Again, the particular values of the absolute probabilities may be found in the Supplementary Material file.) Using the model described here we may also jointly study more pollutants. For instance, we could also include PM 2.5 , nitrogen dioxide, carbon monoxide, among others.
As in Álvarez et al. (2005) and Rodrigues et al. (2015) we could also consider for each pollutant several threshold values and with that estimate the probability that their levels belong to a particular concentration interval. For instance, if we take the values to declare the several levels of emergency alerts in Mexico City related to ozone and PM 10 , i.e., if L = 354 μg/m 3 , which are well above the Mexican environmental thresholds for these pollutants. Hence, we could calculate the probability of having the ozone concentration between 0.095 and 0.154 while having the PM 10 concentration in the interval (214, 354), i.e., the probability of having the ozone concentration above the Mexican environmental threshold, but below the ozone threshold used to declare Phase I of the environmental alert with the PM 10 's indicating that a Phase I could be declared. This would provide information on the possibility of having an environmental alert in the city. Based on this information, environmental authorities could take actions in order to avoid the occurrence of this alert.
One fact that is clear is that the estimated probabilities are not smooth functions of time. However, it is also clear that the empirical (observed) probabilities which are being represented are not smooth either. This could change if more years were considered in the analysis. This would increase the chance that under represented events occur. Just for the sake of comparison we have decided to check if with the 100 simulated chains of our simulated study, smother function would be obtained, for instance, for P(W t =ĩ),ĩ ∈ S 2 , t = 1, 2, . . . , T . However, that did not happen. Nevertheless, in spite of having only twenty five years of observations the estimated results represent well the observed. If smoothness of the functions is of interest, another option to produce smoother curves to model the observed probabilities and hence the probabilities P(W t =ĩ),ĩ ∈ S 2 , t = 1, 2, . . . , T is to use cubic splines.
One may also argue what the gain is in using this type of model when it seems that the empirical results give all the information we need. One positive side is that in the cases where transition events which did not happen during the observational period, but are part of the sample space, would be allowed to be generated in case one wanted to simulate string of events and also if the interest resides in informing the possibility of having such transitions occurring (even though the probabilities might be overestimated).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Leave-one-out cross-validation(LOO)
The LOO method in the case of homogeneous chain is described as follows (Chang 2019). Consider N realisation of a Markov chain with state space S = {1, 2, . . . , J } and let N x = (N x1 , N x2 , . . . , N x J ) be the vector whose ith coordinate N xi records the total number of times in all N realisations of the chains that we have a state x followed by a state i; i, x ∈ S, i.e., xi is the total number of transitions from x to i in the jth chain. Let α = (α 1 , α 2 , . . . , α J ) be the vector of hyperparameters of the Dirichlet prior distributions of the rows of the transition matrix associated with each Markov chain. For a vector x = (x 1 , x 2 xm . In order to describe the LOO in the non-homogeneous case we split into cases according to the order of the chain. The description is made taking into account the context and state spaces of the present work. The selected model (order) is the one producing the lowest LOO.  (1) is as when we have K = 0 and n (i) ml (t) = 1 if at the tth time in the ith chain we have a transition from m to l and it is zero otherwise.
(1) = 1 if the initial state of the ith chain of order K is m and it is zero otherwise and n (i,K ) ml (t) = 1 if the transition at time t in the ith chain is from m to l and it is zero otherwise.

Appendix B. Simulation results
In this section we present the results obtained using the methodology applied in the present context to the case of simulated non-homogeneous Markov chains. Using the estimated transition probabilities in the case of K = 1 we have generated 100 realisations of the non-homogeneous Markov chain. The choice of the order one is made for simplicity. Then, we have applied the MAP and LOO to see if the selected order is the correct one. Results are given in Table 2 as follows.
When we look at the results given in Table 2, we see that the estimated order using either MAP or LOO are the same. Additionally, we have that even with 25 chains (which is the case in the present application) the order is well estimated using either selection method.
Note that one may also generate the initial and transition probabilities using the selected Dirichlet distributions. However, since the present case we already have that in the application we have used those values.
After that, the estimated absolute probabilitiesP(W t =ĩ),ĩ ∈ S 2 , t = 1, 2, . . . , T , are obtained using (7). In Fig. 7 we have these estimated (blue lines) probabilities as well as the observed (dashed black lines) and the 95% credible intervals (green and maroon lines). The numerical values of these probabilities may be obtained upon request.