Márton Ispány
email: ispany.marton@inf.unideb.hu

Valdério A Reisen

Glaura C Franco

Pascal Bondon

Higor H A Cotta

Paulo R P Filho

Faradiba S Serpa

On generalized additive models with dependent time series covariates

Keywords: air pollution, generalized additive model, multicollinearity, principal component analysis, time series, vector autoregressive model

 which is the combination of GAM with the principal component analysis (PCA) and the vector autoregressive (VAR) model. In this paper, some properties of the GAM-PCA-VAR model are discussed theoretically and verified by simulation. A real data set is also analysed with the aim to describe the association between respiratory disease and air pollution concentrations.

Introduction

In the recent literature of time series, there has been an outstanding growth in models proposed for data that do not satisfy the Gaussian assumption. This is mainly the case when the response variable under study is a count series or an integer valued series. Procedures developed to analyse this kind of data comprises, for example, observation driven models, see [START_REF] Benjamin | Generalized autoregressive moving average models[END_REF] and [START_REF] Davis | Observation-driven models for Poisson counts[END_REF], integer valued autoregressive (INAR) processes, see [START_REF] Al-Osh | First-order integer valued autoregressive (INAR(1)) process[END_REF] and [START_REF] Barczy | Additive outliers in INAR(1) models[END_REF], or non-Gaussian state space models, see [START_REF] Durbin | Time series analysis of non-Gaussian observations based on state space models from both classical and Bayesian perspectives[END_REF] and [START_REF] Gamerman | A non-Gaussian family of statespace models with exact marginal likelihood[END_REF].

Particularly in health and environmental studies, where the response variable is typically a count time series, the generalized additive model (GAM) has been widely used to associate the dependent series, such as the number of respiratory or cardiovascular diseases to some pollutant or climate variables, see, for example, [START_REF] Chen | Ambient air pollution and hospital admission in Shanghai, China[END_REF], [START_REF] Nascimento | Association between the concentration of fine particles in the atmosphere and acute respiratory diseases in children[END_REF], [START_REF] Ostro | Air pollution and health effects: A study of medical visits among children in Santiago, Chile. Environ[END_REF], [START_REF] Schwartz | Harvesting and long term exposure effects in the relationship between air pollution and mortality[END_REF], [START_REF] De Souza | Generalized additive models with principal component analysis: an application to time series of respiratory disease and air pollution data[END_REF] and [START_REF] Souza | Principal components and generalized linear modeling in the correlation between hospital admissions and air pollution[END_REF] among others. Therefore, in general, the researches related to the study of the association between pollution and adverse health effects usually consider only one pollutant. This simple model choice may be due to the fact that the pollutants are linearly time correlated variables, see the discussion and references in the recent paper [START_REF] De Souza | Generalized additive models with principal component analysis: an application to time series of respiratory disease and air pollution data[END_REF].

Recently, it has become common practice to use principal component analysis (PCA) in regression models to reduce the dimensionality of an independent set of data, especially the pollutants, which in some instances can include a large number of variables. The PCA is highly indicated to this purpose, as it can handle the multicollinearity problem that can cause biased regression estimates, see, for example, [START_REF] Wang | Analyzing the effects of air pollution and mortality by generalized additive models with robust principal components[END_REF].

Nevertheless, use of PCA in the time series context can bring some misspecifications in the fit of the GAM model, as this technique requires that the data should be independent. This problem arises due to the fact that the principal components are linear combinations of the variables. In this context, as the covariates are time series, the autocorrelation present in the observations are promptly transferred to the principal components, see [START_REF] Zamprogno | PCA in time series with short and long-memory time series[END_REF].

One solution to this issue was recently proposed by [START_REF] De Souza | Generalized additive models with principal component analysis: an application to time series of respiratory disease and air pollution data[END_REF], see, also, [START_REF] Souza | Principal components and generalized linear modeling in the correlation between hospital admissions and air pollution[END_REF], who introduced a model which combines GAM, PCA and the vector autoregressive (VAR) process. The authors suggest to apply the VAR model to the covariates, in order to eliminate the serial correlation and produce white noise processes, which in turn will be used to build the principal components in the PCA. The new variables obtained in the PCA are finally used as covariates in the GAM model, originating the so called GAM-PCA-VAR model. In their work, the authors have focused on presenting the model and showing its superiority compared to the sole use of GAM or the GAM-PCA procedures, but have not deepened on the theoretical properties of the model.

Thus, to cover this gap, this work aims to state and prove some properties of the GAM-PCA-VAR model, as well as to perform some simulation study to check the results for small samples.

The paper is organized as follows. Section 2 presents the main statistical model, GAM-PCA-VAR, addressed here and its related models as GAM, PCA and VAR, in some detail. In Section 3 the theoretical results are proved for the main model. Section 4 discusses the simulation results and Section 5 is devoted to the analysis of a real data set. Section 6 concludes the work.

The GAM-PCA-VAR model

The generalized additive model (GAM), see [START_REF] Hastie | Generalized Additive Models[END_REF] and [START_REF] Wood | Generalized Additive Models: An Introduction with R. 2nd edn[END_REF], with a Poisson marginal distribution is typically used to relate a non-negative integer valued response variable Y with a set of covariates or explanatory variables X 1 , . . . , X p . In GAM the expected value µ = E(Y ) of the response variable depends on the covariates via the formula

g(µ) = β 0 + p i=1 f i (X i ),
where g denotes the link function, β 0 is the intercept parameter and f i 's are functions with a specified parametric form, e.g., they are linear functions f i (x) = β i x, β i ∈ R, i = 1, . . . , p, or non-parametric, e.g., they are simple smoothing functions like splines or moving averages. The unknown parameters β 0 and f i , i = 1, . . . , p can be estimated by various algorithms, e.g., backfitting or restricted maximum likelihood (REML) method. However, if the data observed for variables Y and X i , i = 1, . . . , p, form a time series the observations cannot be considered as a result of independent experiments and the covariates present strong interdependence, e.g., multicollinearity or concurvity, the standard fitting methods result in remarkable bias, see, e.g., [START_REF] Dionisio | A simulation study to quantify the impacts of exposure measurement error on air pollution health risk estimates in copollutant time-series models[END_REF] and [START_REF] De Souza | Generalized additive models with principal component analysis: an application to time series of respiratory disease and air pollution data[END_REF].

Let {Y t } ≡ {Y t } t∈Z be a count time series, i.e., it is composed of non-negative integer valued random variables. We suppose that the explanatory variables form a zero-mean stationary vector time series {X t } ≡ {X t } t∈Z of dimension p, i.e., X t = (X 1t , . . . , X pt ) where denotes the transpose, with the covariance matrix Σ X = E(X t X t ). Let F t denote the σ-algebra which contains the available information up to time t for all t ∈ Z from the point of view of the response variable, e.g., X t is F t-1 -measurable. The GAM-PCA-VAR model is introduced in [START_REF] De Souza | Generalized additive models with principal component analysis: an application to time series of respiratory disease and air pollution data[END_REF] as a probabilistic latent variable model. In this paper, we define this model in a more general form as

Y t | F t-1 ∼ Poi(µ t ), (1)
X t = ΦX t-1 + AZ t (2)
with link

g(µ t ) = β 0 + p i=1 ∞ j=0 f ij (Z i(t-j) ), (3)
where Poi(•) denotes the Poisson distribution, the latent variables {Z t }, Z t = (Z 1t , . . . , Z pt ) , form a zero-mean Gaussian vector white noise process of dimension p with diagonal variance matrix Λ = diag{λ 1 , . . . , λ p }, where λ 1 ≥ λ 2 ≥ . . . ≥ λ p , A is an orthogonal matrix of dimension p×p, Φ is a matrix of dimension p × p, g is a known link function, β 0 denotes the intercept, and f ij 's are unknown functions. For a zero-mean Gaussian vector white noise process {Z t } with covariance matrix Σ we shall use the notation {Z t } ∼ GWN(Σ). See also [START_REF] Brockwell | Time Series: Theory and Methods[END_REF]Definition 11.1.2]. Clearly, for all i, the univariate time series {Z it } ∼ GWN(λ i ), and {Z it } is mutually independent from {Z jt } for all j = i. We assume that all the eigenvalues of Φ are less than 1 in modulus which implies that equation (2) has a unique stationary causal solution. In the case of a Poisson distributed response variable the two widely used link functions are the identity link, g(z) = z, and the canonical logarithmic link, g(z) = log z. The set (β 0 , {f ij }, A, Λ, Φ) forms the parameters of the GAM-PCA-VAR model to be estimated. We remark that in the case of canonical logarithmic link function no additional assumption is needed for the parameters, while in the case of identity link function all the parameters in equation (3), i.e., β 0 and f ij 's, have to be non-negative. It should be also emphasized that the underlying intensity process {µ t } of {Y t } is also a time series with a complex dependence structure, and µ t is F 2), the latent variable can be expressed as Z t = A U t , where U t := X t -ΦX t-1 for all t. Thus, as the first step, the intermediate vector times series {U t } is derived from filtering {X t } by a VAR(1) filter. One can see that {U t } ∼ GWN(Σ U ) where Σ U := AΛA . Then, as the second step, the latent vector time series {Z t } as principal component (PC) vector is derived by principal component transformation of the intermediate vector white noise {U t }. The transformation matrix of the PCA is given by the spectral decomposition of Σ U . Finally, as the third step, the standard GAM with link (3) is fitting for the response time series {Y t } using the latent vector time series {Z t }. The impact of the VAR(1) filter in the first step is to eliminate the serial correlation present in the original covariates. On the other hand, the impact of the PCA in the second step is to eliminate the correlation in the state space of the original covariates. Hence, the result of these two consecutive transformations is the latent vector time series {Z t } whose components, Z it , i = 1, . . . , p, t ∈ Z, are independent Gaussian variables both in space and time. In the case of logarithmic link function, large positive values in a coordinate of the latent variable indicate locally high influence according to this latent factor. On the contrary, large negative values indicate negligible influence on the response, see, for example, [START_REF] Zamprogno | PCA in time series with short and long-memory time series[END_REF]. The order of models in the acronym GAM-PCA-VAR corresponds to these steps starting with the third one and finishing with the first one.

The GAM-PCA-VAR model contains several submodels with particular dependence structure. If Φ = 0 then the VAR equation ( 2) is simplified to a principal component transformation. In this case, we suppose that there is no serial correlation and we only have to handle the correlation in the state space of covariates. We have two transformation steps: PCA and GAM. This kind of models is called GAM-PCA model that is intensively studied nowadays, see, e.g., [START_REF] Roberts | Using supervised principal components analysis to assess multiple pollutant effects[END_REF] and [START_REF] Zhao | A comparison between two GAM models in quantifying relationships of environmental variables with fish richness and diversity indices[END_REF]. Beside the full PCA when all PCs are involved into the GAM, we can fit a restricted PCA model by defining f ij = 0 for all i > r and j ≥ 0 where r < p. In this case, the first rth PCs are applied as covariates in the GAM step. If the matrices in VAR(1) model ( 2) have the following block structures

Φ = Φ q 0 0 0 , A = A q 0 0 I p-q ,
where the eigenvalues of the q × q matrix Φ q are less than one in modulus, A q is an orthogonal matrix of dimension q×q (q ≤ p), and f i1 (z) = β i z with β i ∈ R for i = 1, . . . , r (r ≤ q), f i1 is a general smoothing function for i = q + 1, . . . , p, f ij = 0 otherwise, then we obtain the model that was studied in [START_REF] De Souza | Generalized additive models with principal component analysis: an application to time series of respiratory disease and air pollution data[END_REF] and applied in the data analysis of Section 5. In this model it is supposed that the set of covariates can be partitioned into two sets: (X 1 , . . . , X q ) are normal covariates, e.g., the pollutant variables in the terminology of Section 5, while (X q+1 , . . . , X p ) are so-called confounding variables as trend, seasonality, etc. The normal covariates satisfy a q-dimensional VAR(1) model, however, instead of the all coordinates of the innovation, only its first rth PCs are involved into the GAM taking into consideration that the covariates present strong inter-correlation. Finally, we note that our model can be further generalized by replacing equation ( 2) by the more general VARMA or VARIMA or their seasonal variants (SVARMA or SVARIMA) models.

Since the latent variables {Z t } form a Gaussian vector time series, given a sample (X 1 , Y 1 ), . . . , (X n , Y n ), the log-likelihood can be expressed in an explicit form, see [START_REF] De Souza | Generalized additive models with principal component analysis: an application to time series of respiratory disease and air pollution data[END_REF] for a particular case. Because this log-likelihood is rather complicated a three-stage estimation method is proposed. Firstly, VAR(1) model is fitted to the original covariates by applying standard time series techniques. Secondly, PCA is applied for the residuals defined by Z t = X t -ΦX t-1 , t = 2, . . . , n, where Φ denotes the estimated autoregressive coefficient matrix in the fitted VAR(1) model. Thirdly, GAM model is fitted using the PCs. The approach discussed above is similar to the principal component regression, see, e.g., [START_REF] Jolliffe | Principal Component Analysis[END_REF]Chapter 8], and it can be considered as a three-stage non-linear regression method.

The first two steps of the above proposed parameter estimation method for GAM-PCA-VAR model can be interpreted as consecutive orthogonalizations, firstly in time and then in the state space of covariates. In [START_REF] De Souza | Generalized additive models with principal component analysis: an application to time series of respiratory disease and air pollution data[END_REF]Remark] we argued that the order of VAR filter and PCA can not be interchanged because the orthogonalization in the state space does not eliminate the serial correlation and, as the necessary next step, the orthogonalization in time by VAR filter bring back the inter-correlation between the covariates. In what follows, we demonstrate this phenomena by giving a simple example. Let {X t } be a zero-mean causal VAR(1) process defined by

X t = Ψ X t-1 + W t ,
where {W t } is a zero-mean vector white noise process with variance matrix Σ W . Suppose that the variance matrix Σ X of {X t } is diagonal, i.e., the coordinates of {X t } can be interpreted as PCs after PCA. Then Σ W is not necessarily a diagonal matrix, which implies that a VAR(1) filter may result in an intercorrelated white noise. Namely, consider the following parameters Σ W = AΛA and Ψ = ASA , where Λ and S are diagonal matrices and A is an orthogonal matrix. In other words, we suppose that the orthogonal matrix A in the spectral decomposition of Σ W diagonalizes the autoregressive coefficient matrix as well.

Then, we have, by formula (11.1.13) in [START_REF] Brockwell | Time Series: Theory and Methods[END_REF], that

Σ X = ∞ j=0 Ψ j Σ W (Ψ ) j = ∞ j=0 AS j ΛS j A = Adiag λ i 1 -s 2 i A .
Let σ 2 > max i {λ i } arbitrary and define s i := 1 -λ i /σ 2 for all i. Clearly, Ψ is a causal matrix since all its eigenvalues are less than 1 in modulus and Σ X = σ 2 I, i.e., the coordinates of {X t } are uncorrelated. However, the innovation variance matrix Σ W can be arbitrary proving that the application of VAR filter for a nonintercorrelated vector time series can give inter-correlated vector white noise in its coordinates. Now, we present some particular examples of GAM-PCA-VAR models.

Example 1. One of the simplest GAM-PCA-VAR models is the model with dimension p = 1 and log-linear link function. In this case, there is only one covariate {X t }, and the VAR equation ( 2) is an AR(1) model

X t = φX t-1 + Z t , (4)
where |φ| < 1 which guarantees the existence of a unique stationary causal solution, {Z t } ∼ GWN(λ), λ > 0. We remark that A = 1 in equation ( 2) in order for the model to be identifiable. The link is log-linear expressed as

log µ t = β 0 + β 1 Z t . (5)
The parameter set of this model is (β 0 , β 1 , λ, φ) with parameter space R 2 × R + × (-1, 1). In this model, there is no dimension reduction. Clearly, Z t = X t -φX t-1 , thus the response depends on the covariate through the link

log µ t = γ 0 + γ 1 X t + γ 2 X t-1 , (6)
where there is a one-to-one correspondence between the parameter sets (β 0 , β 1 , φ) and (γ 0 , γ 1 , γ 2 ) defined by the equations γ 0 = β 0 , γ 1 = β 1 and γ 2 = -φβ 1 provided φ = 0. However, if we fit the standard GAM by using the link (6) with covariates X t and X t-1 at time t, we take no count of the interdependence in time series {X t } which can result in biased and inconsistent estimators of the GAM parameters.

Example 2. Define a particular two-dimensional (p = 2) GAM-PCA-VAR model with logarithmic link function in the following way. The two-dimensional covariate vector process {X t }, X t = (X 1t , X 2t ) , satisfies the VAR(1) model

X 1t X 2t = φ 1 0 0 φ 2 X 1(t-1) X 2(t-1) + cos ϕ -sin ϕ sin ϕ cos ϕ Z 1t Z 2t ,
where

|φ 1 | < 1, |φ 2 | < 1 and {Z it } ∼ GWN(λ i ) with λ i > 0, i = 1, 2
, which are independent from each other. Note that the set of two-dimensional orthogonal matrices, A, can be parametrized by an angle parameter ϕ ∈ [0, 2π). We assume that the link is log

µ t = β 0 + β 1 Z 1t .
The parameter set of this model is (

β 0 , β 1 , ϕ, λ 1 , λ 2 , φ 1 , φ 2 ) and the parameter space is R 2 × [0, 2π) × R 2 + × (-1, 1) 2 .
Note that, in this model, there is a PCA step as a dimension reduction since only the first coordinate {Z 1t } of the vector innovation is involved into the GAM as covariate. One can see that the response depends on the covariates through the link

log µ t = γ 0 + γ 1 X 1t + γ 2 X 2t + γ 3 X 1(t-1) + γ 4 X 2(t-1) ,
where γ 0 = β 0 , γ 1 = β 1 cos ϕ, γ 2 = β 1 sin ϕ, γ 3 = -β 1 φ 1 cos ϕ and γ 4 = -β 1 φ 2 sin ϕ. Thus, the intensity process {µ t } depends on all coordinates of X t and X t-1 . Clearly, there is a one-to-one correspondence between the two parameter sets (β 0 , β 1 , ϕ, φ 1 , φ 2 ) and (γ 0 , γ 1 , γ 2 , γ 3 , γ 4 ).

Example 3. A seasonal one-dimensional GAM-PCA-VAR model with linear link function can be defined in the following way. Suppose that the one-dimensional covariate process {X t } satisfies the SAR s (1) model:

X t = φX t-s + Z t ,
where |φ| < 1, {Z t } ∼ GWN(λ) with λ > 0 and s ∈ Z + denotes the seasonal period. The link is linear and is given by

µ t = β 0 + β 1 f (Z t ),
where f : R → R + is a known function and β 0 , β 1 ∈ R + are parameters. The parameter set of this model is (β 0 , β 1 , λ, φ) with parameter space R 3 + × (-1, 1). The response variable depends on the original covariates through the link

µ t = β 0 + β 1 f (X t -φX t-s ).
If the function f is sufficiently smooth we have by approximation f (X t -φX t-s ) ≈ f (X t ) -φf (X t )X t-s , and then

µ t = γ 0 + γ 1 f 1 (X t ) + γ 2 f 2 (X t , X t-s ), (7)
where f 1 , f 2 are known functions and γ 0 = β 0 , γ 1 = β 1 and γ 2 = -β 1 φ. Thus, the response depends on the original covariate and its s-step lagged series through the standard GAM. However, the covariates in equation ( 7) are clearly dependent.

Theoretical results

In this section, we prove some theoretical results for particular classes of GAM-PCA-VAR models. Consider the log-linear model defined by the link

log µ t = β 0 + p i=1 ∞ j=0 β ij Z i(t-j) , (8)
where β 0 , β ij ∈ R, i = 1, . . . , p, j ∈ Z + . The first proposition is about the existence of log-linear GAM-PCA-VAR models.

Proposition 1. Suppose that

σ 2 := p i=1 λ i ∞ j=0 β 2 ij is finite.
Then the GAM-PCA-VAR model with log-linear link (8) has solution {(Y t , X t )} which is a strictly stationary process and E(Y t ) = E(µ t ) = exp(β 0 + σ 2 /2) for all t ∈ Z.

Proof. By conditioning we have that

E(Y t ) = E(E(Y t | F t-1 )) = E(µ t ) = E(exp(log µ t )) = exp(β 0 + σ 2 /2) (9)
is finite since, by equation ( 8), log µ t ∼ N (β 0 , σ

(X t ) = X t-1 , B(Y t ) = Y t-1 for all t ∈ Z.
In the next proposition, we prove that all moments of the log-linear GAM-PCA-VAR model are finite.

Proposition 2. Suppose that σ 2 defined in Proposition 1 is finite. Then all moments of the stochastic process {(Y t , X t )} are finite. In particular, we have, for all t ∈ Z,

Var(Y t ) = exp(2β 0 + σ 2 )(exp(σ 2 ) -1 + exp(-β 0 -σ 2 /2)),
Var(µ t ) = exp(2β 0 + σ 2 )(exp(σ 2 ) -1).

Proof. Let r ∈ N. Define the rth factorial of a non-negative integer k as k [r] := k(k -1) • • • (k -r + 1) and let k [0] := 1. For the rth factorial moment of Y t we have by conditioning that

E(Y [r] t ) = ∞ k=0 k [r] P(Y t = k) = E ∞ k=0 k [r] P(Y t = k | F t-1 ) =E ∞ k=r µ k t (k -r)! e -µt = E(µ r t )
for all t ∈ Z. Similarly to (9), we have that the factorial moments are finite, since

E(Y [r] t ) = E(µ r t ) = E(exp(r log µ t )) = exp{β 0 r + (σr) 2 /2}. (10)
Since the higher order moments can be expressed by the factorial moment via the formula

E(Y r ) = r j=0 S(r, j)E(Y [j] ),
where S(r, j)'s denotes Stirling numbers of the second kind, the finiteness of all higher order moments follows easily. Since {X t } is a Gaussian process all its moments are finite. Finally, the existence of mixed moments follows by the Cauchy-Schwarz inequality.

From Equation (10), we have

Var(µ t ) =E(µ 2 t ) -E 2 (µ t ) = exp(2β 0 + (2σ) 2 /2) -exp(2β 0 + σ 2 ) = exp(2β 0 + σ 2 )(exp(σ 2 ) -1).
Finally, the formula for Var(Y t ) can be derived by

Var(Y t ) = E(Var(Y t | F t-1 )) + Var(E(Y t | F t-1 )) = E(µ t ) + Var(µ t ).
The existence of all moments for the log-linear GAM-PCA-VAR process is to be compared with the same result for the integer valued GARCH, so-called INGARCH, process, see [START_REF] Ferland | Integer-valued GARCH process[END_REF]Proposition 6]. This implies that the log-linear GAM-PCA-VAR process possesses second and higher order structures, e.g., the autocorrelation function, the spectral density function, the cumulants and the higher order spectra exist. Let ρ Y denotes the autocorrelation function of the time series {Y t }. Proposition 3. For the auto-and cross-correlation functions of the GAM-PCA-VAR process {(Y t , X t )} with intensity process {µ t }, we have ρ

Y (h) = c Y ρ(h), ρ µ (h) = c µ ρ(h) and ρ Y µ (h) = c Y µ ρ(h) where ρ(h) := exp   p i=1 λ i ∞ j=0 β i(j+|h|) β ij   -1, h ∈ Z \ {0},
and the constants c Y , c µ , c Y µ are defined by

c Y := (exp(σ 2 )-1+exp(-β 0 -σ 2 /2)) -1 , c µ := (exp(σ 2 )-1) -1 , c Y µ := √ c Y c µ . Moreover, Cov(Y t+h , X t ) = Cov(µ t+h , X t ) = E(Y t+h X t ) = E(µ t+h X t ) = C(h) with C(h) := exp(β 0 + σ 2 /2) × ∞ k=0 Φ k A(λ • β h+k ) if h ≥ 0, ∞ k=0 Φ k-h A(λ • β k ) if h ≤ 0, (11)
where λ := (λ 1 , . . . , λ p ) , β j := (β 1j , . . . , β pj ) , j ∈ Z + , and • denotes the entrywise (Hadamard) product.

Proof. Let h ∈ N. One can see that for the intensity process we have µ t+h = µ

(1)

th µ (2)
th where log µ

(1)

th := β 0 + p i=1 h j=1 β i(h-j) Z i(t+j) , log µ (2) th := p i=1 ∞ j=0 β i(j+h) Z i(t-j) .
Clearly, µ

th is independent of F t-1 and Y t , while µ

(2)

th is F t-1 -measurable. Hence, we have by conditioning that

E(Y t+h Y t ) =E(Y t E(Y t+h | F t+h-1 )) = E(µ t+h Y t ) = E(µ (1) th µ (2) th Y t ) =E(µ (1) th )E(µ (2) th E(Y t | F t-1 )) = E(µ (1) th )E(µ (2) th µ t ) = E(µ t+h µ t ) since µ t is independent of µ (1)
th . This gives the result for h > 0. On the other hand, for all h > 0, again by conditioning,

E(Y t+h µ t ) = E(µ t+h µ t ). Thus Cov(Y t+h , Y t ) = Cov(µ t+h , µ t ) = Cov(Y t+h , µ t ), h ∈ Z \ {0}. Since E(µ t+h µ t ) = E(µ (1)
th µ

(2)

th µ t ) = E(µ (1)
th )E(µ

(2)
th µ t ) similarly to equation [START_REF] Ferland | Integer-valued GARCH process[END_REF] we have

E(µ t+h µ t ) = exp   2β 0 + 1 2 p i=1 λ i   h-1 j=0 β 2 ij + ∞ j=0 (β i(j+h) + β ij ) 2     = exp   p i=1 λ i ∞ j=0 β i(j+h) β ij   E(µ t+h )E(µ t ).
Thus, the first part of the proposition follows by Proposition 2.

Next we prove the formula [START_REF] Hastie | Generalized Additive Models[END_REF] for the cross-correlations of response and covariate variables. Clearly, by conditioning, E(Y t+h X t ) = E(µ t+h X t ) for all h ∈ Z + . On the other hand, for all t ∈ Z, h ∈ Z + , we have X t+h = X

th := h k=1 Φ h-k AZ t+k , X (2)
th := ∞ k=0 Φ h+k AZ t-k .
One can see that X

(1)

th is independent of F t-1 and Y t , while X

th is F t-1measurable. Thus, we have that

E(X t+h Y t ) =E((X (1) th + X (2) th )Y t ) = E(X (1) th )E(Y t ) + E(X (2) th E(Y t | F t-1 )) =E(X (1) th )E(µ t ) + E(X (2) th µ t ) = E(X t+h µ t ). Hence E(Y t+h X t ) = E(µ t+h X t )
for all h ∈ Z and it is enough to compute the cross-correlation between {X t } and {µ t }. Let h ≥ 0. For all ∈ {1, . . . , p}, k ∈ Z + let I h k := {1, . . . , p} × Z + \ ( , k + h) and define the random variables

log ξ th k := β 0 + (i,j)∈I h k β ij Z i(t+h-j) , log η th k := β (k+h) Z (t-k) .
Then µ t+h = ξ th k η th k , where the factors in this decompostion are independent. Since E(µ t+h X t ) = ∞ k=0 Φ k AE(µ t+h Z t-k ) and, using the fact that for Z ∼ N (0, λ) and β ∈ R we have E(Z exp(βZ)) = βλ exp(λβ 2 /2),

E(µ t+h Z (t-k) ) = E(ξ th k η th k Z (t-k) ) = E(ξ th k )E(η th k Z (t-k) ) = E(µ t+h )β (k+h) λ ,
we obtain the formula [START_REF] Hastie | Generalized Additive Models[END_REF]. The proof is similar in the case of h < 0.

Remark 1. It is easy to see that if β ij = β j i for all i, j, then the function ρ is given by ρ(h) = exp(

p i=1 λ i β |h| i /(1 -β 2 i )) -1, h ∈ Z.
If β i 's are all positive then ρ is positive everywhere and we have autocorrelation functions which are similar to what is displayed in Figure 1. For the one-dimensional model in Example 1 we have the cross-correlation function (CCF) C(h) = exp(β 0 + λβ 2 1 /2)λβ 1 φ -h for h ≤ 0 and C(h) = 0 for h > 0. If φ > 0 then, according to positive or negative β 1 , we obtain everywhere positive or negative CCFs. For example, see the CCFs in Figure 2 between the response (Admissions) and pollutants CO, NO 2 that are positive and the CCFs between the response (Admissions) and O 3 , SO 2 that are negative at every lag, respectively.

Consider another widely used link function, the linear one, and define the linear GAM-PCA-VAR model by the link

µ t = β 0 + p i=1 ∞ j=0 β ij f (Z i(t-j) ), (12)
where β 0 , β ij ∈ R + , i = 1, . . . , p, j ∈ Z + are parameters and f : R → R + is a known function, e.g., f (z) = exp(z). Let ϕ(x | λ) denote the probability density function of the normal distribution with mean 0 and variance λ.

Proposition 4. Suppose that, for all i = 1, . . . , p,

∞ j=0 β ij < ∞ and τ i := ∞ -∞ f (x)ϕ(x | λ i )dx < ∞.
Then the GAM-PCA-VAR model with linear link (12) has a strictly stationary solution {(Y t , X t )}. Moreover,

E(Y t ) = E(µ t ) = β 0 + p i=1 τ i ∞ j=0 β ij .
Proof. The proof is similar to the proof of Proposition 1.

Clearly, the assumptions of Proposition 4 do not necessarily garantee the existence of higher order moments of linear GAM-PCA-VAR process. Indeed, the rth order moment E(Y r t ) is finite if and only

∞ -∞ f r (x)ϕ(x | λ i )dx < ∞ for all i where r ≥ 1.

Simulation study

In order to evaluate the effect on the parameter estimation of a GAM model in the presence of temporal correlation in the covariate {X t }, a simulation study was conducted. The data were generated according to the model discussed in Example 1. Three estimation methods were considered: the standard GAM with only one covariate where the estimated parameters were β 0 and β 1 (M1); the standard GAM with two covariates, the original one and its 1-step lagged series, where the estimated parameters were β 0 , β 1 , β 2 and φ = -β 2 /β 1 (M2); the full GAM-PCA-VAR model by the procedure described in Section 2 where all parameters β 0 , β 1 , φ, λ were estimated (M3).

For the model discussed in Example 1 the data were generated under β 0 = 0.2, β 1 = 1, λ = 2 and three scenarios were considered as φ = -0.7, 0.3, 0.9 to model strong negative, small positive and strong positive correlations, respectively. In order to model the impact due to some unobservable variables, e.g., environmental ones in the context of the next section, independent N (0, 0.1) distributed random variables were added to the predictor of log µ t for all t ∈ Z. The sample size n = 1000 and the number of Monte Carlo simulations was equal to 100. The empirical values of mean, bias and mean square error (MSE) are displayed in Table 1. All results were obtained by using R-code. In the case of standard GAM estimation (M1) it can be seen that the estimate of β 1 is heavily affected by the autocorrelation structure present in the covariate, by presenting a negative bias which increases in absolute value as |ϕ| increases. The estimated MSE also increases substantially with |ϕ|. On the other hand, it can also be seen that the fitted standard GAM model tends to severely overestimate β 0 . Contrarily, the estimation methods M2 and M3 work equally well, the estimates of the parameters are very close to the true values with noticeably small MSE. The undoubted advantage of method M3 against M2 is that an AR(1) model is also fitted for the covariate where the innovation variance λ is estimated and which can be applied later in the prediction. In this procedure firstly the covariate variable is predicted by equation ( 4) and then the response variable is predicted by the GAM using the link [START_REF] Chen | Ambient air pollution and hospital admission in Shanghai, China[END_REF].

Application to air pollution data

In this study, the number of hospital admissions (Admissions) for respiratory diseases (RD) as response variable was obtained from the main childrens emergency department in the Vitória Metropolitan Area (called Hospital Infantil Nossa Senhora da Gloria), ES, Brazil. The following atmospheric pollutants as covariates were studied: particulate material (PM 10 ), sulphur dioxide (SO 2 ), nitrogen dioxide (NO 2 ), ozone (O 3 ) and carbon monoxide (CO). For details, e.g., descriptive statistics and basic time series plots, see [START_REF] De Souza | Generalized additive models with principal component analysis: an application to time series of respiratory disease and air pollution data[END_REF]. The data analysed in this section can be obtained from http://wileyonlinelibrary.com/journal/rss-datasets

The graphs of the sampling functions of the autocorrelations and partial autocorrelations in Figure 1 show that the series of the number of hospital admissions for RD possesses seasonal behaviour, which was to be expected for this phenomena. Another characteristic observed in the series was an apparently weak stationarity. Similar graphs for the pollutant series can be found in [START_REF] De Souza | Generalized additive models with principal component analysis: an application to time series of respiratory disease and air pollution data[END_REF]. Figure 2 shows the sample cross-correlation functions (CCF) between the response and pollutant covariates. As we discussed in Remark 1 four CCF's among them present similar behaviour: the impact of pollutants CO and NO 2 is positive while the impact of SO 2 and O 3 are negative to the response variable at every lag. This observation is consistent with the PCA result presented in [START_REF] De Souza | Generalized additive models with principal component analysis: an application to time series of respiratory disease and air pollution data[END_REF], see Table 5, where CO and NO 2 form a joint cluster for PC1. On the other hand, all CCF's possess seasonal behaviour as well. Figure 3 shows the sample cross-correlation functions (CCF) between the response variable and the first three PCs derived from applying PCA for the vector of pollutants. In Section 3.2 of [START_REF] De Souza | Generalized additive models with principal component analysis: an application to time series of respiratory disease and air pollution data[END_REF], see Table 5 there, one can see that the first three components correspond to 83.2% of the total variability. The temporal behaviour of the PCs is also presented in the autocorrelation plots of [START_REF] De Souza | Generalized additive models with principal component analysis: an application to time series of respiratory disease and air pollution data[END_REF]Figure 4]. The autocorrelations and the cross-correlations displayed here presented heavy seasonality as well. On the other hand, the shape of the CCFs for the response and PCs can also be classifed into similar groups to the CCFs in Figure 2. The CCF of PC1 is similar to the one of the PM 10 . The CCF of PC2 displays only negative correlations similar to SO 2 and O 3 , while the CCF of PC3 (Figure 3) displays only positive correlations, see CO and NO 2 in Figure 2.

In order to filter the vigorous seasonality both in the response and pollutant variables, seasonal ARMA filters with a 7-day period were applied. The pollutant vector time series and the one-dimensional response time series were filtered by SVAR 7 (1) and SARMA 7 (1, 1) processes, respectively. The residuals obtained by these filters indicate remaining significant correlations, see the CCFs between these residuals in Figure 4. The significant cross-correlations and their respective lags are presented in Table 2. Clearly, the correlations which belong to the negative lags are spurious. However, the correlations which belong to the positive lags measure the true impact of a covariate. For example, there are significant correlations at lag 2 for pollutants PM 10 , NO 2 and CO equally which could mean that the influence of these pollutants to the response indicates 2 days delay. Contrarily, the influence of the pollutants SO 2 and O 3 presents far delays. Figure 5 shows the sample CCF between the residuals of the response variable and the first three PCs after the filtering. The significant cross-correlations and its respective lags are presented in Table 3. It should be emphasized that there are strong coincidences in the lags between Table 2 and3. For example, the lag 2 in PC1 corresponds to the pollutants PM 10 , NO 2 and CO, the lag 6 in PC1 corresponds to the pollutant CO, while lag 25 in PC1 corresponds to the pollutant O 3 . The lag 12 in PC2 corresponds to the pollutant SO 2 . Finally, the lag 14 corresponds to the pollutant NO 2 and the lag 23 to the pollutants SO 2 and NO 2 . These correspondences are compatible with the clustering derived in [START_REF] De Souza | Generalized additive models with principal component analysis: an application to time series of respiratory disease and air pollution data[END_REF]Table 7]. The fitted GAM-PCA-VAR model with its goodness-of-fit measures are reported in [START_REF] De Souza | Generalized additive models with principal component analysis: an application to time series of respiratory disease and air pollution data[END_REF] as well. We note that in this fitted model f ij = 0 was chosen for all j > 0. In view of the above results the GAM-PCA-VAR model with link log µ t = β 0 + can also be a possible candidate, where I i denotes the set of lags which belong to the significant cross-correlation between the residuals of the response and the ith PC. This model can be fitted by using the procedure described in Section 2.

Conclusions

A hybrid called GAM-PCA-VAR model composed by three statistical tools, the VAR model, PCA and the GAM, with Poisson marginal distribution, was developed in a more general framework than in [START_REF] De Souza | Generalized additive models with principal component analysis: an application to time series of respiratory disease and air pollution data[END_REF]. A three-stage estimation method was proposed and studied by simulation for some examples. Some theoretical properties were also proved. The model was applied to describe the dependence between the number of hospital admissions for respiratory diseases and air pollutant covariates.

An extension of the proposed estimation method for the GAM-PCA-VAR model by a variable selection procedure which ensures that only the significant PCs with their respective lags are involved into the model will be pursed in future works.

Fig. 1 .

 1 Fig. 1. autocorrelation function (ACF) and partial autocorrelation function (PACF) of the response variable.

Fig. 2 .

 2 Fig. 2. Sample cross-correlation function (CCF) of the response and pollutant variables.

Fig. 3 .

 3 Fig. 3. Sample cross-correlation function (CCF) of the response and first three PCs.

Fig. 4 .

 4 Fig. 4. Sample cross-correlation function (CCF) between the response and pollutant variables after the filtering.

Fig. 5 .

 5 Fig. 5. Sample cross-correlation function (CCF) between the response and PCs after the filtering.

  t-1 -measurable for all t ∈ Z. One can see that the time series {X t } of covariates depends on {Z t } by formula X t =

	∞ k=0 Φ k AZ t-k for all t, see [4, Example 11.3.1].
	The dependence of the response time series {Y t } from the explanatory vector
	time series {X t } in the GAM-PCA-VAR model can be described by three trans-
	formation steps. Clearly, by equation (

  2 ), i.e., µ t has a lognormal distribution, and the moment generating function of ξ ∼ N (β 0 , σ 2 ) is given by M ξ (t) := E(exp(tξ)) = exp(β 0 t + (σt) 2 /2). Thus, the non-negative integer valued random variable Y t is finite with probability one for all t ∈ Z. The vector time series {Z t } forms a Gaussian white noise. Hence it is strictly stationary process with backshift operator B(Z t ) = Z t-1 for all t ∈ Z. Since both stochastic processes {Y t } and {X t } depend on {Z t } through time-invariant functionals, we have the strict stationarity of {(Y t , X t )} and B

Table 1 .

 1 Simulation results for model in Example 1

	Estimation method	φ Parameter Mean Bias MSE
	M1: GAM with Xt	-0.7 β0 = 0.2 0.699 0.499 0.253
		β1 = 1	0.507 -0.492 0.244
	M2: GAM with Xt, Xt-1	β0 = 0.2 0.204 0.004 0.001
		β1 = 1	0.999 -0.001 0.0002
		φ = -0.7 -0.7	0 0.0001
	M3: GAM-PCA-VAR	β0 = 0.2 0.205 0.005 0.001
		β1 = 1	0.999 -0.001 0.0002
		φ = -0.7 -0.695 0.004 0.0005
		λ = 2	2.003 0.003 0.008
	M1: GAM with Xt	0.3 β0 = 0.2 0.302 0.102 0.012
		β1 = 1	0.905 -0.095 0.009
	M2: GAM with Xt, Xt-1	β0 = 0.2 0.209 0.009 0.001
		β1 = 1	0.998 -0.002 0.0002
		φ = 0.3	0.3	0 0.0002
	M3: GAM-PCA-VAR	β0 = 0.2 0.209 0.009 0.001
		β1 = 1	0.999 -0.001 0.0002
		φ = 0.3	0.306 0.006 0.0008
		λ = 2	1.995 -0.005 0.009
	M1: GAM with Xt	0.9 β0 = 0.2 1.002 0.802 0.651
		β1 = 1	0.191 -0.809 0.655
	M2: GAM with Xt, Xt-1	β0 = 0.2	0.2	0	0.001
		β1 = 1	1	0 0.0002
		φ = 0.9	0.899 -0.001 0
	M3: GAM-PCA-VAR	β0 = 0.2 0.203 0.003 0.001
		β1 = 1	1	0 0.0002
		φ = 0.9	0.899 -0.001 0.0001
		λ = 2	2.007 0.007 0.0086

Table 2 .

 2 Significant cross-correlations and their respective lags between the response and pollutants after the filtering

			RD×SO2				RD×NO2
	Lag	-19	-14	-6	12	23	-12	2	4	14	22
	Value -0.063 -0.062 -0.042 -0.047 -0.051 -0.044 -0.050 0.048 0.053 -0.044
			RD×PM10		RD×CO	RD×O3
		Lag	2	23	-12	2	6	9	25
		Value -0.044 -0.043 -0.053 -0.048 0.045 0.054 -0.055

Table 3 .

 3 Significant cross-correlations and their respective lags between the response variable RD and PCs after the filtering

			RD×PC1				RD×PC2			RD×PC3
	Lag	-14	-12	2	6	25	-5	-2	5	12	1	14	23
	Value -0.051 -0.046 -0.057 0.046 0.043 -0.048 -0.046 0.048 -0.047 0.042 -0.078 -0.045

Acknowledgments

The authors thank the following agencies for their support: the National Council for Scientific and Technological Development (Conselho Nacional de Desenvolvimento Científico e Tecnológico -CNPq), the Brazilian Federal Agency for the Support and Evaluation of Graduate Education (Coordenação de Aperfeiçoamento de Pessoal de Nível Superior -CAPES), Espírito Santo State Research Foundation (Fundação de Amparo à Pesquisa do Espírito Santo -FAPES) and Minas Gerais State Research Foundation (Fundação de Amparo à Pesquisa do estado de Minas Gerais -FAPEMIG). Márton Ispány was supported by the EFOP-3.6.1-16-2016-00022 project.The project is co-financed by the European Union and the European Social Fund. Pascal Bondon thanks to the Institute for Control and Decision of the Université Paris-Saclay.


Valdério Anselmo Reisen
email: valderioanselmoreisen@gmail.com

Marcio Adriano

Céline Sgrancio

Edson Zambon Lévy-Leduc

Higor Monte

Aranda Henrique

Pascal Cotta

Flávio Bondon

Augusto

Adriano Marcio Sgrancio

Céline Lévy-Leduc

Edson Zambon Monte

Henrique Aranda Cotta

Pascal Bondon

Augusto Flávio

Ziegelmann

Robust factor modelling for high-dimensional time series: An application to air pollution data

Keywords: Factor analysis, Time series, Robustness, Eigenvalues, Reduced rank, Air pollution

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

In the last fifty years, issues related to air pollution have grown into a major problem, specially in developing countries, where the air quality has been degraded as a result of industrialization, population growth, high rates of urbanization and inadequate or non-existent policies to control air pollution. The problems caused by air pollution produce local, regional and global impacts. Among different environmental problems, air pollution is reported to cause the greatest damage to health and loss of quality of life see, for example, WHO [START_REF]Air pollution estimates[END_REF]. The most common health problems caused by air pollution are asthma, rhinitis, burning eyes, fatigue, dry cough, heart and lung diseases and heart failure. The main pollutants are carbon monoxide (CO), sulphur dioxide (SO 2 ), nitrogen oxides (NO x ), ozone (O 3 ) and inhalable particles with diameter smaller than 10 µm (PM 10 ). The papers by Brunekreef and Holgate [START_REF] Brunekreef | Air pollution and health[END_REF], Maynard [START_REF] Maynard | Key airborne pollutants: the impact on health[END_REF], WHO [START_REF]Air quality guidelines: global update 2005[END_REF], Curtis et al. [START_REF] Curtis | Adverse health effects of outdoor air pollutants[END_REF] and Souza et al. [START_REF] Souza | Generalized additive model with principal component analysis: An application to time series of respiratory disease and air pollution data[END_REF] discuss the relationship between these pollutants and health problems. In addition, air pollution contributes to the degradation of the environment, the greenhouse effect among many others problems.

In recent studies related to air pollution, much attention has been paid to mathematical receptor models with the aim to measure and analyse the pollutant concentrations at the source of emission. For this, mathematical and statistical tools are used to provide the identification of the pollutant emission sources from chemical characteristics of the particles on the receiver and the pollutant emission sources see, for example, Seinfeld and Pandis [START_REF] Seinfeld | Atmospheric chemistry and physics: from air pollution to climate change[END_REF]. In the literature, the most studied receptor models are: chemical mass balance (CMB), multivariate analysis, principal component analysis techniques (PCA), factor analysis model (FA), multiple linear regression, cluster analysis and positive matrix factorization (PMF) (Watson et al. [START_REF] Watson | Receptor modeling application framework for particle source apportionment[END_REF]). In particular, the classical factor analysis methodology has been widely used in air pollution analysis specially for the identification of emission sources, the management of monitoring networks, regression analysis, cluster analysis and prediction.

In many practical problems, it is quite common to have observations which accommodate the serial dependence of each component and the interdependence between different components, that is, the data are time-dependent. However, it should be noted that, among the studies that adopted the classical PCA and factor analysis techniques, the time dependence of the data is a commonly neglected feature. A basic assumption of the multivariate statistical tools is that the data are independent in time (see for example Anderson [START_REF] Anderson | An introduction to multivariate statistical analysis[END_REF] and Johnson and Wichern [START_REF] Johnson | Applied multivariate statistical analysis[END_REF]). To deal with autocorrelated data in factorial analysis, Peña and Box [START_REF] Peña | Identifying a simplifying structure in time series[END_REF],

Stock and Watson [START_REF] Stock | Forecasting using principal components from a large number of predictors[END_REF], Lam et al. [START_REF] Lam | Estimation of latent factors for high-dimensional time series[END_REF] and Lam and Yao [START_REF] Lam | Factor modeling for high-dimensional time series: inference for the number of factors[END_REF] studied the factor modelling for multivariate time series from a dimension-reduction point of view. Contrarily to the PCA and factor analysis for independent observations, these papers look for factors which drive the serial dependence of the original time series. Further discussions and additional references can be found in Lam and Yao [START_REF] Lam | Factor modeling for high-dimensional time series: inference for the number of factors[END_REF].

Since factor analysis allows to reduce the order of the estimated model, this technique has been widely used for forecasting. According to Stock and Watson [START_REF] Stock | Forecasting using principal components from a large number of predictors[END_REF], the dimension reduction becomes a central concern for forecasting when the number of candidate predictor series is very large. This issue can make the forecast investigation impractical in a real application, for example in the use of vector autoregressive moving average (VARMA) models with a large number of variables. This high-dimensional problem is simplified by modelling the common dynamics in terms of a relatively small number of unobserved latent factors. Then, forecasting can be carried out in two-step approach: first, a time series of the factors is estimated from the predictors; second, the relationship between the variable to be forecast and the factors is estimated, for example, using a linear regression.

Environmental time series are often of high dimension due to the large number of measurements recorded across many different locations. These data may also present interesting phenomena to be considered from an applied and theoretical point of view. Indeed, the concentration of pollutant may present high peaks, which can be seen as outlying values from an aspect of statistical analysis. Outliers and high dimension data are common in many areas of applied mathematics. Therefore, the methodology proposed here can be widely used in many other areas where the multivariate techniques are the main tools to describe and interpret the data. This is the case of the health science area (Perc [23], Gosak et al. [START_REF] Gosak | The relationship between node degree and dissipation rate in networks of diffusively coupled oscillators and its significance for pancreatic beta cells[END_REF], Souza et al. [START_REF] Souza | Generalized additive model with principal component analysis: An application to time series of respiratory disease and air pollution data[END_REF]), air route network problems (Zhang et al. [START_REF] Zhang | Analysis of flight conflicts in the chinese air route network[END_REF], Lordan et al. [START_REF] Lordan | Robustness of airline route networks[END_REF]), environmental engineering (Zamprogno [START_REF] Zamprogno | PCA applied in time series data with applications to air quality data[END_REF]) and statistical process controls (Vanhatalo and Kulahci [30]), to name a few.

As is well known, outliers can destroy the statistical properties of the estimates such as the sample mean and sample covariance (see, for example, Chang et al. [START_REF] Chang | Estimation of time series parameters in the presence of outliers[END_REF], Tsay [START_REF] Tsay | Outliers, level shifts, and variance changes in time series[END_REF], Chen and Liu [START_REF] Chen | Joint estimation of model parameters and outlier effects in time series[END_REF] and the references therein). Since the parameter estimation is connected with these sample functions, the final estimated time series model can be strongly affected by the outliers. When the series has additive outliers, one way to deal with model estimation is to use robust estimates of these statistics. For a univariate time series, Ma and Genton [START_REF] Ma | Highly robust estimation of dispersion matrices[END_REF] proposed a robust sample autocorrelation function (ACF) based on the robust scale estimate Q n (.) suggested in Rousseeuw and Croux [START_REF] Rousseeuw | Alternatives to the median absolute deviation[END_REF]. This robust ACF estimator was recently studied by Lévy-Leduc et al. [START_REF] Lévy-Leduc | Asymptotic properties of U-processes under long-range dependence[END_REF], Lévy-Leduc et al. [START_REF] Lévy-Leduc | Large sample behaviour of some well-known robust estimators under long-range dependence[END_REF] and Lévy-Leduc et al. [START_REF] Lévy-Leduc | Robust estimation of the scale and the autocovariance function of Gaussian short and long-range dependent processes[END_REF].

This paper considers multivariate time series with additive outliers using the factor analysis technique for dimension reduction. In this context, a robust version of the dimension reduction estimator given in Lam and Yao [START_REF] Lam | Factor modeling for high-dimensional time series: inference for the number of factors[END_REF] is proposed here. Theoretical results are discussed, and the method performance is investigated through Monte Carlo simulations. The proposed methodology is applied to PM 10 concentrations measured at the Air Quality Automatic Monitoring Network (AQAMN), Vitória, Brazil.

The rest of the paper is organized as follows. In Section 2, the model and the estimation methods are presented. Section 3 discusses the asymptotic properties of the robust eigenvalues. Section 4 presents some Monte Carlo experiments. Section 5 considers an application of the proposed methodology and some concluding remarks are provided in Section 6.

Factor model in time series and estimation methods

The factor model and the estimate of the number of factors

Let k be the number of candidate predictor series. Let Z t , t ∈ Z, be a kdimensional zero-mean vector of an observed time series. Let also X t be an unobserved r-dimensional vector of common factors (r ≤ k). It is assumed that Z t is generated by

Z t = PX t + ε t , (1)
where P is an unknown k × r matrix of parameters of rank r, denoted by the factor-loading matrix, and ε t is a k-dimensional white-noise sequence with fullrank covariance matrix Σ ε . When r is small relative to k, the model presented in Equation ( 1) is most useful, since it will result in a multivariate time series model with a reduced dimension and, consequently, can lead to a much simpler multivariate time series model for forecasting.

In the sequel, the following assumption is made.

(A1) X t , t ∈ Z, is a multivariate stationary process and ε t ∼ WN(0, Σ ε ). Moreover, X t and ε t are assumed to be uncorrelated and P P = I r , where I r denotes the r × r identity matrix.

Note that Assumption 1 is to ensure identifiability in Equation (1); see Lam and Yao [START_REF] Lam | Factor modeling for high-dimensional time series: inference for the number of factors[END_REF] and Peña and Box [START_REF] Peña | Identifying a simplifying structure in time series[END_REF] for further details.

It follows from Equation ( 1) and under Assumption (A1) that the covariance matrices of Z t are given by

Γ Z (0) = PΓ X (0)P + Σ ε , (2)
Γ Z (h) = PΓ X (h)P , h ≥ 1, (3)
where

Γ X (h) = E[X t-h X t ] is the covariance matrix of X t .
Based on the observations Z 1 , ..., Z n , the first step in Equation ( 1) is to estimate the number of factors r in order to compute the estimate of the k × r factor loading matrix P. Once P is obtained, the estimator for the factor process and the residuals are respectively given by Xt

= P Z t , (4)
and

ˆ t = (I d -P P )Z t . (5)
For further details on the estimation of P, see Lam and Yao [START_REF] Lam | Factor modeling for high-dimensional time series: inference for the number of factors[END_REF].

Let ΓZ (h) denote the sample covariance matrix of Z t at lag h and let

M = h 0 h=1 ΓZ (h) ΓZ (h) (6)
for a prescribed integer h 0 ≥ 1. Following the lines of Lam and Yao [START_REF] Lam | Factor modeling for high-dimensional time series: inference for the number of factors[END_REF], the estimator for the number of factors r is given by:

r = argmin 1≤i≤R λi+1 / λi , (7)
where r < R < k is a constant, λ1 ≥ . . . ≥ λk are the eigenvalues of M. Under some assumptions, Lam and Yao [START_REF] Lam | Factor modeling for high-dimensional time series: inference for the number of factors[END_REF] derive the asymptotic properties of the above results. According to Lam and Yao [START_REF] Lam | Factor modeling for high-dimensional time series: inference for the number of factors[END_REF], in practice, for example, R = p/2 may be used. Note that, there is no need to extend the test up to p since it is expected that the minimum eigenvalue of M be practically 0, specially for a small n and a large p. It is important to say that when p and n are of the same order, the estimators for eigenvalues are no longer consistent, although the ratio based estimator still can be used, see Lam and Yao [START_REF] Lam | Factor modeling for high-dimensional time series: inference for the number of factors[END_REF]. As previously stated, the aim of this paper is to propose robust estimators of M and r which are based on a robust covariance matrix estimator for Z t . These issues are discussed in the following section.

The robust estimator of the number of factors

Let (Y i ) i≥1 be a stationary Gaussian process. Given the observations Y 1:n = (Y 1 , . . . , Y n ), the estimator of scale proposed by Rousseeuw and Croux [START_REF] Rousseeuw | Alternatives to the median absolute deviation[END_REF] is defined by

Q n (Y 1:n ) = c Y i -Y j ; 1 ≤ i, j ≤ n ( n 2 /4 ) , (8)
where c = 2.21914. Now, consider the following assumption on X t .

(A2) X t = (X 1,t , X 2,t , . . . , X r,t ) is a multivariate stationary zero-mean Gaussian process satisfying

h≥1 |γ X i, j (h)| < ∞, for all i, j ∈ {1, . . . , r}, where γ X i, j (h) = Cov[X i,t , X j,t+h ]
. By using Equations ( 1) and ( 3), (Z t ) is also a multivariate stationary zero-mean Gaussian process satisfying

h≥1 |γ i, j (h)| < ∞, for all i, j ∈ {1, . . . , k}, (9)
where

γ i, j (h) = Cov[Z i,t , Z j,t+h ].
From the estimator Q n defined in Equation ( 8) and from the observations (Z 1 , . . . , Z n ), it is proposed here a robust estimator of γ i, j (h) = Cov(Z i,t , Z j,t+h ) for all i, j in {1, . . . , k} defined as follows

γQ i, j (h) = 1 4 Q 2 n-h (Z i,1:n-h + Z j,h+1:n ) -Q 2 n-h (Z i,1:n-h -Z j,h+1:n ) , (10)
where Z i,1:n-h = (Z i,1 , . . . , Z i,n-h ) and Z j,h+1:n = (Z j,h+1 , . . . , Z j,n ). γQ i, j (h) is the multivariate estimator of the univariate case suggested by Ma and Genton [START_REF] Ma | Highly robust estimation of the autocovariance function[END_REF].

From Equation [START_REF] Johnson | Matrix theory and applications[END_REF], the robust estimator of the covariance matrix of Z t is given by

ΓQ (h) =             γQ 1,1 (h) γQ 1,2 (h) . . . γQ 1,k (h) . . . . . . . . . . . . γQ k,1 (h) γQ k,2 (h) . . . γQ k,k (h)             . ( 11
)
Based on Equations ( 6) and [START_REF] Johnson | Applied multivariate statistical analysis[END_REF], the robust version of the estimator M is suggested here as follows

MQ = h 0 h=1 ΓQ (h) ΓQ (h) . (12)
Therefore, the robust estimator rQ of r is similarly obtained from Equation ( 7) by replacing λi+1 and λi by λQ i+1 and λQ i , respectively, where ( λQ i ) 1≤i≤k are the eigenvalues of MQ .

Theoretical results

This section provides some analytical results to theoretically support the robust approach discussed in the previous section.

Theorem 1. Let h be a fixed positive integer and ΓQ (h

) 1≤i, j≤k = γQ i, j (h) 1≤i, j≤k
, where γQ i, j (h) is defined in Equation [START_REF] Johnson | Matrix theory and applications[END_REF]. Assume that Assumptions (A1) and (A2) hold, then

√ n sup 1≤ j≤k λQ j -λ j = O p (1), as n → ∞,
where ( λQ j ) 1≤ j≤k and (λ j ) 1≤ j≤k denote the eigenvalues of h 0 h=1 ΓQ (h) ΓQ (h) and h 0 h=1 Γ(h)Γ(h) , respectively, where (Γ(h)) 1≤i, j≤k = γ i, j (h) 1≤i, j≤k and h 0 is a fixed integer larger than 1.

The proof of this theorem directly follows from Lemmas 1, 2 and 3 given below and proved in Section 8.

Remark 1. By Theorem 1 and Lam and Yao [12, Proposition 1], it can be seen that the eigenvalues of the robust estimator covariance matrix of Z t have the same rate of convergence as the eigenvalues of the standard estimator of the covariance matrix of Z t .

Lemma 1. Let Ân be a sequence of k × k symmetric matrices and A be a k × k symmetric matrix such that u n ( Ân -A) = O p [START_REF] Anderson | An introduction to multivariate statistical analysis[END_REF], where u n is a sequence of positive numbers tending to infinity as n tends to infinity, then

u n sup 1≤ j≤p |λ j ( Â) -λ j (A)| = O p (1), as n → ∞,
where (λ j ( Â)) 1≤ j≤k and (λ j (A)) 1≤ j≤k are the eigenvalues of Ân and A, respectively. Lemma 2. Let Ân (h) be a sequence of k × k symmetric matrices and A(h) be a k × k symmetric matrix such that u n ( Ân (h) -A(h)) = O p (1), for each fixed h ∈ {1, . . . , h max }, where u n is a sequence of positive numbers tending to infinity as n tends to infinity, then

u n         h max h=1 Ân (h) Ân (h) - h max h=1 A(h)A(h)         = O p (1).
Lemma 3. Let h be a non negative integer and i and j be two integers in {1, . . . , k}. Assume that Assumptions (A1) and (A2) hold, then the robust autocovariance estimator γQ i, j (h) defined in Equation ( 10) satisfies the following central limit theorem

√ n(γ Q i, j (h) -γ i, j (h)) d -→ N(0, σ 2 i, j (h)), as n → ∞,
where

σ 2 i, j (h) = E[ψ(Z i,1 , Z j,1+h ) 2 ] + 2 ≥1 E[ψ(Z i,1 , Z j,1+h )ψ(Z i, +1 , Z j, +1+h )],
where ψ is defined in Equation ( 15).

Monte Carlo studies

This section reports simulation results related to the performance of the proposed methodology previously discussed for finite sample size. In this context, the empirical study considered X t as a VAR(1) model with r = 3. The VAR(1) model was generated with independent white noise vector from N(0, I) and the matrix of the coefficients Φ 1 . This matrix is displayed in Table 1. Note that the Φ coefficients corresponds to a process with no temporal correlation outside the main diagonal. The sample sizes are n = 50, 100, 200, 400, 800 and 1600, and k = 0.2n, 0.5n, 0.8n, and h 0 = 1 (these are also considered in Lam and Yao [START_REF] Lam | Factor modeling for high-dimensional time series: inference for the number of factors[END_REF]). The factor model (Equation ( 1)) was generated as follows: first, all k × r elements of matrix P were generated as independent observations from the uniform distribution on the interval [-1, 1] (see, also, Lam and Yao [START_REF] Lam | Factor modeling for high-dimensional time series: inference for the number of factors[END_REF]). The process ε t in Equation (1) consists of N(0, 1) components independent between each other and across t. The statistical quantities were computed based on 1000 replications. These simulations were ran using the R programming language and the code is available upon request.

[Table 1

about here.]

The main interest in this empirical study is to verify the performance of the statistics r and rQ in the context of VAR(1) models with and without outliers. For this, the estimate of P(r = r) is obtained by computing the relative frequencies of r = r, denoted here as f f req. (r = r), where r is the estimator of r. The rQ estimator was similarly computed.

Table 2 reports the relative frequency estimates ( f f req. (r = 3)) for P(r = 3) of the model in the study. From this, it is observed that the ratio-based estimator of r improves when n is increased. Similar performance of the ratio is also observed when the dimension k increases. These are in accordance with the results given in Table 1 of Lam and Yao [START_REF] Lam | Factor modeling for high-dimensional time series: inference for the number of factors[END_REF], that is, the asymptotic properties of r is corroborated with these finite sample size investigations.

The results related to the alternative method rQ are displayed in Table 3. It shows that this estimator presents similar empirical performance of r, which is an expected result based on the asymptotic theory of both estimators of the ratio (see Remark 1). Therefore, this simple empirical study together with the theoretical results discussed in the previous section give support to using rQ as an alternative method to estimate the number of the factors in the model presented in Equation ( 1).

[Table 2

about here.]

[Table 3 about here.] Now, let X * t , t ∈ Z, be a vector process contaminated by additive outliers defined as follows

X * t = X t + ω • δ t , (13)
where "•" is the Hadamard product (Johnson [START_REF] Johnson | Matrix theory and applications[END_REF]). The vector ω = (ω 1 , ..., ω k ) is a magnitude vector of additive outliers and the random vector δ t = (δ 1,t , ..., δ k,t ) indicates the occurrence of an outlier at time t, in variable i, i = 1, ..., k, such as P(δ i,t = -1) = P(δ i,t = 1) = p/2 and P(δ i,t = 0) = 1p, where E[δ i,t ] = 0 and E[δ 2 i,t ] = Var[δ i,t ] = p. The model described above assumes that X * t and δ t are independent processes. Also, it is assumed that the elements of δ t are not correlated and temporally uncorrelated, i.e., E[δ t δ t ] = Σ δ = diag(p, ..., p) and E[δ t δ t+h ] = 0 for h 0.

Remark 2. δ i,t , i = 1, ..., k, is the product of Bernoulli(p) random variable with Rademacher random variable, the latter equals 1 or -1, both with probability 1/2.

In this empirical investigation, the probability of an outlier occurring at time t is p = 0.05 and, without loss of generality, it is also assumed that ω = [ω 1 , 0, 0] ; that is, X * 1,t , t = 1, ..., n, is the only process in X * t = (X * 1,t , X * 2,t , X * 3,t ) contaminated with outliers and ω 1 = 15.

Table 4 shows the relative frequency estimates, r and rQ , for the dimensional reduction for Model 1, when r = 3 and in the presence of outliers. The standard case, that is, ΓX * and p = 0, is in accordance with the results given by Tables 2 and3 for r and rQ , respectively.

From the cases where there are outliers (column of Table 4 in which p 0), it is clearly perceived that the standard method is destroyed by a single outlier at the time t by presenting a substantial reduction of the correct estimated ratio. This is an expected result since r is based on the standard sample matrix of autocorrelation. Note that, according the simulation plan, the expected number of outliers in the whole series X * 1,t is 5. The conclusions of the empirical study in this section corroborate the Proposition 1 in Cotta et al. [START_REF] Cotta | Robust estimation of covariance and correlation functions of a stationary multivariate process[END_REF], which demonstrates theoretically that one outlier is enough to destroy the properties of the sample autocorrelation function of multivariate processes.

In contrast with the empirical performance of r, rQ keeps almost unchanged with the proportion of the outliers in the present study. The estimated ratios are very close in the contaminated and uncontaminated series. Therefore, the percentage of outliers in this empirical example seems to be, in general, not strong enough to destroy the robustness property of rQ and this indicates that the methodology proposed here may be used when the presence of outliers in the series is uncertain. Other simulation cases with different degrees of contamination present similar conclusions and are available upon request.

[Table 4 about here.]

Application to the pollutant PM 10

This section presents an application of the methodology discussed previously for PM 10 pollutant concentrations measured at the Automatic Quality Monitoring Network (AAQMN) in the Greater Vitória Area (GVR), Espírito Santo, Brazil. GVR is comprised of seven cities with a population of approximately 1.9 million inhabitants in the area of 2,319 km 2 . The AAQMN consists of eight monitoring stations distributed in the cities of GVR; Laranjeiras, Carapina, Camburi, Suá, Vitória (center), Vila Velha (center), Ibes and Cariacica. The pollutant PM 10 , expressed in µg/m 3 and was hourly measured from January 2008 to December 2009, k = 8, though the daily average values (n = 731) are used in this empirical study. This follows the same lines of the application given in Lam and Yao [START_REF] Lam | Factor modeling for high-dimensional time series: inference for the number of factors[END_REF]. Now let Z t , t = 1, ..., 731, be the vector of the PM 10 concentrations, that is, Z t = (Z 1,t , ..., Z 8,t ) ] where Z i,t corresponds to PM 10 concentration at location i.

Figure 1 shows the plots of the PM 10 concentrations for the eight stations. Based on this figure, the series indicated that they present high levels of pollutant concentrations which can be identified, from the statistical point of view, as additive outliers. This is justified by the fact that they produce similar impact on the sample ACF caused by additive outliers, that is, they lead to a reduction of the sample autocorrelation values. Therefore, the high values can be seen as (additive) outliers. Hence, both methods, the robust and non-robust approaches discussed previously, are used here to verify whether or not these high levels make any impact on the factor model estimation.

The sample ACFs displayed in Figures 2 (the classical ACF estimator) and 3 (the robust ACF estimator) show possible seasonal pattern of period s = 7 which is an expected evidence since the data are daily averaged. In terms of the comparison between the sample ACF estimates, as a simple case, the classical sample ACF values at VVCentro station are 0.47, 0.12, 0.15 and 0.13 for lags h = 1, 3, 5, 10, respectively, while the ACF based on the Q n function are 0.54, 0.25, 0.20 and 0.19. This simple case shows that the high levels of PM 10 at VVCentro station indicated an reduction of the sample ACFs values of the classical autocorrelation estimator.

The effect of atypical observations on the estimation of the ACF function is well discussed in Molinares et al. [START_REF] Molinares | Robust estimation in longmemory processes under additive outliers[END_REF] for a single time series. The comparison between the sample ACFs of PM 10 concentrations from the other stations presented similar conclusions.

[Figure 1 From the above discussion, it is expected that the FA estimated model will show different performance for the two methodologies, that is, for the standard and robust ones. The estimates of the number of factors r were computed by performing an eigenanalysis on M and on MQ of Equations ( 6) and ( 12), respectively, with h 0 = 7, to capture the seasonality feature of the data set. The eigenvalues (the scree plot), in decreasing order, and their ratios obtained using ΓZ are shown in Figure 4 (a) and (b), respectively. The corresponding robust versions, i.e., using ΓQ , are shown in Figure 5 (a) and (b), respectively. As can be seen, the robust and the classical ACF estimators suggest r = 2 and r = 1, respectively. This confirms the expected result previously stated. The reduction was not affected when varying the value of h 0 .

[Figure 4 about here.]

[Figure 5 about here.] Figure 6 (a) and (b) display the time series plots of the two components series X1,t and X2,t , respectively, of the estimated factors Xt = ( X1,t , X2,t ) for ΓZ,Q (see Equation ( 4)).

Following similar lines as in the application section of Lam and Yao [START_REF] Lam | Factor modeling for high-dimensional time series: inference for the number of factors[END_REF], it is now addressed the study to verify the percentage of the variability of the pollutant vector Z t explained by P Xt . For this, the PM 10 concentration measured at the Laranjeiras station is used and the original data and the estimated one are displayed in Figures 6 (c) and (d), respectively. From the plots of Figure 6, it is possible to note that the behaviour of the two factors is similar to the ones of the Laranjeiras station, including the high volatility and the periods with peaks of PM 10 concentrations. The estimated PM 10 concentrations (Figure 6(d)) are a linear combination of the estimated factor coefficients with values 0.4753 and 0.8231, for the first and second factor, respectively. This estimated series behaves similarly to the observed values (Figure 6 (c)).

The quantity ||Bu|| 2 /||u|| 2 = 0.0015 where u is the 731 × 1 vector consisting of the PM 10 concentrations of Laranjeiras station index over the same period of time, and B denotes the projection matrix onto the orthogonal complement of the linear space spanned by the two components series X1,t and X2,t . From this, 99,85% of the PM 10 concentrations of Laranjeiras station can be explained by a linear combination X1,t and X2,t . Therefore, these analyses suggest the following factor model for the PM 10 concentrations vector

Z t = p 1 X 1,t + p 2 X 2,t + ε t , (14)
where X 1,t denotes the first factor, X 2,t is the second factor, and ε t is a vector whitenoise process. Finally, for forecasting purpose, Equation ( 14) is a simpler model than a Kmultivariate stationary time series model with dimension k = 8, that is, the h-step ahead forecast for the Z t series is simplified using the formula Ẑ(h)

T +h = P X(h) T +h , where X(h)

T +h is an h-step ahead forecast for {X t }, based on the estimated past values X1 , . . . , XT (see, Lam et al. [START_REF] Lam | Estimation of latent factors for high-dimensional time series[END_REF]).

[Figure 6 about here.]

Conclusions

In this paper, a robust method for high-dimensional time series with additive outliers is proposed. Some theoretical results are discussed and verified through Monte Carlo experiments under different scenarios of outliers contamination. The simulations illustrated the effect of the additive outliers on the reduction of the factor dimension. The empirical investigation showed that the robust method presented better performance compared to the classic procedure of identifying the number of factors being an alternative method when there is any evidence of atypical observations in the multivariate time series data, such as high levels of the pollutants in the pollution area. In addition, the proposed methodology was used to identify pollution behaviour of the pollutant PM 10 , which can be very useful for the management of the air quality network. The results in this paper will hopefully stimulate further research on this topic. By A2 and Equation ( 9), ≥1 |γ (+) i, j ( )| < ∞ and ≥1 |γ (-) i, j ( )| < ∞. The proof of this lemma, thus, follows the same lines as the ones of Lévy-Leduc et al. [START_REF] Lévy-Leduc | Robust estimation of the scale and the autocovariance function of Gaussian short and long-range dependent processes[END_REF]Theorem 2] by replacing X i and X i+h by Z i,t and Z j,t+h , respectively, and the summations on i by summations on t which leads to

√ n -h γQ i, j (h) -γ i, j (h) = 1 √ n -h n-h t=1
ψ(Z i,t , Z j,t+h ) + o P (1), where ψ(x, y) = 1 2 γ i,i (0) + γ j, j (0) + γ i, j (h) + γ j,i (-h) IF

      
x + y γ i,i (0) + γ j, j (0) + γ i, j (h) + γ j,i (-h)

, Q, Φ        - 1 2 γ i,i (0) + γ j, j (0) -γ i, j (h) -γ j,i (-h) IF       
xy γ i,i (0) + γ j, j (0) -γ i, j (h) -γ j,i (-h) , Q, Φ

       , (15)
where IF is defined in Equation ( 20) of Lévy-Leduc et al. [START_REF] Lévy-Leduc | Robust estimation of the scale and the autocovariance function of Gaussian short and long-range dependent processes[END_REF]. By applying Arcones [2, Theorem 4], the result is obtained. Table 1: Φ 1 coefficients matrix of the VAR(1) process.

Φ 1 0.60 0.00 0.00 0.00 -0.50 0.00 0.00 0.00 0.30

  about here.] [Figure 2 about here.] [Figure 3 about here.]

Figure 3 :Figure 4 :Figure 5 :

 345 Figure 3: Robust ACF estimates of the PM 10 pollutant concentrations.

Acknowledgements

Part of the simulation and application results in this paper are in chapters of the Ph.D. thesis of the second author under the supervision Prof. V. A. Reisen. The authors would like to thank CNPq, CAPES and FAPES for their financial support. Part of this paper was revised when Prof. Valdério Reisen was visiting CentraleSupélec (from December 2016 to January 2017 and in July 2018). This author is indebted to CentraleSupélec for its financial support. The authors are grateful to a referee for the time and efforts in providing helpful comments and additional references that have led to clarify and substantially improve the quality of the paper.

Proofs

Proof of Lemma 1. By Weyl's Theorem, see Horn and Johnson [9,p. 239], for all j ∈ {1, . . . , k}, it follows that

By exchanging the role of  and A, for all j ∈ {1, . . . , k}, it follows that

Hence, sup

1≤ j≤k

where X 2 denotes the largest absolute value of the eigenvalues of a matrix X. Since u n ( Ân -A) = O p (1), the result follows.

Proof of Lemma 2. The proof of this lemma directly follows from the application of the continuous mapping theorem; see van der Vaart [29, Theorem 2.3].

Proof of Lemma 3. Observe that the autocovariance of the process (Z i,t + Z j,t+h ) t≥1 at lag is equal to

and that the autocovariance of the process (Z i,t -Z j,t+h ) t≥1 at lag is equal to Relative frequency estimates for P(r = 3). . . . . . . . . . . . . . 28 3 Relative frequency estimates for P( rQ = 3). . . . . . . . . . . . . 29 4 Relative frequency estimates for dimensional reduction, n = 100. . 30