François-Xavier Dupé
email: francois-xavier.dupe@greyc.ensicaen.fr

Luc Brun
email: luc.brun@greyc.ensicaen.fr

Hierarchical bag of paths for kernel based shape classification

Keywords: Shape, Skeleton, Support Vector Machine, Graph Kernel

Graph kernels methods are based on an implicit embedding of graphs within a vector space of large dimension. This implicit embedding allows to apply to graphs methods which where until recently solely reserved to numerical data. Within the shape classification framework, graphs are often produced by a skeletonization step which is sensitive to noise. We propose in this paper to integrate the robustness to structural noise by using a kernel based on a bag of path where each path is associated to a hierarchy encoding successive simplifications of the path. Several experiments prove the robustness and the flexibility of our approach compared to alternative shape classification methods.

Introduction

The skeleton of a 2D shape is defined as the location of the singularities of the signed distance function to the border of the shape. This structure has several interesting properties: it is thin, homotopic to the shape, invariant under rigid transformations of the plane and most importantly it has a natural interpretation as a graph. The representation of a shape by a skeletal (or shock) graph has become popular owing the good properties of this representation inherited from the properties of the skeleton. However, beside all this good properties, the skeletonization is not continuous and small perturbations of the boundary insert structural noise within the graph encoding the shape.

Several graph based methods have been proposed to compute a distance between shapes robust to such a structural noise. Sharvit et al. [START_REF] Sharvit | Symmetry-based indexing of image databases[END_REF] propose a graph matching method based on a graduated assignment algorithm. Siddiqi [START_REF] Siddiqi | Shock graphs and shape matching[END_REF] proposes to transform the shock graph into a tree and then applies a tree matching algorithm. Pellilo [START_REF] Pelillo | Matching hierarchical structures using association graphs[END_REF] uses the same tree representation but transforms the tree matching problem into a maximal clique problem within a specific association graph.

All the above graph methods operate directly on the space of graphs which contains almost no mathematical structure. This lack of mathematical structure forbids the use of basic statistical tools such as the mean or the variance. Graph kernels provide an elegant solution to this problem. Using appropriate kernels, graphs can be mapped either explicitly or implicitly into a vector space whose dot product corresponds to the kernel function. All the "natural" operations on a set of graphs which were not defined in the original graph space are now possible into this transformed vector space. In particular, graph kernels may be combined with the kernelised version of robust classification algorithms such as the Support Vector Machine (SVM).

A Graph kernel used within the shape representation framework should take into account the structural noise induced by the skeletonization process. Bunke [START_REF] Bunke | A family of novel graph kernels for structural pattern recognition[END_REF] proposes to combine edit distance and graph kernels by using a set of n prototype graphs {g 1 , . . . , g n }. Given a graph edit distance d(., .), Bunke associates to each graph g the vector φ(g) = (d(g, g 1 ) . . . , d(g, g n )). The kernel k(g 1 , g 2 ) between the two graphs g 1 and g 2 is then defined as the dot product < φ(g 1 ), φ(g 2 ) >.

Neuhaus [START_REF] Neuhaus | Edit distance based kernel functions for structural pattern classification[END_REF] proposes a similar idea by defining for a prototype graph g 0 , the kernel:

k g0 (g, g ′ ) = 1 2 d 2 (g, g 0 ) + d 2 (g 0 , g ′ ) -d 2 (g, g ′ )
, where d(., .) denotes the graph edit distance. Several graph prototypes may be incorporated by summing or multiplying such kernels. Using both Neuhaus [START_REF] Neuhaus | Edit distance based kernel functions for structural pattern classification[END_REF] and Bunke [START_REF] Bunke | A family of novel graph kernels for structural pattern recognition[END_REF] kernels two close graphs should have close edit distance to the different graph prototypes. The metric induced by such graph kernels is thus relative both to the weights used to define the edit distance and to the graph prototypes. This explicit use of prototype graphs may appear as artificial in some application. Moreover, the definite positive property of these kernels may not in general be guaranteed.

Suard [START_REF] Suard | Kernel on bag of paths for measuring similarity of shapes[END_REF] proposes to use the notion of bag of paths of finite length for shape matching. This method associates to each graph all its paths whose length is lower than a given threshold. The basic idea of this approach is that two close shapes should share a large amount of paths. A kernel between these sets should thus reflect this proximity. However, small perturbations may drastically reduce the number of common paths between two shapes (Section 3.2). Moreover, the straightforward definition of a kernel between set of paths does not lead to a definite positive kernel (Section 2.1).

This paper proposes a new definite positive kernel between set of paths which takes into account the structural noise induced by the skeletonization process. We first present in Section 2 the bag of paths approach for shape similarity. Our contributions to this field are then presented in Section 3. The effectiveness of our method is demonstrated through experiments in Section 4.

Kernels on bag of paths

Let us consider a graph G = (V, E) where V denotes the set of vertices and E ⊂ V × V the set of edges. As mentioned in Section 1 a bag of paths P of length s associated to G contains all the paths of G of length lower than s. We denote by |P | the number of paths inside P . Let us denote by K path a generic kernel between paths. Given two graphs G 1 and G 2 and two paths h 1 ∈ P 1 and h 2 ∈ P 2 of respectively G 1 and G 2 , K path (h 1 , h 2 ) may be interpreted as a measure of similarity between h 1 and h 2 and thus as a local measure of similarity between these two graphs. The aim of a kernel between bags of paths consits to agregate all these local measures between pairs of paths into a global similarity measure between the two graphs.

The max kernel

This first method, proposed by Suard [START_REF] Suard | Kernel on bag of paths for measuring similarity of shapes[END_REF], uses the kernel K path as a measure of similarity and computes for each path h 1 ∈ P 1 the similarity with its closest path in P 2 (max hj ∈P2 K path (h 1 , h j )). A first global measure of similarity between P 1 and P 2 is then defined as:

Kmax (G 1 , G 2 ) = Kmax (P 1 , P 2 ) = 1 |P 1 | hi∈P1 max hj ∈P2 K path (h i , h j ). (1)
The function Kmax (G 1 , G 2 ) is however not symmetric according to G 1 and G 2 . Suard obtains a symmetric function interpreted as a graph kernel by taking the mean of Kmax (G 1 , G 2 ) and Kmax (G 2 , G 1 ):

K max (G 1 , G 2 ) = 1 2 Kmax (G 1 , G 2 ) + Kmax (G 2 , G 1 ) . (2)
This kernel is not positive definite in general. However as shown by Haasdonk [START_REF] Haasdonk | Feature space interpretation of svms with indefinite kernels[END_REF], SVM with indefinite kernels have in some cases a geometrical interpretation as the maximization of distances between convex hulls. Moreover, experiments (section 4, and [START_REF] Suard | Kernel on bag of paths for measuring similarity of shapes[END_REF]) show that this kernel usually leads to valuable results.

The matching kernel

The non definite positiveness of the kernel K max is mainly due to the max operator. Suard [START_REF] Suard | Kernel on bag of paths for measuring similarity of shapes[END_REF] proposes to replace the kernel K path by a kernel which decreases abruptly when the two paths are different. The resulting kernel is defined as:

K matching (G 1 , G 2 ) = K matching (P 1 , P 2 ) = (3) 1 |P 1 | 1 |P 2 | hi∈P1 hj ∈P2 exp -d 2 path (h i , h j ) 2σ 2 .
where d path is the distance associated to the kernel K path and defined by:

d 2 path (h 1 , h 2 ) = K path (h 1 , h 1 ) + K path (h 2 , h 2 ) -2K path (h 1 , h 2 ).
The resulting function defines a definite positive kernel. This kernel relies on the assumption that using a small value of σ, the couple of paths with the smallest distance will predominate the others in equation 3. This kernel may thus lead to erroneous results if the distance are of the same order of magnitude than σ or if several couples of paths have nearly similar distances.

The change detection kernel

Desobry [START_REF] Desobry | An online kernel change detection algorithm[END_REF] proposed a general approach for the comparison of two sets which has straightforward applications in the design of a kernel between bags (sets) of paths. Desobry models the two sets as the observation of two sets of random variables in a feature space and proposes to estimate a distance between the two distributions without explicitly building the pdf of the two sets.

The feature space considered by Desobry is based on the normalised kernel

(K(h, h ′ ) = K path (h, h ′ )/ (K path (h, h)K path (h ′ , h ′ ))
). Using such a kernel we have h 2 K = K(h, h) = 1 for any path. The image in the feature space of our set of paths lies thus on an hypersphere of radius 1 centered at the origin (Fig. 1). Desobry defines a region on this sphere by using a single class ν-SVM. This region corresponds to the density support estimate of the unknown pdf of the set of paths [START_REF] Desobry | An online kernel change detection algorithm[END_REF].

Using Desobry's method, two set of vectors are thus map onto two regions of the unit sphere and the distance between the two regions corresponds to a distance between the two sets. Several kernels based on this mapping have been proposed:

1. Desobry proposed [START_REF] Desobry | An online kernel change detection algorithm[END_REF] to define the distance between the two spherical arcs as a contrast measure defined by:

d 2 Desobry (P 1 , P 2 ) = arccos w 1 K 1,2 w 2 w 1 w 2 arccos ρ 1 w 1 +arccos ρ 2 w 2 .
This distance is connected to the Fisher ratio (see [8, Section IV]). However, the definite positiveness of the Gaussian RBF kernel based on this distance remains to be shown.

2. Suard [START_REF] Suard | Kernel on bag of paths for measuring similarity of shapes[END_REF] proposed the following kernel:

K Suard (G 1 , G 2 ) = K Suard (P 1 , P 2 ) = ρ 1 ρ 2 hi∈P1 hj ∈P2 α 1,i K path (h i , h j )α 2,j with w 1 = (α 1,1 , . . . , α 1,|P1| ) and w 2 = (α 2,1 , . . . , α 2,|P2|
). This kernel is definite positive, but does not correspond to any straightforward geometric interpretation.

Path kernel

All the kernels between bags of paths defined in Section 2 are based on a generic kernel K path between paths. A kernel between two paths [START_REF] Kashima | Marginalized kernel between labeled graphs[END_REF] built by considering each path as a sequence of nodes and a sequence of edges. This kernel denoted K classic is then defined as 0 if both paths have not the same size and as follows otherwise:

h 1 = (v 1 , . . . , v n ) and h ′ = (v ′ 1 , . . . , v ′ p ) is classically
K classic (h, h ′ ) = K v (ϕ(v 1 ), ϕ(v ′ 1 )) |h| i=2 K e (ψ(e vi-1vi ), ψ(e v ′ i-1 v ′ i ))K v (ϕ(v i ), ϕ(v ′ i )) (4)
where ϕ(v) and ψ(e) denote respectively the vectors of features associated to the node v and the edge e. The terms K v and K e denote two kernels between respectively nodes and edge's features. For the sake of simplicity, we have used Gaussian RBF kernels between the attributes of nodes and edges (Section 4).

Hierarchical kernels

Since the main focus of this paper is a new kernel method for shape classification, the construction of skeletal graphs from shapes has been adressed using classical methods. We first build a skeleton using the method proposed by Siddiqi [START_REF] Siddiqi | Shock graphs and shape matching[END_REF].

However the graph we build from the skeleton does not correspond to the shock graph proposed by Siddiqi. Indeed, this graph provides a precise description of the shape but remains sensitive to small perturbations of the boundary. We rather use the construction scheme proposed by Suard [START_REF] Suard | Kernel on bag of paths for measuring similarity of shapes[END_REF] and Ruberto [START_REF] Ruberto | Recognition of shapes by attributed skeletal graphs[END_REF] which consists to select as node all the pixels of the skeleton which correspond to end points or junctions. These nodes are then connected by edges, each edge being associated to one branch of the skeleton. Given a skeletal graph G we valuate each of its edge by an additive weight measure and we consider the maximal spanning tree T of G. The bag of path associated to G is built on the tree T . Note that, the skeletonization being homotopic we have G = T if the 2D shape does not contain any hole.

Bag of path kernel

None of the bag of path kernels proposed by Desobry or Suard (Section 2) is both definite positive and provides a clear geometrical interpretation. We thus propose a new kernel based on the following distance:

d 2 change (P 1 , P 2 ) = arccos w 1 K 1,2 w 2 w 1 w 2 . ( 5
)
This distance corresponds to the angle α between the two mean vectors w 1 and w 2 of each region (Fig. 1). Such an angle may be interpreted as the geodesic distance between two points on the sphere and has thus a clear geometrical interpretation. Based on this distance we use the Gaussian RBF kernel:

K change (G 1 , G 2 ) = K change (P 1 , P 2 ) = exp -d 2 change (P 1 , P 2 ) 2σ 2 . ( 6
)
This kernel is definite positive since the normalized scalar product is positive definite and arccos is bijective on [0, 1]. The Gaussian RBF kernel based on this distance is thus definite positive (see [START_REF] Berg | Harmonic Analysis on Semigroups[END_REF] for further details).

Hierarchical kernel between paths

A mentioned in Section 1, the use of kernels between bags of paths within the shape matching framework relies on the assumption that the graphs associated to two similar shapes share a large amount of similar paths. This assumption is partially false since a small amount of structural noise may have important consequences on the set of paths. Let us for example, consider the small deformation of the square (Fig. 1(b)) represented on Fig. 1(c). This small deformation transforms the central node in Fig. 1(b) into an edge (Fig. 1(c)). Consequently graphs associated to these two shapes only share two paths of length 2 (the ones which connect the two corners on the left and right sides). In the same way, a small perturbation of the boundary of the shape may add branches to the skeleton(Fig. 1(d)). Such additional branches i) split existing edges into two sub edges by adding a node and ii) increase the size of the bag of path either by adding new paths or by adding edges within existing paths. The influence of small perturbations of the shape onto an existing set of paths may thus be modeled by node and edge insertions along these paths. In order to get a path kernel robust against structural noise we associate to each path a sequence of successively reduced paths, thus forming a hierarchy of paths. Our implicit assumption is that, if a path has been elongated by structural noise one of its reduced version should corresponds to the original path.

The reduction of a path is performed either by node removal or edge contraction along the path. Such a set of reduction operations is compatible with the taxinomy of topological transition of the skeleton compiled by Giblin and Kimia [START_REF] Giblin | On the local form and transitions of symmetry sets, medial axes, and shocks[END_REF]. Note that, since all vertices have a degree lower than 2 along the path these operations are well defined. In order to select both the type of operation and the node or the edge to respectively remove or contract we have to associate a weight to each node and edge which reflects its importance according to the considered path and the whole graph. Let us consider a skeletal graph G, its associated maximal spanning tree T and a path h = (v 1 , . . . , v n ) within T . We valuate each operation on h as follows:

Node removal : Let us denote by v i , i ∈ {2, . . . , n -1} the removed node of the path h. The node v i has a degree greater than 2 in T by construction.

Our basic idea consists to valuate the importance of v i by the total weight of the additional branches which justify its existence within the path h. For each neighbor v of v i not equal to v i-1 nor v i+1 we compute the weight W (v) defined as the addition of the weight of the tree rooted on v in T -{e viv } and the weight of e viv . This tree is unique since T is a tree. The weight of the node v i (and the cost of its removal) is then defined as the sum of weight W (v) for all neighbors v of v i (excluding v i-1 and v i ).

After the removal of this node the edges e vi-1vi and e vivi+1 are concatenated into a single edge in the new path h ′ . The weight of this new edge is defined as the sum of the weight of the edges e vi-1vi , e vivi+1 and the weight of the node v i (Fig. 2(a) and (b)). Edge contraction : The cost of an edge contraction is measured by the relevance of the edge which is encoded by its weight. Let us denote by e vivi+1 , i < n the contracted edge of the path h = (v 1 , . . . , v n ). In order to preserve the total weight of the tree after the contraction, the weight of the edge e vivi+1 is equally distributed among the edges of T incident to v i and v i+1 :

∀e ∈ ι(v i ) ∪ ι(v i+1 ) -{e vivi+1 } w ′ (e) = w(e) + w(e vivi+1 ) d(v i ) + d(v i+1 ) -2
where ι(v) and d(v) denote respectively the set of edges incident to v and the cardinal of this set (the vertex's degree). The symbol w(e) denotes the weight of the edge e. For example, the contraction of the edge e 2,3 in Fig. 2(a) corresponds to a cost of w(e 2,3 ) = .5. The contraction of this edge induces the incrementation of the edge's weights w(e 2,1 ),w(e 2,6 ),w(e 3,4 ) by .5/3 ≈ .16

Any additive measure encoding the relevance of a branch of the skeleton may be used as a weight. We choose to use the measure defined by Torsello [START_REF] Torsello | A skeletal measure of 2d shape similarity[END_REF] which associates to each branch of the skeleton (an thus to each edge) the length of the boundaries which contributed to the creation of this branch. Such a measure initially defined for each pixel of the skeleton is trivially additive.

Let us denote by κ the function which applies the cheapest operation on a path. The successive applications of the function κ associate to each path h a sequence of reduced paths (h, κ(h), . . . , κ D (h)) where D denotes the maximal number of reductions. Using K classic for the path comparison, we define the kernel K edit as the mean value of kernels between reduced paths of equal length. Given two paths h and h ′ , this kernel is thus equal to 0 if ||h|-|h ′ || > D. Indeed, in this case the maximal reduction of the longuest path remains longuer than the shortest one. Otherwise, ||h| -|h ′ || ≤ D, and K edit (h, h ′ ) is defined as:

K edit (h, h ′ ) = 1 D + 1 D k=0 D l=0 K classic (κ k (h), κ l (h ′ )) ( 7
)
This kernel is proportional (by a factor D + 1) to a sum of R-convolution kernels [14, Lemma 1] and is thus definite positive.

Since K classic is equal to 0 for paths of different lengths, K edit is indeed equal to a sum of kernels between reduced paths of equal length. For example, given two paths h and h ′ whose respective length is equal to 4 and 3 we have for D = 2:

K edit (h, h ′ ) = 1 3 K classic (κ(h), h ′ ) + K classic (κ 2 (h), κ(h ′ ))
4 Experiments

We used the following features for our experiments: Each node is weighted by its distance to the gravity center of the shape and each edge is assocated to a vector of two features: The first feature corresponds to the edge's weight (section 3.2). The second feature is the angle between the straight line passing through the two nodes of the edge and the principal axis of the shape. These experiments are based on the LEMS [15] database which consists of 99 objects divided into 9 classes. We defined three kernels for these experiments: The kernel K max,classic based on a conjoint use of the kernels K max (equation 2) and K classic (equation 4) has been introduced by Suard [START_REF] Suard | Kernel on bag of paths for measuring similarity of shapes[END_REF]. The kernel K change,classic based on a conjoint use of the kernels K change (equation 6) and K classic allows to evaluate the performances of the kernel K change compared to the kernel K max . Finally, the kernel K new is based on a conjoint used of the two kernels K edit and K change proposed in this paper. The kernel K classic is defined by the two parameters σ edge and σ vertex respectively used by the Gaussian RBF kernels on edges and vertices. The kernel K max,classic does not require additional parameters while K change,classic is based on a ν-SVM and requires thus the parameter ν. σ edge and σ vertex are common to all kernels. The remaining parameters have been been set in order to maximize the performances of each kernel on the experiments below.

Our first experiment compares the distance induced by each kernel k and defined as d 2 (x, x ′ ) = k(x, x)+k(x ′ , x ′ )-2k(x, x ′ ). The mean number of matches for each class is defined as follows: For each shape of the selected class we sort all the shapes of the database according to their distances to the selected shape using an ascending order. The number of good matches of the input shape is then defined as the number of shapes ranked before the first shape which belongs to a different class than the selected one. For example, the 10 nearest neighbors of a hand sorted in an ascending order are represented in Fig. 3(b), the number of good matches of each shape is indicated on the right of the figure. Note that the greater number of good match being obtained for the kernel K new . The mean number of good matches of a class is defined as the mean value of the number of good matches for each shape of the class. The different values represented in Tab. 1(a) represent the mean values of these number of good matches for the classes: hands, tools and dudes (Fig. 3(a)). As indicated by Tab. 1(a), the kernel K max,classic provides stable results but is sensitive to slight perturbations of the shapes as the ones of the class dudes and cannot handle the severe modifications of the hands. The kernel K change,classic leads to roughly similar results on the different classes. Though not presented here, the kernel K matching,classic (equation 3) gives worst results than the others kernels. This result may be explained by the drawbacks of this kernel (Section 2.2). The kernel K new always provides the best results with a good robustness to perturbation on dudes and hands.

Our second experiment evaluates performances of each kernel within a classification framework. To this end, we have trained a SVM on 5 shapes of each of the three classes: dudes, hands, tools on one side and one model of each of the 6 remaining classes on the other side. The SVM margin parameter was selected in order to maximize the number of true positive while having no false positive. Tab. 1(b) shows the number of well classified shapes for each class. The kernel K new gives the best performances especially for the hands where the two missing shapes are the more perturbed ones. The two others kernels present good results and are competitive when shapes are not strongly deformed. This experiment confirms the robustness of our kernel against perturbed shapes.

Conclusion

The bag of path approach is based on a decomposition of the complex graph structure into a set of linear objects (paths). Such an approach benefits of recent advances in both string and vectors kernels. Our graph kernel based on a hierarchy of paths is more stable to small perturbations of the shapes than kernels based solely on a bag of paths. Our notion of path's hierarchy is related to the graph edit distance through the successive rewritings of a path. Our kernel is thus related to the ones introduced by Neuhaus and Bunke.

Fig. 1 .

 1 Fig. 1. Separating two sets using one-class SVM (a). The symbols (w 1 , ρ 1 ) and (w 2 , ρ 2 ) denote the parameters of the two hyperplanes which are represented by dashed lines. Influence of small perturbations on the bag of paths ((b), (c) and (d))

reduction of edge e2, 3 Fig. 2 .

 32 Fig. 2. Different reductions of a path (in gray) within a skeletal tree.

Fig. 3 .

 3 Fig. 3. Five representative shapes of the classes hands, tools and dudes of the LEMS database (a), and (b) the 10 closest shapes from an hand using the distances induced by the kernels K max,classic (1), K change,classic (2) and K new (3).

⋆ This work is performed in close collaboration with the laboratory Cycéron and is supported by the CNRS and the région Basse-Normandie.


Carlo Correa Solci

Valdério Anselmo Reisen

Alessandro Jose Queiroz Sarnaglia

Pascal Bondon

Empirical Study of Robust Estimation Methods for PAR Models with Application to the Air Quality Area

Keywords: Robust estimation, PAR models, Outliers, PM 10 pollutant

. The daily mean particulate matter concentration (PM 10 ) data is used to illustrate the methodologies in a real application, that is, in the Air Quality area.

Introduction

Classical time series analysis generally relies on stationarity assumptions, see, e.g., [START_REF] Fuller | Introduction to Statistical Time Series[END_REF], [START_REF] Priestley | Spectral Analysis and Time Series[END_REF], [START_REF] Brockwell | Time Series: Theory and Methods[END_REF], [START_REF] Shumway | Time Series Analysis and Its Applications -With R Examples[END_REF]. Despite the broad use of stationary tools, in some cases, this requirement is too restrictive. Examples of non-stationary phenomena are unit roots, deterministic trends, heteroskedasticity, among others.

The periodic correlation (PC) or cyclostationarity property, introduced by the seminal paper of [START_REF] Gladyshev | Periodically correlated random sequences[END_REF], deserves special attention due to the fact that this phenomenon is not revealed by usual stationary tools, which may lead to a model misspecification [START_REF] Tiao | Hidden periodic autoregressive-moving average models in time series data[END_REF]. Due to this fact, special methods to identify the presence of PC in time series have been proposed, see, for example, [START_REF] Hurd | Graphical methods for determining the presence of periodic correlation[END_REF] and [START_REF] Bloomfield | Periodic correlation in stratospheric ozone data[END_REF]. PC appears in many areas of application: [START_REF] Gardner | Characterization of cyclostationary random signal processes[END_REF] investigate cyclostationarity in electrical engineering; [START_REF] Lund | Climatological time series with periodic correlation[END_REF] have found PC in climatological time series; [START_REF] Noakes | Forecasting monthly riverflow time series[END_REF] have discovered PC in time series of monthly river flows, among others.

One of the most popular models for PC is the PAR model, which is a generalization of the well-known AR model introduced by [START_REF] Box | Time Series Analysis: Forecasting and Control[END_REF] where the coefficients and orders vary periodically in time. Estimation methods for the PAR model parameters have been studied by many authors among, [START_REF] Basawa | Large sample properties of parameter estimates for periodic ARMA models[END_REF], [START_REF] Anderson | Parameter estimation for periodically stationary time series[END_REF], [START_REF] Sarnaglia | Periodic ARMA models: Application to particulate matter concentrations[END_REF].

Although the PAR model has been applied in several fields, to the best of our knowledge, it is still relatively unexplored in the air quality research field, especially in the context of contaminated data. Among the air pollutants, the particulate matter with diameter smaller than 10 µm (PM 10 ), is recognized for its effects on human health and is one of the most common and important pollution variables collected by an air quality monitoring network, see, for example, [START_REF] Reisen | Modeling and forecasting daily average PM 10 concentrations by a seasonal long-memory model with volatility[END_REF], [START_REF] Souza | Generalized additive model with principal component analysis: An application to time series of respiratory disease and air pollution data[END_REF] and references therein. Note that, the air quality data usually present asymmetric distributions and large peeks of concentrations. Classical estimators such the sample mean, variance and autocovariance functions are affected by these observations. This suggests to consider robust estimators for PAR parameters, like the ones proposed by [START_REF] Shao | Robust estimation for periodic autoregressive time series[END_REF], [START_REF] Sarnaglia | Robust estimation of periodic autoregressive processes in the presence of additive outliers[END_REF][START_REF] Sarnaglia | A robust estimation approach for fitting a PARMA model to real data[END_REF]. One issue of this paper is to fit robust PAR models to PM 10 concentrations.

To the best our knowledge, there are no empirical studies in the literature of cyclostationary processes investigating the behaviour of robust estimators under asymmetric errors and observations which can be identified as atypical or outliers. This paper aims to investigate finite sample properties of such estimators under asymmetric errors and atypical observations through a Monte Carlo study. In addition, a pollutant PM 10 concentration data set is used as an example of application since it may present observations with high levels of pollutant concentrations that may produce sample distributions with heavy tails.

The rest of the paper is organized as follows: Section 2 introduces the well-known PAR model and describes three estimation methods; Section 3 presents and discusses the results of the Monte Carlo experiment; Section 4 illustrates the use of the estimation methodologies with an application to fit a regression model with PAR errors to PM 10 concentrations; Section 5 concludes the paper.

The PAR Model and its Estimation Methods

Let {Y t }, t ∈ Z, be a stochastic process with E(Y 2 t ) < ∞, µ t = E(Y t ) and autocovariance γ t (h) = Cov(Y t , Y t-h ). Process {Y t } has PC or is a periodically stationary process with period S ∈ N + (PS S ), if µ t+S = µ t and γ t+S (h) = γ t (h), t, h ∈ Z, (1)
S being the smallest integer satisfying (1). Now, let t = rS + ν, where r ∈ Z and ν = 1, . . . , S. It follows from (1) that µ rS+ν = µ ν and γ rS+ν (h) = γ ν (h). Therefore, the autocorrelation

ρ t (h) = Corr(Y t , Y t-h ) satisfies ρ rS+ν (h) = ρ ν (h).
In addition, the partial autocorrelation (PACF) defined as

α t (h) = Corr(Y t , Y t-h |Y t-1 , . . . , Y t-h+1 ) t ∈ Z, h ∈ N + ,
see, e.g., [START_REF] Brockwell | Time Series: Theory and Methods[END_REF], is also periodic in time, i.e., α rS+ν (h) = α ν (h). Clearly, the above functions only depend on the period ν and the lag h. When they do not depend on ν, {Y t } is a standard stationary time series in the terminology of Box-Jenkins.For more details, see for example, [START_REF] Mcleod | Diagnostic checking of periodic autoregression models with application[END_REF], [START_REF] Sarnaglia | Robust estimation of periodic autoregressive processes in the presence of additive outliers[END_REF] and references therein.

The standard stationary linear model can be extended to the PS S process {Y rS+ν } via

Y rS+ν = j∈Z ψ j (ν) rS+ν-j ,
where j∈Z |ψ j (ν)| < ∞ for ν = 1, . . . , S. The model is causal when ψ j (ν) = 0 for j < 0. In the same way, the model is invertible when

j≥0 π j (ν)Y rS+ν-j = rS+ν ,
where j≥0 |π j (ν)| < ∞ for ν = 1, . . . , S, see [START_REF] Sarnaglia | Robust estimation of periodic autoregressive processes in the presence of additive outliers[END_REF] and references therein.

The PAR model is a generalization of the well-known AR process and is one of the most used models to fit a PS S time series. The PAR model is given in the following definition.

Definition 1. A zero-mean PS S process {Y rS+ν } follows a PAR(p ν ) model if it satisfies the difference equation

Y rS+ν - pν i=1 φ i (ν)Y rS+ν-i = σ ν rS+ν , (2)
where { t } is a sequence of uncorrelated random variables with E( t ) = 0, E( 2 t ) = 1 and, for each cycle ν = 1, . . . , S, φ ν = (φ 1 (ν), . . . , φ pν (ν)) is the AR coefficient vector with order p ν and σ 2 ν is the error variance.

Conditions to ensure causality of a PAR model can be derived from its vector AR representation, see, e.g., [START_REF] Sarnaglia | Robust estimation of periodic autoregressive processes in the presence of additive outliers[END_REF]. In particular, for a PAR(1) model, the causality condition is

S ν=1 φ 1 (ν) < 1.
(3)

It is assumed here that p 1 = • • • = p S = p and the following notation φ = (φ 1 , . . . , φ S ) = (φ 1 (1), . . . , φ p (1), . . . , φ 1 (S), . . . , φ p (S)) .

The Yule-Walker Estimator (YWE)

Let Y 1 , . . . , Y nS be a sample from the process {Y t }. The estimates of the periodic mean µ ν and autocovariance function γ ν (h) for ν = 1, . . . , S are, respectively, Ȳν = 1 n n-1 r=0 Y rS+ν , and

γν (h) = 1 n n-1 r=0 (Y rS+ν -Ȳν )(Y rS+ν-h -Ȳν-h ),
where Y rS+ν-h is set to zero whenever rS + ν -h < 1 or rS + ν -h > nS. Therefore, the sample ACF is

ρν (h) = γν (h) [γ ν (0)γ ν-h (0)] 1 2
.

The sample PACF can be obtained as in [START_REF] Sakai | Circular lattice filtering using Pagano's method[END_REF] and [START_REF] Shao | Computation and characterization of autocorrelations and partial autocorrelations in periodic ARMA models[END_REF].

The YWE of the PAR coefficients are obtained through the linear equations system

p i=1 φ i (ν)γ ν-i (h -i) = γ ν (h), h = 1, . . . , p, (4)
in which γ ν (h) is replaced by γν (h). The YWE of φ is defined as φ = ( φ 1 , . . . , φ S ) where, for each ν, φν = ( φ1 (ν), . . . , φp (ν)) . Asymptotics results for YWE can be derived under the following assumption.

Assumption 1. {Y t } is a zero-mean causal PAR process and

E(Y 4 t ) < ∞, t ∈ Z.
The following result is due to [START_REF] Sarnaglia | Robust estimation of periodic autoregressive processes in the presence of additive outliers[END_REF]: Under Assumption 1, the YWE estimator φ satisfies

√ n( φ -φ) ; N pS (0, G), n → ∞,
where the ; symbol denotes convergence in distribution.

Other asymptotically equivalent estimators of PAR coefficients are the Least Squares Estimator (LSE) [START_REF] Basawa | Large sample properties of parameter estimates for periodic ARMA models[END_REF] and the Maximum Likelihood Estimator (MLE) [START_REF] Lund | Asymptotics, Nonparametrics, and Time Series, chapter Modeling and Inference for Periodically Correlated Time Series[END_REF]. Due to this equivalence, these estimators will not be considered in this paper.

As well known, the estimators YWE, LSE and MLE, are not resistant in the presence of atypical observations (outliers) see, for example, [START_REF] Shao | Robust estimation for periodic autoregressive time series[END_REF], [START_REF] Sarnaglia | Robust estimation of periodic autoregressive processes in the presence of additive outliers[END_REF]. The lack of robustness of classical estimators has motivated the investigation of robust approaches in the literature. In the following subsections the robust methods introduced by [START_REF] Sarnaglia | Robust estimation of periodic autoregressive processes in the presence of additive outliers[END_REF] and [START_REF] Shao | Robust estimation for periodic autoregressive time series[END_REF] are summarized and the empirical comparison with YWE estimator is presented in the Simulation Section.

2.2

The Robust Yule-Walker Estimator (RYWE) [START_REF] Sarnaglia | Robust estimation of periodic autoregressive processes in the presence of additive outliers[END_REF] have proposed a Robust Yule-Walker Estimator (RYWE) which is based on the autocovariance function proposed in [START_REF] Ma | Highly robust estimation of the autocovariance function[END_REF]. Let y = (y 1 , . . . , y n ) ∈ R n , the robust scale estimator of y proposed by [START_REF] Rousseeuw | Alternatives to the median absolute deviation[END_REF], Q n (y), is defined as the following kth order statistic

Q n (y) = d {|y i -y j | ; 1 ≤ i < j ≤ n} (k) , (5)
where d is a constant factor to ensure Fisher-consistency and k = c 2 ≈ 0.25 n 2 , where c = [n/2]+1 is half of the size n of the vector y. For Gaussian random variables, d = 2.2191. Given a PS S time series, Y 1 , . . . , Y nS , based on (5) [START_REF] Sarnaglia | Robust estimation of periodic autoregressive processes in the presence of additive outliers[END_REF] define the robust sample (periodic) ACV function by

γν (h) = 1 4 [Q 2 n-r+1 (u ν + v ν ) -Q 2 n-r+1 (u ν -v ν )], 0 ≤ h < [(n -1)S + ν], (6)
where

u ν = (Y rS+ν-h , . . . , Y (n-1)S+ν-h ), v ν = (Y rS+ν , . . . , Y (n-1)S+ν ).
Note that, γν (h) does not posses the positive definite property. The RYWE is defined similarly to the YWE. For each ν = 1, . . . , S, replace, in (4), the theoretical ACV, γ ν (h), by its sample robust estimator in Equation 6, γν (h), and solve the resulting linear equations system

p i=1 φ i (ν)γ ν-i (h -i) = γν (h), k = 1, . . . , p, (7)
for φ 1 (ν), . . . , φ p (ν). The RYWE estimator is defined as φ = ( φ 1 , . . . , φ S ) . The white noise variance σ 2 ν can also be robustly estimated by the same argument as in Equation 4with h = 0.

Assumption 2. For any ν = 1, . . . , S, {Y rS+ν ; r ∈ Z} is a mean-zero Gaussian process with strong mixing coefficients α n satisfying: α n ≤ Cn -a , for some a > 1 and C ≥ 1.

The following result is due to [START_REF] Sarnaglia | Robust estimation of periodic autoregressive processes in the presence of additive outliers[END_REF]. Under Assumptions 1 and 2, the RYWE estimator φi (ν) satisfies φi (ν) -φ i (ν) = O p (n -1/2 ), for i = 1, . . . , p, ν = 1, . . . , S.

2.3

The Robust Least Squares Estimator (RLSE) [START_REF] Shao | Robust estimation for periodic autoregressive time series[END_REF] proposes an alternative to the conditional Least Squares Estimator (LSE). The LSE of φ can be defined as the solution of the pS-dimensional estimating equations

S n (φ ν ) = 1 σ ν n-1 r=0 rS+ν Y rS+ν-i = 0, 1 ≤ i ≤ p, 1 ≤ ν ≤ S, (8)
where the error terms are given by rS+ν

= rS+ν (φ ν ) = (Y rS+ν -p i=1 φ i (ν)Y rS+ν-i )/σ ν , 0 ≤ r < n, 1 ≤ ν ≤ S.
Asymptotic properties for the LSE have been studied by [START_REF] Basawa | Large sample properties of parameter estimates for periodic ARMA models[END_REF]. [START_REF] Shao | Robust estimation for periodic autoregressive time series[END_REF] aims to achieve robustness by replacing, in Equation ( 8), rS+ν and Y rS+ν by their robust versions ˘ rS+ν and YrS+ν , respectively, which are defined as

˘ rS+ν = ψ( rS+ν ) (9) and YrS+ν = Y rS+ν , if ˘ rS+ν = rS+ν , p i=1 φ i (ν) YrS+ν-i + σ ν ˘ rS+ν , if ˘ rS+ν = rS+ν . (10)
Therefore, the Robust LSE (RLSE) is the solution of the robustified estimating equations

Sn (φ ν ) = 1 σ ν n-1 r=0 ψ Y rS+ν -p i=1 φ i (ν)Y rS+ν-i σ ν YrS+ν-i = 0, 1 ≤ i ≤ p, 1 ≤ ν ≤ S. (11)
Shao [2008] has considered ψ(•) as the Huber type function, defined by

ψ(x) = ψ c (x) = x, if |x| ≤ c, c sign(x), if |x| > c, ( 12
)
because it is monotonic, which ensures existence and uniqueness of solution to Equation 11. Nevertheless, any odd bounded and differentiable function can be a candidate for the ψ(•) function, including the so-called redescending functions, such as Bisquare, Hampel, Generalized Gauss-Weight, among others. However, they have many zeroes, which may lead to non-optimal solutions.

Assumption 3. The marginal density function f (•) of the error rS+ν in Equation 2 is symmetric about the origin.

The following result has been derived by [START_REF] Shao | Robust estimation for periodic autoregressive time series[END_REF]: Under Assumptions 1 and 3, the RLSE φ satisfies

√ n( φ -φ) ; N pS (0, A), n → ∞,
where the covariance matrix A is given in Equation 14of [START_REF] Shao | Robust estimation for periodic autoregressive time series[END_REF].

In practice, the estimates are obtained using the following iterative procedure starting with an appropriate initial guess value for the RLSE. Suppose φ(l) represents the vector of estimates at the lth iteration. Then, at the (l + 1)th iteration, calculate the residuals

e (l) rS+ν = Y rS+ν - p i=1 φ(l) i (ν)Y rS+ν-i , 1 ≤ ν ≤ S,
where Y t = 0, t ≤ 0, estimate the white noise standard deviation at the period ν, σ ν , by

σ(l) ν = Median e (l) ν , e (l)
S+ν , . . . , e

(l) (n-1)S+ν , 1 ≤ ν ≤ S,
calculate the robust version of e (l)

rS+ν through (9), ȇ(l) rS+ν = ψ(e (l)
rS+ν ), obtain Y (l) rS+ν by (10) with ˘ rS+ν substituted for the robust residual ȇ(l) rS+ν and σ ν replaced with σ(l) ν , and evaluate the solution φ(l+1) of the robustified estimating equations in (11) replacing YrS+ν with Y (l) rS+ν and σ ν with σ(l) ν . Stop the procedure according to some convergence criterion.

Monte Carlo Study

In order to investigate the impact of atypical observations on the estimates obtained from the methods discussed previously, series of periodically stationary processes were generated with and without additive outliers. Let {X rS+ν } be defined as follows

X rS+ν = Y rS+ν + ωV rS+ν (13)
where {Y rS+ν } is a PAR model with S = 4 and coefficients given in Table 1. The parameter values were chosen to have examples of time series models with low (Model 1) and strong (Model 2) correlation dependencies, that is, Model 2 is closer to the non-causality region than Model 1. {V t } is a sequence of independent random variables with P(V t = -1) = P(V t = 1) = ξ/2 and P(V t = 0) = 1-ξ, 0 ≤ ξ < 1, Y t and V s are independent processes for all t, s and ω is the magnitude of the outlier. The sample sizes were taken as small (N = 400) and large (N = 1600), .i.e., n = 100

and n = 400 cycles, respectively, which are common sample sizes in practical situation. The initial value for RLSE was taken as the true parameter vector. The simulation study is divided in these two cases; uncontaminated and contaminated series. The contaminated series were generated from model in (13) with the following specifications: ω = 7 and ξ = 0.01. The effect of the normality departure in the white noise sequence was also studied by generating the random variables t such that t ∼ N (0, 1) and √ 2 t + 1 ∼ χ 2 (1) . In both cases, E( t ) = 0 and E( 2 t ) = 1. For each of these scenarios, 1000 replicates of {Y t } were generated to compute the mean of the empirical Bias and Root Mean Squared Error (RMSE). For the RLSE, according to [START_REF] Shao | Robust estimation for periodic autoregressive time series[END_REF], c = 3.06 was fixed in the Huber function (Equation 12) such that, residuals greater than 3.06 (in absolute value) are regarded as outliers. Other model configurations were also considered in the simulation study such as heavy-tailed distributions, different period lengths, sample sizes and coefficient values and other outlier magnitudes. However, in general, the results leaded to similar conclusions and are not displayed here, but they are available upon request.

Table 1: Parameters of PAR(1) models used in the simulation.

Parameter

Model 1 Model 2 φ 1 (1) 0.9 1.5 φ 1 (2) 0.8 0.8 φ 1 (3) 0.7 1.2 φ 1 (4) 0.6 0.5 λ 0.3024 0.7200

Tables 2 and3 display the Bias and RMSE for Models 1 and 2, respectively. For illustration purpose, the empirical distributions of

√ n( φi (ν) -φ i (ν)), √ n( φi (ν) -φ i (ν)) and √ n( φi (ν) -φ i (ν))
for Model 1, with n = 400, under the uncontaminated and contaminated Gaussian scenarios, are presented in Figures 1 and2, respectively. The same plots for the uncontaminated case with asymmetric errors are displayed in Figure 3. For Model 1 (Table 2), under the uncontaminated Gaussian case, the reduction of the Bias and RMSE when increasing the sample size suggests that all estimators are consistent. It is observed that, in this scenario, YWE and RLSE present better results. The findings for the asymmetric uncontaminated case show that the RYWE is extremely affected by skewness of the data, presenting a persistent Bias which does not seem to vanish by increasing the sample size. The results for heavy tailed errors are similar, which show that the RYWE method is very sensitive to departures from normality. This indicates that the normality requirement in Assumption 2 is crucial to ensure asymptotic properties of the RYWE. The YWE and RLSE do not seem to be affected by non-Gaussian errors (both asymmetric and heavy tailed). This gives empirical evidence that the technical requirement of symmetric errors (Assumption 3) may be over restrictive to ensure asymptotic normality of the RLSE. As expected, atypical observations increase the Bias and RMSE of the YWE. Under normally distributed errors, both RYWE and RLSE show robustness with Bias and RMSE almost unchanged with the presence of outliers. For asymmetric errors, the RLSE is the only one which presents good performance.

The conclusions for Model 2 (Table 3) are similar to the previous case. It is worth noting that in this stronger dependence scenario, there is a overall reduction of the Bias and RMSE. Another remarkable fact is that, in this case, the RYWE does not seem to be strongly affected by asymmetric errors in both uncontaminated and contaminated scenarios.

From Figure 1, it can be seen that the empirical distributions are virtually the same and they have shape very close to the N(0,1) distribution. This corroborates the asymptotic results of the standardized estimators even for a small sample size. However, the scale of the RYWE distribution is slightly greater than of those of the YWE and RLSE. Figure 2 shows the robustness to outliers of RYWE and RLSE methods under Gaussian errors, while the distribution of YWE is shifted to the left due to the well-known memory loss property. Its scale is also increased as a result of the contamination. Figure 3 illustrates the prominent shift to the right of the RYWE distribution The PM 10 data set corresponds to daily average concentrations from January 1st, 2014 to December 29th, 2015 which kept the sample size multiple of the natural choice to the period length S = 7. Due to skewness and some evidences of time varying variance, the natural logarithm transformation (log) was used and the plot of the log(PM 10 ) is displayed in Figure 4. From this figure, one can see large peeks of PM 10 concentration which may be viewed here as outliers and, as mentioned previously, these high levels can provoke serious damage to some statistics, such as the mean and the standard deviation and, therefore, may affect the sample correlation structure of the series, causing misleading results. The existence of any outlier's effect will be discussed in the estimation parameter model (next subsection). It can also be seen the presence of sinusoidal deterministic trends. Analysis of the periodogram (Figure 5) corroborates to this result and indicates that the frequency 2/N , corresponding approximately to a yearly cycle, has a large contribution to the overall variance of the data. The high frequency peaks of the periodogram correspond to weekly periodicity and, according to the daily periodic boxplots displayed in Figure 6, they can be explained by a level decrease in the weekends. This is an expected finding due to fact that the traffic and civil construction decrease in the region in the weekend days. The above preliminary analysis of the series suggests that a deterministic trend must be firstly removed from log(PM 10 ) before further analysis and this is discussed in the next subsection in which a linear model with errors following a PAR model is fitted to the series.

Estimated Model

According to the previous statistical analysis of the log(PM 10 ) series, the following model is suggested here to fit the data log(PM 10,t ) = µ + α 1 sat t + α 2 sun t + β 1 cos t +β 2 t + Y t ;

(14)

Y t = pt i=1 φ i (t)Y t-i + σ t t , (15)
q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0.1 0.2 0.3 0.4 0.5 0.00 0.10 0.20 0.30

Frequency

Figure 5: Periodogram of the log(PM 10 ) time series. q q q q q q q q q Wed Thu Fri Sat Sun Mon Tue 2.5 3.0 3.5 4.0 with the sinusoidal covariate: cos t = cos( 2πt 365.25 ), t = 1, . . . , 728; the linear term t; and a "day of the week" factor with the levels: Week (the reference level); Saturday (represented by the dummy variable sat t which takes value 1 for Saturdays); and Sunday (sun t which takes value 1 for Sundays) and S = 7. The above model means that in the business days the regular level of log(PM 10 ) is µ, on Saturdays it suffers an increase of α 1 and on Sundays it is increased by α 2 and it has a long-run cyclic trend, represented by cos t and a linear term t. The terms sin t and cos t turned out to be statistically insignificant and they were removed from the model.

The model in Equations 14 and 15 will be fitted based on following two steps procedure: (1) the linear model in ( 14) will be estimated through the ordinary least squares procedure; and (2) the PAR model in (15) will be fitted to the residuals of the linear model in step (1), where the AR orders p 1 , . . . , p S will be identified through the Schwartz Information Criterion (BIC) proposed by [START_REF] Schwarz | Estimating the dimension of a model[END_REF] and adapted to the periodic scenario by [START_REF] Mcleod | Diagnostic checking of periodic autoregression models with application[END_REF].

At the first step, the linear model in Equation 14was fitted and the estimated coefficients are displayed in Table 4. As expected, there were negative effects of Saturday and Sunday, which led to a decrease of log(PM 10 ) levels during the weekends. The BIC criterion was used to identify the order of the model (see [START_REF] Sarnaglia | Robust estimation of periodic autoregressive processes in the presence of additive outliers[END_REF] for more details) and the results are displayed in Table 5. In order to keep consistency with the simulation study, c = 3.06 was fixed in the Huber function (Equation 12). Note that the PAR Finally, for all models, the residuals have not passed the Jarque-Bera normality test [START_REF] Jarque | Efficient tests for normality, homoscedasticity and serial independence of regression residuals[END_REF], presenting p-values < 0.05 which is an expected result due to the skewness revealed in the data.

Conclusions

This paper reviews different estimation methodologies for PAR models. More specifically, the method considered are: the so-called YWE [START_REF] Mcleod | Diagnostic checking of periodic autoregression models with application[END_REF], the RYWE [START_REF] Sarnaglia | Robust estimation of periodic autoregressive processes in the presence of additive outliers[END_REF] and the RLSE [START_REF] Shao | Robust estimation for periodic autoregressive time series[END_REF]. The finite sample performance of these methods was compared through a Monte Carlo experiment. The performance of RLSE is remarkably good under uncontaminated and contaminated scenarios, even under asymmetric errors, which violates Assumption 3. The RYWE is quite resistant to outliers, however it has a poor performance under asymmetric errors, mainly under weak correlation scenarios. As expected, YWE empirical distribution is resistant to departures from normality, however this estimator is completely affected by the presence of outliers. In order to illustrate the methodologies considered in this paper, the daily mean PM 10 concentrations collected at the air quality monitoring station, located at Enseada do Suá, ES, Brazil, was considered as an application. The estimation and modelling results revealed outlier effects on the estimates.

Figure 1 :Figure 2 :

 12 Figure 1: Empirical Distributions of √ n( φi (ν) -φ i (ν)) (blue lines), √ n( φi (ν) -φ i (ν)) (green lines) and the √ n( φi (ν) -φ i (ν)) (red lines)for Model 1 with ω = 0, n = 400 and normal errors.

Figure 3 :Figure 4 :

 34 Figure 3: Empirical Distributions of √ n( φi (ν) -φ i (ν)) (blue lines), √ n( φi (ν) -φ i (ν)) (green lines) and the √ n( φi (ν) -φ i (ν)) (red lines)for Model 1 with ω = 0, n = 400 and asymmetric errors.

Figure 6 :

 6 Figure 6: Daily box-plots of the log(PM 10 ) time series.

Figure 7 :Figure 8 :Figure

 78 Figure 7: ACF of the residuals of the YWE fit.

Table 2 :

 2 Bias and RMSE for Model 1 and outliers with probability ξ = 0.01.

Table 3 :

 3 Bias and RMSE for Model 2 and outliers with probability ξ = 0.01.[START_REF] Reisen | An M-estimator for the long-memory parameter[END_REF] and references therein. The data set considered in this paper is the pollutant Particulate Matter with diameter smaller than 10 µm (PM 10 ), measured hourly, in µg/m 3 , collected at the station located in Enseada do Suá area.

					YWE	RYWE	RLSE
	ω	t	n	φ 1 (ν)	Bias RMSE	Bias RMSE	Bias RMSE
				1.5	-0.009 0.055	-0.009 0.100	0.000	0.055
			100	0.8 1.2	-0.004 0.033 -0.006 0.040	-0.005 0.050 -0.008 0.066	-0.004 0.034 -0.006 0.040
		N (0, 1)		0.5 1.5	-0.008 0.033 -0.003 0.026	-0.008 0.043 -0.003 0.037	-0.008 0.033 -0.001 0.026
			400	0.8 1.2	-0.001 0.016 -0.002 0.019	-0.001 0.021 -0.002 0.027	-0.001 0.016 -0.002 0.020
	0			0.5 1.5	-0.002 0.015 -0.009 0.055	-0.003 0.018 0.025 0.098	-0.003 0.015 0.001 0.043
			100	0.8 1.2	-0.007 0.036 -0.007 0.040	0.011 0.019	0.048 0.069	-0.005 0.028 -0.005 0.033
		χ 2 (1) -1 √ 2		0.5 1.5	-0.009 0.033 -0.002 0.026	-0.004 0.039 0.032 0.050	-0.006 0.027 0.000 0.021
			400	0.8 1.2	-0.001 0.016 -0.002 0.019	0.016 0.023	0.026 0.035	-0.001 0.013 -0.002 0.015
				0.5	-0.003 0.015	0.003	0.017	-0.002 0.012
				1.5	-0.174 0.246	-0.023 0.122	-0.014 0.061
			100	0.8 1.2	-0.047 0.076 -0.079 0.125	-0.011 0.058 -0.013 0.079	-0.014 0.040 -0.013 0.048
		N (0, 1)		0.5 1.5	-0.029 0.056 -0.167 0.189	-0.010 0.046 -0.020 0.052	-0.013 0.038 -0.013 0.032
			400	0.8 1.2	-0.042 0.050 -0.080 0.093	-0.007 0.025 -0.015 0.038	-0.007 0.019 -0.008 0.022
	7			0.5 1.5	-0.023 0.032 -0.180 0.253	-0.006 0.021 0.022 0.120	-0.007 0.018 -0.013 0.054
			100	0.8 1.2	-0.049 0.081 -0.082 0.126	0.011 0.016	0.057 0.083	-0.008 0.032 -0.009 0.036
		χ 2 (1) -1 √ 2		0.5 1.5	-0.029 0.053 -0.172 0.193	-0.003 0.041 0.022 0.054	-0.010 0.028 -0.010 0.025
			400	0.8 1.2	-0.041 0.050 -0.083 0.097	0.015 0.012	0.028 0.039	-0.005 0.015 -0.007 0.018
				0.5	-0.024 0.033	0.003	0.020	-0.005 0.014
	Reisen et al. [2018],

Table 4 :

 4 Estimated coefficients of the linear model.

	Parameter	µ	α 1	α 2	β 1	β 2
	Estimate	3.2650 -0.0921 -0.2579 0.0640 0.0003

Acknowledgements

The results in this paper are part of the Master thesis of the first author in the PPGEA-UFES under supervision of the second and third authors. The authors gratefully acknowledge partial financial support from FAPES/ES, CAPES/Brazil and CNPq/Brazil. This paper was revised when Prof. Valderio Reisen was visiting CentraleSupélec (01 to 03/2018). This author is indebted to CentraleSupélec for its financial support. The authors are grateful to the Editor and the referee for the time and efforts in providing very constructive and helpful comments that have led to clarify and substantially improve the quality of the paper.

4 An Application to the Air Quality Area (the PM 10 Data)

The application is based on a data set (air pollutant variables) collected at Automatic Air Quality Monitoring Network (RAMQAr) in the Great Vitória Region GVR-ES, Brazil, which is composed by nine monitoring stations placed in strategic locations and accounts for the measuring of several atmospheric pollutants and meteorological variables in the area. GVR is comprised of seven cities with a population of approximately 1.9 million inhabitants in an area of 2319 km 2 . The region is situated along the South Atlantic coast of Brazil (latitude 20°19 15 S, longitude 40°20 10 W) and has a tropical humid climate, with average temperatures ranging from 24 °C to 30 °C. The estimates of the PAR coefficients provided by YWE, RYWE and RLSE methods are given in Table 6. Based on these results, it is clear the presence of Periodic Correlation in the data, since the AR coefficients and orders are not constant over the seasons. In general, the methods selected different orders and presented quite different coefficient estimates. This indicate that the high levels of the pollutant were stronger enough to provoke changes in the parameter estimates, that is, this reveals that the high levels of the pollutant PM 10 presented the effects of additive outliers according to the discussion presented in the Simulation Section. 1 0.614 0.482 0.630 0.529 0.595 0.361 0.474 2 0.000 0.000 -0.107 -0.143 -0.258 0.000 0.000 3 0.000 0.000 0.451 0.000 0.000 0.000 0.000 4 0.000 0.000 -0.374 0.000 0.000 0.000 0.000 φi (ν) 1 0.661 0.513 0.479 0.376 0.638 0.293 0.522 2 0.000 0.000 -0.018 0.000 -0.167 0.000 0.000 3 0.000 0.000 0.271 0.000 0.000 0.000 0.000 4 0.000 0.000 -0.240 0.000 0.000 0.000 0.000

The fitting performance will be accessed through the in-sample Root Mean Square Error (RMSE), Mean Absolute Percentage Error (MAPE), symmetric MAPE (sMAPE) and Median of Absolute Deviation (MAD). The RMSE and the MAD are well-known and, for a discussion of MAPE and sMAPE quantities see [START_REF] Flores | A pragmatic view of accuracy measurement in forecasting[END_REF]. The results are presented in Table 7. As can be seen, the RLSE are RYWE very competitive by presenting very similar results and they are slightly smaller than YWE method. This may corroborate the previous discussion related the effect of high level concentrations on the model estimation. It can be seen that all the models were able to fully explain the correlation structure of the data, despite the eventual outliers effect. Based on the ACF of the residuals, the three estimation methods are comparable since all the estimated residuals look like a white noise process.