Bartolomeu Zamprogno

Valdério A Reisen

Pascal Bondon

Neyval Higor Cotta

Costa Reis Júnior

Higor H Aranda Cotta

Neyval C Reis

Principal component analysis with autocorrelated data

Keywords: Principal component analysis, autocorrelation, cross-correlation, eigenvalue, air pollution

Introduction

PCA is one of the most widely used multivariate techniques to reduce the dimension of a data set while keeping most of the variability of the data. To clarify how important this technique is, [START_REF] Richman | A principal component analysis of sulphur concentrations in the western United States[END_REF] has shown that between 1983 and 1985 over 60 PCA applications, or similar techniques have appeared in the meteorological/climatological journals. More recently, between 1999 and 2000, 53 of the 215 articles of the International Journal of Climatology have applied PCA. This represents 25%, a rate not achieved by any other statistical technique Jollife (2002, page 71).

The use of PCA goes beyond reducing the dimension of data. For example, [START_REF] Karar | Source apportionment of PM 10 at residential and industrial sites of an urban region of Kolkata, India[END_REF] have used PCA as a grouping tool of pollution sources, and [START_REF] Romero | Daily rainfall affinity areas in Mediterranean Spain[END_REF], [START_REF] White | Climate regionalization and rotation of principal components[END_REF] and [START_REF] Cohen | Classification of 500 mb height anomalies using obliquely rotated principal components[END_REF] have applied PCA to identify homogeneous sub-regions of climatic stations in a large geographical area. Besides the use of PCA as a classification tool, several studies have used the technique to extenuate the multicollinearity in a regression analysis context and to detect outliers, see e.g., [START_REF] Liu | Simulation of the daily average PM 10 concentrations at Ta-Liao with Box-Jenkins time series models and multivariate analysis[END_REF], Wang and Pham (2011a), [START_REF] Souza | Principal components and generalized linear modeling in the correlation between hospital admissions and air pollution[END_REF], [START_REF] 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] Reisen | Robust factor modeling for high-dimensional time series: an application to air pollution data[END_REF]. PCA has also been used as a step procedure in other multivariate techniques such as factor analysis, canonical correlation analysis, and discriminant analysis, see e.g., Jollife (2002, Chapter 9). For example, in the financial area, [START_REF] Matteson | Dynamic orthogonal components for multivariate time series[END_REF] have proposed a PCA based approach to modeling the conditional mean vector and conditional covariance matrix of a stationary multivariate autoregressive and conditionally heteroscedastic time series. [START_REF] Hu | Principal volatility component analysis[END_REF] have extended the idea of PCA to principal volatility component analysis with a focus on the dynamic dependence of volatility.

In the statistical control process area, [START_REF] Vanhatalo | Impact of autocorrelation on principal components and their use in statistical process control[END_REF] have illustrated the impact of the autocorrelation on the descriptive ability of PCA and on the monitoring process control. [START_REF] Vanhatalo | On the structure of dynamic principal component analysis used in statistical process monitoring[END_REF] have proposed a driven method to determine the maximum number of lags in dynamic PCA in multivariate time series analysis and a method for determining the number of principal components (PCs) to retain. In the high-dimension setting, [START_REF] Hellton | The impact of measurement error on principal component analysis[END_REF] have addressed the problem of the impact of measurement error on PCA.

In the domain of air quality monitoring, the identification of pollution sources using PCA has been considered by many authors. For example, in the network management context, Pires et al. (2008a,b) have used PCA with monitored pollutant concentrations to manage the monitoring network of the metropolitan area of Porto (Portugal) to reduce costs. The authors have proposed to select only one station among those belonging to a same cluster and having similar concentrations behaviours. They have concluded that six stations instead of ten are sufficient to measure the level of concentration of sulphur dioxide (SO 2 ), and no more than two stations are required for monitoring the PM less than 10 µm in diameter (PM 10 ). [START_REF] Lu | Performance assessment quality monitoring networks using principal component analysis and cluster analysis[END_REF] have evaluated the performance of PCA and cluster analysis for the management of the local air quality monitoring network of Hong Kong (China) with the aim to identify city areas with similar air pollution behaviours and to locate emission sources. They have found that the monitoring stations could be grouped into different classes based on air pollution behaviours.

One of the usual assumptions of PCA is that the data are independent in time. Nevertheless, PCA has been widely used with time series which are time-correlated, without justification. For example, the pollution data considered in the above cited papers are time-dependent. Not taking into account the time-dependent structure of the data may lead to misleading analysis and interpretations. It is essential to recognize that neglecting the required data assumption when using standard statistical methods like PCA may produce biased estimates and spurious results see e.g., [START_REF] Vanhatalo | Impact of autocorrelation on principal components and their use in statistical process control[END_REF].

The effect of time-correlation on model estimation using PCA is also one of the main contribution of [START_REF] Souza | Generalized additive models with principal component analysis: an application to time series of respiratory disease and air pollution data[END_REF], where the multicollinearity issue when using pollutants as covariates in the generalized additive model is solved using PCA, and where it is suggested to use a multivariate time series model to remove the temporal correlation of the covariates. Following similar lines, [START_REF] Melo | Correlação entre percepção do incômodo e exposição ao material particulado presente na atmosfera e sedimentado[END_REF] and Melo et al. (to appear) have considered PCA in a logistic regression model to quantify the association between the pollutants and perceived annoyance. The methodologies proposed in these three papers were mainly based on the theoretical results discussed in [START_REF] Zamprogno | O uso e interpretação de análise de componentes principais, em séries temporais, com enfoque no gerenciamento da qualidade do ar[END_REF]. Wang and Pham (2011b) have also considered PCA in the regression model to quantify the relationship between morbidity and pollutants; however, the temporal correlation of the variables was ignored by the authors.

The purpose of this paper is to generalize the use of PCA, mainly developed for independent observations, to multivariate time series. The effect of different correlation structures of multivariate stationary processes on the interpretation and inference of the PCs is illustrated. The study is justified empirically and theoretically, and a real data set of pollutant concentrations is considered as an example of application. Due to the serial correlation in the data, the PCs are shown to be autocorrelated and cross-correlated. Thus, this paper suggests to pre-whiten the data with a linear model to attenuate the time-correlation before applying PCA. This whitening technique has been considered by some authors in the econometric area, but without discussing the consequence of neglecting the temporal correlation. For example, [START_REF] Matteson | Dynamic orthogonal components for multivariate time series[END_REF] and [START_REF] Hu | Principal volatility component analysis[END_REF] applied vector autoregressive (VAR) models to remove the serial correlation of time series of stock returns before carrying out PCA of the residuals.

The manuscript is structured as follows: Section 2 considers the time series model and theoretical properties of PCA with autocorrelated data. Monte Carlo simulations are addressed in Section 3. Section 4 discusses the real data application and Section 5 concludes the paper.

PCA with time series data

Let X t = [X 1t , . . . , X kt ] , t ∈ Z, be a k-dimensional linear process defined by

X t = µ + ∞ j=0 Ψ j ε t-j , (1)
where

µ ∈ R k , ε t = [ε 1t , . . . , ε kt ] is a vector white noise process such that E(ε t ) = 0 and Γ ε (h) = Cov(ε t , ε t+h ) = E(ε t ε t+h ) = Σ ε if h = 0, 0 if h = 0, (2)
Σ ε is a nonsingular matrix, and the Ψ j 's are k×k matrices of real coefficients satisfying Ψ 0 = I, I being the identity matrix, and ∞ j=0 tr(Ψ j Σ ε Ψ j ) < ∞, where tr(A) denotes the trace of a square matrix A. It follows from (1) and (2) that X t is a second-order stationary process with mean µ and covariance matrix

Γ X (h) = Cov(X t , X t+h ) = E((X t -µ)(X t+h -µ) ) = ∞ j=0 Ψ j Σ ε Ψ j+h , (3)
for all h ≥ 0. In the following, it is assumed without loss of generality that µ = 0.

In the analysis of a multivariate data set, PCA looks for linear combinations of the components capturing the highest percentage of variation of the data. This technique depends exclusively on the covariance or the correlation matrix of the data, see, e.g., [START_REF] Jollife | Principal component analysis[END_REF]. PCA is well suited for time-independent observations since it explains only the contemporaneous correlation of the data and does not take into account the time-correlation. Specifically, PCA calculates the characteristic roots and vectors of Γ X (0). Let λ 1 ≥ . . . ≥ λ k ≥ 0 be the non necessarily distinct eigenvalues of Γ X (0) with corresponding orthonormal (with respect to the usual inner product) eigenvectors p 1 , . . . , p k (p i p i = 1 and p i p j = 0 when i = j). Then Γ X (0)p i = λ i p i for i = 1, . . . , k, and P Γ X (0)P = Λ where P is the k × k matrix whose ith column is p i and Λ is the k × k diagonal matrix whose ith diagonal element is λ i , i.e., Λ = diag(λ 1 , . . . , λ k ), see e.g., Banerjee and Roy (2014, Theorem 11.27). Equivalently, Γ X (0) admits the so-called spectral decomposition

Γ X (0) = P ΛP = k i=1 λ i p i p i .
(4)

The PC vector process is given by Y t = P X t , i.e., Y t = [Y 1t , . . . , Y kt ] where Y it = p i X t for i = 1, . . . , k. The following proposition summarizes some properties of the covariance of Y t .

Proposition 2.1. Let X t be defined by (1), λ 1 ≥ . . . ≥ λ k ≥ 0 be the eigenvalues of Γ X (0) with corresponding orthonormal eigenvectors p 1 , . . . , p k , and Y it = p i X t be the ith PC for i = 1, . . . , k. Then,

a) Var(Y it ) = p i Γ X (0)p i = λ i , b) Cov(Y it , Y jt ) = p i Γ X (0)p j = 0 when i = j, c) Cov(Y it , Y j(t+h) ) = p i Cov(X t , X t+h )p j = p i Γ X (h)p j for i, j = 1, . . . , k and h = 0.
Proof. a) and b) follow directly from (4), and c) results from Y it = p i X t .

Remark 1. Propositions 2.1(a),(b) appear in [START_REF] Anderson | An Introduction to Multivariate Statistical Analysis[END_REF] and are the particular cases of an uncorrelated process, i.e, when X t = ε t in (1). Proposition 2.1(c) shows that the autocovariances (i = j) and the cross-covariances (i = j) of the PCs are non-zero. This induces some issues discussed below in descriptive and inferential procedures of PCA in the case of time series.

Remark 2. If some eigenvalues are equal, the corresponding eigenvectors and PCs are not uniquely defined. Nevertheless, the vector space generated by these eigenvectors is unique, see e.g., Harville (1997, pages 537-538).

Remark 3. The properties of PCs discussed here are based on the eigenvalues and eigenvectors of the covariance matrix of X t . These properties are still valid for the eigenvalues and eigenvectors obtained from the autocorrelation matrix. It is quite common in practice to compute PCA based on the autocorrelation matrix, especially when the unities and the variances of the variables are different. However, this is not the case in the application problem discussed in this paper. One of the advantages of using sample PCs based on the covariance matrix is that the statistical inferences for the population PCs are easier than those of the sample PCs obtained from the correlation matrix. For a discussion of the advantages and disadvantages of using sample PCA of the covariance matrix instead of the autocorrelation matrix, see Jollife (2002, Chapter 2).

Remark 4. Let X t be defined by (1). It follows from (3) that

tr(Γ X (0)) = tr(Σ ε ) + tr( ∞ j=1 Ψ j Σ ε Ψ j ).
(5)

Let A n = n j=1 Ψ j Σ ε Ψ j , and A ij n be the (i, j)th element of A n for 1 ≤ i, j ≤ k. Then, tr( lim n→∞ A n ) = k i=1 ( lim n→∞ A ii n ) = lim n→∞ k i=1 A ii n = lim n→∞ tr(A n ). (6)
Since A n is a nonnegative definite matrix, tr(A n ) ≥ 0. Then, lim n→∞ tr(A n ) ≥ 0, and we deduce from ( 5) and ( 6) that tr(Γ X (0

)) ≥ tr(Σ ε ). Now, tr(Γ X (0)) = tr(P ΛP ) = tr(Λ) = tr(Γ Y (0)) = k i=1 λ i .
Therefore, the PCs of X t present more variability than the ones of ε t . This can lead to a wrong use of PCA technique if the time-correlation of X t is ignored.

A parametric class of models satisfying ( 1) is the k-dimensional vector seasonal autoregressive moving average (VSARMA) process with non-seasonal orders p and q, seasonal orders P and Q, and season s ∈ N -{0}. This process is defined by the difference equation

φ(B)Φ(B s )X t = θ(B)Θ(B s )ε t , (7)
where ε t is a vector white noise with E(ε t ) = 0 and Γ ε (h) given by ( 2), and B is the backward operator, i.e., BX t = X t-1 for any process X t . The matrix-valued polynomials φ(•), θ(•), Φ(•) and Θ(•) given by

φ(z) = I -φ 1 z -• • • -φ p z p , θ(z) = I + θ 1 z + • • • + θ q z q , Φ(z) = I -Φ 1 z -• • • -Φ P z P , Θ(z) = I + Θ 1 z + • • • + Θ Q z Q ,
satisfy that det(φ(z)Φ(z s )) = 0 and det(θ(z)Θ(z s )) = 0 for all z ∈ C such that |z| ≤ 1. These two conditions are known as the causality and invertibility properties, respectively. Additional conditions have to be imposed in order to obtain an identifiable model, see e.g. Brockwell and Davis (1991, page 431) and Reinsel (1997, section 2.3). In ( 7), the matrix parameters φ i 's, θ i 's, Φ i 's and Θ i 's are unknown and have to be estimated from the observed data X 1 , . . . , X n .

The VSARMA process has a short-memory correlation structure in the sense that the sequence of matrices Γ X (h) for h ∈ Z is summable. The vector seasonal autoregressive fractionally integrated moving average (VSARFIMA) process is a linear process defined by an extension of the difference equation ( 7). This process has a long-memory behaviour in the sense that the matrices Γ X (h) are only square summable, see [START_REF] Chung | Sample means, sample autocovariances, and linear regression of stationary multivariate long memory processes[END_REF]. A VSARFIMA model is used in Section 4.

As mentioned in Remark 4, when X t is time-correlated, the PCs of X t have larger variances than the ones of ε t . One way to mitigate this effect is to apply to X t a multivariate linear filter, such as the VSARMA filter before applying PCA. In this context, PCA is applied to ε t in place of X t in (7).

The VAR(1) model is the particular case of ( 7) where X t satisfies the difference equation X t = ΦX t-1 + ε t with Φ a matrix parameter. This model is widely used in modelling multivariate time series. Proposition 2.2 illustrates the effect of temporal correlation on the PCs Y t when X t is a VAR(1) process. This result can be extended to more general processes. For example, this is well-known that the VAR(p) model can be written as a VAR(1) process see e.g., Lutkepohl (2005, page 15) and Hamilton (1994, page 259).

Proposition 2.2. Let X t be a stationary VAR(1) process. Then Γ X (h) = Γ X (0)(Φ h ) and Γ Y (h) = ΛP (Φ h ) P for all h ≥ 0.

Proof. It follows from Brockwell and Davis (1991, Example 11.3.1)

that Ψ j = Φ j in (1). Then (3) implies that Γ X (h) = Γ X (0)(Φ h ) for all h ≥ 0. Since Y t = P X t , Γ Y (h) = P Γ X (h)P = P Γ X (0)(Φ h ) P = P P ΛP (Φ h ) P = ΛP (Φ h ) P for all h ≥ 0.
Remark 5. Consider the particular VAR(1) process where Φ = diag(φ 1 , . . . , φ k ) with

|φ i | < 1 for i = 1, . . . , k. Then, it results from (3) that the (i, j)th element of Γ X (h), Γ ij X (h), is given by Γ ij X (h) = ∞ l=0 φ l i Σ ij ε φ l+h j = φ h j /(1 -φ i φ j )Σ ij ε , (8)
for all h ≥ 0. Therefore,

tr(Γ Y (0)) = tr(Γ X (0)) = tr(Σ ε ) + k i=1 φ 2 i /(1 -φ 2 i )Σ ii ε . (9)
It follows from ( 9) that the variability of the PCs of X t increases as |φ i | increases, and may be much larger than the one of the PCs of ε t . Furthermore, since Γ Y (h) = ΛP (Φ h ) P , its (i, j)th element, Γ ij Y (h), is given by

Γ ij Y (h) = λ i k l=1 φ h l p li p lj , (10)
for all h ≥ 0, where

p i = [p 1i , . . . , p ki ] . Suppose that φ i = φ for i = 1, . . . , k. Then, Γ X (h) = φ h Γ X (0) = φ h /(1 -φ 2 )Σ ε for all h ≥ 0.
The eigenvectors of Γ X (h) and Σ ε are the same, while the eigenvalues of Γ X (h) are the ones of Σ ε multiplied by φ h /(1 -φ 2 ). We have Γ Y (h) = ΛP (Φ h ) P = φ h Λ . Then, when Σ ε is not diagonal, the components of X t are cross-correlated, while the components of the PCs are not, for all h ≥ 0. Observe, using (10), that the components of the PCs are generally cross-correlated when the parameters φ i 's are not all equal.

The VMA(1) model is the particular case of (7) where X t satisfies the difference equation X t = ε t + Θε t-1 with Θ a matrix parameter. Proposition 2.3 gives the expressions of Γ X (h) and Γ Y (h) when X t is a VMA(1) process. As for the VAR(1) model, this result can be extended to more complicated processes.

Proposition 2.3. Let X t be a VMA(1) process where all the eigenvalues of Θ are less than one in modulus. Then

Γ X (h) =      Σ ε + ΘΣ ε Θ if h = 0, Σ ε Θ if h = 1, 0 if h > 1, and Γ Y (h) =      Λ if h = 0, P Σ ε Θ P if h = 1, 0 if h > 1.
Proof. The expression of Γ X (h) follows from the difference equation X t = ε t + Θε t-1 .

On the other hand, Γ Y (h) = P Γ X (h)P for h ∈ Z.

Remark 6. Consider the particular VMA(1) process where Θ = diag(θ 1 , . . . , θ k ) with

|θ i | < 1 for i = 1, . . . , k. We deduce from Proposition 2.3 that Γ ij X (h) =      (1 + θ i θ j )Σ ij ε if h = 0, θ j Σ ij ε if h = 1, 0 if h > 1. Therefore, tr(Γ Y (0)) = tr(Γ X (0)) = tr(Σ ε ) + k i=1 θ 2 i Σ ii ε . (11) It follows from (11) that tr(Σ ε ) ≤ tr(Γ Y (0)) ≤ 2 tr(Σ ε ). If, θ i = θ for i = 1, . . . , k, Σ ε , Γ X (0)
and Γ X (1) have the same eigenvectors, while the eigenvalues of Γ X (0) and Γ X (1) are the ones of Σ ε multiplied by 1 + θ 2 and θ, respectively. Furthermore, in this case, we deduce from Proposition 2.3 that Γ Y (1) = θP Σ ε P = θ/(1 + θ 2 )P Γ X (0)P = θ/(1 + θ 2 )Λ. Then, the components of the PCs are not cross-correlated for all h ≥ 0.

In practice, Γ X (0) is unknown and must be estimated from a set of observations X 1 , . . . , X n of X t . The sample estimate of Γ X (0) is

ΓX (0) = 1 n n t=1 X t X t , (12)
ΓX ( 0) is symmetric and non-negative definite with spectral decomposition

ΓX (0) = BLB , (13)
where L = diag(l 1 , . . . , l k ), l 1 ≥ . . . ≥ l k ≥ 0 are the eigenvalues of ΓX (0), and B is an orthonormal matrix whose ith column b i is an eigenvector associated to l i for i = 1, . . . , k. Each eigenvalue l i is an estimate of λ i . Suppose that the eigenvalues of Γ X (0) are distinct, i.e., λ 1 > . . . > λ k . In this case, P is unique in (4). Let D = √ n(L -Λ) and G = √ n(B -P ). Under additional assumptions, Taniguchi and Krishnaiah (1987, Theorem 1) have shown that for model (1), the joint distribution of D and G converges as n tends to infinity. If X t is Gaussian, then the limiting joint distribution of D and G is normal with D and G independent and the diagonal elements of D are independent.

A major concern about using PCA is how many PCs should be selected. Several criteria have been proposed in the literature such as the eigenvalues plot of [START_REF] Jollife | Principal component analysis[END_REF] and the mean eigenvalue test of [START_REF] Perez-Neto | How many principal components? Stopping rules for determining the number of non-trivial axes revisited[END_REF]. Assume that the random variables X t are mutually independent and identically distributed with finite moments and λ 1 > . . . > λ k > 0. Fujikoshi (1980, Theorem 1) has generalized Anderson (2003, Theorem 13.5.1) to non Gaussian data and has shown that √ n(l i -λ i ) has the limiting normal distribution N(0, 2λ 2 i + κ i 4 ), where κ i 4 is the fourth-order cumulant of the ith component X it of X t for all i = 1, . . . , k. Therefore, an asymptotic confidence interval (ACI) of significance level α for λ i is given by

l i - 2l 2 i + κi 4 n z α 2 ≤ λ i ≤ l i + 2l 2 i + κi 4 n z α 2 , ( 14
)
where κi 4 is the sample estimate of κ i 4 , F (z α 2 ) = 1 -α 2 and F is the cumulative distribution function of the N(0, 1) random variable. Now, let andc

τ m = (λ 1 +• • •+λ m )/(λ 1 + • • •+λ k ) be
R m = (l 1 +• • •+l m )/(l 1 +• • •+l k ) be an estimate of τ m . Fujikoshi (1980, Theorem 3) implies that √ n(R m -τ m ) has the limiting normal distribution N(0, k i=1 T 2 i (2λ 2 i +κ i 4 )), where T i = (c i -τ m )/(λ 1 + • • • + λ k ),
i = 1 for i = 1, . . . , m, c i = 0 for i = m + 1, . . . , k. Therefore, an ACI of significance level α for τ m is R m - k i=1 T 2 i (2l 2 i + κi 4 ) n z α 2 ≤ τ m ≤ R m + k i=1 T 2 i (2l 2 i + κi 4 ) n z α 2 , ( 15
) where Ti = (c i -R m )/(l 1 + • • • + l k ).

Numerical experiments

This section presents finite sample size studies to illustrate the effect of time-correlation on the eigenvalues of Γ X (0) and on the interpretation of PCA. For this purpose, we consider VAR(1) processes with different correlation structures. The calculus and simulations were coded with R Core Team (2019, Version 3.6.2) and are available upon request. Let X t = ΦX t-1 + ε t , where matrix Σ ε is given by

Σ ε =     10 0 0 0 0 5 0 0 0 0 3 0 0 0 0 1     (16)
and the matrix parameters Φ are displayed in Table 1. The correlation structures of X t depend on Φ. The white noise model X t = ε t is denoted by Model 1. Since the matrices of the parameters of Models 2 and 3 are diagonal and Σ ε is also diagonal, it follows from (8) that the covariance matrices of these models are also diagonal and have the same eigenvectors which correspond to the natural basis of R 4 . Since all φ i 's are equal in Model 2, the eigenvalues of Γ X (0) in Models 1 and 2 are proportional, which is not the case in Models 1 and 3. Contrarily to the three first models, Models 4 and 5 present cross-correlations between the components of X t at different lags h. According to Proposition 2.2, Γ X (h) = Γ X (0)(Φ h ) for all h ≥ 0. Therefore, if the entries of Γ X (0) are nonnegative, large positive entries of Φ implies large positive cross-covariances. In this sense, Model 5 presents stronger cross-covariances than Model 4. These correlation structures may seriously affect the analysis and interpretation of the PCA. In particular, a significant impact occurs when using Models 4 and 5, which have large positive degrees of the correlations. These issues are discussed as follows.

The covariance matrices Γ X (0) of the VAR(1) Models 2 to 5 are displayed in Table 2. As expected, Model 5 displays the largest covariances. For each model, it can be seen that tr(Γ X (0)) ≥ tr(Σ ε ), as mentioned in Remark 4. Table 1.: Matrix parameters Φ of VAR(1) in Models 2 to 5.

Model 2

Model 3 0.3 0.0 0.0 0.0 0.8 0.0 0.0 0.0 0.0 0.3 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0 0.3 0.0 0.0 0.0 0.3 0.0 0.0 0.0 0.0 0.3 0.0 0.0 0.0 0.3 Model 4 Model 5 0.3 0.0 0.1 0.1 0.3 0.5 0.7 0.4 0.0 0.3 0.0 0.0 0.0 0.3 0.0 0.0 0.2 0.0 0.3 0.0 0.0 0.0 0.3 0.0 0.0 0.1 0.0 0.3 0.8 0.6 0.0 0.3 Table 2.: Covariance matrices Γ X (0) of the VAR(1) Models 2 to 5.

Model 2

Model 3 10.99 0.00 0.00 0.00 27.78 0.00 0.00 0.00 0.00 5.49 0.00 0.00 0.00 6.67 0.00 0.00 0.00 0.00 3.30 0.00 0.00 0.00 3.30 0.00 0.00 0.00 0.00 1.10 0.00 0.00 0.00 1.01 Model 4 Model 5 11.11 0.01 0.88 0.04 29.29 1.09 0.79 25.43 0.01 5.49 0.00 0.18 1.09 5.49 0.00 1.37 0.88 0.00 3.90 0.00 0.79 0.00 3.30 0.21 0.04 0.18 0.00 1.17 25.43 1.37 0.21 38.98 Table 3 shows, for each VAR(1) model, the eigenvalues λ i 's of Γ X (0) with their respective percentage of variability λ i /(λ 1 + • • • + λ 4 ). As expected, Models 1 and 2 display the same percentages since the λ i 's are proportional. Model 3 presents more variability than Models 1 and 2 because λ 1 is much larger than the other eigenvalues. Since the parameters Φ of Models 2 and 4 are close, the associated eigenvalues of Γ X (0) and their percentages of variability are similar. A very distorted case of the percentages is observed between Model 2 and Model 5. The large positive cross-covariance in Model 5 drastically increases the variability of the eigenvalues of Γ X (0), and the first PC captures almost all the variability. This is a problem of high practical relevance, for example in the context of reducing the data dimension. Now, more general VAR(1) models are considered in the study. The matrix Σ ε Table 3.: Eigenvalues of Γ X (0) of the VAR(1) Models 1 to 5 with their percentages of variability.

Model 4. Note that some autoregressive parameters are negative, which implies that the models may produce negative autocorrelations. These negative correlations may lead to different impacts on the inferential analysis compared to the previous cases. The covariance matrices Γ X (0) of the VAR(1) Models 7 to 10 are Table 4.: Matrix parameters Φ of VAR(1) Models 7 to 10.

λ 1 λ 2 λ 3 λ 4 % λ 1 % λ 2 % λ 3 %
Model 7 Model 8 0.2 0.0 0.0 0.0 -0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0 -0.5 0.0 0.0 0.0 -0.1 0.0 0.0 0.0 0.0 -0.3 0.0 0.0 0.0 0.9 Model 9

Model 10 0.4 0.1 0.3 0.1 0.6 0.3 0.6 0.03 0.0 0.8 0.4 0.0 -0.1 0.2 -0.1 0.2 0.2 0.0 0.3 0.0 0.1 -0.8 0.4 0.5 0.0 0.0 0.6 -0.4 0.2 0.0 0.1 -0.5 presented in Table 5. For each model, we have tr(Γ X (0)) ≥ tr(Σ ε ) in agreement with Remark 4. The trace of Γ X (0) represents the total variability of the PCs of X t and increases from Model 6 to Model 10.

Table 6 shows the eigenvalues λ i 's of the matrices Γ X (0) for each VAR(1) model and their respective percentage of variability λ i /(λ 1 + • • • + λ 4 ). Comparing with Table 3, it can be seen that the cross-covariances in Models 7 to 10 do not have drastic effects in the interpretation of PCA compared to Model 6. On the contrary, the percentages of variability are very stable across different correlation structures.

Samples of size n = 1000 of Models 6, 8, 9 and 10 with Gaussian innovations, were generated and the sample autocorrelation and cross-correlation functions (ACF and CCF) of the PCs were computed. The number of replications was 500. The mean of some of these quantities are displayed in Figures 1 and2 for Models 6 and 8 and Models 9 and 10, respectively. 240.04 193.95 104.52 81.18 575.20 44.82 183.86 120.35 193.95 399.10 110.20 Figure 1a) shows that the PCs are neither autocorrelated nor cross-correlated in the case of a white noise. Figure 1b) shows that the PCs may be autocorrelated and cross-correlated when the matrix parameter Φ is diagonal but the diagonal elements are not all equal. These features become more clear for Models 9 and 10. Indeed, Figure 2 shows that the full correlation structure of the data is transferred to the PCs in the case of general matrices Φ and Σ ε . These empirical evidences corroborate and illustrate Proposition 2.1.

The numerical experiments discussed in this section confirm that time-correlations in the vector X t have impacts on PCA. Therefore, it is necessary to introduce procedures that allow the use of PCA with multivariate time-correlated data. This paper suggest to pre-processing the data with a multivariate linear filter in order to whiten the data before applying PCA. This is explored in the application Section. Note that transforming the data with linear filters to attenuate the temporal structure in multivariate techniques has been also addressed in the recent work of [START_REF] Jaimungal | Consistent functional PCA for financial timeseries[END_REF], [START_REF] Greenaway-Mcgrevy | Estimating the number of common factors in serially dependent appoximate factor models[END_REF] and [START_REF] Hu | Principal volatility component analysis[END_REF].

Application to PM 10 data

PCA is used here to identify cities areas with similar PM 10 concentrations, without ignoring the time dependence of the data. We investigate whether or not the temporal correlation of the variables affects PCA and its interpretation. In general, this issue is not addressed in applied works, see e.g., Pires et al. (2008a,b). All the results in this section were also obtained with R Core Team (2019, Version 3.6.2).

The data set was collected at the automatic air quality monitoring network (AAQMN) in the Greater Vitória Region (GVR) in Brazil. The eight monitoring stations are located at urban sites of four cities in the GVR. Additionally to PM 10 concentrations, the AAQMN monitors the total suspended particles (TSP), ozone (O 3 ), nitrogen oxides (NO x ), carbon monoxide (CO), hydrocarbons (HC) and meteorological variables. The PM 10 concentrations (µg/m 3 ) were measured in eight stations, from January 2005 to December 2009. The daily averages at the eight stations constitute the time series X t which are plotted in Figure 3. The sample ACF of each component of X t are plotted in Figure 4. This figure shows a strong weekly seasonal behaviour which is expected with daily pollution data. In addition, the sample autocorrelations are positive and decrease slowly, which is typical of a long memory seasonal time series.

=           
0.27 -0.13 0.17 -0.06 -0.03 0.13 -0.01 0.02 0.02 -0.05 0.10 0.01 -0.05 -0.01 0.08 0.12 0.08 -0.05 0.07 -0.04 0.07 0.07 -0.05 0.09 0.18 -0.07 0.04 0.06 0.01 0.02 0.01 0.06 0.09 -0.01 0.04 0.01 0.04 -0.03 0.07 0.09 0.08 -0.01 0.09 -0.02 -0.06 0.09 0.00 0.08 0.06 0.02 0.02 -0.05 0.02 -0.03 0.09 0.06 0.04 0.00 0.06 0.00 -0.08 0.06 0.05 0.06

          
.

Apart from the first diagonal element, all the coefficients of Φ are quite small, which indicates that the fractional filtering giving Ẑt extracts almost all the temporal correlation of X t . Figure 5 plots the sample ACF of each component of the residual εt = Ẑt -Φ Ẑt-1 and clearly shows that these components are white noises. Note that, even if the nondiagonal autoregressive parameters are very small, they should not be ignored in the use of PCA. Now, we investigate the temporal correlation effect in the analysis and interpretation of PCA applied to PM 10 data. The sample estimate ΓX (0) of Γ X (0) is given by ( 12) and its spectral decomposition is (13). Let Γε (0) = (1/n) n t=1 εt ε t with the spectral decomposition Γε (0) = CM C , where M = diag(m 1 , . . . , m k ), m 1 ≥ . . . ≥ m k ≥ 0 are the eigenvalues of Γε (0), and C is an orthonormal matrix whose ith column c i is an eigenvector associated to m i for i = 1, . . . , k.

As addressed in Remark 3, in practice, it is more common to compute the PCs based on the eigenvectors and eigenvalues derived from standardized variables, i.e, from the correlation matrix. This is the cases when the components of the vector X t have distinct units and very different variances. The PM 10 concentrations are measured with the same unit and have similar standard deviations; the minimum and maximum standard deviations are 7.4 µg/m 3 (Cariacica) and 13.12 µg/m 3 (Laranjeiras). In addition, the percentages of cumulative variation explained by the PCs obtained from the autocorrelation and autocovariance matrices are very close. For example, the cumulative percentages for the first three PCs are 60.40, 72.24, 82.94, and 61.21 70.89 78.75 in the first and second cases, respectively. Thus, there do not seem to be any noticeable differences between the PCs from the sample correlation and the autocovariance matrices.

In Table 8, the four columns corresponding to the PCA of ΓX (0) display the eigenvectors b i 's, the eigenvalues l i 's, the proportions l

i /(l 1 + • • • + l 8 )'s and the cumulative proportions (l 1 +• • •+l i )/(l 1 +• • •+l 8 )'s for i = 1, . . . , 4.
The four columns corresponding to the PCA of Γε (0) display the eigenvectors c i 's, the eigenvalues m i 's, the proportions

m i /(m 1 + • • • + m 8 )'s and the cumulative proportions (m 1 + • • • + m i )/(m 1 + • • • + m 8
)'s for i = 1, . . . , 4. For both PCA, the main part of the variability is captured by the first PC, namely 61% for the PCA of ΓX (0) and 57% for the PCA of Γε (0). The proportions for the other PCs are quite similar for both PCA. To group the monitoring stations in classes, we select for each PC the stations with the highest factor loading in absolute value. The coefficients in bold are larger than 0.37 in absolute value. Selecting these coefficients, we retain the class CL1 : VixCentro, Ibes and Cariacica for the 1st PC of ΓX (0), the class CL2 : Laranjeiras and Carapina for the 2nd PC of ΓX (0), the class CL3 : VVCentro for the 3rd PC of ΓX (0), the class CL4 : Camburi and Sua for the 4th PC of ΓX (0), and the class CL1 : Sua, VixCentro and Ibes for the 1st PC of Γε (0), the class CL2 : Laranjeiras, Carapina and Cariacica for the 2nd PC of Γε (0), and the classes CL3 : CL4 : Camburi and VVCentro for the 3rd and the 4th PC of Γε (0), respectively. Note that four PCs are necessary in the PCA of ΓX (0) to encompass the eight stations, while three PCs are enough in the PCA of Γε (0).

Figure 6 shows the average daily profile of daily average PM 10 concentrations at the monitoring stations, grouped by the correspondent PC/CL category. Similar profiles of PM 10 concentrations are observed in all sites belonging to the same PC/CL category. However, it is clear that the associations PC/CL obtained with Γε (0) are better balanced and discriminate the data more clearly.

Following the same approach as Pires et al. (2008a,b), the number of monitoring stations that should be maintained among the eight corresponds to the maximum number of selected PCs. Based on the PCs of ΓX (0), the four stations Ibes, Laranjeiras, VVCentro and Camburi are maintained, while the analysis of the PCs of Γε (0) leads q q q q q q q 0 10 20 30 40 50 60 PC1/CL1 of Γ ^X(0) PM10 (µg/m 3 ) q q q q q q q q q q q q q q Mon Tues Wed Thurs Fri Sat Sun VixCentro Ibes Cariacica q q q q q q q 0 10 20 30 40 50 60 PC2/CL2 of Γ ^X(0) PM10 (µg/m 3 ) q q q q q q q Mon Tues Wed Thurs Fri Sat Sun Laranjeiras Carapina q q q q q q q 0 10 20 Mon Tues Wed Thurs Fri Sat Sun VVCentro q q q q q q q 0 10 20 30 40 50 60 PC4/CL4 of Γ ^X(0) PM10 (µg/m 3 ) q q q q q q q Mon Tues Wed Thurs Fri Sat Sun Camburi Suá q q q q q q q 0 10 20 30 40 50 60 PC1/CL1 of Γ ^ε ^(0) PM10 (µg/m 3 ) q q q q q q q q q q q q q q Mon Tues Wed Thurs Fri Sat Sun Suá VixCentro Ibes q q q q q q q 0 10 20 30 40 50 60 PC2/CL2 of Γ ^ε ^(0) PM10 (µg/m 3 ) q q q q q q q q q q q q q q Mon Tues Wed Thurs Fri Sat Sun Laranjeiras Carapina Cariacica q q q q q q q 0 10 20 30 40 50 60 PC3/CL3 of Γ ^ε ^(0) PM10 (µg/m 3 ) q q q q q q q Mon Tues Wed Thurs Fri Sat Sun Camburi VVCentro to retain only the three stations, Ibes, Laranjeiras and Camburi. The equipment of the others stations may be moved to alternative areas of interest to cover a larger area of the GVR. Figure 7 plots the sample ACF of the PCs of original and filtered PM 10 concentrations. Figure 7a) shows that the PCs are autocorrelated in the case of a correlated time series. Since the filtered time series εt is almost a white noise, the autocorrelations in Figure 7b) are very small.

Conclusion

This paper has investigated the effect of time-correlation on the PCA technique. It was shown that the PCs are generally cross-correlated and present more variability compared to the case of time uncorrelated data. Explicit calculations have illustrated the effect of time-correlation on the PCs when the data follow a VAR(1) and a VMA(1) model. The theoretical results were illustrated empirically through Monte Carlo simulations. It was found that large positive cross-covariance radically increases the variability of the PCs, and the first PC captures almost all the variability. Therefore, when the data are strongly time-correlated, it is recommended to apply a linear filter for whitening the data before PCA. The proposed methodology was applied to PM 10 concentrations to identify redundant air quality measurements. The PCA of the filtered data is more parsimonious and leads to retaining fewer monitoring stations.

  the fraction of the variance explained by the first m PCs, where 1 ≤ m < k, and

Figure 1 .

 1 Figure 1.: ACF and CCF plots of some sample PCs of Models 6 and 8.

Figure 2 .

 2 Figure 2.: ACF and CCF plots of some sample PCs of Models 9 and 10.

  Figure 4.: Sample ACF of the daily average of the PM 10 concentrations.

Figure 5 .

 5 Figure 5.: Sample ACF of the residuals of the fitted VSARFIMA model to PM 10 concentrations.

Figure 6 .

 6 Figure 6.: Average daily profile of PM 10 concentrations grouped by the PC/CL category.

Figure 7 .

 7 Figure 7.: Sample ACF of the PCs of original and filtered PM 10 concentrations.

Table 5 .

 5 : Covariance matrices Γ X (0) of the VAR(1) Models 7 to 10.

		Model 7			Model 8
	132.29 30.00	42.73	54.39 169.33 24.00	49.47	44.29
	30.00	58.00	33.00	70.00	24.00	77.33	31.43 116.67
	42.73	33.00	85.33	89.23	49.47	31.43	64.65	53.70
	54.39	70.00	89.23 337.25 44.29 116.67 53.70 477.78
		Model 9			Model 10

Table 8 .

 8 : PCA of original and filtered PM 10 concentrations.

	Station	1	PCA of ΓX (0) 2 3	4	1	PCA of Γε (0) 2 3	4
	Laranjeiras	-0.3002	0.7193	-0.1756	0.1460	-0.3067	0.7090	-0.0529	0.1606
	Carapina	-0.3554	-0.4004	0.2628	0.1750	-0.3536	-0.5233	0.0368	0.0669
	Camburi	-0.3472	0.1700	0.0502	0.7019	-0.3166	0.0560	0.7079	0.5055
	Sua	-0.3632	0.2163	0.0406	-0.6118	-0.3722	0.2283	-0.3546	-0.1360
	VixCentro	-0.3864	-0.2265	-0.1026	-0.1629	-0.3856	-0.0222	-0.2168	-0.2125
	Ibes	-0.3869	0.1787	0.2359	-0.2271	-0.3935	0.0625	-0.1563	0.1426
	VVCentro	-0.3055	-0.2942	-0.8391	0.0141	-0.3222	-0.0087	-0.4764	-0.7571
	Cariacica	-0.3721	-0.2766	0.3542	0.0507	-0.3669	-0.4044	-0.2652	0.2383
	Eigenvalue	4.8971	0.7744	0.6282	0.4973	4.5586	0.7462	0.6412	0.6050
	Proportion	61.22	9.68	7.85	6.22	56.98	9.32	8.01	7.56
	Cumulative	61.22	70.90	78.75	84.97	56.98	66.30	74.31	81.87

Acknowledgements

The authors are indebted to an anonymous reviewer for providing insightful comments that led to improve substantially the manuscript. Some results of this paper appear in the PhD thesis of Zamprogno (2013) under the supervision of Prof. V. A. Reisen in the Environmental Engineering PhD program at UFES (PPGEA), Brazil. The authors would like to thank CNPq, CAPES and FAPES for their financial support. Part of this paper was revised when Prof. V. A. Reisen was visiting CentraleSupélec (12/2018 to 01/2019). This author is indebted to CentraleSup ĺec for its financial support. This research was also partially supported by DATAIA Convergence Institute as part of the "Programme d'Investissement d'Avenir, (ANR17-CONV-0003) operated by Cen-traleSupélec", and by the iCODE Institute, research project of the IDEX Paris-Saclay, and by the Hadamard Mathematics LabEx (LMH) through the grant number ANR-11-LABX-0056-LMH in the "Programme d'Investissement d'Avenir".

We fit a VSARFIMA model with season s = 7 to X t . The estimator suggested by [START_REF] Reisen | A semiparametric approach to estimate two seasonal fractional parameters in the SARFIMA model[END_REF] is used to estimate the fractional parameters at the lung-run (d) and at the seasonal period s=7 (D) with the bandwidth m = n 0.5 . The estimates ( d, D) and their standard deviation (σ( d), σ( D)) are displayed in Table 7. We see that these fractional parameters are significant for each station. For each i = 1, . . . , 8, we build the series Ẑit = (1 -B) di (1 -B s ) Di X it and we fit a VSARMA model (7) to Ẑt . Following the standard methodology, we choose the orders (p, q, P, Q) with an information criterion, namely the bias-corrected Akaike information criterion (AICC), see Brockwell and Davis (1991, Section 9.2). This criterion selects a