Raphaël Raynaud
email: raphael.raynaud@ens.fr

Emmanuel Dormy
email: dormy@phys.ens.fr

Intermittency in spherical Couette dynamos

We investigate dynamo action in three-dimensional numerical simulations of turbulent spherical Couette flows. Close to the onset of dynamo action, the magnetic field exhibits an intermittent behavior, characterized by a series of short bursts of the magnetic energy separated by low-energy phases. We show that this behavior corresponds to the so-called on-off intermittency. This behavior is here reported for dynamo action with realistic boundary conditions. We investigate the role of magnetic boundary conditions in this phenomenon.

I. INTRODUCTION

First suggested by Joseph Larmor in 1919, dynamo action, i.e. the self-amplification of a magnetic field by the flow of an electrically conducting fluid, is considered to be the main mechanism for the generation of magnetic fields in the universe for a variety of systems, including planets, stars, and galaxies [START_REF]Mathematical aspects of natural dynamos[END_REF]. Dynamo action is an instability by which a conducting fluid transfers part of its kinetic energy to magnetic energy.

In experiments, it is rather difficult to achieve a regime of self-excited dynamo action. The low value of the magnetic Prandtl number of liquid metals requires the injection of a sufficiently high mechanical power, and thus generates turbulent flows, before reaching the dynamo threshold. Dynamo action was first observed experimentally only in 2001, in Karlsruhe [START_REF] Stieglitz | [END_REF] and Riga [3], and then in 2007 with a von Kármán swirling flow of liquid sodium [4].

In parallel with these approaches, numerical simulations have been carried out to model either laboratory experiments or astrophysical systems, for which the spherical geometry is relevant. We investigate spherical Couette flow and focus on the characteristics of the magnetic field close to the dynamo onset. We observe a series of short bursts of the magnetic energy separated by low-energy phases. This intermittent behavior, also known as on-off intermittency or blowout bifurcation, is usually interpreted as the effect of a multiplicative noise acting on a bifurcating system [5,6].

On-off intermittency has so far never been observed in dynamo experiments, except in the case of an externally amplified magnetic field [7]. In contrast, it has been reported in a small number of numerical simulations [8][9][10], all relying on a flow in a periodic geometry produced by a periodic analytic forcing. Here we investigate the influence of a realistic choice of boundary conditions on this phenomenon.

II. GOVERNING EQUATIONS

The spherical Couette flow geometry consists of two concentric spheres in differential rotation: the outer sphere, of radius r o , is rotating around the vertical axis e z with an angular velocity Ω, and the solid inner sphere, of radius r i , is rotating at velocity Ω + ∆Ω around an axis that can make an angle θ with e z . The aspect ratio χ = r i /r o is set to 0.35 to mimic that of Earth's liquid core. The spherical shell in between the two spheres is filled with an incompressible conducting fluid of kinematic viscosity ν, electrical conductivity σ, and density ρ. Its magnetic per-meability µ 0 is that of vacuum. The magnetic diffusivity η is defined as η = 1/(µ 0 σ).

We describe the problem in the reference frame rotating with the outer sphere. This introduces two extra terms in the governing equations: the Coriolis force and the centrifugal acceleration. The latter can be rewritten in the form 1 2 ∇ Ω 2 s 2 , where s denotes the distance to the axis of rotation. This term is a gradient and can be added to the pressure term which acts as a Lagrange multiplier to enforce the solenoidal condition on the velocity field. To establish the set of equations for this system, we rely on the same non-dimensional form as in [11]: the velocity u is scaled by r i ∆Ω, the magnetic field B by ρµ 0 r i r o (Ω + ∆Ω) ∆Ω, and the length scale by r o . The Navier-Stokes equation governing the fluid velocity u then takes the form

∂ u ∂t + (u • ∇) u + 2 ERe (e z × u) = - 1 Re ∇Π+ 1 Re ∇ 2 u + 1 Re 1 E + Re χ (∇ × B) × B , (1)
and the induction equation for the magnetic field B,

∂ B ∂t = ∇ × (u × B) + 1 Rm ∇ 2 B . (2)
Both fields are solenoidal

∇ • u = 0 , ∇ • B = 0 . (3)
The dimensionless parameters are the Ekman number E = ν/ Ωr 2 o , the Reynolds number Re = (r o r i ∆Ω) /ν, the magnetic Prandtl number P m = ν/η, and the magnetic Reynolds number Rm = ReP m . The potential Π includes all gradient terms (the pressure term as well as the centrifugal effect introduced above). The Reynolds number varies with the rotation rate of the inner sphere, while the Ekman number is inversely proportional to the rotation rate of the outer sphere. When the latter is at rest, the Ekman number tends toward infinity and the Coriolis term in the Navier-Stokes equation vanishes. In our simulations, the Ekman number is set to 10 -3 . This moderate value yields a moderate computing time.

We impose no slip boundary conditions for the velocity field on both spheres. Magnetic boundary conditions are of three types. The first one can only be applied to the inner sphere, as it implies a meshing of the bounding solid domain. The inner sphere can be a conductor with the same electric and magnetic properties as the fluid. In that case the magnetic diffusion equation is discretized and solved in the solid conductor (we refer to this set of boundary conditions as "conducting"). The outer sphere as well as the inner sphere can be electrical insulators. In that case the magnetic field is continuous across the boundary and matches a potential field, decaying away from the boundary. The spherical harmonic expansion allows an explicit and local expression for these boundary conditions (we refer to these boundary conditions as "insulating"). In addition, the use of high-magnetic-permeability boundary conditions may enhance dynamo action [12]. Therefore, we also used boundary conditions which enforce the magnetic field to be normal to the boundary. This is equivalent to assuming that the medium on the other side of the boundary has an infinitely larger permeability (we refer to these boundary conditions as "ferromagnetic"). The different configurations investigated in this study are summarized in Table I. We integrated our system with parody [13], a parallel code which has been benchmarked against other international codes. The vector fields are transformed into scalars using the poloidal-toroidal decomposition. This expansion on a solenoidal basis enforces the constraints (3). The equations are then discretized in the radial direction with a finite-difference scheme on a stretched grid. On each concentric sphere, variables are expanded using a spherical harmonic basis (i.e., generalized Legendre polynomials in latitude and a Fourier basis in longitude). The coefficients of the expansion are identified with their degree l and order m. The simulations were performed using from 150 to 216 points in the radial direction, and the spherical harmonic decomposition is truncated at (l max , m max ) = (70, 20). We observe for both spectra a decrease of more than two orders of magnitude over the range of l and m. This provides an empirical validation of convergence. We checked on a few critical cases that the results are not affected when the resolution is increased to l max = 100 .

Let us define the non-dimensional kinetic and magnetic energy densities as

E k = 1 V s Vs u 2 dx , (4)
E b = 1 V s 1 Re 1 E + Re χ Vs B 2 dx , (5)
in which the unit of energy density is ρ (r i ∆Ω) 2 .

In the above expressions, V s refers to the volume of the spherical shell. In addition, we also investigate the symmetry of the flow and the symmetry of the magnetic field with respect to the equatorial plane. To that end, we define the contributions to the energy densities corresponding to the symmetric and antisymmetric components of the velocity (respectively E kS and E kA ) and magnetic field (respectively E bS and E bA ). The symmetric and antisymmetric contributions to the kinetic energy density respectively correspond to the flows

u S (x, y, z) = 1 2 [u(x, y, z) + u(x, y, -z)] , (6)
u A (x, y, z) = 1 2 [u(x, y, z) -u(x, y, -z)] . (7)
In contrast, the symmetries are reversed for the magnetic field. This comes from the fact that the magnetic field is a pseudovector (i.e, the curl of a vector). Then,

B S (x, y, z) = 1 2 [B(x, y, z) -B(x, y, -z)] , (8)
B A (x, y, z) = 1 2 [B(x, y, z) + B(x, y, -z)] . (9
) According to our definition, the dipolar component is symmetric.

III. DIRECT NUMERICAL SIMULATIONS

As shown by [11], contra-rotation is more efficient than co-rotation for dynamo action. In order to introduce more control over the system, we let the angle θ between the axes of rotation of both spheres take any value in [0, π]. Contrary to our expectations, we do not significantly lower the dynamo threshold with the inclination of the rotation axis of the inner sphere. In fact, for θ = π/2, the fluid is mainly in co-rotation with the outer sphere, and dragged only by a thin layer on the inner sphere, which is not sufficient to trigger dynamo action. In our parameter regime, the best configuration seems to remain θ = π, when the two spheres are in contra-rotation. We therefore keep this parameter fixed in the rest of the study.

A. Role of boundary conditions

Let us first investigate the dynamo transition in this setup at fixed magnetic Prandtl number P m = 0.2, using the Reynolds number as the controlling parameter. With a conducting inner sphere and an insulating outer sphere, we find a critical magnetic Reynolds number Rm c ∈ [291.0; 292.0], which is in good agreement with [11]. Above the dynamo onset, the magnetic field displays an intermittent behavior characterized by series of short bursts of the magnetic energy separated by low energy phases (see Fig. 1). When the distance to the threshold increases, bursts become more and more frequent and eventually intermittency disappears.

Changing the boundary conditions generally leads to different thresholds for dynamo action. Using ferromagnetic boundary conditions, we find a critical magnetic Reynolds number Rm I). At lower Reynolds numbers (top), we see in linear scale the characteristic intermittent bursting. Intermittency gradually disappears at higher Reynolds numbers and the field reaches saturation.

involves larger numerical resolutions. In order to maintain the hydrodynamic Reynolds number at values which involve a moderate resolution, we therefore had to increase the magnetic Prandtl number from 0.2 to 0.4 . We then obtain the dynamo onset for Rm c ∈ [530.0; 534.8] .

We emphasize that we observe the same intermittent regime with all the above choices of boundary conditions as long as the magnetic Reynolds number is close enough to the onset of the instability.

For all boundary conditions, we observe that the dominant mode is predominantly of quadrupolar symmetry [the larger poloidal and toroidal modes are the (l = 2, m = 0) and (l = 1, m = 0) modes, respectively]. For these Reynolds numbers, the flow is predominantly equatorially symmetric (E kA ≪ E kS ).

B. Increasing the magnetic Prandtl number

Having assessed that the intermittent behavior of the magnetic field near onset could be observed with three different sets of boundary conditions, we restrict here our attention to simulations with ferromagnetic boundary conditions.

Figure 2 presents the results we obtain at P m = 2 . Close to the threshold, the magnetic field still exhibits intermittency, but the nature of the process has significantly changed. There is now a clear distinction between two different regimes: phases of dynamo activity separated by phases of pure exponential decay. Both seem to alternate randomly. When the dynamo is active, the magnetic field still displays a quadrupolar symmetry. In contrast, we observe the emergence of an axial dipole during decaying phases. The change of the global symmetry of the field coincides with the change of slope in the decaying phases [see Fig. 2 (bottom) and Fig. 3]. This change of slope is associated with a slower decay of the dipolar component over the quadrupolar mode.

IV. DISCUSSION

A. Canonical model for on-off intermittency

The simplest model that exhibits on-off intermittency is [14]

Ẋ = [a + ζ(t)] X -X 3 , ( 10
)
where a is the distance to the threshold, and ζ a Gaussian white noise of zero mean value and amplitude D, defined as

ζ(t)ζ(t ′ ) = Dδ(t -t ′ )
where indicates the average over realizations (ensemble average). In the absence of noise, the system undergoes a supercritical pitchfork bifurcation at a = 0 . If a is sufficiently small, the fluctuations lead to on-off intermittency, with bursts (a+ζ > 0) followed by decays (a+ζ < 0). During the off phases, one can neglect nonlinearities and write Ẏ = a+ζ(t), with Y = ln(X). Thus, ln(X) should follow a random walk, with a small positive bias. Since solutions of Eq. ( 10) mimic solutions of the magnetohydrodynamics equations we observe in Fig. 1, we further investigate some properties of the model. (i) Equation (10) leads to a stationary probability density function (PDF) of the form [START_REF] Stratonovitch | Topics in the Theory of Random Noise[END_REF] P

(X) ∝ X (2a/D)-1 e -X 2 /D , (11)
which diverges at the origin for 0 s = (2a/D) -1 < 1 . (ii) In addition, all the moments of X must follow a linear scaling with a. (iii) Finally, another characteristic of this model is that the distribution of the duration of the off phases T off follows a power law behavior, P (T off ) ∼ T -α off , with α = -3/2. To compare these predictions to our results, we rely as in [10] on the magnetic energy density as a global measure of the magnetic field strength.

B. Predictions and results

Figure 4 shows the PDFs of the magnetic energy for a set of simulations at different Reynolds numbers. At low Rm, the PDF is characterized by a linear scaling on a log-log plot. The cutoff at low energies is not predicted by the theory, which considers the limit E b → 0 . For Rm > 310 , the magnetic energy fluctuates around a mean value and the PDF no longer scales as a power law. We see in Fig. 5 that the coefficient s is proportional to the distance to the threshold. Examples of the fit of the exponent s are presented in Fig. 6. The values of the coefficient are mainly affected by the range over which the data are fitted. Thus we select a range as large as possible. We then randomly sample this range with half-size subintervals. We then compute the mean slope and its standard deviation (represented in Fig. 5 with error bars).

We then investigate the linearity of the moments. Figure 7 shows our results for the first and second moments of the magnetic energy. We see that the mean magnetic energy grows linearly as a function of the magnetic Reynolds number. The second moments seem to follow the same linear trend, but only at the lower values of the magnetic Reynolds number. Deviations at larger values of Rm are expected, as this description is only valid in the limit Rm → Rm c . The duration t I of the time series used to compute these values ranges from 3.2 10 3 to 1.4 10 4 U.T. (values are presented in Table II). These integration times are quite significant for a fully three dimensional set of partial differential equations but are necessarily short compared to the ones usually used with simplified models such as Eq. (10). To quantify the uncertainty associated with the moment values, we sampled the integration time with sub-intervals randomly chosen. We then computed the moments on the full interval (symbols in Fig. 7) and the standard deviation on the sub-intervals (reported as error bars). The subintervals can be set from t I /4 to t I /10 without affecting these estimates. Finally, we also tested the distribution of the duration time of the off phases. A definitive validation would require longer simulations, in order to have a significant number of off phases. For this reason, we can not rely on the simulations immediately above the threshold. Despite these short-comings, an illustrative case is presented in Fig. 8. Numerical values are given in Table III. We define an off phase by a magnetic energy below a threshold value (indicated here by the horizontal black line). The PDF (b) decay is close to a power law with the expected -3/2 value for the exponent (dashed line).

To conclude, we emphasize that the predictions of the model are consistent with the threedimensional simulations, and thus confirm the on-off hypothesis for the observed intermittency at low magnetic Prandtl number.

C. Simulations at higher magnetic Prandtl number

The simulations we performed at P m = 2 exhibit a peculiar behavior of the magnetic field. This can be better understood by ex- amining the dynamics of the flow. Indeed, we also carried out purely hydrodynamic simulations at Re = 1200 and observed intermittent transitions between two states. This kind of intermittent behavior of the flow was not reported in [11], but has been observed experimentally [START_REF] Zimmerman | [END_REF]. One state is characterized by larger fluctuations of the energy as we can see in Fig. 9. In addition, the analysis of the energy spectra reveals that the m = 3 modes dominate over the m = 2 modes during the "laminar" phases, whereas both are of the same order during the "turbulent" phases. Duration of the "turbulent" phases tends to increase gradually with the increase of the Reynolds number, so that the intermittent behavior of the flow eventually disappears and is thus no longer present in the simulations at higher Reynolds number in which we have identified on-off intermittency. Dynamo action is inhibited during the "laminar" phases (when the m = 3 modes dominate), which highlights the mechanism which leads to the peculiar behavior of the magnetic field, as we can see in Fig. 10. In contrast, "turbulent" phases favor dynamo action, and one must wait a change in the flow to see the restart of dynamo action after a phase of decay. Moreover, in a full magnetohydrodynamics simulation, we can artificially suppress the m = 3 modes of the velocity field by setting them equal to zero at each time step. We check that it is sufficient to suppress intermittency of the flow. Then, we observe that the phases of exponential decay are also suppressed and the dynamo is no longer intermittent.

V. CONCLUSION

Despite the fact that on-off intermittency has so far never been observed in dynamo experiments, we showed that the phenomenon can appear in numerical simulations of dynamo action using realistic boundary conditions. We identified in several cases the predicted behavior of the PDF of the magnetic energy, linear scaling of the moments, and distribution of the duration of the off phases. In addition, we tested these properties for three different boundary conditions (conducting inner core with insulating outer sphere, insulating or ferromagnetic spheres). Finally, we pointed out a different kind of intermittency due to hydrodynamic transitions that appears at lower Reynolds numbers.

To explain the absence of on-off intermittency in the experiments, several reasons have already been invoked [17]. One explanation could be the imperfectness of the bifurcation (due for instance to Earth's ambient magnetic field). Since it has been shown that low-frequency noise controls on-off intermittency [14], another possible explanation could be that the low-frequency fluctuations are too small. However, the lack of experimental observations of on-off intermittency remains an open question and needs further investigations.

FIG. 1 .

 1 FIG. 1. Time evolution of the magnetic energy in linear (left) and log scale (right) for increasing Reynolds numbers at Pm = 0.2, using a conducting inner sphere (B.C.1 in TableI). At lower Reynolds numbers (top), we see in linear scale the characteristic intermittent bursting. Intermittency gradually disappears at higher Reynolds numbers and the field reaches saturation.

FIG. 2 .

 2 FIG.2. Time evolution of the magnetic energy in linear and log scales for Pm = 2 and Re = 1200, using ferromagnetic boundary conditions. Instead of bursts, we now observe phases of dynamo activity which seem to randomly alternate with phases of exponential decay. The latter are no longer chaotic and are instead characterized by two different decay rates.

FIG. 3 .

 3 FIG.3. (Color online) Time evolution of the symmetric (dashed line) and antisymmetric (solid line) part of the magnetic energy. We focus on one of the decaying phases presented in Fig.2. The decaying phase is characterized by a change of the dominant symmetry, as we can see on the visualizations of the magnetic field lines. The color insets respectively correspond to snapshots in the quadrupolar phase (left) and dipolar phase (right).

FIG. 4 .

 4 FIG. 4. Estimates of the probability density functions. Statistics are done from time series of the magnetic energy obtained with ferromagnetic boundary conditions, for Pm = 0.2.

FIG. 5 .

 5 FIG. 5. Fit of the coefficient s = (2a/D) -1, taking into account the linear domain of the PDFs in the intermittent regime only. Statistics are done from time series of the magnetic energy obtained with ferromagnetic boundary conditions, for Pm = 0.2.

FIG. 6 .

 6 FIG. 6. (Color online) Examples of the fit (dashed lines) of the probability density functions (solid and dotted lines). Statistics are done from time series of the magnetic energy obtained with ferromagnetic boundary conditions, for Pm = 0.2.

4 FIG. 7 .

 47 FIG. 7. (Color online)Moments of the magnetic energy as function of the magnetic Reynolds number. Black squares represent the mean. The dashed line fits these data points with an error about 10% on the slope coefficient. The second moments (red circles) have been rescaled to match with the mean at Rm = 300 .

FIG. 8 .

 8 FIG. 8.Distribution of the duration of the off phases for a simulation performed at Pm = 0.2 and Re = 1470, with conducting boundary conditions. We define an off phase by a magnetic energy below a threshold value (indicated here by the horizontal black line). The PDF (b) decay is close to a power law with the expected -3/2 value for the exponent (dashed line).

FIG. 10 .

 10 FIG.10. (Color online) Time evolution at Re = 1200 and Pm = 2 of the total magnetic energy E b (dashed green line) and the kinetic energies E2 (dotted black line) and E3 (solid red line) for the m = 2 and m = 3 modes, respectively. When the latter becomes larger than the former, the dynamo is no longer sustained and the magnetic energy exponentially decays.

TABLE .

 . III. Estimate of the exponent α of the PDF of the duration of the off phases for different threshold values and different ranges over which the fit is done. The standard error on the estimate of α is about 1%. Range values correspond to log 10 T off [x axis in Fig. 8(b)].

ACKNOWLEDGMENTS

The numerical simulations have been carried out at CEMAG, CINES, and MESOPSL. We thank S. Fauve, F. Pétrélis, and M. Schrinner for fruitful discussions and comments. We are most grateful to C. Gissinger for technical assistance.


Mor Ndongo

Abdou Kâ Diongue

Simplice Dossou-Gbété

Seasonal Fractional ARIMA model with BL-GARCH type innovations

Keywords: Long memory, Seasonality, Volatility clustering, Leverage effects, CSS method, Monte Carlo experiments

In this paper, we introduce the class of seasonal ARFIMA models with bilinear GARCH (BL-GARCH) type innovations that are capable of capturing simultaneously four key properties of non-linear time series: long range dependence, seasonality, volatility clustering and leverage effects. Stationarity and invertibility conditions are derived and conditional sum of squares (CSS) estimation of the model is also considered. Under some assumptions, we show that the resulting estimators are consistent and asymptotically normally distributed. Monte carlo simulation results are presented to evaluate the small-sample performance of the CSS method for various models.

Introduction

When dealing with empirical time series arising from diverse fields of applications, we are confronted with the phenomenon of long memory or long range dependence. A popular way to analyze a long memory time series is to use fractionally integrated autoregressive moving average (FARIMA) processes introduced by Granger and Joyeux [START_REF] Granger | An intoduction to long memory time series and fractional differencing[END_REF] and Hosking [START_REF] Hosking | Fractional Differencing[END_REF]. Furthermore, most of time series in real life may have a persistent periodic behavior, in addition to long term structure. Unfortunately, the FARIMA model does not allow to take into account a periodic or cyclical behavior. Thus, models dedicated to take seasonal or cyclical components with long memory have been developed. Recent contributions related to the seasonal FARIMA model are Porter-Hudak [START_REF] Porter-Hudak | An application of the seasonal fractionally differenced model to the monetary aggegrates[END_REF], Hassler [START_REF] Hassler | Mis)specification of long memory in seasonal time series[END_REF], Arteche and Robinson [START_REF] Arteche | Semiparametric inference in seasonal and cyclical long memory processes[END_REF] and Reisen et al. [START_REF] Reisen | Estimation of seasonal fractionally integrated processes[END_REF]. All of the aforementioned works assume that the conditional variance of time series is a constant over time. However, non-constant variance in non-linear time series is a challenging modelling exercise, considered among many other things by Tong [START_REF] Tong | Non-linear Time Series: A Dynamical System Approach[END_REF]. In particular, the stylized fact that the volatility of financial time series is non-constant has been long recognized in the literature (see e.g. Bollerslev [START_REF] Bollerslev | Generalized Autoregressive Conditional Heteroskedasticity[END_REF], Bollerslev, Engle and Woodridge [START_REF] Bollerslev | A Capital As-set Pricing Model With Time-Varying Covariance[END_REF], and Weiss [START_REF] Weiss | ARMA Models With ARCH Errors[END_REF]). Thus, the methodology for modeling time series with long memory behavior has been extended to long memory time series with a time-varying conditional variance. See for instance, Ling and Li [START_REF] Ling | On Fractionally Integrated Autoregressive Moving-Average Time series Models with Conditional Heteroscedasticity[END_REF] who developed the ARFIMA model with GARCH type innovations, and recently Reisen et al [START_REF] Reisen | Modeling and forecasting daily average PM 10 concentrations by a seasonal ARFIMA model with volatility[END_REF] examine the daily average PM 10 concentrations using a seasonal ARFIMA model with GARCH errors. Given that, these models cannot allow to capture asymmetries or non-linearities (Black [3]), we introduce in this paper a new class of seasonal ARFIMA with BL-GARCH type innovations. This approach allows to model simultaneously: long memory, seasonality, volatility clustering and leverage effects, often observed in financial or economics time series. The stationarity and invertibility of the proposed model are analyzed and, the conditional sum of squares (CSS) estimation of the model parameters is considered, and we study its asymptotic properties. Also, the small-sample performance of the CSS approach is illustrated, using Monte Carlo simulations.

The article is organized as follows. In Section 2, the class of SARFIMA-BL-GARCH model is introduced with some important properties concerning the positivity of the conditional variance, as well as, for stationarity and invertibility. In Section 3, the conditional sum of squares estimation is presented, and we derive its asymptotic properties. Finally, Section 4 calibrates the performance of the estimation procedure through Monte Carlo simulations, while Section ?? provides concluding remarks.

Model and probabilistic properties

In this section, we introduce the model, we will work with, and we give conditions of existence and invertibility.

Model

Let B be the back shift operator satisfying BY t = Y t-1 for any process (Y t ) t∈Z , and s ∈ N * the seasonal period, then the polynomial of non-seasonal orders p and q, seasonal orders P and Q are respectively defined by:

φ(B) = 1 -φ 1 B -φ 2 B 2 -• • • -φ p B p θ(B) = 1 + θ 1 B + θ 2 B 2 + • • • + θ q B q Φ(B s ) = 1 -Φ s B s -Φ 2s B 2s -• • • -Φ Ps B Ps Θ(B s ) = 1 + Θ s B s + Θ 2s B 2s + • • • + Θ Qs B Qs .
It is assumed that these polynomials have no common zeros and satisfy the conditions Φ(z s )φ(z) 0 and Θ(z s )θ(z) 0 for |z| = 1. Furthermore, in the above equations, (Φ i ) 1≤i≤P , (φ j ) 1≤ j≤p , (Θ k ) 1≤k≤Q and (θ l ) 1≤l≤q are unknown parameters.

A zero-mean process (X t ) t∈Z is said a seasonal ARFIMA model with bilinear GARCH type errors, denoted hereafter as SARFIMA(p, d, q) × (P, D, Q) s -BL-GARCH(r, m), if the following equation is satisfied

φ(B)Φ(B s )(1 -B) d (1 -B s ) D X t = θ(B)Θ(B s )ε t , (2.1)
where the long-memory parameters d and D are fractional parameters at the zero (or long-run) and seasonal frequencies, respectively . The process (ε t ) t∈Z in equation (2.1) is a BL-GARCH(r, m) model defined by:

ε t = h t Z t (2.2) h 2 t = a 0 + r i=1 a i ε 2 t-i + m j=1 b j h 2 t-j + r * k=1 c k ε t-k h t-k , (2.3)
where a 0 , a i , b j and c k are constant parameters, r, m and r * are non-negative integers with r * = min(r, m), h 2 t is the conditional variance of the process (ε t ) t∈Z given the σ-fields I t-1 generated by the past information {ε t-1 , ε t-2 , . . .}. In above, the process (Z t ) t∈Z is a sequence of independent identically normally distributed random variables with mean 0 and variance 1. This model has the advantage of being characterized by a more flexible parametric structure. In this model, leverage effects are explained by the interactions between past observations and volatilities. More precisely, for c k < 0 a positive quantity is added to

h 2 t if ε t-k < 0 while the same quantity is subtracted if ε t-k > 0. If c k = 0 for all k, the model in (2.1)-(2.
3) is the seasonal fractional ARIMA model with generalized autoregressive conditional heteroscedasticity disturbances (SARFIMA(p, d, q) × (P, D, Q) s -GARCH(r, m)) introduced by Reisen et al [START_REF] Reisen | Modeling and forecasting daily average PM 10 concentrations by a seasonal ARFIMA model with volatility[END_REF]. In addition, if D = P = Q = 0, we obtain the fractionally integrated process with GARCH erros (ARFIMA(p, d, q)-GARCH(r, m)) proposed by Baillie et al [START_REF] Baillie | Analyzing industrialized countries inflation by the fractionally integrated ARFIMA-GARCH model[END_REF] and Ling and Li [START_REF] Ling | On Fractionally Integrated Autoregressive Moving-Average Time series Models with Conditional Heteroscedasticity[END_REF].

According to Giraitis and Leipus [START_REF] Giraitis | A generalized fractionally differencing approach in long memory modelling[END_REF] or Reisen et al. [START_REF] Reisen | Estimation of seasonal fractionally integrated processes[END_REF], one can easily show that

(I -B) d (I -B s ) D = [ s 2 ] j=0 (1 -e iλ j B)(1 -e -iλ j B) d j = [ s 2 ] j=0 (1 -2 cos(λ j )B + B 2 ) d j , (2.4)
with

d 0 = d+D 2 , d j = D, for j = 1, . . . , [ s 2 ] -1, d [ s 2 ] = D 2
, and λ j = 2π j s , for j = 0, . . . , [ s 2 ], where [x] is the greatest integer small than or equal to x, and i is the complex number such that i 2 = -1. Otherwise, by means of the expansion we have:

[ s 2 ] j=0 (1 -2cosλ j B + B 2 ) d j = +∞ j=0 ψ j (d, ν)B j ,
where the coefficients ψ j (d, ν) j≥0 are given by:

ψ j (d, ν) = 0≤l 0 ,••• ,l [ s 2 ] ≤ j, l 0 +•••+l [ s 2 ] = j C l 0 (d 0 , ν 0 ) • • • C l [ s 2 ] d [ s 2 ] , ν [ s 2 ] , (2.5)
and where d = (d 0 , . . . ,

d [ s 2 ] ), ν = (ν 0 , . . . , ν [ s 2 ]
) with ν j = cos(λ j ), for j = 0, . . . , [ s 2 ]. The weights (C l (d i , ν i )) l∈Z are the Gegenbauer polynomials and they can be computed using the following recursion formula:

           C 0 (d i , ν i ) = 1 C 1 (d i , ν i ) = 2d i ν i C j (d i , ν i ) = 2ν i d i -1 j + 1 C j-1 (d i , ν i ) -2 d i -1 j + 1 C j-2 (d i , ν i ) , ∀ j > 1.

Probabilistic properties

We now specify some probabilistic properties of the model 2.1 -2.3. First, we give the conditions for which the conditional variance h 2 t , defined in (2.3), is non-negative. It is important in practice for estimation theory (using quasi-maximum likelihood methods) that a model as in (2.2) and (2.3) does not generate negative conditional variance h 2 t in sample, since the log quasi-likelihood involves a term in log(h 2 t ), which explodes to -∞ as h 2 t approaches 0, and is ill-defined for h 2 t ≤ 0. The second set of conditions concerns the second-order stationary of the BL-GARCH process and the stationary for a seasonal fractional ARIMA model. The BL-GARCH model with gaussian distributions was recently proposed by Storti and Vitale [START_REF] Storti | BL-GARCH models and asymmetries in volatility[END_REF]. This study was extended by Diongue et al [START_REF] Diongue | BL-GARCH models with elliptical distributed innovations[END_REF], considering elliptical distributions, and they give new probabilistic results concerning the stationarity of the process and the moments.

For the positivity of the conditional variance, Storti and Vitale [START_REF] Storti | BL-GARCH models and asymmetries in volatility[END_REF] show that, for a 0 > 0 and s i > 0 (i = 1, . . . , max(r, m)r * ), with s i = a r * +i if r > m, or s i = a r * +i if r < m, a sufficient conditions for h 2 t > 0 is given by:

c 2 i < 4a i b i , for i = 1, . . . , r * . (2.6)
They give also the condition of second-order stationary for the BL-GARCH model, restricting to the case r = m. Here, we extend this result to the more general case in which r m. The result is stated in the following corollary.

Corollary 2.1. The BL-GARCH(r, m) process, defined by equations (2.2) and (2.3), is stationary in widesense if and only if all the roots B i of the polynomial

Π (B) = 1 - r * i=1 π i B i - max(r,m) i=r * +1 α i B i , (2.7)
with B being the Backward operator, π i = a i + b i , for i = 1, . . . , r * , lie outside the unit circle. Note that in (2.7), we set

α i = a i if r > m, α i = b i if r < m.
In other words, the condition can be expressed as follows:

Π (B) = 0 for |B i | > 1, i = 1, . . . , max(r, m). (2.8)
Proof. The proof of Corollary 2.1 follows mainly the lines of the proof of Theorem 2 in Storti and Vitale [START_REF] Storti | BL-GARCH models and asymmetries in volatility[END_REF].

Concerning the seasonal fractional ARIMA model, Reisen et al. [START_REF] Reisen | Estimation of seasonal fractionally integrated processes[END_REF] (see also Giraitis and Leipus [START_REF] Giraitis | A generalized fractionally differencing approach in long memory modelling[END_REF]) show that it is stationary and invertible if and only if the following conditions hold:

|D + d| < 1 2 and |D| < 1 2 .
(2.9)

Therefore, using the model properties in Reisen et al. [START_REF] Reisen | Estimation of seasonal fractionally integrated processes[END_REF] and the conditions listed in Storti and Vitale [START_REF] Storti | BL-GARCH models and asymmetries in volatility[END_REF], the following theorem is established for the SARFIMA(p, d, q) × (P, D, Q) s -BL-GARCH(r, m) model.

Theorem 2.2. Let (X t ) t∈Z be generated by equations (2.1), (2.2) and (2.3). We assume that the conditions (2.6) and (2.8) hold. We suppose that the polynomials Φ(z s )φ(z) and Θ(z s )θ(z) have no common zeros, and the long memory parameters d and D satisfy condition (2.9). Then, the following statements hold: a) if Φ(z s )φ(z) 0, for |z| = 1, then (X t ) t∈Z is second-order stationary and has the following representation

X t = ∞ j=0 ψ j (d, ν) Θ(z s )θ(z) Φ(z s )φ(z) ε t-j . (2.10)
Hence, (X t ) is strictly stationary and ergodic.

b) If Θ(z s )θ(z) 0, for |z| ≤ 1, then (X t ) is invertible, and that is, (ε t ) can be written as

ε t = ∞ j=0 π j (d, ν) Φ(z s )φ(z) Θ(z s )θ(z) X t-j , (2.11)
where the weights π j (d, ν)

j≥0 are such that π j (d, ν) = ψ j (-d, ν), with the coefficients ψ j (d, ν) j≥0
given by equation (2.5).

Proof. The proof of this theorem is given in the appendix.

Under stationarity and invertibility conditions, MA(∞) and AR(∞) representations are respectively:

X t = ∞ j=0 c j ε t-j and ε t = ∞ j=0 c j X t-j , (2.12)
where the coefficients (c j ) j≥0 and (c j ) j≥0 are determined respectively by:

Θ(z s )θ(z) ∞ j=0 ψ j (d, ν)z j = Φ(z s )φ(z) ∞ j=0 c j z j , |z| ≤ 1 (2.13)
and

Φ(z s )φ(z) ∞ j=0 π j (d, ν)z j = Θ(z s )θ(z) ∞ j=0 c j z j , |z| ≤ 1. (2.14)
In the particular case where P = Q = 0, it is easy to verify that the coefficients (c j ) j≥0 and (c j ) j≥0 can be computed using the following recursion formula:

c 0 = 1 and c j = ψ j (d, ν) + min( j,p) i=1 φ i c j-i - min( j,q) i=1 θ i ψ j-i (d, ν) , ∀ j ≥ 1 (2.15) and c0 = 1 and c j = π j (d, ν) - min( j,p) i=1 φ i π j-i (d, ν) + min( j,q) i=1 θ i c j-i , ∀ j ≥ 1.
(2.16)

Conditional Sum of Squares estimation

Chung [START_REF] Chung | Estimating a generalized long memory process[END_REF] proposed a method based on maximization of the CSS function. In this section, we define the CSS method in the estimation of the SARFIMA(p, d, q) × (P, D, Q) s -BL-GARCH(r, m) process defined by equations (2.1), (2.2) and (2.3), and we establish the asymptotic properties of the CSS estimator.

Definition of the estimator

Suppose that X 1 , . . . , X n are generated by the model (2.1), (2.2) and (2.3). Denote by γ = (d, D, φ 1 , . . . , φ p , θ 1 , . . . , θ q , Φ 1 , . . . ,

Φ P , Θ 1 , . . . , Θ Q ) T , δ = (a 0 , a 1 , . . . , a r , b 1 , . . . , b m , c 1 , . . . , c r * ) T ,
and ω = (γ T , δ T ) T the parameter vector to be estimated. We assume that the parameter ω satisfies the stationary conditions given in Theorem 2.2, and ω 0 = (γ T 0 , δ T 0 ) T is the true value of ω and is in the interior of the compact set Λ ⊆ R p+q+P+Q+r+m+r * +3 . Under the assumption of normality of the standardized innovation (Z t ) t∈Z , the likelihood function is equal to

L(ω) = L(ω | Z 1 , . . . , Z n ) = n t=1 g(Z t , ω),
where g(Z t , ω) is the conditional Gaussian distribution function. Thus, the conditional sum of squares estimator ωn of ω in Λ maximizes the conditional log-likelihood on I 0 , (ignoring the constant)

L(ω) = 1 n n t=1 l t , l t = - 1 2 log(h 2 t ) - ε 2 t 2h 2 t .
(3.17)

Asymptotic properties

The following theorem provides some results related to the asymptotic properties of the SARFIMA-BL-GARCH process estimators, obtained by the CSS approach.

Theorem 3.1. Suppose that (X t ) t∈Z are generated by equations (2.1), (2.2) and (2.3). Assume that the hypothesis of Theorem 2.2 holds. Then, the conditional sum of squares estimator ωn of ω, obtained by maximizing the conditional log-likelihood function (3.17

Ω 0 = diag(Ω γ 0 , Ω δ 0 )
, Ω γ 0 and Ω δ 0 are values of Ω γ and Ω δ at ω = ω 0 , with

Ω γ = E 1 h 2 t ∂ε t ∂γ ∂ε t ∂γ T + 1 2h 4 t ∂h 2 t ∂γ ∂h 2 t ∂γ T and Ω δ = E 1 2h 4 t ∂h 2 t ∂δ ∂h 2 t ∂δ T .
(c) Further, the information matrix Ω γ and Ω δ can be estimated consistently by:

Ωγ = 1 n n t=1

Proof. The proof is given in Appendix

It is important to not that ε t , h t all depend on the theoretically infinite past history of (X t ) or (ε t ). We choose the presample estimates of h t and ε 2 t to be n t=1 ε 2 t /n. As noted in Ling and Li [START_REF] Ling | On Fractionally Integrated Autoregressive Moving-Average Time series Models with Conditional Heteroscedasticity[END_REF], this will not affect asymptotic efficiency and other asymptotic properties (see also Bollerslev [4] and Weiss [START_REF] Weiss | ARMA Models With ARCH Errors[END_REF]).

Monte Carlo Simulation

In this section, we study the finite sample performance of the CSS method described previously, to estimate the parameters of SARFIMA(p, d, q) × (P, D, Q) s -BL-GARCH(r, m) process. The simulation results give the average values, the root mean square error (RMSE) and the mean absolute error (MAE) of the estimation procedure based on 1000 replications. All calculations were carried out using an R programming environment (see [START_REF] Development | R: A language and environment for statistical computing[END_REF]) on a Pentium (R) Dual-Core CPU, 2.60 GHz (2 CPUs) computer.

Because there is no known technique for generating an exact Seasonal ARFIMA model with BL-GARCH type innovations, we adapt the method developed by Stoev and Taqqu [START_REF] Stoev | Simulation methods for linear fractional stable motion and FARIMA using the Fast Fourier Transform[END_REF] for ARFIMA time series with stable innovations. Thus, we approximate the path X t , t = 1, . . . , n by the truncation moving average

X t = M j=0 c j ε t-j , (4.18)
where M is the truncation parameter and fixed in this study to 10 000, and the non random constants (c j ) j≥0 are determined by equation (2.13). The sum in (4.18) is computed using the Fast Fourier Transform (FFT) algorithm, as suggested by Stoev and Taqqu [START_REF] Stoev | Simulation methods for linear fractional stable motion and FARIMA using the Fast Fourier Transform[END_REF]. The BL-GARCH innovations (ε t ) are generated using the method developed in Diongue et al. [START_REF] Diongue | BL-GARCH models with elliptical distributed innovations[END_REF]. The simulation algorithm generates n + 500 observations for each series, saving only the last n. This operation is performed in order to avoid dependence on initial values.

Thus, for fixed value of the seasonal period (s = 4), we carry out an experiment of 1000 samples for the four processes summarized in Table 1. For each model, we use four different sample sizes (n = 100, n = 300, n = 500 and n = 1000). The simulation results are presented in Tables 2345. In these tables, the sample sizes are given in the first columns. The estimations of the parameters are given in the next columns, the root mean square error (RMSE) is given in the row below (in parentheses) and the mean absolute error (MAE) is given under the row of the RMSE (in brackets).

• In the first experiment, a SARFIMA(0, d, 0)×(0, D, 0)-BL-GARCH(1, 1) model is considered (i.e. the model without short memory parameters). The simulation results are summarized in Table 2. Results reveal that parameter estimates are satisfactory, even for small sample sizes (n = 100, n = 300), in the sense that the RMSE and also the MAE are small. We can also remark that the impact of the sample size n on the estimation method. Indeed, when the sample size increases (n = 1000), the results improve significantly too. Figure 1 displays the corresponding boxplots2 and shows the relative scatter of the 1000 estimates. The vertical axis in the figure indicates the deviation from the nominal value of the parameters. This figure confirms the previous results since we observe that confidence intervals are small, and shows also that the impact of the sample sizes on this dispersion.

• In the second experiment, we consider the preceding model (Model 1) with short memory part. Consequently, this experiment is designed to examine the relative performances of the estimators when there are long-memory and short memory components simultaneously. The results are presented in Tables 345. We observe that, when short memory components are introduced in the model, the estimation of long memory parameters is disturbed. Indeed, the RMSE and MAE obtained in Tables 3, 4 and 5 are larger than those presented in Table 2. This phenomenon is already observed in the literature (e.g. Boutahar et al. [START_REF] Boutahar | Estimation Methods of the Long Memory Parameter: Monte Carlo and Application[END_REF], Diongue and Guégan [START_REF] Diongue | Estimation of k-Factor GIGARCH Process: A Monte Carlo Study[END_REF]). However, we observe a significant improvement when the sample size becomes large (n = 500, n = 1000). We remark also that the estimators of the parameters for the BL-GARCH errors seem not to be affected by the presence of the AR and MA components, particularly for larger sample sizes (n = 500, n = 1000). 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 d n=100 n=300 n=500 n=1000

q q q q q q q q q q q q 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 D n=100 n=300 n=500 n=1000 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 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.2 0.4 0.6 a0 n=100 n=300 n=500 n=1000 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 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.2 0.4 0.6 a1 n=100 n=300 n=500 n=1000 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 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.2 0.4 0.6 0.8 1.0 b1 n=100 n=300 n=500 n=1000 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 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.2 0.4 0.6 c1 n=100 n=300 n=500 n=1000

Empirical study

The daily returns of the kenyan stock market index are used for the empirical study in this article to illustrate both the usefulness and the applicability of the proposed model and estimation method.

The raw series has a sample size of 6181 observations, measured form March 1st of 1991 to January 31st of 2008. Table ?? gives the summary statistics of the returns of kenyan stock index for the full sample.

For modeling purpose, the time series is divided into two parts; learning and prediction sets. The 5816 observations from March 1st of 1991 until January 31st of 2007 are considered as learning set and the remaining 365 observations are considered for the prediction study.

Conclusion

In this article, we have developed the family of seasonal ARFIMA process with BL-GARCH type errors, which should prove useful in many fields of time series analysis. It is much more flexible in the simultaneous modeling of long-memory behaviour, seasonal components, volatility clustering and leverage effects; often encountered in financial data. Under some assumptions, the model is shown to be stationary and invertible. The parameter estimation problem is addressed to the conditional sum of squares procedure proposed by Chung [START_REF] Chung | Estimating a generalized long memory process[END_REF], and under some regularly conditions we establish the asymptotic properties of the resulting estimates. Finite sample behaviours of this method were studied using Monte Carlo simulations. It indicate that the approach can yield asymptotic efficient estimates. Consequently, it is reasonable method to deal with seasonal long-range dependent data containing volatility clustering and leverage effects.

Since the results from the model presentation and the estimation methodology are encouraging, it will be interesting to examine the empirical application of the SARFIMA-BL-GARCH model in financial data. We raise also the question of the "best" estimation procedure. This issue appears to be further investigated.

X t = Θ(B s )θ(B) Φ(B s )φ(B) (1 -B) -d (1 -B s ) -D ε t = Ψ 1 (B)ε t ,
where

Ψ 1 (z) = Θ(z s )θ(z) Φ(z s )φ(z) (1 -z) -d (1 -z s ) -D . Now the power series expansion (1 -z) -d (1 -z s ) -D converges
for all |z| ≤ 1, when |D + d| < 1/2 and |D| < 1/2 (see Reisen et al. [START_REF] Reisen | Estimation of seasonal fractionally integrated processes[END_REF] for more details). The series Θ(z s )θ(z) converge also for all z, since it is a polynomial, and that 1/Φ(z s )φ(z) converges for all |z| ≤ 1 when all the roots of the equation Φ(z s )φ(z) = 0 lie outside the circle |z| = 1. Thus, the power series expansion of Ψ 1 (z) converges for all |z| ≤ 1, and so (X t ) t exists with representation (2.10). From Corollary 2.1, the process (ε t ) t is second-order stationary and hence E(ε 2 t ) < ∞. Similar to the proof of Theorem 2.2 in Ling and Li [START_REF] Ling | On Fractionally Integrated Autoregressive Moving-Average Time series Models with Conditional Heteroscedasticity[END_REF], we can show that (X t ) t is also second-order stationary. Indeed, we have

E X 2 t = E                ∞ k=0 c k ε t-k        2         (6.19) = ∞ k=0 c 2 k E ε 2 t-k + k 1 ,k 2 =0, k 1 <k 2 c k 1 c k 2 E ε t-k 1 ε t-k 2 , (6.20)
where the coefficients (c k ) k≥0 are determined by equation (2.13). By representation (2.2), we have (ε t ) t is a measurable function of i.i.d. random variables Z t 's with mean zero and variance 1. Hence, the second term of equation (6.20) is equal to zero. Thus, we have:

E X 2 t = ∞ k=0 c 2 k E ε 2 t-k =        ∞ k=0 c 2 k        2 E ε 2 t . (6.21)
Under the conditions (2.9), we have c 2 k < ∞, and we know that the right side of equation (6.21) is finite, and so (X t ) t is second-order stationary. If (ε t ) t is a function of i.i.d random variables, hence so is (X t ) t . Consequently, (X t ) t is strictly stationary and ergodic.

b) The proof is similar to a) except that the conditions are required on the convergence of

Ψ 2 (z) = Φ(z s )φ(z) Θ(z s )θ(z) (1 -z) d (1 -z s ) D , for all |z| ≤ 1.
A similar argument as in a) can be used to show that the power series expansion of Ψ 2 (z) is convergent for all |z| ≤ 1, and so (X t ) t is invertible and we know that (2.11) Proof. The proof of Lemma 6.1 follows mainly the lines of the proof of Theorem 3.1. in Ling and Li [START_REF] Ling | On Fractionally Integrated Autoregressive Moving-Average Time series Models with Conditional Heteroscedasticity[END_REF].

To prove the consistency and asymptotic normality of the CSS estimator, we will check the conditions listed by Basawa et al. [START_REF] Basawa | Asymptotic Properties of Maximum Likelihood Estimators for Stochastic Processes[END_REF]. They have analyzed the asymptotic properties of ML estimator in process with dependent observations, by given a set of sufficient conditions. Thus to prove condition (a), we need to check the following three conditions:

(1) -n -1 ∂l t (ω 0 ) ∂ω P -→ 0;

(2) There exists a nonrandom matrix M(ω 0 ) > 0 such that for all ε > 0,

P -n -1 ∂ 2 l t (ω 0 ) ∂ω∂ω T ≥ M(ω 0 ) > 1 -ε, for all n ≥ n 1 (ε);
(3) There exists a constant M < ∞ such that

E ∂ 3 l t (ω) ∂ω i ∂ω j ∂ω k < M, for all ω ∈ Λ,
where ω i is the ith component of ω.

Consider the first-order derivatives of the function l t with respect to the parameters. Thus, we have for each t:

∂l t ∂γ = 1 2h 2 t ε 2 t h 2 t -1 ∂h 2 t ∂γ - ε t h 2 t ∂ε t ∂γ and ∂l t ∂δ = 1 2h 2 t ε 2 t h 2 t -1 ∂h 2 t ∂δ . (6.22)
These derivatives involve that of ε t and h 2 t , which we now specify as follows:

∂ε t ∂d = - ∞ k=1 1 k ε t-k ; ∂ε t ∂D = - ∞ k=1 1 k ε t-ks ; ∂ε t ∂φ j = -φ -1 (B)ε t-j , j = 1, . . . , p; ∂ε t ∂θ j = -θ -1 (B)ε t-j , j = 1, . . . , q; ∂ε t ∂Φ j = -Φ -1 (B s )ε t-js , j = 1, . . . , P; ∂ε t ∂Θ j = -Θ -1 (B s )ε t-js , j = 1, . . . , Q; ∂h 2 t ∂γ = 2 r i=1 a i ε t-i ∂ε t-i ∂γ + m j=1 b j ∂h 2 t-j ∂γ + r * k=1 c k ε t-k ∂h t-k ∂γ + r * k=1 c k h t-k ∂ε t-k ∂γ ,
and similarly, we have

∂h 2 t ∂δ = εt + m j=1 b j ∂h 2 t-j ∂δ + r * k=1 c k ε t-k ∂h t-k ∂δ , where εt = (1, ε 2 t-1 , . . . , ε 2 t-r , h 2 t-1 , . . . , h 2 t-m , ε 2 t-1 h 2 t-1 , . . . , ε 2 t-r * h 2 t-r * ) T .
Using the expression in equation (6.22), we have E ∂l t ∂ω ω=ω 0 = 0. According the assumptions of Theorem 2.2 and Corollary 2.1, the process (X t ) t∈Z and (ε t ) t∈Z are strictly stationary and ergodic. Hence, applying the ergodic theorem, we obtain condition (1). To check condition (2), we need to compute the derivatives of order two as follows:

∂ 2 l t ∂γ∂γ T = - 1 h 2 t ∂ε t ∂γ ∂ε t ∂γ T - 1 2h 4 t ∂h 2 t ∂γ ∂h 2 t ∂γ T ε 2 t h 2 t + ε 2 t h 2 t -1 ∂ ∂γ 1 2h 2 t ∂h 2 t ∂γ T - 2ε t h 2 t ∂ε t ∂γ ∂h 2 t ∂γ T + ε t h 2 t ∂ 2 ε t ∂γ∂γ T , and
∂ 2 l t ∂δ∂δ T = - 1 2h 4 t ∂h 2 t ∂δ ∂h 2 t ∂δ T ε 2 t h 2 t + 1 2 ε 2 t h 2 t -1 ∂ ∂δ 1 h 2 t ∂h 2 t ∂δ T .
Similarly, we find (∂ 2 l t )/(∂γ∂δ T ) and (∂ 2 l t )/(∂δ∂γ T ).

According to Lemma 6.1, the matrix Ω 0 is positive definite, and hence for any constant vector C 0,

lim n→+∞ P n -1 C T ∂ 2 l t (ω 0 ) ∂ω ∂ω T C + C T Ω 0 C ≤ ∆ = 1.
In other words, we have ∀ ε > 0, there exists n 1 = n 1 (ε) such that ∀ n ≥ n 1 :

P n -1 C T ∂ 2 l t (ω 0 ) ∂ω ∂ω T C + C T Ω 0 C ≤ ∆ -1 < ε =⇒ P n -1 C T ∂ 2 l t (ω 0 ) ∂ω ∂ω T C + C T Ω 0 C ≤ ∆ > 1 -ε 14
=⇒ P n -1 C T ∂ 2 l t (ω 0 ) ∂ω ∂ω T C + C T Ω 0 C ≤ C T Ω 0 C 2 > 1 -ε =⇒ P n -1 C T ∂ 2 l t (ω 0 ) ∂ω ∂ω T C ≤ - C T Ω 0 C 2 > 1 -ε =⇒ P -n -1 C T ∂ 2 l t (ω 0 ) ∂ω ∂ω T C ≥ C T Ω 0 C 2 > 1 -ε
So the condition ( 2) is satisfied. To prove condition (3), we remark that:

∂ 3 ε t ∂d 3 = ln 3 (1 -B)ε t =        - ∞ k=1 1 k B k        3 ε t . =⇒ E ∂ 3 ε t ∂d 3 2 = E         k 1 ,k 2 ,k 3 =1 1 k 1 k 2 k 3 ε t-k 1 -k 2 -k 3         2 = k 1 ,k 2 ,k 3 =1 1 k 2 1 k 2 2 k 2 3 E ε 2 t-k 1 -k 2 -k 3 = π 2 6 3 E ε 2 t < ∞
Similarly, we can show that E ∂ 3 ε t ∂ω 3 i 2 < ∞, where ω i is the ith component of ω. Thus, differentiating ∂ 2 l t /∂ω∂ω T and using the preceding discussion, we obtain condition (3).

Finally, using the first-order Taylor expansion of the expression ∂ L(ω) ∂ ω at ω = ω 0 , and the conditions (1), ( 2) and (3), we obtain the existence and the consistency of the CSS estimator ωn .

To prove condition (b) of Theorem 3.1, we remark at first, by Lemma 6.1, we have: where η 0 is an arbitrary constant vector with η T 0 η 0. We have

E η T 0 ∂l t (ω 0 ) ∂ω = η T 0 E ∂l t (ω 0 ) ∂ω = 0.
Then η T 0 ∂l t (ω 0 ) ∂ω t are a sequence of martingale differences and S n is the sum of martingale difference.

Hence S n is a martingale and we have 1 n E(S 2 n ) = η T 0 Ω 0 η > 0. From the stationarity and the ergodicity of (X t ) and (ε t ), we have: Using the first-order Taylor expansion of the expression ∂ L(ω) ∂ ω at ω = ω 0 , and the conditions (6.23), (6.24) and (6.25), we have:

1 n E S 2
√ n ( ωn -ω 0 ) D -→ N(0, Ω -1 0 ), as n -→ +∞. The condition (c) of Theorem 3.1 is directly obtained from Lemma 6.1.

  ), has the following properties (a) The CSS estimator ωn exists and satisfies: ∂L(ω) ∂ω = 0 and ωn P -→ ω 0 , as n -→ +∞. (b) For such a sequence, √ n ( ωn -ω 0 ) D -→ N(0, Ω -1 0 ), as n -→ +∞, where P -→ and D -→ denote respectively convergence in probability and in distribution. Furthermore,

Figure 1 :

 1 Figure 1: Boxplot of estimates of d = 0.1, D = 0.3, a 0 = 0.01, a 1 = 0.09, b 1 = 0.9 and c 1 = 0.15 in Model 1, with sample size n = 100, n = 300, n = 500 and n = 1000, based on 1000 replications.

  2 l t (ω 0 ) ∂ω ∂ω T C a.s -→ -C T Ω 0 C. For any given C, let 0 < ∆(C) < C T Ω 0 C/2. Then ∀ ε > 0, lim n→+∞ P n -1 C T ∂ 2 l t (ω 0 ) ∂ω ∂ω T C + C T Ω 0 C > ε = 0Particularly, we have for ε = ∆:

  (0, Ω 0 ), as n -→ +∞.(6.25)

Table 1 :

 1 Data generating processes (DGPs).

	DGPs	d	D	φ	θ	a 0	a 1	b 1	c 1
	Model 1	0.10	0.30	-	-	0.01	0.09	0.9	0.15
	Model 2	0.10	0.30	0.50	-	0.01	0.09	0.9	0.15
	Model 3	0.10	0.30	-	0.30	0.01	0.09	0.9	0.15
	Model 4	0.10	0.30	0.50	0.30	0.01	0.09	0.9	0.15

Table 2 :

 2 CSS estimation for Model 1.

	sizes	d	D	a 0	a 1	b 1	c 1
	100	0.0980	0.3137	0.0617	0.0825	0.8156	0.1844
		(0.0732)	(0.0893)	(0.1022)	(0.0945)	(0.2286)	(0.1187)
		[0.0605]	[0.0727]	[0.0557]	[0.0750]	[0.1433]	[0.0903]
	300	0.1013	0.3158	0.0190	0.0857	0.8883	0.1552
		(0.0503)	(0.0520)	(0.0246)	(0.0490)	(0.0721)	(0.0529)
		[0.0408]	[0.0419]	[0.0119]	[0.0378]	[0.0462]	[0.0415]
	500	0.1019	0.3125	0.0140	0.0897	0.8926	0.1541
		(0.0387)	(0.0413)	(0.0102)	(0.0335)	(0.0383)	(0.0409)
		[0.0305]	[0.0326]	[0.0060]	[0.0261]	[0.0289]	[0.0320]
	1000	0.0996	0.3097	0.01148	0.0896	0.8972	0.1499
		(0.0264)	(0.0280)	(0.0042)	(0.0198)	(0.0211)	(0.0251)
		[0.0213]	[0.0223]	[0.0030]	[0.0158]	[0.0164]	[0.0202]

Table 3 :

 3 CSS estimation for Model 2.

	sizes	d	D	φ	a 0	a 1	b 1	c 1
	100	0.1005	0.3116	0.4907	0.0763	0.0763	0.8114	0.1789
		(0.1068) (0.0991) (0.1414) (0.2127) (0.0964) (0.2406) (0.1214)
		[0.0906] [0.0789] [0.1100] [0.0694] [0.0763] [0.1476] [0.0917]
	300	0.0942	0.3083	0.4984	0.0211	0.0837	0.8825	0.1500
		(0.0870) (0.0548) (0.1002) (0.0316) (0.0467) (0.0932) (0.0531)
		[0.0767] [0.0431] [0.0842] [0.0134] [0.0368] [0.0519] [0.0409]
	500	0.0957	0.3074	0.4989	0.0135	0.0876	0.8957	0.1511
		(0.0771) (0.0413) (0.0856) (0.0096) (0.0331) (0.0381) (0.0375)
		[0.0671] [0.0326] [0.0719] (0.0056) [0.0256] [0.0278] [0.0296]
	1000	0.1004	0.3058	0.4973	0.0113	0.0901	0.8978	0.1529
		(0.0634) (0.0284) (0.07138) (0.0043) (0.0205) (0.0216) (0.0261)
		[0.0526] [0.0229] [0.0584] [0.0031] [0.0159] [0.0166] [0.0207]
					8

Table 4 :

 4 CSS estimation for Model 3.

	sizes	d	D	θ	a 0	a 1	b 1	c 1
	100	0.1265	0.3116	0.2945	0.0760	0.0778	0.8061	0.1748
		(0.1050) (0.0973) (0.1366) (0.2014) (0.0987) (0.2481) (0.1191)
		[0.0863] [0.0779] [0.1090] [0.0692] [0.0773] [0.1514] [0.0886]
	300	0.1268	0.3064	0.3011	0.0199	0.0826	0.8882	0.1508
		(0.0792) (0.0518) (0.0844) (0.0271) (0.0466) (0.0778) (0.0522)
		[0.0658] [0.0413] [0.0686] [0.0123] [0.0363] [0.0470] [0.0400]
	500	0.1350	0.3061	0.3118	0.0135	0.0871	0.8961	0.1481
		(0.0738) (0.0409) (0.0758) (0.0086) (0.0325) (0.0357) (0.0366)
		[0.0612] [0.0323] [0.0612] [0.0056] [0.0257] [0.0270] [0.0289]
	1000	0.1381	0.3000	0.3181	0.0110	0.0896	0.8972	0.1446
		(0.0615) (0.0292) (0.0564) (0.0044) (0.0199) (0.0221) (0.0258)
		[0.0509] [0.0228] [0.0458] [0.0030] [0.0158] [0.0170] [0.0207]

Table 5 :

 5 CSS estimation for Model 4.

	sizes	d	D	φ	θ	a 0	a 1	b 1	c 1
	100	0.0904	0.3137	0.4753	0.2626	0.0790	0.0789	0.7974	0.1830
		(0.1081) (0.1003) (0.3194) (0.3437) (0.2122) (0.1002) (0.2625) (0.1277)
		[0.0931] [0.0805] [0.2336] [0.2599] [0.0721] [0.0785] [0.1603] [0.0947]
	300	0.0841	0.3107	0.4650	0.2453	0.0199	0.0810	0.8901	0.1555
		(0.0840) (0.0527) (0.2680) (0.2717) (0.0251) (0.0465) (0.0734) (0.0527)
		[0.0735] [0.0415] [0.1867] [0.1935] [0.0122] [0.0365] [0.0460] [0.0411]
	500	0.0880	0.3116	0.4728	0.2588	0.0139	0.0862	0.8963	0.1526
		(0.0775) (0.0421) (0.2332) (0.2302) (0.0090) (0.0320) (0.0355) (0.0378)
		[0.0672] [0.0334] [0.1603] [0.1602] [0.0058] [0.0255] [0.0269] [0.0297]
	1000	0.0924	0.3053	0.4560	0.2493	0.0116	0.0896	0.8964	0.1485
		(0.0679) (0.0303) (0.2216) (0.2016) (0.0047) (0.0197) (0.0224) (0.0262)
		[0.0582] [0.0238] [0.1459] [0.1347] [0.0032] [0.0155] [0.0170] [0.0211]

  holds.

		∂ε t ∂γ	∂ε t ∂γ T +	1 2h 4 t	∂h 2 t ∂γ	∂h 2 t ∂γ T	and	Ω δ = E	1 2h 4 t	∂h 2 t ∂δ	∂h 2 t ∂δ T ,
	where	a.s

-→ denotes convergence almost surely.

Boxplot of estimates in Model 2, Model

and Model

are omitted in the paper but are available to authors upon request.

Appendix : Proofs of Theorems

Proof of Theorem 2.2

The proof of this theorem follows mainly the proof of Theorem 2.3 in Ling and Li [START_REF] Ling | On Fractionally Integrated Autoregressive Moving-Average Time series Models with Conditional Heteroscedasticity[END_REF]. a) The process in (2.1) can be written as Proof of Theorem 3.1 The proof of Theorem 3.1 uses the following lemma. Lemma 6.1. Suppose that (X t ) t∈Z are generated by equations (2.1), (2.2) and (2.3), and the assumptions of Theorem 3.1 hold. Then


Mor Ndongo

Kâ Abdou

Aliou Diongue

Ouagnina Diop

Hili

Estimation of Fractional ARIMA Process with Stable Innovations: A Monte Carlo Study

Keywords: Stable distributions, Long memory, Whittle estimates, CSS method, MHD estimator, Markov chains Monte Carlo

.

Introduction

Modelling time series using fractionally integrated processes with finite variance has been the purpose of many investigations. See for example, [START_REF] Granger | An introduction to long memory time series and fractional differencing[END_REF] and [START_REF] Hosking | Fractional Differencing[END_REF] who introduced the Gaussian FARIMA(p, d, q) model. However, most of time series in real life may display high variability, in addition to long-range dependence structure. To take into account these stylized facts, the theory of αstable FARIMA(p, d, q) process was introduced and developed by [START_REF] Kokoszka | Fractional ARIMA with stable innovations[END_REF]. Examples where such models appears to be appropriate have been found by [START_REF] Stanislavsky | FARIMA modeling of solar flare activity from empirical time series of soft X-ray solar emission[END_REF], who described solar flare activity, and [START_REF] Agnieszka | Robust Wavelet-Domain Estimation of the Fractional Difference Parameter in Heavy-Tailed Time Series: An Empirical Study[END_REF] or [START_REF] Ilow | Forecasting network traffic using FARIMA models with heavy-tailed innovations[END_REF] who analyzed network traffic data.

The parameter estimation is a crucial question for the use of these parametric models. The usual response to such difficulties is that the density functions of these processes cannot be written in a closed form in the sense that it is not expressible in terms of known elementary functions. As a result, Maximum Likelihood estimation is often very difficult (for detailed discussion, see [START_REF] Calder | Inferences for linear processes with stable noise[END_REF]). The estimation via alternative methods such as Quasi-Maximum Likelihood (QML) method and Generalized Method of Moments (GMM) also presents difficulties. For example, QML is infeasible because the variance of the error term may be infinite. For GMM care must be taken when choosing moment conditions because the stable distribution does not have a finite absolute moment of order higher than the tail index α.

To solve the parameter estimation problems of the stable FARIMA(p, d, q) process, [START_REF] Kokoszka | Parameter estimation for infinite variance fractional ARIMA[END_REF] have used a variant of Whittle's method. They show the consistency of the resulting estimators and derive its asymptotic distribution. Otherwise, [START_REF] Diongue | Estimation of k-Factor GIGARCH Process: A Monte Carlo Study[END_REF] have developed the theory of seasonal fractional ARIMA model with infinite variance innovations which is a direct generalization of stable FARIMA model. To achieve the parameter estimation of this model, [START_REF] Ndongo | Estimation of long-memory parameters for seasonal fractional ARIMA with stable innovations[END_REF] considered the semiparametric method studied in Reisen et al. (2006) and two Whittle approaches: the classical Whittle method (e.g. [START_REF] Embrechts | Modelling Extremal Events New York[END_REF], [START_REF] Mikosch | Parameter estimation for ARMA models with infinite variance innovations[END_REF], [START_REF] Kokoszka | Parameter estimation for infinite variance fractional ARIMA[END_REF]) and a method based on a Markov chains Monte Carlo (MCMC) procedure. The authors studied their performance using Monte Carlo simulations and show that the best results are obtained by the MCMC-Whittle method.

The main objective of this article is to propose two parametric estimation methodologies for α-stable FARIMA process. For instance, we consider the Minimum Hellinger Distance (MHD) estimator introduced by Beran (1977) and the Conditional Sum of Squares (CSS) approach developed by [START_REF] Chung | Estimating a generalized long memory process[END_REF]. The advantage of using these procedures is that they do not require any prior knowledge about the driving noise distribution. Monte Carlo simulations are performed to study the finite sample properties of the estimators. For comparison purpose, we consider also the MCMC Whittle method developed by [START_REF] Ndongo | Estimation of long-memory parameters for seasonal fractional ARIMA with stable innovations[END_REF].

The paper is organized as follows. The next section briefly reviews the α-stable FARIMA(p, d, q) model and its probabilistic properties. Section 3 presents the different parameter estimation techniques. Section 5 calibrates the performance of the estimation procedures through Monte Carlo experiments while Section 6 provides concluding remarks.

The FARIMA(p, d, q)-SαS model

In order to model a time series with long memory behaviour and infinite variance components, [START_REF] Kokoszka | Fractional ARIMA with stable innovations[END_REF] proposed a class of models, the Fractional Autoregressive Integrated Moving Average (FARIMA, in short) time series with symmetric α-stable (SαS) innovations. These models, denoted hereafter by FARIMA(p, d, q)-SαS, where p, q ∈ N, d a fractional number and α ∈ (1, 2], are an infinite variance counterpart to the Gaussian FARIMA model introduced by [START_REF] Granger | An introduction to long memory time series and fractional differencing[END_REF] and [START_REF] Hosking | Fractional Differencing[END_REF]. In this Section, we examine the definition and the basic characteristics of an stable FARIMA model, and we refer to [START_REF] Kokoszka | Fractional ARIMA with stable innovations[END_REF] or [START_REF] Stoev | Simulation methods for linear fractional stable motion and FARIMA using the Fast Fourier Transform[END_REF] for a full description. Thus, we assume that:

0 < d < 1 - 1 α , and 1 < α ≤ 2. (2.1)
We define a centered fractional ARIMA process (Y t ) t∈Z with SαS innovations, as the stationary solution to the back-shift operator equation:

Φ(B)Y t = Θ(B)(1 -B) -d Z t , t ∈ Z, (2.2)
where d is (fractional) differencing degree, Φ and Θ are autoregressive and moving average polynomials of degrees p, q, respectively:

Φ(z) = 1 -φ 1 z -φ 2 z 2 -• • • -φ p z p , and Θ(z) = 1 + θ 1 z + θ 2 z 2 + • • • + θ q z q . (2.3)
It is assumed that these polynomials have no common zeros, and have all its roots outside the unit disk. The innovations (Z t ) t∈Z in (2.2) are independently and identically distributed SαS random variables, with 1 < α ≤ 2. The tail exponent α determines the rate at which the tails of the distribution tape off. When α = 2, the Gaussian distribution is obtained and when α < 2, the variance is infinite. On the other hand, the probability density function of a symmetric standard α-stable random variable can be expressed as:

f (x; α) = 1 π +∞ 0 exp (-|t| α ) cos (tx) dt.
(2.4) Formula (2.4) does not have closed form expression, except in the three cases: Levy, Cauchy and Gaussian distributions corresponding respectively to the cases: α = 1/2, α = 1 and α = 2. However, it can be numerically integrated. All the relevant properties of stable distributions can be found in [START_REF] Samorodnitsky | Stable non-Gaussian random processes : Stochastic models with infinite variance[END_REF].

The following stationarity and invertibility conditions of model (2.2) are established and proved in [START_REF] Kokoszka | Fractional ARIMA with stable innovations[END_REF].

• Under condition (2.1), the process (Y t ) t∈Z is stationary and have a unique moving average representation

Y t = ∞ j=0 c j Z t-j .
(2.5)

The coefficients (c j ) j≥0 in (2.5) satisfy the following relation

c 0 = 1 and c j = b j (d) + min( j,p) i=1 φ i c j-i + min( j,q) i=1 θ i b j-i (d), j ≥ 1, (2.6)
where the weights (b j (d)) j≥0 are given by:

b 0 (d) = 1 and b j (d) = Γ( j + d) Γ(d)Γ( j + 1) , j = 1, 2, . . . . (2.7)
• In addition, if α > 1 and |d| < 1 -1 α the solution is invertible and we have:

Z t = ∞ j=0 c j Y t-j , (2.8)
where the coefficients (c j ) j≥0 are defined by: c0 = 1 and c j = b j (-d) -

min( j,p) i=1 φ i b j-i (-d) + min( j,q) i=1 θ i c j-i , j ≥ 1.
(2.9)

• The power transfer function of the process (Y t ) t∈Z is

g(λ) = |Θ(e iλ )| 2 |Φ(e iλ )| 2 |1 -e iλ | -2d , for -π ≤ λ ≤ π.
(2.10)

Estimation methods

In this section, we discuss briefly the application of three methods to the estimation of fractional ARIMA model with α-stable distributions. The first one is the Markov Chains Monte Carlo (MCMC) Whittle approach, the second deals with the Conditional Sum of Squares (CSS) procedure while the third is based to the Minimum Hellinger Distance Estimator (MHDE).

Let Y 1 , . . . , Y T be an observed finite sequence generated by a stationary and invertible symmetric αstable FARIMA process defined by equation (2.2). Assume that the tail index parameter α is known and denote by ψ = d, φ 1 , . . . , φ p , θ 1 , . . . , θ q the (p + q + 1)-dimensional vector of interest, where φ 1 , . . . , φ p and θ 1 , . . . , θ q are the coefficients of the autoregressive and moving average polynomials respectively, and d is the long memory parameter. Assume that ψ 0 is the true value of ψ and is in the interior of the compact set Ψ ⊆ R p+q+1 .

The MCMC Whittle method

This estimation method is developed by [START_REF] Ndongo | Estimation of long-memory parameters for seasonal fractional ARIMA with stable innovations[END_REF] to estimate the long memory parameters for seasonal fractional ARIMA model with SαS innovations, which is a direct generalization of FARIMA-SαS model. It is an alternative to the classical Whittle approach used by Klüppelberg andMikosch (1993, 1994), [START_REF] Mikosch | Parameter estimation for ARMA models with infinite variance innovations[END_REF] and [START_REF] Kokoszka | Parameter estimation for infinite variance fractional ARIMA[END_REF]. The procedure is based on approximation of the Whittle likelihood function using Bayesian approach. In this work, this estimator is denoted by ψ W and can be obtained by minimizing the following likelihood function:

L W (ψ) = 1 N N j=1 1 g(λ j , ψ) , (3.11)
where g(λ, ψ) is the power transfer function of the process defined by equation (2.10). The sequences λ 1 , . . . , λ N are generated using a Metropolis-Hastings algorithm. A detailed revision of this procedure is found in Section 3.2 of [START_REF] Ndongo | Estimation of long-memory parameters for seasonal fractional ARIMA with stable innovations[END_REF].

The CSS method

The CSS approach is studied by [START_REF] Chung | Small sample bias in conditional sum-of-squares estimators of fractionally integrated ARMA models[END_REF] to estimate the parameters of Gaussian FARIMA model. As remarked by [START_REF] Chung | Estimating a generalized long memory process[END_REF], the advantage of using the CSS estimation is that the normality assumption is not essential. Hence, this method can overcome difficulties arising from the lack of closed expression for the likelihood. Consequently, it seems to be a better alternative to the estimation of stable FARIMA process. Thus, the conditional sum of squares estimator ψ CS S of ψ 0 is the value of ψ which minimizes

S (ψ) = T t=1 Z t (ψ) 2 .
However, given the observation Y 1 , . . . , Y T the innovation (Z t ) t=1,...,T cannot be directly computed, since an infinite sample would be needed. Nevertheless, they may be estimated by:

Z t (ψ) = t-1 j=0 c j (ψ)Y t-j , t = 1, . . . , T, (3.12)
where the coefficients (c j (ψ)) j≥0 can be determined by equation (2.9).

Minimum Hellinger Distance Estimator

The Minimum Hellinger Distance (MHD in short) method is initially introduced by Beran (1977) for parametric models, in the i.i.d case. The procedure based on dependent data is studied in [START_REF] Hili | On the estimation of nonlinear time series models[END_REF]. Recently, Kamagate and Hili (2012) use the MHD method to estimate the parameters of Gaussian FARIMA process. The basic idea for this approach is to minimize the Hellinger distance between the probability density function of the innovations of the process and a parametrized random function. In this section, we investigate the MHD procedure to estimate the parameters of stable FARIMA model. To do so, let's denote by f (.) = f ψ 0 (.) the probability density function of the innovation (Z t ) t∈Z . In order to construct a MHD estimator of ψ 0 , we need a nonparametric estimator of the density function f ψ 0 (.), which can be obtained using the nonparametric kernel estimator of form:

f ψ, T (x) = 1 T h T T t=1 K       x -Z t (ψ) h T       , (3.13)
where Z t (ψ) is an estimate of Z t (ψ) for t = 1, . . . , T , given by equation (3.12), (h T ) is a sequence of bandwidths and K is a kernel function on R.

The Hellinger distance between the density function f ψ 0 and the nonparametric kernel estimator f ψ,T is defined as the L 2 -norm of the difference square root of these two functions, i.e.

HD T f ψ,T , f ψ 0 = f 1 2 ψ,T -f 1 2 ψ 0 2 = R f 1 2 ψ,T (x) -f 1 2 ψ 0 (x) 2 dx 1/2 . (3.14)
The minimum Hellinger distance estimator of ψ 0 is the value ψ MHD that minimizes (3.14), with respect to ψ ∈ Ψ, that is,

ψ MHD = arg min ψ∈Ψ HD T f ψ,T , f ψ 0 .
To implement the MHD estimator, we take the kernel function K to be the density of N(0, 1) distribution.

For the choice of the bandwidth h T , we consider the bandwidth proposed in Wu and Mielniczuk (2002) for long-range dependence sequences, that is:

h T = T ν L(T ), (3.15)
where -1 < ν < 0 and L(.) is a slowly varying function at ∞. However, it is important to recognize that when using the bandwidth in (3.15), an additional choice needs to be made, which is that the optimal value of ν. The simulation in this paper (Section 5.1) show that the MHD method is sensitive to the choice of ν. Thus, a data-driving selection method of ν would be necessary. Ideally, ν should be selected so that it minimizes the mean-square error (MSE) of the MHD estimator. Precisely, let ψ (ν) be the MHD estimator for a given ν in ] -1, 0[. Then, we look for ν opt such that:

ν opt = arg min ν 1 R R j=1 ψ (ν) j -ψ 0 2 , (3.16)
where R is the number of replications, ψ 0 is the nominal value of ψ, and ψ (ν) j is the estimate for sample j.

4. Asymptotic properties of the MHD estimator

Monte Carlo simulations

In this section, we provide numerical assessment of the estimation techniques discussed in Section 3, using Monte Carlo experiments. To appreciate the precision of the estimators, we evaluate the mean, the root mean square error (RMSE) and the mean absolute error (MAE) of the estimation procedures. The α-stable FARIMA models are generated using the Fast Fourier Transform algorithm described in [START_REF] Stoev | Simulation methods for linear fractional stable motion and FARIMA using the Fast Fourier Transform[END_REF]. All calculations are carried out using an R programming environment.

Monte Carlo study comparing CSS, MHDE and Whittle

As data generating processes, we use three α-stable FARIMA models listed in Table 1, with α = 1.2. For each model, we generate R = 1500 replications of lengths T = 250, T = 500 and T = 1500. For the MHD estimator, we set L(T ) = log(T ) and we use several values of ν in the interval ] -1, 0[ with a step of 0.02, and we determine the optimal value of ν suitable for each model in Table 1, via the criteria given in (3.16). Thus, the MHD estimates obtained with this optimal value of ν are compared to the CSS and the MCMC Whittle procedures. All the results are presented in Tables 234. In these tables, the sample sizes are given in the first columns, the estimations of the parameters in the next columns, the root mean square error (RMSE) in the row below (in parentheses) and the mean absolute error (MAE) under the row of the RMSE (in brackets). • In the first experiment, we focus our attention to the estimation of long memory parameter. In this design, we consider Model 1 in Table 1. To show the impact of the bandwidth on the MHD estimation, we select some values of ν = {-0.9, -0.5, -0.3, -0.2, -0.1} and compare it with the CSS and Whittle methods in Table 2. Tables for other values of ν are available upon request. We observe that all these three methods perform well, since the RMSE and MAE of the estimates are in general small. We note that the MHD estimator can be worse or better than the CSS and Whittle estimators depending on the values of ν considered. Indeed, the MHD procedure is the best only for ν = -0.3, ν = -0.2 and ν = -0.1. The simulation results show the impact of the bandwidth on the MHD estimator. In terms of RMSE and MAE, the optimal bandwidth is obtained for ν = -0.1. However, the gain from ν = -0.1 to ν ≥ -0.2 is marginal, as it can be seen in Figure 1 which represent the RMSE of the estimates dMHDE for a given value of ν and sample size T = 250. We remark also that the CSS method performs better than the MCMC Whittle procedure. The Monte Carlo experiments indicate also the impact of the sample sizes on these estimation methods. Indeed, when the the sample size increases, the results improve significantly too. To check the asymptotic distribution of d, we plot in Figure 2, the boxplot of the 1500 estimates. The vertical axis indicates the deviation from the nominal value of d. From this figure, we observe that the scatter of the CSS and Whittle methods is larger than those in the MHD estimator. However, this dispersion improves when the sample size T increases. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 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.00 0.05 0.10 0.15 0.20 0.25 0.30 deviation from nominal value • Now consider the previous model in which we introduce a short memory part (i.e. Model 2 and Model 3 in Table 1). This experiment is designed to examine the relative performance in the estimation techniques, when there are long memory and short memory components simultaneously. As in the previous experiment, we choose at first the optimal bandwidth for the MHD estimator, and compare the results to CSS and Whittle procedures. Thus using the criteria in (3.15), we take ν = -0.1 and ν = -0.2 for Model 2 and Model 3, respectively. The results are gathered in Tables 3 and4. They show that, when short memory components are introduced in Model 1, the estimation of the long memory parameter is disturbed. Indeed, the RMSE and MAE obtained in Tables 3 and4 are larger than those presented in Table 2. This phenomenon is already observed in the literature (see for example [START_REF] Boutahar | Estimation Methods of the Long Memory Parameter: Monte Carlo Analysis and Application[END_REF] or Diongue and Gégan (2008)). However, we can note a significant improvement when the sample size becomes larger (T = 500, T = 1500), only for the MMHD and CSS methods. Concerning the MHD estimator, we see that the averages of dMHDE are closed to the true values, moreover it is always better than the CSS method, and it dominates the Whittle procedure. Hence, the MHD estimator shows more resistance to the presence of the short memory components than the CSS and Whittle approaches. For the Whittle estimator, there is a wide different between the true value of d and the average dW , and the results still remain bad, even for larger sample sizes (T = 500, T = 1500). Consequently, the AR and MA components can be a source of the bias of the Whittle estimator. It is important to note that the results obtained for the CSS method are acceptable.

Parameters

MHDE CSS Whittle d d d d d d d sizes ν = -0.9 ν = -0.5 ν = -0.3 ν = -0.2 ν = -0.

Robustness of the MHD estimator

Here, we investigate the robustness of the MHD estimator. The most commonly accepted mathematical method for treating contaminated data comes from to Huber (1964[START_REF] Huber | Robust Statistics[END_REF]; it is this model upon which much of the robustness approach is based. Broadly speaking, the contaminated model addresses outliers by positing that the observed data belongs to a mixture of two distributions. The mixture is of the form:

f (λ) ψ, T (r) = (1 -λ) f ψ, T (r) + λ U [0, 1] (r), for ψ ∈ Ψ, where 0 < λ < 1, f (λ)
ψ, T (r) is the contaminated density function and U [0, 1] (r) is the uniform density function on the interval [0, 1]. Then, we next replaced f ψ, T (.) by f (λ) ψ, T (.) in equation (3.14) and we calculated the MHD estimator using several values of λ. In this experiment, we consider an FARIMA(0, d, 0)-SαS model, with α = 1.6, d = 0.2. We set the sample size to T = 250 and we choose λ = {0.1, 0.2, 0.25, 0.3, 0.4, 0.45, 0.5}. The results are summarized in Table 4. Noticeably, the MHD estimator provides good estimates of d, when we use a contaminated density. Indeed, we can see that the averages of d are closed to the true value of d for all value of λ. Thus, the MHD method is robust to the perturbation. This properties of the MHD estimate is also pointed in Beran (1977).

Conclusion

In this article, we have dealt with a special class of long-memory model with stable innovations. Three parameter estimation techniques for FARIMA(p, d, q)-SαS process have been considered. Our empirical Monte Carlo investigation indicates that the two estimators proposed in this work, namely the MHD and CSS methods perform well even for small sample sizes. Consequently, they are reasonable methods for dealing with long-range dependent data with infinite variance components. The comparison of the different estimators showed that the MHD estimator is the best one, in terms of RMSE and MAE. The simulation results carried out the fact that the estimator of the long memory parameter is disturbed when there is the presence of short memory components, for all three methods. In this context, the results improve when the sample size increases, only for the MHD and CSS procedures, and it still remain bad for the Whittle approach. This paper focuses on the estimation of long memory time series with stable innovations, via Monte Carlo experiments. It will be interesting to investigate the theoretical asymptotic properties for these estimation methods. These problems will be examined on a forthcoming paper.

Figure 2 :

 2 Figure 2: Boxplot of estimates of d in Model 1 under MHD, CSS and MCMC Whittle estimators with sample sizes T = 250, T = 500 and T = 1500, based on 1500 replications.

Table 1 :

 1 Data generating processes for fixed α = 1.2.

		d	φ	θ
	Model 1	0.15	-	-
	Model 2	0.15	0.60	-
	Model 3	0.15	-	0.40

Table 2 :

 2 Monte Carlo study comparing CSS, MHDE and Whittle for FARIMA(0, d, 0)-SαS model, with d = 0.15, for α = 1.2.

	0.05
	0.04
	0.03
	0.02
	-1.0	-0.8	-0.6	-0.4	-0.2	0.0

RMSE

Figure 1: RMSE of estimates of d in Model 1, under MHD estimator and T = 250.

  On the other hand, this estimator has the advantage of being relatively simple to implement, compared to the MHD estimator which requires lot of times in order to minimize the objective function (3.14).

		MHDE		CSS		Whittle
	Sizes	d	φ	d	φ	d	φ
	250	0.144	0.627	0.119	0.640	0.034	0.656
		(0.039) (0.086)	(0.071) (0.094)	(0.124) (0.116)
		[0.024] [0.042]	[0.048] [0.062]	[0.116] [0.101]
	500	0.149	0.617	0.126	0.630	0.050	0.637
		(0.031) (0.077)	(0.061) (0.081)	(0.109) (0.099)
		[0.018] [0.031]	[0.040] [0.050]	[0.100] [0.082]
	1500	0.152	0.604	0.139	0.611	0.061	0.622
		(0.017) (0.040)	(0.038) (0.043)	(0.095) (0.074)
		[0.010] [0.015]	[0.025] [0.027]	[0.089] [0.058]

Table 3 :

 3

		MHDE		CSS		Whittle
	Sizes	d	θ	d	θ	d	θ
	250	0.154	0.409	0.142	0.416	0.063	0.414
		(0.014) (0.050)	(0.033) (0.064)	(0.090) (0.064)
		[0.012] [0.022]	[0.022] [0.039]	[0.086] [0.046]
	500	0.154	0.407	0.147	0.406	0.069	0.405
		(0.012) (0.042)	(0.024) (0.046)	(0.082) (0.047)
		[0.009] [0.017]	[0.017] [0.025]	[0.080] [0.033]
	1500	0.152	0.403	0.150	0.4011	0.073	0.402
		(0.007) (0.038)	(0.015) (0.026)	(0.077) (0.033)
		[0.004] [0.008]	[0.010] [0.013]	[0.076] [0.023]

Monte Carlo study comparing CSS, MHDE (ν = -0.1) and Whittle for FARIMA(1, d, 0)-SαS model, with d = 0.15, φ = 0.6, for α = 1.2.

Table 4 :

 4 Monte Carlo study comparing CSS, MHDE (ν = -0.2) and Whittle for FARIMA(0, d, 1)-SαS model, with d = 0.15, θ = 0.4, for α = 1.2.

Table 5 :

 5 Robustness of the MHDE (ν = -0.1) with d 0 = 0.20, for α = 1.6 and T = 250