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

Zambon Monte

Glaura Da Conceição Franco

Adriano Márcio Sgrancio

Fábio Alexander Fajardo Molinares

Pascal Bondon

Flávio Augusto Ziegelmann

Bovas Abraham

Robust estimation of fractional seasonal processes: modeling and forecasting daily average SO 2 concentrations

Keywords: Forecasting, Robust periodogram, Large peak, Outliers, SO 2 pollutant, Long-memory

This paper deals with the estimation of seasonal long-memory time series models in the presence of 'outliers'. It is long known that the presence of outliers can lead to undesirable effects on the statistical estimation methods, for example, substantially impacting the sample autocorrelations. Thus, the aim of this work is to propose a semiparametric robust estimator for the fractional parameters in the seasonal autoregressive fractionally integrated moving average (SARFIMA) model, through the use of a robust periodogram at both very low and seasonal frequencies. The model and some theories related to the estimation method are discussed. It is shown by simulations that the robust methodology behaves like the classical one to estimate the long-memory parameters if there are no outliers (no contamination). On the other hand, in the contaminated scenario (presence of outliers), the standard methodology leads to misleading results while the proposed method is unaffected. The methodology is applied to model and forecast

Introduction

Sulfur dioxide (SO 2 ) is part of a group of highly reactive gases known as "oxides of sulfur". The main sources of SO 2 emissions into the atmosphere come from fossil fuels (coal and petroleum products), power and other industrial plants. Other emission sources of SO 2 are as follows: industrial processes such as steel and mining, burning fuels containing high sulfur level due to transportation vehicles such as locomotives and large ships, pulp industry, and natural sources such as volcanic emissions (Andersson et al. [START_REF] Andersson | Composition and evolution of volcanic aerosol from eruptions of Kasatochi, Sarychev and Eyjafjallajükull in 2008-2010 based on Caribic observations[END_REF]). High levels of pollution emissions can result in peaks of SO 2 concentration in the atmosphere.

The high SO 2 concentrations are associated with several effects on population health, as for instance the increase of chronic bronchitis (Kanaroglou et al. [START_REF] Kanaroglou | Estimation of sulfur dioxide air pollution concentrations with a spatial autoregressive model[END_REF]). The exposure to SO 2 , even for short periods of time (ranging from 5 minutes to 24 hours), may cause respiratory problems, including bronchoconstrictor and increased asthma symptoms. Consequently, these can cause an increase in hospital admissions of the population at high risk, such as elderly, children and asthmatics.

Several statistical methods have been employed to study the impacts of SO 2 pollution on the environment, such as: principal component and cluster analysis (Cheng and Lam [START_REF] Cheng | Synoptic typing and its application to the assessment of climatic impact on concentrations of sulfur dioxide and nitrogen oxides in Hong Kong[END_REF]); artificial neural networks for predicting sulfur dioxide concentrations (Saral and Ertürk [START_REF] Saral | Prediction of ground level SO 2 concentration using artificial neural networks[END_REF]) and other pollutants (Brunelli et al. [START_REF] Brunelli | Two-days ahead prediction of daily maximum concentrations of SO 2 , O 3 , PM 10 , NO 2 , CO in the urban area of Palermo, Italy[END_REF]); multiple regression models to investigate the influence of emission sources and meteorological conditions on SO 2 pollution (Luvsan et al. [START_REF] Luvsan | The influence of emission sources and meteorological conditions on SO 2 pollution in Mongolia[END_REF]); among others. In terms of forecasting, statistical models based on multiple regression and time series tools, such as the autoregressive integrated moving average (ARIMA) model have been widely used.

Many time series data exhibit a seasonality pattern. For pollutants, seasonal variation is often associated with the changes in meteorological parameters. Furthermore, the reduction of emission levels during weekends can also induce seasonality in the series. Thus, it is very important to consider statistical tools which take into account the seasonality effect.

More recently, several authors have studied time series with long-range dependence (often named long memory property). In time domain analysis, the long-memory dependence is generally characterized by a slow decay of the autocorrelation function. Granger and Joyeux [START_REF] Granger | An introduction to long-memory times series models and fractional differencing[END_REF] and Hosking [START_REF] Hosking | Fractional differencing[END_REF] proposed the autoregressive fractionally integrated moving average (ARFIMA) model to describe the long-memory feature of a time series. The ARFIMA process is an extension of the ARIMA, where the parameter of integration d assumes fractional values. The methods proposed for estimating the parameters of the ARFIMA model are classified either as parametric or semiparametric. Parametric methods consist of simultaneous estimation of model parameters, usually by maximum likelihood. On the other hand, semiparametric estimation is performed in two steps: firstly the parameter d is estimated, and, in a second step, using the estimated d from the first step, the autoregressive and moving average parameters are estimated. The most popular semiparametric estimator was proposed in Geweke and Porter-Hudak [START_REF] Geweke | The estimation and application of long memory times series model[END_REF]. For a recent review on this subject, see Palma [START_REF] Palma | Long-Memory time series: theory and methods[END_REF], for instance. Modifications of this estimator have been developed by Reisen [START_REF] Reisen | Estimation of the fractional difference parameter in the Arima(p, d, q) model using the smoothed periodogram[END_REF], Lobato and Robinson [START_REF] Lobato | Averaged periodogram estimation of long memory[END_REF], Delgado and Robinson [START_REF] Delgado | Optimal spectral bandwidth for long memory[END_REF], Velasco [START_REF] Velasco | Non-Gaussian log-periodogram regression[END_REF], among others. More specifically, Sena et al. [START_REF] Sena | Correlated error in the parameters estimation of the ARFIMA model: a simulated study[END_REF], Molinares et al. [START_REF] Molinares | Robust estimation in longmemory processes under additive outliers[END_REF] and Reisen et al. [START_REF] Reisen | An m-estimator for the longmemory parameter[END_REF] investigated the estimator under various model specifications, such as the presence of non-gaussian errors and outliers.

Due to the important features of the ARFIMA model, it has often been used to analyze environmental problems. As an example, Iglesias et al. [START_REF] Iglesias | Data analysis using regression model with missing observations and long memory[END_REF] adopted an ARFIMA model as the underlying generating process of pollutant concentrations which present long memory and missing values. These phenomena are widely found in time series of different areas of interest.

An extension of the ARFIMA to handle time series with seasonality is the SARFIMA model, which has already been addressed in the literature by several authors. Porter-Hudak [START_REF] Porter-Hudak | An application of the seasonal fractionally differenced model to the monetary aggregates[END_REF] applied the seasonal fractionally differenced model to study the monetary aggregates of the United States. Reisen et al. [START_REF] Reisen | Modeling and forecasting daily average PM 10 concentrations by a seasonal long-memory model with volatility[END_REF] estimated the fractional and seasonal parameters of SARFIMA models by means of a semiparametric procedure, considering a nonconstant conditional error variance. Reisen et al. [START_REF] Reisen | A semiparametric approach to estimate two seasonal fractional parameters in the Sarfima model[END_REF] studied the properties of the SARFIMA model when the data exhibits one and two seasonal periods and short-memory components. Ye et al. [START_REF] Ye | Improving estimation of the fractionally differencing parameter in the Sarfima model using tapered periodogram[END_REF] proposed a new method to estimate the fractional difference parameter in the SARFIMA model using tapered periodogram. Hsu and Tsai [START_REF] Hsu | Semiparametric estimation for seasonal long-memory time series using generalized exponential models[END_REF] proposed a semiparametric estimation method for seasonal long-memory time series, by considering a generalized exponential model in the frequency domain. Chan and Tsai [START_REF] Chan | Inference of seasonal long-memory aggregate time series[END_REF] have studied the autocorrelation structure and the spectral density function of SARFIMA processes in aggregated discrete-time process, using the Whittle estimator. Tsai at al. [START_REF] Tsai | Inference of seasonal long-memory time series with measurement error[END_REF] evaluated the properties of the Whittle likelihood estimation in SARFIMA models in the presence of measurement errors.

Furthermore, environmental time series often contain unusual observations influenced by external events that can cause changes in their dynamics. These observations transient or permanent changes are known in the literature as outliers (sometimes referred to as aberrant or atypical observations) and, depending on their nature, their effects on inference can be substantial. For example, the presence of outliers can increase the estimated variance of the stochastic process, which implies a decrease in the estimated autocorrelations and consequent loss of information about the autocorrelation structure of the process see, for example, Chan [START_REF] Chan | A note on time series model specification in the presence outliers[END_REF][START_REF] Chan | Outliers and financial time series modelling: a cautionary note[END_REF]. Although long-memory models in the presence of outliers has recently been a subject of much interest to some researchers, especially in the areas of economics and finance (Beran [5]; Franses et al. [START_REF] Franses | Long memory and level shifts: reanalyzing inflation rates[END_REF]; Tolvi [START_REF] Tolvi | Long memory and outliers in stock market returns[END_REF]), very few papers have been devoted to this topic in the environmental field. One of the contributions of this paper is to fill this gap. To this purpose, it concatenates longmemory modeling (with more than one fractional parameter), robust estimation and pollution data.

Based on the above discussion, the main goal of this paper is to propose a robust semiparametric estimator for both non-seasonal and seasonal long-memory parameters in SARFIMA models based on a robust autocovariance function estimator. Some model and estimation properties are discussed and a small sample size investigation is conducted to show the method's performance under contaminated and uncontaminated SARFIMA models. In addition, the method is applied to model and forecast SO 2 concentrations, measured at the air quality automatic monitoring network (AQAMN) of the Greater Vitória Region (GVR), ES, Brazil. This series exhibits seasonality, long memory behaviour and high levels (or large peaks) of SO 2 concentrations. These observations may produce sample densities with heavy tails and it can strongly influence sample functions such as the standard mean, the covariance and the periodogram. Since the estimation of time series models is connected with these sample functions, the final estimated model can be strongly affected by the large peaks of concentration.

This paper is organized as follows. Section 2 introduces the model, discusses its properties and summarizes its estimation methods. Section 3 presents a simulation study and Section 4 deals with the analysis, modeling and forecasting of SO 2 concentrations. Some conclusions are drawn in Section 5.

Model definitions and a robust parameter estimation method

Model definition and properties

Let X t ≡ {X t } t∈Z be a zero-mean process with infinite MA representation X t = Ψ(B) t where

Ψ(B) = ψ(B)∇ α 1 ,...,α L ω 1 ,...,ω L (B) = ψ(B) L ι=1 (1 -2 cos ω ι .B + B 2 ) α ι , (1)
{ t } t∈Z is a white noise with E( t ) = 0; Var( t ) = σ 2 , ω ι ∈ (-π, π]; B is the backward shift operator satisfying B k Y t = Y t-k for any process {Y t } t∈Z , α ι = -d ι , d ι ∈ R (d ι > -1
), (d ι 0) and d ι is defined as the fractionally differencing parameter. The following assumptions are made for the power expansion of Ψ(B):

(A1) The function ψ(z) is analytic inside and on the unit circle, |z| ≤ 1 for all z C, that is, the coefficients

ψ j , j = 0, 1, ..., of ψ(z) satisfy ∞ j=0 ψ 2 j < ∞. (A2) |d ι | < 1 2 ∀ ι = 1, ..., L.
Under assumptions A1 and A2, X t is a stationary process with spectral density of the form:

f X (ω) = g ψ (ω)|ω| -2d L ι=1 ξ ι j=1 |ω -ω ι j | -2d ι , (2)
where ω ∈ (-π, π] and g ψ (ω) is a continuous function, bounded above and away from zero and ω ι j 0 are poles for j = 1, ..., ξ ι , ι = 1, ..., L. The autocovariance function of X t behaves like γ X (h) ∼ K j 2d ι -1 cos( jω) as h → ∞ and K is a constant that does not depend on h. See, for example, Giraitis and Leipus [START_REF] Giraitis | A generalized fractionally diferencing approach in long-memory modeling[END_REF], Palma [START_REF] Palma | Long-Memory time series: theory and methods[END_REF], Arteche [START_REF] Arteche | Semiparametric robust test on seasonal or cyclical long-memory time series[END_REF], Arteche and Robinson [START_REF] Arteche | Semiparametric inference in seasonal and cyclical long-memory processes[END_REF], Reisen et al. [START_REF] Reisen | A semiparametric approach to estimate two seasonal fractional parameters in the Sarfima model[END_REF] and references therein.

For suitable choices of the fractionally differencing parameters d ι , ι = 1, ..., L, X t may have a finite number of zeros or singularities of order d 1 ,...,d L on the unit circle which allows the modeling of long and short memory data containing seasonal periodicities.

A stationary long-memory time series with memory parameter ϕ ∈ (-0.5, 0.5) has autocovariance function γ(h) and spectral density f (ω) satisfying, respectively,

(A3) γ(h) ∼ h 2ϕ-1 L 1 (h), as h -→ ∞, (3)
(A4)

f (ω) ∼ |ω| -2ϕ L 2 (ω), as ω -→ 0. ( 4
)
where L 1 and L 2 are slowly varying functions, L 1 at infinity and L 2 at zero (see, for example, Taqqu [START_REF] Taqqu | Fractional brownian motion and long-range dependence[END_REF]). When 0 < ϕ < 1/2, the autocovariance is not absolutely summable and the spectral density becomes unbounded at zero frequency. The properties given by Assumptions 3 and 4 can also be extended for any frequency ω ∈ (-π, π], as here discussed (see Remark 2 for the spectral density function) and in Reisen et al. [START_REF] Reisen | A semiparametric approach to estimate two seasonal fractional parameters in the Sarfima model[END_REF] among others. Now, let ψ(B) =

θ p (B)
φ q (B) where θ q (z) and φ p (z) are polynomials with order p and q, respectively. It is assumed that these polynomials have no common roots and satisfy the conditions φ p (z), θ q (z) 0, for all z C, such that |z| ≤ 1. Then, X t becomes the ARUMA(p, d 1 ,..., d L , q) model introduced by Giraitis and Leipus [START_REF] Giraitis | A generalized fractionally diferencing approach in long-memory modeling[END_REF]. Let now d ι , ι = 1, ..., L, satisfying the following assumption, (A5)

|d ι | <        1/2, 0 < ω < π, 1/4, ω = 0, π.
Under the conditions on φ p (z) and θ q (z) and Assumption 5, Theorem 2 in Giraitis and Leipus [START_REF] Giraitis | A generalized fractionally diferencing approach in long-memory modeling[END_REF] states that the process X t is causal, invertible and has the unique stationary solution

X t = θ q (B) φ p (B) ∇ -d 1 ,...,-d L ω 1 ,...,ω L t = ∞ j=0 ψ j ∇ -d 1 ,...,-d L ω 1 ,...,ω L t-j . (5)
Additionally, the authors show that the spectral density is given by

f X (ω) = |θ q (e -iω )| 2 φ p (e -iω )| -2 f ∇ (ω) (6)
where

f ∇ (ω) = σ 2 2π |∇ -d 1 ,...,-d L (e -iω )| 2 and the autocovariance is γ X (h) = L j=1 |θ q (e -iω j )| 2 φ p (e -iω j )| -2 a j | h | 2d J -1 (cos hω j + o(1)) as h → ∞, (7)
where a j , j = 1, ..., L, are constants that depend on d j , that is, they do not depend on h.

Remark 1.

If there exists at least one d ι > 1/2, the process in Equation ( 5) is non-stationary and therefore the spectral representation given in Equation ( 6) does not exist. Nevertheless, the model is still adequate to adjust time series with longmemory and seasonality. There is a large amount of papers in the literature related to the estimation of non-stationary ARFIMA models (see, for example, Hurvich and Ray [START_REF] Hurvich | Estimation of the memory parameter for nonstationary or noninvertible fractionally integrated processes[END_REF] and Olbermann et al. [START_REF] Olbermann | Invariance of the first difference in ARFIMA models[END_REF]).

The SARFIMA model is a particular case of the ARUMA process, that is, it is intrinsically related to the model and assumptions described before. As previously discussed, the SARFIMA model has been widely studied, theoretically and empirically, and applied to a variety of real data set (see, for example, Ray [START_REF] Ray | Long-range forecasting of ibm product revenues using a seasonal fractionally differenced arma model[END_REF], Marques [START_REF] Marques | Empirical aspects of the whittle-based maximum likelihood method in jointly estimating seasonal and non-seasonal fractional integration parameters[END_REF], Hassler [START_REF] Hassler | Mis)specification of long memory in seasonal time[END_REF], Reisen et al. [START_REF] Reisen | Estimating seasonal long-memory processes: A monte carlo study[END_REF][START_REF] Reisen | Estimation of seasonal fractionally integrated processes[END_REF], Arteche and Robinson [START_REF] Arteche | Semiparametric inference in seasonal and cyclical long memory processes[END_REF], Palma and Chan [START_REF] Palma | Efficient estimation of seasonal long-range-dependent processes[END_REF], among others). Reisen et al. [START_REF] Reisen | A semiparametric approach to estimate two seasonal fractional parameters in the Sarfima model[END_REF] studied the SARFIMA process which encompasses two seasonal fractional and short memory parameters by deriving the model and the asymptotic properties of the semiparametric ordinary least square estimator (OLS). Under some conditions, they showed that the fractional OLS estimators are asymptotically Normally distributed. In their study, the theoretical properties were investigated for finite sample sizes under different scenarios. In this direction, hereafter, the SARFIMA model considered is a particular case of the one discussed by Reisen et al. [START_REF] Reisen | A semiparametric approach to estimate two seasonal fractional parameters in the Sarfima model[END_REF], that is, the SARFIMA model with two fractional parameters (at zero and seasonal frequencies), but in the context of robustness of the OLS long-memory estimators.

Let now

∇ d = (1 -B) d (1 -B s ) D , (8)
where d = (d, D) is the memory parameter vector; d and D are the fractional parameters at the zero (or long-run) and seasonal frequencies, respectively, satisfying Assumption 5 and s ∈ N * = N -{0} is the seasonal length. D = 0 implies that the process does not have seasonal poles. The process X t ≡ {X t } t∈Z is now defined as a zero-mean SARFIMA process satisfying

X t = ψ(B)∇ -d t = θ q (B) φ p (B) ∇ -d t , (9)
where the process t , the fractional vector d = (d, D) and the coefficients φ p (B) and θ q (B) satisfy the conditions previously stated.

In addition, the fractional filters are

(1 -B k ) x = ∞ j=0 x j -B k j , k = 1, s, and x = d, D, where x j = Γ(x + 1) Γ( j + 1)Γ(x -j + 1)
, and Γ(•) is the well-known Gamma function.

Remark 2. If |d + D| < 1/2, |d| < 1/2 and |D| < 1/2,
X t is a stationary and invertible process. These conditions are derived by Assumption 5. At seasonal frequencies ω s ∈ (0, π], the spectral density becomes unbounded and behaves as

f (ω + ω s ) ∼ C 1 |ω| -2D , ω → 0, (10)
and at the zero frequency,

f (ω) ∼ C 2 |ω| -2(d+D) , ω → 0, (11)
where 0 < C 1 , C 2 < ∞ (see Proposition 1 in Reisen et al. [START_REF] Reisen | A semiparametric approach to estimate two seasonal fractional parameters in the Sarfima model[END_REF]). From the above, it can be seen that the process has spectral density with poles at zero and seasonal frequencies.

If, in addition to long-memory and periodicity features, X t also presents outliers, it is necessary to build a model estimation method that encompasses this type of characteristic. To achieve this, the fractional parameter estimation procedure here suggested concatenates the methods given in Reisen et al. [START_REF] Reisen | A semiparametric approach to estimate two seasonal fractional parameters in the Sarfima model[END_REF], Lévy-Leduc et al. [START_REF] Lévy-Leduc | Robust estimation of the scale and the autocovariance function of Gaussian short and long-range dependent processes[END_REF][START_REF] Lévy-Leduc | Large sample behaviour of some well-known robust estimators under long-range dependence[END_REF][START_REF] Lévy-Leduc | Asymptotic properties of U-processes under long-range dependence[END_REF] and Molinares et al. [START_REF] Molinares | Robust estimation in longmemory processes under additive outliers[END_REF] to estimate the fractional parameters of the model presented in Equation ( 9). The robust estimation method is discussed in the next sub-section.

A robust semiparametric estimator for the vector d in the SARFIMA model

Let {X 1 , . . . , X n } be a sample from the process in Equation ( 9). The periodogram function of X t is given by

I n,X (ω j ) = 1 2πn n t=1 X t e iω j t 2 , (12)
where ω j = 2π j n , j = 1, . . . , ( n 2 -1), are the Fourier frequencies. Since X t is a stationary process, I n,X (ω j ) can be also written as follows:

I n,X (ω j ) = 1 2π n-1 k=-(n-1) γX (k)e -iω j k , (13)
where γX (k) is the sample autocovariance function of {X 1 , . . . , X n }.

An alternative robust spectral estimator for the ARFIMA process was proposed in Molinares et al. [START_REF] Molinares | Robust estimation in longmemory processes under additive outliers[END_REF] which replaces the classical sample autocovariance γX (k) in Equation ( 13) by the robust autocovariance function given in Ma and Genton [START_REF] Ma | Highly robust estimation of the autocovariance function[END_REF]. The main asymptotic results of the robust autocorrelation function (ACF) for long-memory processes are discussed in Lévy-Leduc et al. [START_REF] Lévy-Leduc | Robust estimation of the scale and the autocovariance function of Gaussian short and long-range dependent processes[END_REF]. Here, as previously mentioned, the robust ACF estimator (Ma and Genton [START_REF] Ma | Highly robust estimation of the autocovariance function[END_REF]) will be used to obtain the estimates of the parameters d and D in Model [START_REF] Chan | A note on time series model specification in the presence outliers[END_REF]. Therefore, the estimation approach here proposed is an extension to the methods discussed by Reisen et al. [START_REF] Reisen | A semiparametric approach to estimate two seasonal fractional parameters in the Sarfima model[END_REF].

For a sample ν = (ν 1 , ν 2 , . . . , ν n ) , Rousseeuw and Croux [START_REF] Rousseeuw | Explicit scale estimators with high breakdown point[END_REF][START_REF] Rousseeuw | Alternatives to the median absolute deviation[END_REF] suggested the scale robust estimator Q n (.), which is based on the τth order statistic of n 2 distances {|ν j -ν k |, j < k}, and can be written as

Q n (ν) = c × {|ν j -ν k |; j < k} (τ) , ( 14
)
where c is a constant used to guarantee consistency (c = 2.2191 for the normal distribution) and τ =

( n 2 )+2 4 + 1.
Based on the scale robust estimator Q n (ν), Ma and Genton [START_REF] Ma | Highly robust estimation of the autocovariance function[END_REF] proposed the following robust sample autocovariance function

γ Q n ,y (h) = 1 4 Q 2 n -h,y (u + v) -Q 2 n -h,y (u -v) , (15)
where u and v are vectors containing the initial nh and the final nh observations, respectively, of a time series sample Y 1 , ..., Y n of a process Y t∈Z with absolutely summable autocovariance function, that is, ∞ h=0 |γ y (h)| < ∞ where γ y (h) is the autocovariance of Y t at the lag h .

As pointed out by the authors, the robust estimator of the autocorrelation function could be obtained by dividing γ Q n ,y (h) in Equation ( 15) by the product of Q n ,y (u) and Q n ,y (v). However, this would not be a natural autocorrelation estimator because it would not be bounded between -1 and 1. Thus, a highly robust autocorrelation estimator can be computed by

ρ Q n ,y (h) = Q 2 n -h,y (u + v) -Q 2 n -h,y (u -v) Q 2 n -h,y (u + v) + Q 2 n -h,y (u -v) . ( 16
)
It can be shown that | ρ Q n ,y (h)| ≤ 1 for all h. As previously mentioned, the asymptotic properties of γ Q n ,y (h) when the time series Y t has short and long-memory properties are discussed in Lévy-Leduc et al. [START_REF] Lévy-Leduc | Robust estimation of the scale and the autocovariance function of Gaussian short and long-range dependent processes[END_REF][START_REF] Lévy-Leduc | Large sample behaviour of some well-known robust estimators under long-range dependence[END_REF][START_REF] Lévy-Leduc | Asymptotic properties of U-processes under long-range dependence[END_REF]. Here, some of their results are addressed in the following remarks.

Remark 3. Under the assumption that Y t∈Z follows a Gaussian long-memory process with memory parameter 0 < ϕ < 1/2 (see Assumptions 3 and 4 ), Lévy-Leduc et al. [START_REF] Lévy-Leduc | Robust estimation of the scale and the autocovariance function of Gaussian short and long-range dependent processes[END_REF] showed asymptotic results for Q n ,y (.) and γ Q n ,y (h). In particular, the authors demonstrated that for 1/4 < ϕ < 1/2, the robust autocovariance estimator γ Q n ,y (h) has the same asymptotic behavior as the classical autocovariance estimator γy (h). In this case, there is no loss of efficiency. For 0 < ϕ < 1/4, γ Q n ,y (h) has the same rate of convergence of γy (h), but with different variances. The standard Gaussian ARFIMA(p, ϕ, q) is one particular case of the model discussed by the authors. Remark 4. Note that, although the Gaussian distribution is required to obtain the asymptotic distribution properties of Q n ,y (.) and γ Q n ,y (h), the finite sample size investigation given in Lévy-Leduc et al. [START_REF] Lévy-Leduc | Robust estimation of the scale and the autocovariance function of Gaussian short and long-range dependent processes[END_REF] showed that the estimation method also performs well under non Gaussian observations, that is, the robust autocovariance estimator does not seem to be affected by the skewness of the data. Additionally, Molinares et al. [START_REF] Molinares | Robust estimation in longmemory processes under additive outliers[END_REF], Lévy-Leduc et al. [START_REF] Lévy-Leduc | Robust estimation of the scale and the autocovariance function of Gaussian short and long-range dependent processes[END_REF][START_REF] Lévy-Leduc | Large sample behaviour of some well-known robust estimators under long-range dependence[END_REF] and Sarnaglia et al. [START_REF] Sarnaglia | Robust estimation of periodic autoregressive processes in the presence of additives outliers[END_REF], using finite sample size, discuss the finite sample robustness property of γ Q n ,y (h) under different scenarios of the data contaminated with additive outliers, that is, observations which can produce skewness in the series. Their investigations strongly suggest the use of γ Q n,y (h) in this context. A thorough search of the relevant literature on this topic indicated that, under skewness distribution, the asymptotic property of γ Q n,y (h) is very difficult to be obtained and it is still an open problem.

Remark 5. From Equation [START_REF] Brunelli | Two-days ahead prediction of daily maximum concentrations of SO 2 , O 3 , PM 10 , NO 2 , CO in the urban area of Palermo, Italy[END_REF], it can be seen that the autocovariances of the process X t , given in Equation ( 9), show an asymptotic slow decay typical of Equation (3) but with oscillations that depend on the frequency ω, as shown in Proposition 1 of Reisen et al. [START_REF] Reisen | A semiparametric approach to estimate two seasonal fractional parameters in the Sarfima model[END_REF]. Additionally, from Remark 2, the process has spectral density with poles at zero and seasonal frequencies. Under the assumption that the innovations t are from a Gaussian distribution, the SARFIMA model in Equation ( 9) satisfies the model conditions given in Lévy-Leduc et al. [START_REF] Lévy-Leduc | Robust estimation of the scale and the autocovariance function of Gaussian short and long-range dependent processes[END_REF] (see Assumptions 3 and 4) and, therefore, the asymptotic results for the robust autocovariance discussed in Lévy-Leduc et al. [START_REF] Lévy-Leduc | Robust estimation of the scale and the autocovariance function of Gaussian short and long-range dependent processes[END_REF] for the ARFIMA model are also valid for samples from the SARFIMA process (see, Remark 3). The finite sample properties discussed in Remark 4 are also expected for the model here studied (see the finite sample size investigation analysis in Subsection 3).

Based on the previous discussion and on Molinares et al. [START_REF] Molinares | Robust estimation in longmemory processes under additive outliers[END_REF], for the sample X 1 , ..., X n of Model ( 9), a robust spectral estimator can be computed as follows

I Q n ,X (ω) = 1 2π |h|<n κ(h) γ Q n,X (h) cos(hω), (17)
where κ(h) is defined as

κ(h) =        1, |h| ≤ ξ, 0, |h| > ξ. (18)
κ(h) is a particular case of the lag window functions used in classical spectral theory to obtain a consistent spectral estimator, and ξ is the truncation point which is a function of n,

say ξ = G(n), where G(n) must satisfy G(n) → ∞, n → ∞, with G(n) n → 0. G(n) is usually chosen to be G(n) = n β
, where 0 < β < 1 (see, e.g., Priestley [39,). In addition, the robust ACF estimator given in Equation ( 16) does not have the same finite-sample properties as the classical one. For large h, the number of observations in the calculation of γ Q n,X (h) is very small and, consequently, ρ Q n,X (h) becomes very unstable. Therefore, the bandwidth ξ will avoid these undesirable covariance estimates in the calculation of the estimator given in Equation [START_REF] Giraitis | A generalized fractionally diferencing approach in long-memory modeling[END_REF]. See Assumption 3 in Molinares et al. [START_REF] Molinares | Robust estimation in longmemory processes under additive outliers[END_REF] for more discussion on the choice of G(n). Note that, similar to the classical spectral estimation theories, other lag window functions can be used to obtain a robust spectral estimator, as discussed in Molinares et al. [START_REF] Molinares | Robust estimation in longmemory processes under additive outliers[END_REF].

To estimate d = (d, D) for the SARFIMA process in Equation 9, Reisen et al. [START_REF] Reisen | A semiparametric approach to estimate two seasonal fractional parameters in the Sarfima model[END_REF] suggest the ordinary least squares estimator (OLS) dCL = ( dCL , DCL ) computed from the approximated multiple linear regression equation

log I n,X (ω k j ) a 0 -D log 2 sin sω k j 2 2 -d log 2 sin ω k j 2 2 + u k j , (19)
where a 0 is a constant and

u k j = log I n,X (ω k j ) f X (ω k j ) -E log I n,X (ω k j ) f X (ω k j ) .
The frequencies ω k j , k = 0, 1, ..., [ s 2 ], 1 ≤ j ≤ M, are defined as

ω k j = 2πk s + 2π j n , k = 0, 2πk s ± 2π j n , k = 1, ..., [ s 2 ] -1. ( 20
)
ω [ s 2 ] j =        2π[ s 2 ]
s -2π j n , s even;

2π[ s 2 ]
s ± 2π j n , s odd. In the above equations, M = M(n) is the bandwidth that has to satisfy

(A6) M n α log M + 1 M → 0 , as n → ∞ ,
for some α > 0. One suggestion is to use an α that avoids overlapping frequencies, for example, M < n-1 2s (see, also, Remark 7). Under some assumptions which also includes Gaussian innovations, Reisen et al. [START_REF] Reisen | A semiparametric approach to estimate two seasonal fractional parameters in the Sarfima model[END_REF] establish that

√ M( dCL -d) → N W -1 b, π 2 6 W -1 , (21)
where b and W are a vector and a 2 × 2 matrix of constants, respectively. See, Theorems 1 and 2 in Reisen et al. [START_REF] Reisen | A semiparametric approach to estimate two seasonal fractional parameters in the Sarfima model[END_REF].

In this work, a robust OLS estimator dR = ( dR , DR ) is proposed replacing I n,X (ω k j ) in Equation ( 19) by I Q n,X (ω k j ), given in Equation [START_REF] Giraitis | A generalized fractionally diferencing approach in long-memory modeling[END_REF]. The choice of the bandwidths M = M(n) will be based on Assumption 6 under restrictions of Equation [START_REF] Hosking | Fractional differencing[END_REF]. Remark 6. Note that, although the asymptotic distribution of dR is still an open problem, the consistency and distribution properties of the estimators, finite sam-ple properties of the estimators of the long-memory parameter and the sample autocovariance functions in the ARFIMA model are investigated by Molinares et al. [START_REF] Molinares | Robust estimation in longmemory processes under additive outliers[END_REF] and Lévy-Leduc et al. [START_REF] Lévy-Leduc | Robust estimation of the scale and the autocovariance function of Gaussian short and long-range dependent processes[END_REF], respectively. These works together with the simulation results displayed in the following section for the SARFIMA process support the use of the robust function (Equation 17) as an alternative spectral estimator to obtain robust fractional parameter estimates in a real time series that presents long-memory and additive outliers features. The dR estimator is implemented in R-project [START_REF]R: a language and enviroment for statistical computing[END_REF] and the code can be obtained upon request.

Remark 7. Note that, if AR and/or MA coefficients are introduced in the model these short-memory parameters lead to biased estimates and the bias (positive or negative) will depended on the size of the bandwidth M and on the values of the AR/MA parameters. This issue is well-documented in Reisen et al. [START_REF] Reisen | A semiparametric approach to estimate two seasonal fractional parameters in the Sarfima model[END_REF], Reisen [START_REF] Reisen | Estimation of the fractional difference parameter in the Arima(p, d, q) model using the smoothed periodogram[END_REF] and Hurvich and Ray [START_REF] Hurvich | Estimation of the memory parameter for nonstationary or noninvertible fractionally integrated processes[END_REF]. See, also, Hassler [START_REF] Hassler | Mis)specification of long memory in seasonal time[END_REF] for the flexible SARFIMA model.

A simulation study

In this section, the finite sample performance of the OLS robust fractional parameter estimator ( dR = ( dR , DR ) ) is investigated through Monte Carlo experiments for SARFIMA models, with p = q = 0 and s = 12, and with i.i.d innovations from a N(0, 1) distribution. For comparison purposes, the classical OLS fractional estimator ( dCL = ( dCL , DCL ) ), which is based on the classical periodogram (Reisen et al. [START_REF] Reisen | A semiparametric approach to estimate two seasonal fractional parameters in the Sarfima model[END_REF]) is also considered in the simulation study. Since the simulated models here do not have AR and/or MA components, the OLS estimates were computed based on the bandwidth M = [(n -2s)/2s] in order to avoid overlapping frequencies (see, also, Assumption 6 and Remark 7). For the robust spectral estimator (Equations ( 17) and ( 18)), the bandwidth was G(n) = n 0.7 (see, Molinares et al. [START_REF] Molinares | Robust estimation in longmemory processes under additive outliers[END_REF]). The sample size is n = 1000 and the mean and standard deviation were calculated over 1000 replications. Other scenarios with smaller sample sizes (for example, n = 500, s = 3) were also investigated resulting in similar conclusions (they are available upon request).

Figure 1 shows the box-plots with the results of both estimators for series generated without outliers (no contamination). It can be seen that, in general, both methods perform similarly, that is, under the scenario of a non-contaminated time series both estimation methods lead to comparable results, with estimates close to the real values of d and D. However, the robust estimates have slightly smaller variation compared with the classical ones (see, also, Molinares et al. [START_REF] Molinares | Robust estimation in longmemory processes under additive outliers[END_REF] and Lévy-Leduc et al. [START_REF] Lévy-Leduc | Robust estimation of the scale and the autocovariance function of Gaussian short and long-range dependent processes[END_REF] for ARFIMA models).

Figures 2 and3 present the estimated densities of the classical and robust estimators for SARFIMA models with d = D = 0.1 and d = 0.3 and D = 0.1, respectively. Although, as previously stated, the asymptotic distribution of dR was not yet proved, it can be seen that the standardized estimates are fairly close to the density of the N(0, 1) distribution. Therefore, this simple simulation study can give some support for the theoretical discussion presented in the previous remarks and the use of the proposed estimation method in the context of this paper, that is, in time series with long-memory, seasonality and outliers. The discussion of the performance of the robust autocovariance estimate under non-Gaussian distribution made in Remark 4 is also valid for the parameter estimation. Figure 4 displays the densities of the estimates when the innovations follow a Student's t-distribution with 3 degrees of freedom. The performance of the estimates was similar to the Gaussian case. See, also, Lévy-Leduc et al. [START_REF] Lévy-Leduc | Robust estimation of the scale and the autocovariance function of gaussian short and long-range dependent processes[END_REF] (Section 3.3) for non-Gaussian observations.

To investigate the robustness property of dR = ( dR , DR ) in a time series under outliers, the SARFIMA model contaminated with additive outliers was simulated according to the model structure given in Lévy-Leduc et al. [START_REF] Lévy-Leduc | Robust estimation of the scale and the autocovariance function of gaussian short and long-range dependent processes[END_REF] and Molinares et al. [START_REF] Molinares | Robust estimation in longmemory processes under additive outliers[END_REF] and this is summarized below. Let now Y t be defined by

Y t = X t + W t , (22)
where the parameter represents the magnitude of an additive outlier and W t is a random variable with probability distribution

P (W t = -1) = P (W t = 1) = δ/2 and P (W t = 0) = 1 -δ ,
where E[W t ] = 0 and E[W 2 t ] = Var(W t ) = δ. Note that Equation ( 22) is based on the parametric models proposed by Fox [START_REF] Fox | Outliers in time series[END_REF]. W t is the product of Bernoulli (δ) and Rademacher random variables; the latter equals 1 or -1, both with probability 1/2. X t and W t are independent random variables.

Figure 5 presents the box-plots for the classical and robust estimators considering the series with outliers ( = 15 and δ = 0.05). As can be seen from the box-plots, the classical estimator is clearly affected by additive outliers while the robust one keeps almost the same picture of the non-contaminated scenario. This simple investigation leads to analogous conclusions given in Lévy-Leduc et al. [START_REF] Lévy-Leduc | Robust estimation of the scale and the autocovariance function of gaussian short and long-range dependent processes[END_REF] and Molinares et al. [START_REF] Molinares | Robust estimation in longmemory processes under additive outliers[END_REF], when the generated process follows an ARFIMA model, that is, the classical OLS fractional estimator is completely influenced by the outliers while, in general, the robust one is not. In the non-Gaussian series, the methods displayed a similar performance and the results are available upon request.

An application to SO 2 pollutant

The daily average SO 2 concentration is expressed in µg/m 3 and was measured at the Air Quality Automatic Monitoring Network (AQAMN) of Cariacica, which belongs to the Metropolitan area of the Great Vitória Region (GVR) -ES -Brazil. GVR is comprised of seven cities with a population of about 1.7 million inhabitants in an area of 2,331 km 2 . The region is situated along the South Atlantic coast of Brazil (latitude 20 • 19S, longitude 40 • 20W) and has a warm tropical climate, with average temperatures ranging from 24 • C (Celsius) to 30 • C.

The raw series has a sample size of 1826 daily observations (Figure 6), measured from January 1st 2005 to December 31st 2009. The maximum concentration is generally observed in the winter months (southern hemisphere) from July to September. It can be seen that the series has some high concentration of SO 2 in different points in time. As previously mentioned, these large peaks can be viewed here as outliers, since their values can provoke serious damage to the statistical functions, such as the mean and the standard deviation and, therefore, may affect the correlation structure of the series, leading to misleading results.

Since the variability in the series did not seem to be stable, a log transformation was used (Z t = log(X t )). Moreover, the series was divided into two parts: learning and prediction sets. The 1626 observations from January 1st 2005 to June 14th 2009 were considered as learning set and the remaining 200 observations were kept for the forecasting study. The potential effect of the large concentrations on the sample summary functions is now addressed The sample autocorrelation (ACF), the partial autocorrelation (PACF) and the periodogram functions of Z t are shown in Figure 7. These plots indicate possible seasonal behavior with a period equal to seven, which is an expected result since the data are daily mean levels.

The robust ACF, PACF and periodogram functions are displayed in Figures 8(a), 8(b) and 8(c), respectively, to compare them with the classical ones in order to examine whether there is any effect of the large peaks on these functions. As can be seen from these figures, the high levels of concentration reduce the size of the classical ACF and PACF functions while increase the peaks of the periodogram; that is, the classical periodogram across the frequencies close to zero are much higher than the robust ones. In particular, the sample ACFs are smaller than the robust ones, for instance, for lags h = 1, 3, 5, 10, the robust autocorrelations were equal to 0.72, 0.49, 0.48 and 0.41, while the classical autocorrelations were equal to 0.62, 0.42, 0.41 and 0.35, respectively. Such a behavior was theoretically justified in Corollary 1 of Molinares et al. [START_REF] Molinares | Robust estimation in longmemory processes under additive outliers[END_REF]. Note that the long memory property of the series is well observed by looking at the periodogram plots in Figures 7(c The robust SARFIMA modeling strategy follows similar steps suggested in Hosking [START_REF] Hosking | Fractional differencing[END_REF] and investigated by Reisen [START_REF] Reisen | Estimation of the fractional difference parameter in the Arima(p, d, q) model using the smoothed periodogram[END_REF] and Reisen and Lopes [START_REF] Reisen | Some simulations and applications of forecasting long-memory time-series models[END_REF]. In semiparametric procedures, the estimation of the model parameters is performed in two steps: firstly, the parameter vector d R is estimated based on the procedures presented in Section 2.2. Secondly, the truncated filter (1-B) dR (1-B s ) DR is used to filter the observations. This new series is used to estimate the autoregressive and moving average parameters. The fitted models and their accuracy are discussed in the next subsections.

Adjusted models

The robust and classical estimates of the parameter vector d are displayed in Table 1 for different bandwidths M, which corresponds to 0.4 < α < 0.55, as discussed in Section 2.2. The values in brackets correspond to the standard deviations (s.d.). To compute the robust estimates, the bandwidth in Equation [START_REF] Granger | An introduction to long-memory times series models and fractional differencing[END_REF] was fixed as G(n) = n 0.7 . Note that the estimates in Tables 1 do not satisfy all the stationary conditions, that is, the values do not satisfy |d + D| < 0.5 (see Remarks 1 and 2), but the model still has the mean reverting property in the sense that its cumulative impulse response weights sum to a finite number. In addition, the plot in Figure 6 suggests a phenomenon of mixture of stationary and nonstationary blocks, which may imply that 0.5 < | d + D| < 1.0, although the individual estimates are in the stationary region. However, this possible nonstationarity behavior of the series can be reduced to a stationary time series by differencing the series with a long-memory filter, as is discussed in what follows. As mentioned in Remark 7, in the semiparametric approach, the choice of estimates depends on the size of the bandwidth. For example, large M leads to more bias in the fractional estimators, when there are short-run components in the model. As expected all estimates in Table 1 using the classical periodogram were lower than those adopting the robust periodogram. This is due to the presence of observations with large peaks in the SO 2 concentrations, as previously discussed. The robust estimates at zero frequency are very stable across the values of M. This is an indication that, if there is a short-memory parameter in the estimated model, it is not large enough to make an impact on the estimates of d (see, also, the ACF plots in Figure 9). Therefore, M = 33 was chosen and thus, using the robust periodogram, dR = 0.4510 and DR = 0.2610. Note that the robust estimates have smaller s.d. than the standard ones, which is an expected result (see, for example, Molinares at al. [START_REF] Molinares | Robust estimation in longmemory processes under additive outliers[END_REF] and Lévy-Leduc et al. [START_REF] Lévy-Leduc | Robust estimation of the scale and the autocovariance function of Gaussian short and long-range dependent processes[END_REF] for ARFIMA models). Now, using the fractional differencing parameter estimates, showed Table 1 with M = 33, the series ηt = (1 -B) dR (1 -B s ) DR Z t , t = 1, ..., 1626, was obtained. Figure 9 shows the sample robust ACF and robust PACF functions of ηt . Based on the sample robust ACF and PACF functions displayed in Figure 9, some models were considered to fit ηt and the one which presented the smallest AIC, equal to 1,570.21, was the MA(1) model. Table 2 shows the estimated model. The MA(1) estimate was computed by the Hannan-Rissanen algorithm (see, Brockwell and Davis [START_REF] Brockwell | Time series: theory and methods[END_REF]), in which the classical autocovariance was replaced by the robust one. Although the MA(1) component is adding only a small contribution, the robust SARFIMA(0, d R , 1) × (0, D R , 0) 7 was chosen for the SO 2 average data.

Model adequacy is now addressed (Table 3). The Box-Pierce and Ljung-Box statistics (robust tests) demonstrated that the sample residuals are not timecorrelated. In addition, the results indicated that the residuals are not normally distributed, which was an expected result since the original data is skewed to the right. Other classical residual plots were also analyzed and they led to similar conclusions. These plots are available upon request. An alternative way to demonstrate the usefulness, quality and forecasting performance of the proposed model is to compare it to the standard SARFIMA models. The SARFIMA(0, d CL , 1) × (0, D CL , 0) 7 model (with AIC equal to 1,640.03 and M = 33) was chosen to fit the raw data among other candidates and thus, using the classical periodogram, dCL = 0.4301 and DCL = 0.1925. The residuals of standard SARFIMA model were right-skewed and uncorrelated. The forecast performance of the robust SARFIMA(0, d R , 1) × (0, D R , 0) 7 and the standard SARFIMA(0, d CL , 1) × (0, D CL , 0) 7 models is discussed in the next subsection.

Forecasting investigation

The objective of this section is to verify whether the standard SARFIMA model may lead to less accurate forecasts compared to the robust estimation. As stated before, the observations from June 15th 2009 to December 31st 2009 were discarded from the modeling stage (200 observations) to be used in the out-ofsample forecast study. Forecasts one to ten steps ahead are considered.

To measure the accuracy of the forecasts one to ten steps ahead, the criteria used was the prediction mean square error (PMSE) and the values are displayed in Table 4. From this table, it can be seen that the robust SARFIMA model yields more accurate forecasts than the standard SARFIMA model, especially for long-term forecasts. In addition, the Diebold-Mariano test [START_REF] Diebold | Comparing predictive accuracy[END_REF] showed that the forecast of the two models are significantly different and confirmed the superiority forecasting performance of the robust SARFIMA model. For example, for h = 1 and 10, the statistical Diebold-Mariano test gave p-values of 0.4 × 10 -2 and 0.035, respectively, showing rejection of the null hypothesis in favor of the robust SARFIMA forecasting performance. Figure 10 presents a visual analysis of the one-step-ahead forecast values of the robust SARFIMA model, that is, from June 15th 2009 to December 31st 2009. It indicated a reasonably good performance of the model and estimation method proposed here to fit the data. The robust SARFIMA model is able to capture the dynamics of the SO 2 series for one-step-ahead forecasts. These results corroborate the better performance of the fitted robust SARFIMA model to forecast the SO 2 concentration over the standard SARFIMA model, especially for longer lead times.

Concluding remarks

This article considered a robust SARFIMA model, with two fractional parameters, in the presence of additive outliers. To estimate the fractional parameters d and D, this paper proposed an estimation procedure which is based on the methods suggested in Reisen et al. [START_REF] Reisen | A semiparametric approach to estimate two seasonal fractional parameters in the Sarfima model[END_REF], Lévy-Leduc et al. [START_REF] Lévy-Leduc | Robust estimation of the scale and the autocovariance function of Gaussian short and long-range dependent processes[END_REF] and Molinares et al. [START_REF] Molinares | Robust estimation in longmemory processes under additive outliers[END_REF]. Simulation studies demonstrated that the estimation method works well when the series is contaminated with additive outliers. As an example of application, the daily average SO 2 concentration was analysed to show the usefulness of the proposed methodology. This series presented seasonality, long memory phenomena and some occasional large concentrations. As mentioned in the introduction, the high levels of SO 2 are regarded as potential outliers. Based on the robust estimator of the fractional parameters, the SARFIMA(0, d R , 1) × (0, D R , 0) 7 model was used to fit the real data. The results suggested that the residuals were uncorrelated and not normally distributed. The standard estimation method of the SARFIMA model was also considered in the forecasting performance. The PMSE indicated that the robust SARFIMA model had a better accuracy, especially to forecast for long lead times. The robust method is an attractive procedure for estimating the parameters of the SARFIMA model with long memory, seasonality and additive outliers and it can be easily used in application areas. The results in this paper will hopefully stimulate further research on using robust estimation methods and long-memory models to represent and forecast environmental time series.

  Estimates of D (classical DCL and robust DR ).

Figure 1 :

 1 Figure 1: Box-plots of the estimates dCL , dR , DCL and DR for the SARFIMA model with p = q = 0, d = D = 0.1 and s = 12. Non-contaminated series.

  Empirical densities of N(0, 1), dCL and dR .

  Empirical densities of N(0, 1), DCL and DR .

Figure 2 :

 2 Figure 2: Empirical densities of the N(0, 1) and of the standardized estimates for the SARFIMA with p = q = 0, d = D = 0.1 and s = 12.

  Empirical densities of N(0, 1), dCL and dR .

  Empirical densities of N(0, 1), DCL and DR .

Figure 3 :

 3 Figure 3: Empirical densities of the N(0, 1) and of the standardized estimates for the SARFIMA model with p = q = 0, d = 0.3, D = 0.1 and s = 12.

  Empirical densities of N(0, 1), dCL and dR .

  Empirical densities of N(0, 1), DCL and DR .

Figure 4 :

 4 Figure 4: Empirical densities of the the standardized estimates for the SARFIMA model with innovations generated from a t-student distribution (3 d.f.) and p = q = 0, d = 0.1, D = 0.1 and s = 12.

  Estimates of D.

Figure 5 :

 5 Figure 5: Box-plots of the estimates dCL , DCL , dR and DR for the SARFIMA model with p = q = 0, d = D = 0.1 and s = 12, for series with outliers.

3 )Figure 6 :

 36 Figure 6: Daily SO 2 concentrations (µg/m 3 ) from 2005/01/01 to 2009/12/31.

  ) and 8(c). Both plots indicate high values for the frequencies close to zero.Additionally, in Figures7(d) and 8(d) the log-periodogram is plotted against the log of the frequencies, for classical and robust cases, respectively. The figures

  Log(per.) versus log(freq.).

Figure 7 :

 7 Figure 7: (a) ACF, (b) PACF and (c) periodogram of Z t and (d) log(I n (ω j )) versus log(ω j ).

  Log(rob-per.) versus log(freq.).

Figure 8 :

 8 Figure 8: (a) Robust ACF, (b) robust PACF and (c) robust periodogram of Z t and (d) log(I Q n (ω j )) versus log(ω j ) (d).

Figure 9 :

 9 Figure 9: Robust ACF and robust PACF of ηt .

3 )Figure 10 :

 310 Figure 10: Forecasted values by the robust SARFIMA model and SO 2 concentrations (µg/m 3 ) from June 15th 2009 to December 31st 2009, one-step-ahead.

Table 1 :

 1 Estimates of d R and D R for different bandwidths (M) for Z t , using the robust and classical periodogram

			Robust estimates
	M	dR	sd( dR )	DR	sd( DR )
	25 0.4706 (0.0465) 0.2934 (0.0304)
	28 0.4675 (0.0423) 0.2723 (0.0285)
	33 0.4510 (0.0398) 0.2610 (0.0280)
	38 0.4565 (0.0359) 0.2391 (0.0262)
	43 0.4534 (0.0337) 0.2202 (0.0254)
	49 0.4511 (0.0307) 0.2100 (0.0240)
	57 0.4539 (0.0281) 0.1724 (0.0229)
			Classical estimates
	M	dCL	sd( dCL )	DCL	sd( DCL )
	25 0.4295 (0.0578) 0.2214 (0.0378)
	28 0.4303 (0.0552) 0.2080 (0.0372)
	33 0.4301 (0.0500) 0.1925 (0.0351)
	38 0.4221 (0.0481) 0.1923 (0.0351)
	43 0.4099 (0.0451) 0.1952 (0.0341)
	49 0.3983 (0.0420) 0.1661 (0.0329)
	57 0.3876 (0.0386) 0.1440 (0.0315)

Table 2 :

 2 Adjusted robust SARFIMA model for SO 2 concentrations

	Parameter Estimate	sd
	d R	0.4511 0.0398
	D R	0.2610 0.0228
	θ 1	0.0225 0.0145

Table 3 :

 3 Tests for normality* and non correlation (robust tests)**

	Shapiro-Wilk* Jarque-Bera* Box-Pierce** Ljung-Box**
	<0.0001	<0.0001	0.8403	0.8404

Note: the p-values correspond to the robust Box-Pierce and robust Ljung-Box test statistics with lag 1. Other lags were tested and they presented similar conclusions.

Table 4 :

 4 PMSE of the fitted models values to the SO 2 concentration

	Horizon St. SARFIMA (A) Rob. SARFIMA (B) [(A/B)-1]*100
	1	0.0875	0.0857	2.09%
	2	0.1088	0.1034	5.27%
	3	0.1159	0.1079	7.44%
	4	0.1209	0.1126	7.37%
	5	0.1241	0.1140	8.86%
	6	0.1255	0.1156	8.60%
	7	0.1289	0.1161	11.01%
	8	0.1392	0.1264	10.15%
	9	0.1380	0.1240	11.25%
	10	0.1372	0.1236	10.99%

Note: St. = standard; and, Rob. = robust.

Acknowledgements

The authors gratefully acknowledge partial financial support from Fapes/ES, Fapemig/MG and CNPq/Brazil. Part of this paper was revised when Prof. Valde-rio Reisen was visiting the CentraleSupélec (12/2016 to 01/2017). This author is indebted to CentraleSupélec for its financial support. Bovas Abraham was partially supported by the Natural and Engineering Sciences Council of Canada. The authors are grateful to the Editor and referees for the time and efforts in providing very constructive and helpful comments that have led to clarify and substantially improve the quality of th e paper.