Youssef Esstafa
email: youssef.esstafa@univ-lemans.fr

Generalized least squares estimation of fractional autoregressive models

Keywords: AMS 2000 subject classifications: Primary 62M10, secondary 91B84 Fractional Brownian motion, Fractional Gaussian noise, FrAR models, Generalized least squares estimator, Consistency, Asymptotic normality

In this paper we derive the asymptotic properties of the generalized least squares estimator (GLSE) of autoregressive models endowed with fractional Gaussian noise (the socalled fractional autoregressive models). We establish the consistency and the asymptotic normality of the GLSE. Some simulation studies are presented to corroborate our theoretical work.

Introduction

The error term in the classical autoregressive models is usually assumed to be a strong white noise (namely an independent and identically distributed (iid for short) sequence of random variables with mean 0 and common variance), a martingale difference or a stationary process with a bounded spectral density (see [START_REF] Anderson | Strong consistency of least squares estimates in normal linear regression[END_REF], [START_REF] Anderson | Strong consistency of least squares estimates in dynamic models[END_REF], [START_REF] Lai | Strong consistency of least-squares estimates in regression models[END_REF], [START_REF] Lai | Strong consistency of least squares estimates in multiple regression[END_REF], [START_REF] Lai | Strong consistency of least squares estimates in multiple regression ii[END_REF] and [START_REF] Solo | Strong consistency of least squares estimators in regression with correlated disturbances[END_REF]). In the late 1990s, Francq, C. and Zakoïan, J.-M. have released the assumption of the independence on the noise and have considered the problem of estimating ARMA models with uncorrelated but non-independent error term, the so-called weak ARMA models (see [START_REF] Francq | Estimating linear representations of nonlinear processes[END_REF]). This extension allowed ARMA models to be less restrictive and to cover a large class of general nonlinear models. Nevertheless, the uncorrelatedness of the noise assumed in these classes of models considerably restricts their application to short-memory stochastic processes.

Yet, in many scientific disciplines and applied fields such as hydrology, climatology, economics, finance, to name a few, time series exhibit long-range dependence (see for instance [START_REF] Granger | An introduction to long-memory time series models and fractional differencing[END_REF], [START_REF] Fox | Large-sample properties of parameter estimates for strongly dependent stationary Gaussian time series[END_REF], [START_REF] Dahlhaus | Efficient parameter estimation for self-similar processes[END_REF], [START_REF] Hosking | Fractional differencing[END_REF], [START_REF] Beran | Long-memory processes[END_REF], [START_REF] Palma | Long-memory time series[END_REF], among others), that is dependence between distant observations is so strong that classical or weak ARMA models are unable to identify and then to take into consideration.

The fractional Gaussian noise (see [START_REF] Mandelbrot | Fractional Brownian motions, fractional noises and applications[END_REF]) and the fractional autoregressive integrated moving average (FARIMA, for short) models (see [START_REF] Granger | An introduction to long-memory time series models and fractional differencing[END_REF] and [START_REF] Hosking | Fractional differencing[END_REF]) are widely used to model the long memory phenomenon. We shall consider in this work the problem of generalized least squares estimation of autoregressive models endowed with fractional Gaussian noise. More explicitly, we shall study the asymptotic properties of the generalized least squares estimator of the parameters of the centered stationary fractional autoregressive process of order 1 (FrAR(1)) (X t ) t∈Z defined, for all t ∈ Z, by

X t = a 0 X t-1 + H t , (1)
where ( H t ) t∈Z is a fractional Gaussian noise and where the autoregressive parameter a 0 and the Hurst index or Hurst parameter H associated with the fractional Gaussian noise ( H t ) t∈Z are assumed to satisfy the following standard assumption:

(A0):

a 0 ∈ (-1, 1) and H ∈ (0, 1).

The condition a 0 ∈ (-1, 1) ensures the causality and then the second-order stationarity of the process (X t ) t∈Z . In fact, for AR(p) models a causal weak-sense stationary solution is obtained once the roots of the corresponding autoregressive polynomial lie outside the closed unit disk.

The process ( H t ) t∈Z is the increment process of the fractional Brownian motion (B H (t)) t∈R of Hurst index H, namely, for any t ∈ Z, it holds

H t = B H (t + 1) -B H (t).
Based on the definition and the properties of the fractional Brownian motion (B H (t)) t∈R , the process ( H t ) t∈Z is the unique continuous stationary Gaussian centered process with autocovariance function defined for k ∈ Z by

γ H (k) = 1 2 |k -1| 2H -2|k| 2H + |k + 1| 2H . (2)
A Taylor expansion of the function : x → (1 -x) 2H -2 + (1 + x) 2H around 0 implies that for large k,

γ H (k) = 1 2 k 2H (1/k) = H(2H -1)k 2H-2 + o(k 2H-2 ).
The last result implies that γ H (k) → 0 as k → ∞. So the process ( H t ) t∈Z is mixing (see [START_REF] Ito | On the ergodicity of a certain stationary process[END_REF]). Mixing is a stronger property than ergodicity, this result is known as the Khinchin's theorem (see [START_REF] Khinchin | Mathematical foundations of statistical mechanics dover publications[END_REF]). Hence lim k→∞ γ H (k) = 0 implies that the process ( H t ) t∈Z is ergodic. When 1/2 < H < 1, the series of the autocovariances of ( H t ) t∈Z is absolutely divergent and the process ( H t ) t∈Z exhibits long-range dependence.

In the literature, many approaches have been proposed to estimate the parameters of strongly dependent Gaussian processes. We cite for example the R/S method, periodogram estimation and maximum likelihood estimation. [START_REF] Mandelbrot | Limit theorems on the self-normalized range for weakly and strongly dependent processes[END_REF] and [START_REF] Mandelbrot | Robust R/S Analysis of Long Run Serial Correlation[END_REF] have established the theoretical properties of the R/S estimates. [START_REF] Mohr | Modeling data as a fractional Gaussian noise[END_REF], [START_REF] Graf | Long-range correlations and estimation of the self-similarity parameter[END_REF], [START_REF] Geweke | The estimation and application of long memory time series models[END_REF], [START_REF] Fox | Large-sample properties of parameter estimates for strongly dependent stationary Gaussian time series[END_REF] and [START_REF] Dahlhaus | Efficient parameter estimation for self-similar processes[END_REF] have considered the problem of periodogram estimation. [START_REF] Mcleod | Preservation of the rescaled adjusted range: 1. a reassessment of the hurst phenomenon[END_REF] have discussed computational considerations involved in the application of maximum likelihood estimation. [START_REF] Sweeting | Uniform asymptotic normality of the maximum likelihood estimator[END_REF] has proposed assumptions to obtain the consistency and the asymptotic normality of maximum likelihood estimator. [START_REF] Dahlhaus | Efficient parameter estimation for self-similar processes[END_REF] has established the asymptotic properties of the exact maximum likelihood estimation of the parameters of strongly dependent Gaussian processes.

In the particular case of autoregressive models endowed with a strongly dependent noise, [START_REF] Brouste | Asymptotic properties of the MLE for the autoregressive process coefficients under stationary gaussian noise[END_REF] have proposed a new approach based on a transformation of the initial model to establish the asymptotic properties of the maximum likelihood estimator. See also [START_REF] Yajima | On estimation of a regression model with long-memory stationary errors[END_REF] for the case of a regression problem with nonstochastic regressors.

We study in this paper the asymptotic properties of the generalized least squares estimator of autoregressive models endowed with fractional Gaussian noise. We propose a simple and natural transformation of the initial model and we establish the convergence and asymptotic normality of the generalized least squares estimator based on this transformation. These properties are obtained using only the standard hypothesis (A0). We propose simple proofs of the asymptotic properties of the generalized least squares estimator of the autoregressive parameters. We use the closed form of the estimator to provide very intuitive proofs. We also extend the support of the memory parameter by considering that 0 < H < 1. Contrary therefore to [START_REF] Dahlhaus | Efficient parameter estimation for self-similar processes[END_REF] and [START_REF] Brouste | Asymptotic properties of the MLE for the autoregressive process coefficients under stationary gaussian noise[END_REF], we take into consideration not only long memory dynamics but also short memory fluctuations. We also show numerically that our estimator is preferable (in terms of quadratic risk) to the Whittle estimator proposed by [START_REF] Dahlhaus | Efficient parameter estimation for self-similar processes[END_REF] when the sample size is relatively moderate.

The paper is organized as follows. In Section 2, we introduce our model and we study the asymptotic properties of our estimator. In Section 3, we show that our results obtained in the framework of fractional autoregressive models of order 1 can easily be extended to the case of autoregressive models of order p ≥ 1 induced by fractional Gaussian noise. We present in Section 4 some simulation studies to illustrate our theoretical results. We also offer comparisons in terms of the quadratic risk of our estimator with that of Whittle (see for example [START_REF] Dahlhaus | Efficient parameter estimation for self-similar processes[END_REF]). Proofs of our two main results are set out to Section 5. All figures are collected in the last section of the paper.

In all this work, x T designates the transpose of the vector x. We shall also use the matrix norm defined by

A = sup x R k 2 ≤1 Ax R k 1 = ρ 1/2 (A T A), when A is a R k 1 ×k 2 matrix, x 2 R k 2 = x T
x is the Euclidean norm of the vector x ∈ R k 2 , and ρ(•) denotes the spectral radius.

Generalized least squares estimation of FrAR models of order 1

In this section we present the parametrization that is used in the sequel and we state the asymptotic properties of the GLSE of FrAR(1) when the Hurst exponent H is known. An application of these asymptotic properties can be found in Ben [START_REF] Ben Hariz | Fast and Asymptotically-efficient estimation in a Fractional autoregressive process[END_REF] where a two steps procedure is performed when H is unknown.

Statement of the problem and notations

First, note that the ordinary least squares estimator of the parameter a 0 in (1) is biased. To see that, let â(OLS) n be the standard least squares estimator of a 0 . It is well known that

â(OLS) n = n t=1 X t X t-1 n t=1 X 2 t-1
.

Thanks to (1), one has

â(OLS) n -a 0 = n t=1 H t X t-1 n t=1 X 2 t-1
.

The ergodicity of ( H t ) t∈Z and the stationarity of (X t ) t∈Z imply that

â(OLS) n -a 0 a.s. ---→ n→∞ E H 2 X 1 Var (X 1 )
.

From ( 2) and ( 16), one can easily obtain that

E H 2 X 1 = a 0 -2 2 + (a 0 -1) 2 2 j≥2 a j-2 0 j 2H .
Since the function H ∈ (0, 1) -→ j≥2 a j-2 0 j 2H is bijective, we deduce that E[ H 2 X 1 ] = 0 if and only if H = 1/2. This leads to the standard case of an autoregressive model where the noise is an iid sequence with standard normal distribution.

Due to the bias problem of the ordinary least squares estimator (problem caused by the correlation of the noise), we consider here the generalized least squares estimation procedure. This technique consists in transforming Model (1) to another autoregressive model, with the same unknown autoregressive parameter a 0 , endowed with a standard white noise and consider then the ordinary least squares estimation on the new model.

From (1), one can write

   X 2 . . . X n+1    = a 0    X 1 . . . X n    +    H 2 . . . H n+1    . (3)
The vector ( H 2 , . . . , H n+1 ) T is a centered normal random vector with covariance matrix

Ω n,H := [γ H (j -i)] 1≤i,j≤n = [Cov( H i , H j )]
= P n,H D 1/2 n,H P T n,H where D 1/2 n,H = diag((λ (n) 1,H ) 1/2 , . . . , (λ (n) n,H ) 1/2 ) and the λ (n) i,H are the eigenvalues of Ω n,H . The matrix Ω 1/2 n,H is the square root of Ω n,H since it satisfies Ω 1/2 n,H Ω 1/2 n,H = Ω n,H . In view of this last remark, Ω -1/2 n,H ( H 2 , . . . , H n+1 )
T is a standard normal random vector. Thus, a rather natural and simple transformation that can be proposed is the following:

Y n,H = a 0 Z n,H + U n,H , (4)
where

Y n,H = Ω -1/2 n,H    X 2 . . . X n+1    , Z n,H = Ω -1/2 n,H    X 1 . . . X n    and U n,H = Ω -1/2 n,H    H 2 . . . H n+1    . (5)
The random variable ân is called generalized least squares estimator if it satisfies, almost surely, ân = argmin a∈(-1,1)

1 n U n,H (a) 2 R n , (6)
where

U n,H (a) = Ω -1/2 n,H    e H 2 (a) . . . e H n+1 (a)    (7)
and

e H t (a) = X t -aX t-1 (8) with e H t (a) = X t = 0 if t ≤ 0. Note that U n,H (a 0 ) = U n,H .
The resolution of the optimization problem presented in (6) leads to an explicit expression of ân . One can easily prove that

ân = Z T n,H Y n,H Z n,H 2 R n . (9)
Moreover, in view of this last equality and Equation ( 4), one can deduce that

ân -a 0 = Z T n,H U n,H (a 0 ) Z n,H 2 R n . ( 10
)

Asymptotic properties

The asymptotic properties of the generalized least squares estimator of the fractional autoregressive model of order 1 are stated in the following two theorems.

Theorem 1. (Consistency). We assume that (X t ) t∈Z satisfies (1). Under (A0), we have

ân P ---→ n→∞ a 0 ,
where (â n ) n≥1 is a sequence of generalized least squares estimators as defined in (9).

The proof of this theorem is given in Subsection 5.2.

Theorem 2. (Asymptotic normality). We assume that (X t ) t∈Z satisfies (1). Under (A0), the sequence ( √ n(â n -a 0 )) n≥1 has a limiting centered normal distribution with variance 1 -a 2 0 . The proof of this theorem is set out to Subsection 5.3.

Generalized least squares estimation of FrAR models of order p

The results obtained in the previous section can be extended to fractional autoregressive models of order p. To be more precise, suppose that the process (X t ) t∈Z admits an autoregressive representation of the form

X t = a 01 X t-1 + • • • + a 0p X t-p + H t , (11)
where ( H t ) t∈Z is a fractional Gaussian noise with 0 < H < 1. To ensure the stability of the process (X t ) t∈Z in this case, we suppose that its reverse characteristic polynomial has no root in and on the complex unit circle, i.e. we assume that

1 -a 01 z -• • • -a 0p z p = 0 for |z| ≤ 1. ( 12
)
In view of (11), it can readily be shown that

   X 2 . . . X n+1    =    X 1 • • • X 2-p . . . . . . X n • • • X n+1-p       a 01 . . . a 0p    +    H 2 . . . H n+1    . ( 13
)
We readopt the same transformation as the one used in the previous section to obtain a vector representation similar to the one introduced in (4). Multiplying all the terms of Equation ( 13) by the appropriate matrix

Ω -1/2
n,H , one has

Y n,H = Z n,H A 0 + U n,H (A 0 ),
where Y n,H and U n,H (•) are defined in ( 5) and ( 7) and where

Z n,H = Ω -1/2 n,H    X 1 • • • X 2-p . . . . . . X n • • • X n+1-p    and A 0 =    a 01 . . . a 0p    .
Similarly as before, we again show that the generalized least squares estimator A n of A 0 takes a closed form and it is explicitly given by:

A n = argmin A 1 n U n,H (A) 2 R n = Z T n,H Z n,H -1 Z T n,H Y n,H = A 0 + Z T n,H Z n,H -1 Z T n,H U n,H (A 0 ) ,
where the minimum is taken over all vectors A = (a

1 , • • • , a p ) T of R p-1 × R * such that (12) is satisfied.
The following theorem gives the asymptotic properties of the generalized least squares estimator of autoregressive models of order p induced by fractional Gaussian noise. The proof of this theorem is similar to the one of Theorems 1 and 2 detailed in Section 5.

Theorem 3. We assume that (X t ) t∈Z satisfies (11). Suppose that 0 < H < 1 and that

A 0 = (a 01 , • • • , a 0p ) T is such that Condition (12) is verified. Then, the sequence ( √ n( A n -A 0 )) n≥1 has a limiting centered normal distribution with variance Γ -1 (A 0 ), where Γ (A 0 ) is the limit in probability of Z T n,H Z n,H /n.

Numerical illustrations

In this section, we investigate the finite sample properties of the asymptotic results that we introduced in this work. For that, we use Monte Carlo experiments. The numerical illustrations of this section are made with the open source statistical software R (see R Development Core Team, 2017) or (see http://cran.r-project.org/).

Simulation studies

We first study numerically the behavior of the generalized least squares estimator for fractional autoregressive models of the form

X t = aX t-1 + H t , (14)
where the unknown autoregressive parameter is taken as a = 0.3, and where ( H t ) t∈Z is a fractional Gaussian noise with Hurst index H = 0.7.

We simulate N = 2, 000 independent trajectories of sizes n = 30, n = 500 and n = 2, 000 of Model ( 14).

Figures 1, 2 and 3 illustrate the distribution of the generalized least squares estimator ân of the autoregressive parameter a. We can notice that these numerical results are consistent with the convergence established in Theorem 2.

In Figures 4, 5 and6, we corroborate the convergence of the estimator 1 -â2 n of the asymptotic variance of the generalized least squares estimator ân given in Theorem 2. We can remark that even for a small sample size the estimator is very close to the true value of the parameter. This suggests that our estimator is not greedy in terms of sample size. This phenomenon is confirmed by Figures 89101112. For different values of n and for this same set of parameters (i.e. a = 0.3 and H = 0.7), we compare in Figure 7 the quadratic risk of the generalized least squares estimator ân (in black) with that of the Whittle estimator considered for example by [START_REF] Dahlhaus | Efficient parameter estimation for self-similar processes[END_REF] (in green) based on N = 2, 000 independent replications for each value of n. It is clear in this case that when the sample size is relatively small our estimator makes less error than the Whittle estimator. It is also worth noting that this closed form of our estimator makes it preferable even when the sample size is large.

To have a more global view on the behavior of our estimator when the sample size is relatively small, we consider the cases where a = -0.9, -0.5, 0, 0.5, 0.9 in Model (14). For each value of a, we vary H over the interval [0.01, 0.99] with a step of 0.01. For each value of the autoregressive parameter a and the memory parameter H, we simulate N = 1, 000 independent trajectories of (14). We compare in Figures 89101112, the logarithm of the quadratic risk of the generalized least squares estimator ân (in black) with that of the Whittle estimator (in green) of N = 1, 000 replications of Model ( 14) with different values of a and H. In each figure, four sample sizes were considered. From left to right and top to bottom, the size takes the values n = 50, 100, 150, 200. Clearly, for the different values of a and n, our estimator shows a certain stability unlike that of Whittle where the quadratic risk is characterized by a certain instability especially when H is close to 0 and 1. This observation highlights the fact that the generalized least squares estimator ân gives very satisfactory results in the different memory scenarios presented by the process. Note that for all the illustrations given in Figures 89101112, H sweeps the entire memory interval (0, 1). The Whittle estimator seems to be less precise for several values of a when the process is short memory (i.e. when H ∈ (0, 1/2]).

Proofs

In all our proofs, K is a positive constant that may vary from line to line.

Preliminary results

In this subsection, we state some results that will be useful for the proofs of the asymptotic properties of the generalized least squares estimator introduced in Section 2. We use (A0) to write Z n,H as a function of U n,H (a 0 ). This expression will provide very intuitive proofs of the asymptotic properties of the GLSE. We also give the asymptotic behavior of the components of the inverse variance matrix of the fractional Gaussian noise ( H t ) t∈Z based on the behavior in a neighbourhood of zero of the spectral density of this process.

5.1.1. The expression of Z n,H as a function of U n,H (a 0 )
We recall that, when |a| < 1, the polynomial a(z) = 1 -az is invertible and we have

a -1 (z) = j≥0 a j z j . ( 15
)
This implies that under (A0) one may write the process (X t ) t∈Z as a linear combination of the infinite past ( H t-j ) j≥0 . More explicitly, we have

X t = (1 -a 0 L) -1 H t = j≥0 a j 0 H t-j , (16)
where L stands for the back-shift operator.

In view of the definition of the process (e H t (a)) t∈Z introduced in ( 8) and ( 15), one can express X t when 1 ≤ t ≤ n + 1 as a function of the finite past (e H t-j (a 0 )) 0≤j≤t-1 . Indeed, Equation (8) and Assumption (A0) imply that

X t = (1 -a 0 L) -1 e H t (a 0 ) = t-1 j=0 a j 0 e H t-j (a 0 ), ∀t ∈ {1, . . . , n + 1} . (17)
From ( 17), we deduce that the vector (X 1 , . . . , X n ) T can be rewritten in the form

     X 1 X 2 . . . X n      =       1-1 j=0 a j 0 e H 1-j (a 0 ) 2-1 j=0 a j 0 e H 2-j (a 0 ) . . . n-1 j=0 a j 0 e H n-j (a 0 )       =      e H 1 (a 0 ) e H 2 (a 0 ) + a 0 e H 1 (a 0 ) . . . e H n (a 0 ) + • • • + a n-1 0 e H 1 (a 0 )      =      e H 1 (a 0 ) e H 2 (a 0 ) . . . e H n (a 0 )      + a 0      e H 0 (a 0 ) e H 1 (a 0 ) . . . e H n-1 (a 0 )      + • • • + a n-1 0      e H 2-n (a 0 ) e H 3-n (a 0 ) . . . e H 1 (a 0 )      = n-1 j=0 a j 0      e H 1-j (a 0 ) e H 2-j (a 0 ) . . . e H n-j (a 0 )      = n-1 j=0 a j 0      L j+1      e H 2 (a 0 ) e H 3 (a 0 ) . . . e H n+1 (a 0 )           , ( 18
)
where the shift operator L j+1 acts on all the components of the vector (e H 2 (a 0 ), . . . , e H n+1 (a 0 )) T , i.e.

L j+1    e H 2 (a 0 ) . . . e H n+1 (a 0 )    =    e H 1-j (a 0 ) . . . e H n-j (a 0 )    .
In view of ( 18), it follows that the vector Z n,H defined in ( 7) can be expressed as a function of U n,H (a 0 ):

Z n,H = Ω -1/2 n,H      X 1 X 2 . . . X n      = n-1 j=0 a j 0      L j+1 Ω -1/2 n,H      e H 2 (a 0 ) e H 3 (a 0 ) . . . e H n+1 (a 0 )           = n j=1 a j-1 0 L j U n,H (a 0 ) . (19)
5.1.2. Asymptotic properties of the components of the inverse covariance matrix of ( H t ) 1≤t≤n

By [START_REF] Fox | Large-sample properties of parameter estimates for strongly dependent stationary Gaussian time series[END_REF], [START_REF] Sinai | Self-similar probability distributions[END_REF], the spectral density f H (•) of the fractional Gaussian noise

( H t ) t∈Z is defined for all ω ∈ [-π, π] by: f H (ω) = 1 2π h∈Z γ H (h)e -iωh = C H (1 -cos(ω)) j∈Z 1 |ω + 2jπ| 1+2H = C H (1 -cos(ω)) |ω| -1-2H + C H (1 -cos(ω)) j∈Z j =0 1 |ω + 2jπ| 1+2H ,
where C H = 2Γ (2H + 1) sin(πH) and Γ (•) denotes the Gamma function.

As ω → 0, we have

f H (ω) ∼ C H 2 |ω| 1-2H . ( 20
)
Proposition 1. Under Assumption (A0), the following three assertions hold true.

1. The covariance matrix Ω n,H of the centered Gaussian random vector ( H 1 , . . . , H n ) T is a symmetric Toeplitz positive definite matrix. 2. For large k, the behavior of the autocovariance γ H (k) is given by:

(Ω n,H ) 1,k = γ H (k) = H(2H -1)k 2H-2 + o(k 2H-2 ). ( 21
)
3. The inverse autocovariances of ( H t ) 1≤t≤n satisfy

Ω -1 n,H j,k ∼ 1 k -j 2H Γ (2H)e iπH C H π 2 when |k -j| → ∞. ( 22
)
Proof. For the first point, observe that Ω n,H can be considered as a self-adjoint endomorphism of the vector space C n . Note now that for any

C ∈ C n \ {0 C n }, the spectral representation of Ω n,H implies that C T Ω n,H C = n j,k=1 C k (Ω n,H ) k,j C j = n j,k=1 C k π -π f H (ω)e i(j-k)ω dω C j = π -π n j=1 C j e ijω 2 f H (ω)dω.
The function ω → n j=1 C j e ijω is an analytic function on the disc and is not identically zero. Since f H (ω) > 0 for ω = 0 the last integral is positive. Thus, we conclude that the symmetric Toeplitz matrix Ω n,H is positive definite.

In view of (2), a Taylor expansion of the function

: x → (1 -x) 2H -2 + (1 + x) 2H around 0 implies that for large k, (Ω n,H ) 1,k = γ H (k) = 1 2 k 2H (1/k) = H(2H -1)k 2H-2 + o(k 2H-2 ).
This gives the second point of the proposition.

The components of the inverse matrix of Ω n,H can also be expressed as a function of the spectral density of ( H t ) t∈Z (see [START_REF] Fox | Large-sample properties of parameter estimates for strongly dependent stationary Gaussian time series[END_REF]). More precisely, we have for all j, k = 1, . . . , n,

Ω -1 n,H j,k = 1 (2π) 2 π -π 1 f H (ω) e i(k-j)ω dω. ( 23
)
Using the parity of the function f H (•), observe that when j = k we have

Ω -1 n,H j,j = 1 (2π) 2 π -π 1 f H (ω) dω = 1 2π 2 π 0 1 f H (ω) dω. ( 24
)
From equivalence in (20), one has, when ω → 0,

1 f H (ω) = 2 C H |ω| 2H-1 + o 2 C H |ω| 2H-1 .
This implies that for κ > 0 there exists δ κ > 0 such that for any ω ∈] -δ κ , δ κ [ we have

(1 -κ) 2 C H |ω| 2H-1 ≤ 1 f H (ω) ≤ (1 + κ) 2 C H |ω| 2H-1 . ( 25
)
Thanks to ( 24), ( 25) and the boundedness of the function ω → 1/f H (ω) on ]δ κ , π], one can show that there exists a positive constant K such that for any j = 1, . . . , n,

Ω -1 n,H j,j ≤ 1 + κ C H π 2 δκ 0 ω 2H-1 dω + 1 2π 2 π δκ 1 f H (ω) dω ≤ δ 2H κ (1 + κ) 2HC H π 2 + π -δ κ 2π 2 sup ω∈]δκ,π] 1 f H (ω) ≤ K . ( 26
)
When j = k (suppose without loss of generality that k > j), one can use Lemma 6.3.2 (page 442) of [START_REF] Ablowitz | Complex Variables: Introduction and Applications[END_REF] which is an analog version of Watson's Lemma (see for instance [START_REF] Watson | The harmonic functions associated with the parabolic cylinder[END_REF] for the original proof or [START_REF] Miller | Applied asymptotic analysis[END_REF] for a more recent development) to obtain the asymptotic equivalent of the (j, k)-th component of Ω -1 n,H when |k -j| → ∞ stated in the third point of the proposition.

Remark 2. From (21), for any ρ > 0 there exists N ρ ∈ N such that for any k ≥ N ρ we have

H {(2H -1) -ρ |2H -1|} k 2H-2 ≤ (Ω n,H ) 1,k ≤ H {ρ |2H -1| + (2H -1)} k 2H-2 . Since -1 < 2H -1 < 1 we obtain that (Ω n,H ) 1,k ≤ H {ρ |2H -1| + 1} k 2H-2 . ( 27
)
Observe that for k = 1, . . . , N ρ , we have in view of (2) that the autocovariances (Ω n,H ) 1,k are bounded, i.e. there exists K H,ρ > 0 such that

(Ω n,H ) 1,k ≤ K H,ρ . Since k 2H-2 ≤ 1 for k = 1, . . . , N ρ , it follows that (Ω n,H ) 1,k ≤ K H,ρ k 2H-2 , k = 1, . . . , N ρ . ( 28
)
Consequently, we obtain from ( 27) and (28) that for any k ≥ 1,

(Ω n,H ) 1,k ≤ K k 2H-2 , ( 29
)
where

K = max(H {ρ |2H -1| + 1} , K H,ρ ).
A similar calculation can be done to show that there exists a positive constant K such that for any j, k = 1, 2, . . . , we have

Ω -1 n,H j,k ≤ K 1 k -j 2H .
(30)

Proof of Theorem 1

The proof of convergence in probability of (â n ) n≥1 to a 0 will be done in two steps. We first show that ( Z n,H 2 R n /n) n≥1 converges in mean square to 1/(1 -a 2 0 ) and we prove in a second time that (Z T n,H U n,H (a 0 )/n) n≥1 converges in mean square to 0.

5.2.1. Proof of the convergence in mean square of ( Z n,H 2 R n /n) n≥1 to 1/(1 -a 2 0 )
To establish the mean square convergence of (

Z n,H 2 R n /n) n≥1 to 1/(1 -a 2 0 )
, we consider the following lemmas and remark.

Lemma 1. Under Assumption (A0), we have

E 1 n Z n,H 2 R n ---→ n→∞ 1 1 -a 2 0 .
Proof. In view of ( 19), one has

E 1 n Z n,H 2 R n = 1 n n j 1 =1 n j 2 =1 a j 1 +j 2 -2 0 E L j 1 U T n,H (a 0 ) L j 2 U n,H (a 0 ) = 1 n n j 1 =1 n j 2 =1 a j 1 +j 2 -2 0 n k=1 E L j 1 (U n,H (a 0 )) k L j 2 (U n,H (a 0 )) k = 1 n n j 1 =1 n j 2 =1 a j 1 +j 2 -2 0 n k=1 n r 1 =1 n r 2 =1 Ω -1/2 n,H k,r 1 Ω -1/2 n,H k,r 2 E e H r 1 +1-j 1 (a 0 )e H r 2 +1-j 2 (a 0 ) = 1 n n j 1 =1 n j 2 =1 a j 1 +j 2 -2 0 n r 1 =j 1 n r 2 =j 2 Ω -1 n,H r 1 ,r 2 (Ω n,H ) r 1 +1-j 1 ,r 2 +1-j 2 .
The stationarity of ( H t ) t∈Z implies that

• for j 2 < j 1 : (Ω n,H ) r 1 +1-j 1 ,r 2 +1-j 2 = (Ω n,H ) r 2 ,r 1 +j 2 -j 1 ,
• for j 2 ≥ j 1 :

(Ω n,H ) r 1 +1-j 1 ,r 2 +1-j 2 = (Ω n,H ) r 1 ,r 2 +j 1 -j 2 .
Thus, we obtain

E 1 n Z n,H 2 R n = 1 n n j 1 =1 j 1 -1 j 2 =1 a j 1 +j 2 -2 0 n r 1 =j 1 n r 2 =j 2 Ω -1 n,H r 1 ,r 2 (Ω n,H ) r 2 ,r 1 +j 2 -j 1 + 1 n n j 1 =1 n j 2 =j 1 a j 1 +j 2 -2 0 n r 1 =j 1 n r 2 =j 2 Ω -1 n,H r 1 ,r 2 (Ω n,H ) r 1 ,r 2 +j 1 -j 2 .
This implies that

E 1 n Z n,H 2 R n = 1 n n j 1 =1 j 1 -1 j 2 =1 a j 1 +j 2 -2 0 n r 1 =j 1 (I n ) r 1 ,r 1 +j 2 -j 1 - 1 n n j 1 =1 j 1 -1 j 2 =1 a j 1 +j 2 -2 0 n r 1 =j 1 j 2 -1 r 2 =1 Ω -1 n,H r 1 ,r 2 (Ω n,H ) r 2 ,r 1 +j 2 -j 1 + 1 n n j 1 =1 n j 2 =j 1 a j 1 +j 2 -2 0 n r 2 =j 2 (I n ) r 2 ,r 2 +j 1 -j 2 - 1 n n j 1 =1 n j 2 =j 1 a j 1 +j 2 -2 0 j 1 -1 r 1 =1 n r 2 =j 2 Ω -1 n,H r 1 ,r 2 (Ω n,H ) r 1 ,r 2 +j 1 -j 2 .
We have

1 n n j 1 =1 j 1 -1 j 2 =1 a j 1 +j 2 -2 0 n r 1 =j 1 (I n ) r 1 ,r 1 +j 2 -j 1 = 0 and 1 n n j 1 =1 n j 2 =j 1 a j 1 +j 2 -2 0 n r 2 =j 2 (I n ) r 2 ,r 2 +j 1 -j 2 = 1 n n j=1 a 2j-2 0 (n -j + 1),
thus we may write

E 1 n Z n,H 2 R n = 1 n n j=1 a 2j-2 0 (n -j + 1) -W n,H , (31)
where

W n,H = 1 n n j 1 =1 j 1 -1 j 2 =1 a j 1 +j 2 -2 0 n r 1 =j 1 j 2 -1 r 2 =1 Ω -1 n,H r 1 ,r 2 (Ω n,H ) r 2 ,r 1 +j 2 -j 1 + 1 n n j 1 =1 n j 2 =j 1 a j 1 +j 2 -2 0 j 1 -1 r 1 =1 n r 2 =j 2 Ω -1 n,H r 1 ,r 2 (Ω n,H ) r 1 ,r 2 +j 1 -j 2 . ( 32
)
We shall prove that W n,H ---→ n→∞ 0. Consider the first term on the right side of Equation (32) (the second term can be treated in a similar way) and note that r 2 -r 1 ≤ j 1 -r 1 -1 and r 2 -r 1 + j 1 -j 2 ≤ j 1 -r 1 -1. By the properties of the components of the matrix Ω n,H and its inverse given in ( 29) and ( 30), there exists a positive constant K such that

1 n n j 1 =1 j 1 -1 j 2 =1 a j 1 +j 2 -2 0 n r 1 =j 1 j 2 -1 r 2 =1 Ω -1 n,H r 1 ,r 2 (Ω n,H ) r 2 ,r 1 +j 2 -j 1 ≤ K n n j 1 =1 j 1 -1 j 2 =1 |a 0 | j 1 +j 2 -2 n r 1 =j 1 j 2 -1 r 2 =1 1 (r 2 -r 1 ) 2H 1 (r 2 -r 1 + j 1 -j 2 ) 2-2H ≤ K n n j 1 =1 j 1 -1 j 2 =1 |a 0 | j 1 +j 2 -2 n r 1 =j 1 j 2 -1 r 2 =1 1 (j 1 -r 1 -1) 2H 1 (j 1 -r 1 -1) 2-2H = K n n j 1 =1 n r 1 =j 1 |a 0 | j 1 -1 (j 1 -r 1 -1) 2 j 1 -1 j 2 =1 (j 2 -1) |a 0 | j 2 -1 . Since j 1 -1 j 2 =1 (j 2 -1) |a 0 | j 2 -1 = 1 -|a 0 | j 1 -1 {(j 1 -1)(1 -|a 0 |) + 1} (1 -|a 0 |) 2 ≤ 1 (1 -|a 0 |) 2 ,
we deduce that

1 n n j 1 =1 j 1 -1 j 2 =1 a j 1 +j 2 -2 0 n r 1 =j 1 j 2 -1 r 2 =1 Ω -1 n,H r 1 ,r 2 (Ω n,H ) r 2 ,r 1 +j 2 -j 1 ≤ K n(1 -|a 0 |) 2 n j 1 =1 n r 1 =j 1 |a 0 | j 1 -1 (j 1 -r 1 -1) 2 .

Now we consider the change of indices r

1 = j 1 + k -1 to obtain that 1 n n j 1 =1 j 1 -1 j 2 =1 a j 1 +j 2 -2 0 n r 1 =j 1 j 2 -1 r 2 =1 Ω -1 n,H r 1 ,r 2 (Ω n,H ) r 2 ,r 1 +j 2 -j 1 ≤ K n(1 -|a 0 |) 2 n j 1 =1 |a 0 | j 1 -1 k≥1 1 k 2 ---→ n→∞ 0.
The second term on the right side of Equation ( 32) can be treated in the same way. The proof of the convergence of the sequence (W n,H ) n≥1 towards 0 is then completed. Finally, from (31) we conclude that

E 1 n Z n,H 2 R n = n + 1 n 1 -a 2n 0 1 -a 2 0 - 1 n 1 -a 2n 0 1 -a 2 0 - 1 n 1 -a 2n-2 0 n(1 -a 2 0 ) + a 2 0 (1 -a 2 0 ) 2 -W n,H ---→ n→∞ 1 1 -a 2 0 ,
and the lemma is proved.

Remark 3. Denote by Σ n,X ,H the variance matrix of the random vector (X 1 , . . . , X n ) T , i.e. for any 1 ≤ i, j ≤ n, (Σ n,X ,H ) i,j = Cov(X i , X j ), and let Z

(1) n,H be the random vector defined by Z

(1)

n,H = Σ -1/2 n,X ,H (X 1 , . . . , X n ) T . From (7), one has 1 n Z n,H 2 R n = 1 n Z T n,H Z n,H = 1 n    X 1 . . . X n    T Ω -1 n,H    X 1 . . . X n    = 1 n (Z (1) n,H ) T Σ 1/2 n,X ,H Ω -1 n,H Σ 1/2 n,X ,H Z (1)
n,H .

The covariance matrix Σ 1/2 n,X ,H Ω -1 n,H Σ 1/2
n,X ,H is a symmetric real matrix. So there exists an orthogonal matrix

Q n,H such that Λ n,H = Q T n,H Σ 1/2 n,X ,H Ω -1 n,H Σ 1/2 n,X ,H Q n,H is a diagonal matrix. Hence, the random variable Z n,H 2 R n /n can be rewritten as 1 n Z n,H 2 R n = 1 n (Q T n,H Z (1) n,H ) T Λ n,H Q T n,H Z (1) n,H = 1 n n j=1 η (n) j,H Q T n,H Z (1) n,H j 2 , ( 33
)
where η

(n) 1,H , . . . , η (n)
n,H are the eigenvalues of the symmetric real matrix

Σ 1/2 n,X ,H Ω -1 n,H Σ 1/2 n,X ,H . The expression of Z n,H 2 R n /n given in (33) has the advantage of involving the standard normal random vector Q T n,H Z (1)
n,H . This will considerably simplify the calculations in the following.

We show in the lemma below the boundedness of the eigenvalues η

(n) 1,H , . . . , η (n)
n,H of the matrix

Σ 1/2 n,X ,H Ω -1 n,H Σ 1/2
n,X ,H uniformly in n. Lemma 2. Under Assumption (A0) and for any 1 ≤ j ≤ n, we have

1 (1 + |a 0 |) 2 ≤ η (n) j,H ≤ 1 (1 -|a 0 |) 2 .
Proof. Let f X (•) be the spectral density of (X t ) t∈Z and denote by f H (•) the spectral density of the stationary process ( H t ) t∈Z . One can express f X (•) as a function of f H (•). Indeed, using ( 16) and the stationarity of (X t ) t∈Z and ( H t ) t∈Z we have for any ω ∈ R,

f X (ω) = 1 2π h∈Z γ X (h)e -iωh = 1 2π h∈Z Cov(X t , X t-h )e -iωh = 1 2π h∈Z Cov   j∈N a j 0 H t-j , k∈N a k 0 H t-h-k   e -iωh = 1 2π h∈Z j,k∈N a j+k 0 γ H (h + k -j)e -iω(h+k-j) e -iω(j-k) = f H (ω)   j∈N a j 0 e -iωj     j∈N a j 0 e -iωj   = f H (ω) j∈N a j 0 e -iωj 2 ,
where γ X (•) denotes the autocovariance function of the process (X t ) t∈Z and z designates the complex conjugate of the complex number z. Furthermore, Assumption (A0) implies that the modulus of the complex number a 0 e -iω is strictly less than 1. Hence, we deduce that the spectral density of the stationary process (X t ) t∈Z is related to f H (•) as follows:

f X (ω) = f H (ω) 1 1 -a 0 e -iω 2 = f H (ω) 1 -2a 0 cos(ω) + a 2 0 , ∀ω ∈ R. ( 34
)
From this last equation and taking into consideration the fact that

(1-|a 0 |) 2 ≤ 1-2a 0 cos(ω)+a 2 0 ≤ (1 + |a 0 |) 2 , one obtains that for any ω ∈ R, (1 -|a 0 |) 2 f X (ω) ≤ f H (ω) ≤ (1 + |a 0 |) 2 f X (ω). ( 35
)
Note now that for any vector V of R n×1 , we have

V T Ω n,H -(1 + |a 0 |) 2 Σ n,X ,H V = n k=1 n j=1 V j Ω n,H -(1 + |a 0 |) 2 Σ n,X ,H j,k V k = n k,j=1 V j π -π f H (ω) -(1 + |a 0 |) 2 f X (ω) e iω(k-j) dωV k = π -π f H (ω) -(1 + |a 0 |) 2 f X (ω) V T Υ n (ω)V dω, (36)
where

Υ n (ω) = V n (ω)V * n (ω) with V n (ω)
= (e iω , e 2iω , . . . , e niω ) T . The notation V * designates the conjugate transpose of the vector V .

Furthermore, the real V T Υ n (ω)V is a non-negative number. In fact,

V T Υ n (ω)V = V T V n (ω)V * n (ω)V = (V * n (ω)V ) * (V * n (ω)V ) = |V * n (ω)V | 2 .
Then, from ( 35) and ( 36), one has

V T (Ω n,H -(1 + |a 0 |) 2 Σ n,X ,H )V ≤ 0. (37)
Now, let ζ be an element of the spectrum of

Ω -1/2 n,H Σ n,X ,H Ω -1/2 n,H . Then, there exists C 1 ∈ R n×1 \ {0 R n×1 } such that Ω -1/2 n,H Σ n,X ,H Ω -1/2 n,H C 1 = ζC 1 . Hence, it follows that C T 1 Ω -1/2 n,H Σ n,X ,H Ω -1/2 n,H C 1 = ζ C 1 2 R n . Taking C 2 = Ω -1/2
n,H C 1 , we obtain from this last equation that

C T 2 Σ n,X ,H C 2 = ζ Ω 1/2 n,H C 2 2 R n . Using (37), we deduce that C T 2 Ω n,H C 2 ≤ (1 + |a 0 |) 2 ζ Ω 1/2 n,H C 2 2 R n . Since C 1 = Ω 1/2
n,H C 2 = 0 R n×1 as an eigenvector, we conclude that

ζ ≥ C T 2 Ω n,H C 2 Ω 1/2 n,H C 2 2 R n 1 (1 + |a 0 |) 2 = 1 (1 + |a 0 |) 2 . ( 38
)
A similar calculation can be done to prove that, for any vector

V of R n×1 , V T (Ω n,H -(1 -|a 0 |) 2 Σ n,X ,H )V ≥ 0.
Arguing as before, one obtains that

ζ ≤ C T 2 Ω n,H C 2 Ω 1/2 n,H C 2 2 R n 1 (1 -|a 0 |) 2 = 1 (1 -|a 0 |) 2 . Since for any α ∈ R, det Σ 1/2 n,X ,H Ω -1 n,H Σ 1/2 n,X ,H -αI n = det Σ 1/2 n,X ,H Ω -1 n,H Σ n,X ,H -αI n Σ -1/2 n,X ,H = det Ω -1 n,H Σ n,X ,H -αI n = det Ω -1/2 n,H Ω -1/2 n,H Σ n,X ,H Ω -1/2 n,H -αI n Ω 1/2 n,H = det Ω -1/2 n,H Σ n,X ,H Ω -1/2 n,H -αI n , the spectrum of the matrix Σ 1/2 n,X ,H Ω -1 n,H Σ 1/2 n,X ,H is equal to the set of the eigenvalues of Ω -1/2 n,H Σ n,X ,H Ω -1/2
n,H , the lemma is then proved.

Lemma 3. Under Assumption (A0), we have

1 n Z n,H 2 R n L 2 ---→ n→∞ 1 1 -a 2 0 .
Proof. In view of Lemma 2, Remark 3 and (33), we have

Var 1 n Z n,H 2 R n = 1 n 2 n j=1 η (n) j,H 2 Var Q T n,H Z (1) n,H j 2 = 2 n 2 n j=1 η (n) j,H 2 ≤ 2 n(1 -|a 0 |) 4 ---→ n→∞ 0. Now, observe that E 1 n Z n,H 2 R n - 1 1 -a 2 0 2 = Var 1 n Z n,H 2 R n + E 1 n Z n,H 2 R n - 1 1 -a 2 0 2 .
Thanks to Lemma 1 and the convergence stated just before, the expected result is obtained.

Proof of the convergence in mean square of (Z T

n,H U n,H (a 0 )/n) n≥1 to 0 Lemma 4. Under Assumption (A0), we have

1 n Z T n,H U n,H (a 0 ) L 2 ---→ n→∞ 0.
Proof. Due to the fact that e H t (a 0 ) is centered for any t, we notice that for fixed (i, j, k) the joint cumulant of e H 0 (a 0 ), e H i (a 0 ), e H j (a 0 ) and e H k (a 0 ) is given by:

cum(e H 0 (a 0 ), e H i (a 0 ), e H j (a 0 ), e H k (a 0 )) = E e H 0 (a 0 )e H i (a 0 )e H j (a 0 )e H k (a 0 ) -E e H 0 (a 0 )e H i (a 0 ) E e H j (a 0 )e H k (a 0 ) -E e H 0 (a 0 )e H j (a 0 ) E e H i (a 0 )e H k (a 0 ) -E e H 0 (a 0 )e H k (a 0 ) E e H i (a 0 )e H j (a 0 ) .
From ( 7) and ( 19), the positive definiteness of Ω n,H and the fact that the process (e H t (a 0 )) t∈Z defined in ( 8) is centered, we obtain that

E 1 n Z T n,H U n,H (a 0 ) 2 = 1 n 2 n j 1 =1 n j 2 =1 a j 1 +j 2 -2 0 E L j 1 U n,H (a 0 ) T U n,H (a 0 ) L j 2 U n,H (a 0 ) T U n,H (a 0 ) = 1 n 2 n j 1 ,j 2 =1 a j 1 +j 2 -2 0 n k 1 ,k 2 =1 n r 1 ,...,r 4 =1 Ω -1/2 n,H k 1 ,r 1 Ω -1/2 n,H k 1 ,r 2 Ω -1/2 n,H k 2 ,r 3 Ω -1/2 n,H k 2 ,r 4 × E e H r 1 +1-j 1 (a 0 )e H r 2 +1 (a 0 )e H r 3 +1-j 2 (a 0 )e H r 4 +1 (a 0 ) = ∆ 1,n,H + ∆ 2,n,H + ∆ 3,n,H + ∆ 4,n,H , (39)
where

∆ 1,n,H = 1 n 2 n j 1 ,j 2 =1 a j 1 +j 2 -2 0 n r 1 =j 1 n r 2 =1 n r 3 =j 2 n r 4 =1 Ω -1 n,H r 1 ,r 2 Ω -1
n,H r 3 ,r 4

× cum e H r 1 +1-j 1 (a 0 ), e H r 2 +1 (a 0 ), e H r 3 +1-j 2 (a 0 ), e H r 4 +1 (a 0 ) ,

∆ 2,n,H = 1 n 2 n j 1 ,j 2 =1 a j 1 +j 2 -2 0 n r 1 =j 1 n r 2 =1 n r 3 =j 2 n r 4 =1 Ω -1 n,H r 1 ,r 2 Ω -1 n,H r 3 ,r 4 × E e H r 1 +1-j 1 (a 0 )e H r 2 +1 (a 0 ) E e H r 3 +1-j 2 (a 0 )e H r 4 +1 (a 0 ) , ∆ 3,n,H = 1 n 2 n j 1 ,j 2 =1 a j 1 +j 2 -2 0 n r 1 =j 1 n r 2 =1 n r 3 =j 2 n r 4 =1 Ω -1 n,H r 1 ,r 2 Ω -1 n,H r 3 ,r 4 × E e H r 1 +1-j 1 (a 0 )e H r 3 +1-j 2 (a 0 ) E e H r 2 +1 (a 0 )e H r 4 +1 (a 0 ) and ∆ 4,n,H = 1 n 2 n j 1 ,j 2 =1 a j 1 +j 2 -2 0 n r 1 =j 1 n r 2 =1 n r 3 =j 2 n r 4 =1 Ω -1 n,H r 1 ,r 2 Ω -1 n,H r 3 ,r 4
× E e H r 1 +1-j 1 (a 0 )e H r 4 +1 (a 0 ) E e H r 2 +1 (a 0 )e H r 3 +1-j 2 (a 0 ) .

Since the vector ( H 1 , . . . , H n+1 ) T is Gaussian random vector, any joint cumulant involving three (or more) elements of ( H 1 , . . . , H n+1 ) T is zero. The random variables e H r 1 +1-j 1 (a 0 ), e H r 2 +1 (a 0 ), e H r 3 +1-j 2 (a 0 ) and e H r 4 +1 (a 0 ) in ( 39) are elements of the Gaussian vector ( H 1 , . . . , H n+1 ) T so cum e H r 1 +1-j 1 (a 0 ), e H r 2 +1 (a 0 ), e H r 3 +1-j 2 (a 0 ), e H r 4 +1 (a 0 ) = 0.

Thus, we have ∆ 1,n,H = 0.

Let Ω n+1,H be the covariance matrix of the random vector ( H 1 , . . . , H n+1 ) T . Thanks to the stationarity of the process ( H t ) t∈Z one has (Ω n+1,H ) r 3 +1-j 2 ,r 4 +1 = (Ω n,H ) r 3 -j 2 ,r 4 when j 2 , r 4 = 1, . . . , n and j 2 + 1 ≤ r 3 ≤ n. Thus ∆ 2,n,H can be rewritten as:

∆ 2,n,H = 1 n 2 n j 1 ,j 2 =1 a j 1 +j 2 -2 0 n r 1 =j 1 n r 2 =1 n r 3 =j 2 n r 4 =1 Ω -1 n,H r 1 ,r 2 Ω -1 n,H r 3 ,r 4 (Ω n+1,H ) r 1 +1-j 1 ,r 2 +1 (Ω n+1,H ) r 3 +1-j 2 ,r 4 +1 = 1 n 2 n j 1 ,j 2 =1 a j 1 +j 2 -2 0 n r 1 =j 1 n r 2 =1 n r 3 =j 2 +1 n r 4 =1 Ω -1 n,H r 1 ,r 2 Ω -1 n,H r 3 ,r 4 (Ω n+1,H ) r 1 +1-j 1 ,r 2 +1 (Ω n,H ) r 3 -j 2 ,r 4 + 1 n 2 n j 1 ,j 2 =1 a j 1 +j 2 -2 0 n r 1 =j 1 n r 2 =1 n r 4 =1 Ω -1 n,H r 1 ,r 2 Ω -1 n,H j 2 ,r 4 (Ω n+1,H ) r 1 +1-j 1 ,r 2 +1 (Ω n+1,H ) 1,r 4 +1 .
This implies that

∆ 2,n,H = 1 n 2 n j 1 ,j 2 =1 a j 1 +j 2 -2 0 n r 1 =j 1 n r 2 =1 n r 3 =j 2 +1 Ω -1 n,H r 1 ,r 2 (Ω n+1,H ) r 1 +1-j 1 ,r 2 +1 (I n ) r 3 ,r 3 -j 2 + 1 n 2 n j 1 ,j 2 =1 a j 1 +j 2 -2 0 n r 1 =j 1 +1 n r 2 =1 n r 4 =1 Ω -1 n,H r 1 ,r 2 Ω -1 n,H j 2 ,r 4 (Ω n,H ) r 1 -j 1 ,r 2 (Ω n+1,H ) 1,r 4 +1 + 1 n 2 n j 1 ,j 2 =1 a j 1 +j 2 -2 0 n r 2 =1 n r 4 =1 Ω -1 n,H j 1 ,r 2 Ω -1 n,H j 2 ,r 4 (Ω n+1,H ) 1,r 2 +1 (Ω n+1,H ) 1,r 4 +1 = 1 n 2 n j 1 ,j 2 =1 a j 1 +j 2 -2 0 n r 1 =j 1 n r 2 =1 n r 3 =j 2 +1 Ω -1 n,H r 1 ,r 2 (Ω n+1,H ) r 1 +1-j 1 ,r 2 +1 (I n ) r 3 ,r 3 -j 2 + 1 n 2 n j 1 ,j 2 =1 a j 1 +j 2 -2 0 n r 1 =j 1 +1 n r 4 =1 Ω -1 n,H j 2 ,r 4 (Ω n+1,H ) 1,r 4 +1 (I n ) r 1 ,r 1 -j 1 + 1 n 2 n j 1 ,j 2 =1 a j 1 +j 2 -2 0 n r 2 =1 n r 4 =1 Ω -1 n,H j 1 ,r 2 Ω -1 n,H j 2 ,r 4 (Ω n+1,H ) 1,r 2 +1 (Ω n+1,H ) 1,r 4 +1 . Since j 1 , j 2 = 1, . . . , n, (I n ) r 1 ,r 1 -j 1 = (I n ) r 3 ,r 3 -j 2 = 0, we obtain that ∆ 2,n,H = 1 n 2 n j 1 ,j 2 =1 a j 1 +j 2 -2 0 n r 2 =1 n r 4 =1 Ω -1 n,H j 1 ,r 2 Ω -1 n,H j 2 ,r 4 (Ω n+1,H ) 1,r 2 +1 (Ω n+1,H ) 1,r 4 +1 =   1 n n j=1 a j-1 0 n r =1 Ω -1 n,H j,r (Ω n+1,H ) 1,r +1   2 . ( 40
)
Using a similar approach as that given in the proof of Lemma 1, one can confirm that the term in the right side of Equation ( 40) is equal to the square of the expected value of the random variable Z T n,H U n,H (a 0 )/n. Similar calculations to those of the proof of Lemma 1 imply that

E 1 n Z T n,H U n,H (a 0 ) ---→ n→∞ 0.
Let us be more precise. From ( 19), we have

E 1 n Z T n,H U n,H (a 0 ) = 1 n n j=1 a j-1 0 n k=1 E L j (U n,H (a 0 )) k (U n,H (a 0 )) k = 1 n n j=1 a j-1 0 n k=1 n r 1 ,r 2 =1 Ω -1/2 n,H k,r 1 Ω -1/2 n,H k,r 2 E e H r 1 +1-j (a 0 )e H r 2 +1 (a 0 ) (41) = 1 n n j=1 a j-1 0 n r 1 =j n r 2 =1 Ω -1 n,H r 1 ,r 2 (Ω n+1,H ) r 1 +1-j,r 2 +1 . ( 42
)
As before, we proceed by cases to deduce that

E 1 n Z T n,H U n,H (a 0 ) = 1 n n j=1 a j-1 0 n r 1 =j+1 n r 2 =1 Ω -1 n,H r 1 ,r 2 (Ω n,H ) r 1 -j,r 2 + 1 n n j=1 a j-1 0 n r 2 =1 Ω -1 n,H j,r 2 (Ω n+1,H ) 1,r 2 +1 = 1 n n j=1 a j-1 0 n r 1 =j+1 (I n ) r 1 ,r 1 -j + 1 n n j=1 a j-1 0 n r 2 =1 Ω -1 n,H j,r 2 (Ω n+1,H ) 1,r 2 +1 = 1 n n j=1 a j-1 0 n r 2 =1 Ω -1 n,H j,r 2 (Ω n+1,H ) 1,r 2 +1 = 1 n n j=1 a j-1 0 Ω -1 n,H j,j (Ω n+1,H ) 1,j+1 + 1 n n j=1 a j-1 0 1≤r 2 ≤n r 2 =j Ω -1 n,H j,r 2 (Ω n+1,H ) 1,r 2 +1 .
Finally, in view of ( 29) and ( 30) we conclude that

E 1 n Z T n,H U n,H (a 0 ) ≤ K    1 n n j=1 |a 0 | j-1 + 1 n n j=1 |a 0 | j-1 1≤r ≤n r =j 1 (r -j) 2H 1 r 2-2H    . (43)
Using the fact that the autoregressive parameter a 0 is assumed to belong to (-1, 1), Cesàro's Lemma implies that the first term in the right hand side of (43) tends to 0.

Proceeding by a simple calculation and using Cesàro's Lemma, one can also prove that the second term in the right hand side of (43) converges to 0. More explicitly,

1 n n j=1 |a 0 | j-1 1≤r ≤n r =j 1 (r -j) 2H 1 r 2-2H = 1 n n r =1 1 r 2-2H 1≤j≤n j =r |a 0 | j-1 1 (r -j) 2H = 1 n n r =1 1 r 2-2H |a 0 | r -1 1-r ≤j≤n-r j =0 |a 0 | j 1 j 2H = 1 n n r =2 1 r 2-2H |a 0 | r -1 -1 j=1-r |a 0 | j 1 j 2H + 1 n n r =1 1 r 2-2H |a 0 | r -1 n-r j=1 |a 0 | j 1 j 2H ≤ 1 n n r =1 1 r 1-2H |a 0 | r -1 max 1 |a 0 | , 1 |a 0 | r r 2H + 1 n n r =1 1 r 2-2H |a 0 | r 1 -|a 0 | n-r 1 -|a 0 | ---→ n→∞ 0.
Therefore, we deduce that

∆ 2,n,H ---→ n→∞ 0.
It is obvious that by exchanging the roles of j 1 and j 2 and of r 2 and r 4 in the expression of ∆ 4,n,H , we obtain that ∆ 4,n,H = ∆ 2,n,H . Then, the sequence (∆ 4,n,H ) n≥1 also converges to 0.

It remains to show that (∆ 3,n,H ) n≥1 converges to 0 to obtain the convergence in mean square of (Z T n,H U n,H (a 0 )/n) n≥1 to 0 and thus the convergence in probability of the generalized least squares estimator ân defined in (9) to the true parameter a 0 .

We use the stationarity of the process (e H t (a 0 )) t∈Z and we follow the same calculation procedures developed to obtain the explicit expression of ∆ 2,n,H given in ( 40). This leads to

∆ 3,n,H = 1 n 2 n j 1 ,j 2 =1 a j 1 +j 2 -2 0 n r 1 =j 1 n r 2 =1 n r 3 =j 2 n r 4 =1 Ω -1 n,H r 1 ,r 2 Ω -1 n,H r 3 ,r 4 (Ω n,H ) r 1 +1-j 1 ,r 3 +1-j 2 (Ω n+1,H ) r 2 +1,r 4 +1 = 1 n 2 n j 1 ,j 2 =1 a j 1 +j 2 -2 0 n r 1 =j 1 n r 2 =1 n r 3 =j 2 n r 4 =1 Ω -1 n,H r 1 ,r 2 Ω -1 n,H r 3 ,r 4 (Ω n,H ) r 1 +1-j 1 ,r 3 +1-j 2 (Ω n,H ) r 2 ,r 4 = 1 n 2 n j 1 ,j 2 =1 a j 1 +j 2 -2 0 n r 1 =j 1 n r 3 =j 2 n r 4 =1 Ω -1 n,H r 3 ,r 4 (Ω n,H ) r 1 +1-j 1 ,r 3 +1-j 2 (I n ) r 1 ,r 4 = 1 n 2 n j 1 ,j 2 =1 a j 1 +j 2 -2 0 n r 1 =j 1 n r 3 =j 2 Ω -1 n,H r 3 ,r 1 (Ω n,H ) r 1 +1-j 1 ,r 3 +1-j 2 = 1 n 2 n j 1 =1 j 1 j 2 =1 a j 1 +j 2 -2 0 n r 1 =j 1 n r 3 =j 2 Ω -1 n,H r 3 ,r 1 (Ω n,H ) r 1 +1-j 1 ,r 3 +1-j 2 + 1 n 2 n j 1 =1 n j 2 =j 1 +1 a j 1 +j 2 -2 0 n r 1 =j 1 n r 3 =j 2 Ω -1 n,H r 3 ,r 1 (Ω n,H ) r 1 +1-j 1 ,r 3 +1-j 2 .
The stationarity of ( H t ) t∈Z implies that • for j 2 ≤ j 1 :

(Ω n,H ) r 1 +1-j 1 ,r 3 +1-j 2 = (Ω n,H ) r 1 +j 2 -j 1 ,r 3 ,

• for j 2 > j 1 :

(Ω n,H ) r 1 +1-j 1 ,r 3 +1-j 2 = (Ω n,H ) r 1 ,r 3 +j 1 -j 2 .
It follows that

∆ 3,n,H = 1 n 2 n j 1 =1 j 1 j 2 =1 a j 1 +j 2 -2 0 n r 1 =j 1 n r 3 =j 2 Ω -1 n,H r 3 ,r 1 (Ω n,H ) r 1 +j 2 -j 1 ,r 3 + 1 n 2 n j 1 =1 n j 2 =j 1 +1 a j 1 +j 2 -2 0 n r 1 =j 1 n r 3 =j 2 Ω -1 n,H r 3 ,r 1 (Ω n,H ) r 1 ,r 3 +j 1 -j 2 = 1 n 2 n j 1 =1 j 1 j 2 =1 a j 1 +j 2 -2 0 n r 1 =j 1 n r 3 =1 Ω -1 n,H r 3 ,r 1 (Ω n,H ) r 1 +j 2 -j 1 ,r 3 - 1 n 2 n j 1 =1 j 1 j 2 =1 a j 1 +j 2 -2 0 n r 1 =j 1 j 2 -1 r 3 =1 Ω -1 n,H r 3 ,r 1 (Ω n,H ) r 1 +j 2 -j 1 ,r 3 + 1 n 2 n j 1 =1 n j 2 =j 1 +1 a j 1 +j 2 -2 0 n r 1 =1 n r 3 =j 2 Ω -1 n,H r 3 ,r 1 (Ω n,H ) r 1 ,r 3 +j 1 -j 2 - 1 n 2 n j 1 =1 n j 2 =j 1 +1 a j 1 +j 2 -2 0 j 1 -1 r 1 =1 n r 3 =j 2 Ω -1 n,H r 3 ,r 1 (Ω n,H ) r 1 ,r 3 +j 1 -j 2 . (44)
Thus we have

∆ 3,n,H = 1 n 2 n j 1 =1 (n -j 1 + 1)a 2j 1 -2 0 - 1 n 2 n j 1 =1 j 1 j 2 =1 a j 1 +j 2 -2 0 n r 1 =j 1 j 2 -1 r 3 =1 Ω -1 n,H r 3 ,r 1 (Ω n,H ) r 1 +j 2 -j 1 ,r 3 - 1 n 2 n j 1 =1 n j 2 =j 1 +1 a j 1 +j 2 -2 0 j 1 -1 r 1 =1 n r 3 =j 2 Ω -1 n,H r 3 ,r 1 (Ω n,H ) r 1 ,r 3 +j 1 -j 2 = 1 n E 1 n Z n,H 2 R n . (45)
Now, using Lemma 1 we deduce that

∆ 3,n,H ---→ n→∞ 0.
The lemma is then proved.

The convergences established in Lemmas 3 and 4 imply that the joint vector

(Z T n,H U n,H (a 0 ), Z n,H 2 R n )/n converges in probability to (0, 1/(1 -a 2 0 ))
. Let now g be the continuous function defined on R × R * by g (x, y ) = xy -1 , it is obvious that the vector (0, 1/(1 -a 2 0 )) belongs to the set of continuity points of the function g . The continuous mapping theorem yields

ân -a 0 = g Z T n,H U n,H (a 0 )/n, Z n,H 2 R n /n P ---→ n→∞ g 0, 1/(1 -a 2 0 ) = 0.
The proof of Theorem 1 is then complete.

Remark 4. The following equivalent formulation of Equality (41)

E 1 √ n Z T n,H U n,H (a 0 ) = 1 √ n n j=1 a j-1 0 n r 1 =j n r 2 =1 Ω -1 n,H r 1 ,r 2 E e H r 1 +1-j (a 0 )e H r 2 +1 (a 0 ) (46)
will be helpful in the proof of Theorem 2.

Proof of Theorem 2

We establish in this section the proof of the asymptotic normality of the generalized least squares estimator ân given in (9) under the causality of (X t ) t∈Z (i.e. when the true parameter a 0 is assumed to belong to (-1, 1)) and in the case where H is known.

Recall that the normal distribution admits moments of all orders and is completely determined by its moments. The expected convergence in distribution stated in Theorem 2 can be obtained using the method of moments (see Theorem 30.2. of [START_REF] Billingsley | Probability and measure[END_REF]).

In our case, Theorem 2 will be proved as soon as for any r ≥ 1

lim n→∞ E √ n(â n -a 0 ) r = E [G r 1 ] , (47)
where G 1 is normally distributed with mean 0 and variance 1 -a 2 0 . In view of the explicit expression of âna 0 (see ( 10)), the random variable √ n(â n -a 0 ) can be rewritten in the form:

√ n(â n -a 0 ) = 1 √ n Z T n,H U n,H (a 0 ) 1 n Z n,H 2 R n . ( 48
)
From Lemma 3, the random variable

Z n,H 2 R n /n converges in mean square to 1/(1 -a 2 0 ). This implies that the sequence ( Z n,H 2 R n /n) n≥1 converges in probability to 1/(1 -a 2 0
). The idea of the proof is to show using the method of moments that Z T n,H U n,H (a 0 )/ √ n has a limiting centered normal distribution with variance 1/(1 -a 2 0 ). The proof will be concluded by applying Slutsky's theorem to

Z n,H 2 R n /n and Z T n,H U n,H (a 0 )/ √ n. The asymptotic distribution of the sequence (Z T n,H U n,H (a 0 )/ √ n) n≥1
is given in the following lemma.

Lemma 5. Under Assumption (A0), the sequence of random variables (Z T n,H U n,H (a 0 )/ √ n) n≥1 has a limiting centered normal distribution with variance 1/(1 -a 2 0 ). Proof. Using the existence of moments of all orders of Z T n,H U n,H (a 0 )/ √ n and in view of Theorem 30.2. of [START_REF] Billingsley | Probability and measure[END_REF], the proof of Lemma 5 is obtained once

lim n→∞ E 1 √ n Z T n,H U n,H (a 0 ) r = E [G r 2 ] , for all r ≥ 1, (49)
where G 2 is normally distributed with mean 0 and variance 1/(1 -a 2 0 ). To prove (49), in a first step we consider an odd order r and in a second step r will be assumed to be even.

Recall that for all r ≥ 1, we have

E [G r 2 ] = 0 if r is odd 1/(1 -a 2 0 ) r /2 (r -1)!! if r is even, (50)
where

r !! = r (r -2)(r -4) • • • 1.
We start from (19), we use the fact that the random vector (e H 1 (a 0 ), . . . , e H n+1 (a 0 )) T is a centered Gaussian random vector and we exploit the expression of moments in terms of cumulants to obtain that

E 1 √ n Z T n,H U n,H (a 0 ) r = 1 n r /2 n j 1 ,j 2 ,...,jr =1 a r i=1 (j i -1) 0 n k 1 =j 1 n k 2 =1 n k 3 =j 2 n k 4 =1 • • • n k 2r -1 =jr n k 2r =1 r i=1 Ω -1 n,H k 2i-1 ,k 2i ×E r i=1 e H k 2i-1 +1-j i (a 0 )e H k 2i +1 (a 0 ) = 1 n r /2 n j 1 ,j 2 ,...,jr =1 a r i=1 (j i -1) 0 n k 1 =j 1 n k 2 =1 n k 3 =j 2 n k 4 =1 • • • n k 2r -1 =jr n k 2r =1 r i=1 Ω -1 n,H k 2i-1 ,k 2i × π B∈π cum e H i (a 0 ) : i ∈ B , (51)
where π runs through the list of all unordered pair partitions of {k 1 + 1 -j 1 , k 2 + 1, k 3 + 1 -j 2 , k 4 + 1, . . . , k 2r -1 + 1 -j r , k 2r + 1} and B runs through the list of all blocks of the partition π (see [START_REF] Brillinger | Time series: data analysis and theory[END_REF]).

Remark 5. In the above expression, we only kept unordered pair partitions of

{k 1 + 1 -j 1 , k 2 + 1, k 3 + 1 -j 2 , k 4 + 1, . . . , k 2r -1 + 1 -j r , k 2r + 1}
. This is due to the fact that the random vector (e H 1 (a 0 ), . . . , e H n+1 (a 0 )) T is a centered Gaussian vector. Any joint cumulant involving three (or more) elements of (e H 1 (a 0 ), . . . , e H n+1 (a 0 )) T is zero. Moreover, since the joint cumulant of just one random variable is its expectation and the process (e H t (a 0 )) t∈Z is centered, the first cumulant of any element of (e H 1 (a 0 ), . . . , e H n+1 (a 0 )) T vanishes. Thanks to the above remark, we may write (51) as

E 1 √ n Z T n,H U n,H (a 0 ) r = 1 n r /2 n j 1 ,j 2 ,...,jr =1 a r i=1 (j i -1) 0 n k 1 =j 1 n k 2 =1 n k 3 =j 2 n k 4 =1 • • • n k 2r -1 =jr n k 2r =1 r i=1 Ω -1 n,H k 2i-1 ,k 2i × π B∈π E e H i (a 0 ) : i ∈ B . ( 52
)
We strength the fact that when r = 1, (52) is ( 46) and for r = 2 we retrieve an analogous formula to (39).

The following remark will help us to develop (51) more explicitly.

Remark 6.

• In Equation (51), the exchange of the roles of indices with even sub-indices does not modify the term in the right hand side. It is the same when we exchange the roles of indices with odd sub-indices. For example, for the indices k 1 and k 3 it holds that

1 n r /2 n j 1 ,j 2 ,...,jr =1 a r i=1 (j i -1) 0 n k 1 =j 1 n k 2 =1 n k 3 =j 2 n k 4 =1 • • • n k 2r -1 =jr n k 2r =1 r i=1 Ω -1 n,H k 2i-1 ,k 2i × π B∈π cum e H i (a 0 ) : i ∈ B = 1 n r /2 n j 1 ,j 2 ,...,jr =1 a r i=1 (j i -1) 0 n k 1 =j 2 n k 2 =1 n k 3 =j 1 n k 4 =1 • • • n k 2r -1 =jr n k 2r =1 r i=1 Ω -1 n,H k 2i-1 ,k 2i × π B∈π cum e H i (a 0 ) : i ∈ B ,
• Again in Equation (51) when r is odd, every pair partition π of the set of indices {k 1 + 1j 1 , k 2 + 1, k 3 + 1 -j 2 , k 4 + 1, . . . , k 2r -1 + 1 -j r , k 2r + 1} contains at least one block of the type {k 2i-1 + 1 -j i , k i + 1}. To see that, let us consider the case when r = 3. The third moment of Z T n,H U n,H (a 0 )/ √ n is given by:

E 1 √ n Z T n,H U n,H (a 0 ) 3 = 1 n 3/2 n j 1 ,j 2 ,j 3 =1 a j 1 +j 2 +j 3 -3 0 n r 1 =j 1 n r 2 =1 n r 3 =j 2 n r 4 =1 n r 5 =j 3 n r 6 =1 Ω -1 n,H r 1 ,r 2 Ω -1 n,H r 3 ,r 4 Ω -1
n,H r 5 ,r 6

× E e H r 1 +1-j 1 (a 0 )e H r 2 +1 (a 0 )e H r 3 +1-j 2 (a 0 )e H r 4 +1 (a 0 )e H r 5 +1-j 3 (a 0 )e H r 6 +1 (a 0 ) = 1 n 3/2 n j 1 ,j 2 ,j 3 =1 a j 1 +j 2 +j 3 -3

0 n r 1 =j 1 n r 2 =1 n r 3 =j 2 n r 4 =1 n r 5 =j 3 n r 6 =1 Ω -1 n,H r 1 ,r 2 Ω -1 n,H r 3 ,r 4 Ω -1 n,H r 5 ,r 6
× E e H r 1 +1-j 1 (a 0 )e H r 2 +1 (a 0 ) E e H r 3 +1-j 2 (a 0 )e H r 4 +1 (a 0 ) E e H r 5 +1-j 3 (a 0 )e H r 6 +1 (a 0 )

+ E e H r 1 +1-j 1 (a 0 )e H r 2 +1 (a 0 ) E e H r 3 +1-j 2 (a 0 )e H r 5 +1-j 3 (a 0 ) E e H r 4 +1 (a 0 )e H r 6 +1 (a 0 )

+ E e H r 1 +1-j 1 (a 0 )e H r 2 +1 (a 0 ) E e H r 3 +1-j 2 (a 0 )e H r 6 +1 (a 0 ) E e H r 4 +1 (a 0 )e H r 5 +1-j 3 (a 0 )

+ E e H r 1 +1-j 1 (a 0 )e H r 3 +1-j 2 (a 0 ) E e H r 2 +1 (a 0 )e H r 4 +1 (a 0 ) E e H r 5 +1-j 3 (a 0 )e H r 6 +1 (a 0 )

+ E e H r 1 +1-j 1 (a 0 )e H r 3 +1-j 2 (a 0 ) E e H r 2 +1 (a 0 )e H r 5 +1-j 3 (a 0 ) E e H r 4 +1 (a 0 )e H r 6 +1 (a 0 )

+ E e H r 1 +1-j 1 (a 0 )e H r 3 +1-j 2 (a 0 ) E e H r 2 +1 (a 0 )e H r 6 +1 (a 0 ) E e H r 4 +1 (a 0 )e H r 5 +1-j 3 (a 0 )

+ E e H r 1 +1-j 1 (a 0 )e H r 4 +1 (a 0 ) E e H r 2 +1 (a 0 )e H r 3 +1-j 2 (a 0 ) E e H r 5 +1-j 3 (a 0 )e H r 6 +1 (a 0 )

+ E e H r 1 +1-j 1 (a 0 )e H r 4 +1 (a 0 ) E e H r 2 +1 (a 0 )e H r 5 +1-j 3 (a 0 ) E e H r 3 +1-j 2 (a 0 )e H r 6 +1 (a 0 )

+ E e H r 1 +1-j 1 (a 0 )e H r 4 +1 (a 0 ) E e H r 2 +1 (a 0 )e H r 6 +1 (a 0 ) E e H r 3 +1-j 2 (a 0 )e H r 5 +1-j 3 (a 0 )

+ E e H r 1 +1-j 1 (a 0 )e H r 5 +1-j 3 (a 0 ) E e H r 2 +1 (a 0 )e H r 3 +1-j 2 (a 0 ) E e H r 4 +1 (a 0 )e H r 6 +1 (a 0 )

+ E e H r 1 +1-j 1 (a 0 )e H r 5 +1-j 3 (a 0 ) E e H r 2 +1 (a 0 )e H r 4 +1 (a 0 ) E e H r 3 +1-j 2 (a 0 )e H r 6 +1 (a 0 )

+ E e H r 1 +1-j 1 (a 0 )e H r 5 +1-j 3 (a 0 ) E e H r 2 +1 (a 0 )e H r 6 +1 (a 0 ) E e H r 3 +1-j 2 (a 0 )e H r 4 +1 (a 0 )

+ E e H r 1 +1-j 1 (a 0 )e H r 6 +1 (a 0 ) E e H r 2 +1 (a 0 )e H r 3 +1-j 2 (a 0 ) E e H r 4 +1 (a 0 )e H r 5 +1-j 3 (a 0 )

+ E e H r 1 +1-j 1 (a 0 )e H r 6 +1 (a 0 ) E e H r 2 +1 (a 0 )e H r 4 +1 (a 0 ) E e H r 3 +1-j 2 (a 0 )e H r 5 +1-j 3 (a 0 )

+ E e H r 1 +1-j 1 (a 0 )e H r 6 +1 (a 0 ) E e H r 2 +1 (a 0 )e H r 5 +1-j 3 (a 0 ) E e H r 3 +1-j 2 (a 0 )e H r 4 +1 (a 0 ) .

We can use this last remark and Equations ( 46) and ( 39) to show that

E 1 √ n Z T n,H U n,H (a 0 ) 3 =   1 √ n n j 1 =1 a j 1 -1 0 n r 1 =j 1 n r 2 =1 Ω -1 n,H r 1 ,r 2 E e H r 1 +1-j 1 (a 0 )e H r 2 +1 (a 0 )   × 1 n n j 2 ,j 3 =1 a j 2 +j 3 -2 0 n r 3 =j 2 n r 4 =1 n r 5 =j 3 n r 6 =1 Ω -1 n,H r 3 ,r 4 Ω -1 n,H r 5 ,r 6 × E e H r 3 +1-j 2 (a 0 )e H r 4 +1 (a 0 ) E e H r 5 +1-j 3 (a 0 )e H r 6 +1 (a 0 ) + E e H r 3 +1-j 2 (a 0 )e H r 5 +1-j 3 (a 0 ) E e H r 4 +1 (a 0 )e H r 6 +1 (a 0 ) + E e H r 3 +1-j 2 (a 0 )e H r 6 +1 (a 0 ) E e H r 4 +1 (a 0 )e H r 5 +1-j 3 (a 0 ) = E 1 √ n Z T n,H U n,H (a 0 ) E 1 √ n Z T n,H U n,H (a 0 ) 2 .
In view of (39), (42) and Remarks 5 and 6, one can rewrite the r -th moment of the random variable Z T n,H U n,H (a 0 )/ √ n when r is odd in the form:

E 1 √ n Z T n,H U n,H (a 0 ) r = E 1 √ n Z T n,H U n,H (a 0 ) E 1 √ n Z T n,H U n,H (a 0 ) r -1 . (53)
The calculation given in the proof of Lemma 4 to show the convergence to 0 of the sequence (E[Z T n,H U n,H (a 0 )/n]) n≥1 remains valid when we multiply the general term of the last sequence by √ n and when we use the fractional version of Cesàro's Lemma 1 to conclude. In fact, from (43) we have

E 1 √ n Z T n,H U n,H (a 0 ) ≤ K    1 √ n n j=1 |a 0 | j-1 + 1 √ n n j=1 |a 0 | j-1 1≤r ≤n r =j 1 (r -j) 2H 1 r 2-2H    . ( 54
)
1 Recall that the fractional version of Cesàro's Lemma states that for (ht ) t≥0 a sequence of positive real numbers, κ > 0 and c ≥ 0 we have

lim t→∞ ht t 1-κ = |κ| c ⇒ lim n→∞ 1 n κ n t=0 ht = c.

Since

√ j |a 0 | j-1 → 0 when j → ∞, fractional version of Cesàro's Lemma implies that the first term in the right hand side of (54) tends to 0 when n goes to infinity. Based on the calculation given in the proof of Lemma 4 and using the fractional version of Cesàro's Lemma, the second term in the right hand side of (54) verifies

1 √ n n j=1 |a 0 | j-1 1≤r ≤n r =j 1 (r -j) 2H 1 r 2-2H ≤ 1 √ n n r =1 1 r 1-2H |a 0 | r -1 max 1 |a 0 | , 1 |a 0 | r r 2H + 1 √ n n r =1 1 r 2-2H |a 0 | r 1 -|a 0 | n-r 1 -|a 0 | ---→ n→∞ 0.
Then the sequence (E[Z T n,H U n,H (a 0 )/

√ n]) n≥1 converges to 0. Furthermore, Equation ( 45) and Lemma 1 imply that

lim n→∞ E 1 √ n Z T n,H U n,H (a 0 ) 2 = lim n→∞ E 1 n Z n,H 2 R n = 1 1 -a 2 0 .
Using (53) we deduce that

lim n→∞ E 1 √ n Z T n,H U n,H (a 0 ) 3 = 0. (55)
When r = 4, the set {k 1 + 1 -j 1 , k 2 + 1, k 3 + 1 -j 2 , k 4 + 1, . . . , k 2r -1 + 1 -j r , k 2r + 1} can be partitioned by pairs in two different ways. The first is to consider partitions that contain at least one block of type {k 2i-1 + 1 -j i , k i + 1} and the second is to take partitions that do not contain any such block. In the first case, we use the same remarks introduced to obtain (53) and we factorize by the expectation of Z T n,H U n,H (a 0 )/ √ n. In the second case, we use the result stated in Lemma 1. More precisely, we have

E 1 √ n Z T n,H U n,H (a 0 ) 4 = 3 E 1 n Z n,H 2 R n 2 + E 1 √ n Z T n,H U n,H (a 0 ) E 1 √ n Z T n,H U n,H (a 0 ) 3 .
Therefore, in view of Lemma 1 and (55) we have

lim n→∞ E 1 √ n Z T n,H U n,H (a 0 ) 4 = 3 (1 -a 2 0 ) 2 .
We can generalize this result by an immediate recurrence on r to obtain

E 1 √ n Z T n,H U n,H (a 0 ) r =        E 1 √ n Z T n,H U n,H (a 0 ) E 1 √ n Z T n,H U n,H (a 0 ) r -1 if r is odd (r -1)!! E 1 n Z n,H 2 R n r /2 + E 1 √ n Z T n,H U n,H (a 0 ) E 1 √ n Z T n,H U n,H (a 0 ) r -1 if r is even.
This implies that the r -th moment of the random variable Z

T n,H U n,H (a 0 )/ √ n converges to    0 if r is odd (r -1)!! 1 1-a 2 0 r /2 if r is even, ( 56
)
and this is exactly the result recalled in (50). The lemma is then proved.

The convergences established in Lemmas 3 and 5 imply (by Slutsky's theorem) that the vector

(Z T n,H U n,H (a 0 )/ √ n, Z n,H 2 R n /n) converges in distribution to (G 2 , 1/(1 -a 2 0 ))
, where G 2 is the random variable defined in the proof of Lemma 5. Continuous mapping theorem allows to complete the proof of Theorem 2. -

Figures

-3 -2 -1 0 1 2 3 -4 -3 -2 -1 0 1 2 3 Normal Q-Q Plot of n (a ^n -a0) Distribution of n (a ^n -a0) Density -4 -3 -2 -1 0 1 2 3 0.0 0.1 0.2 0.3 0.4 -3 -2 -1 0 1 2 3 -0.6 -0.2 0.2 0.4 Normal Q-Q Plot of
3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 Normal Q-Q Plot of n (a ^n -a0) Distribution of n (a ^n -a0) Density -3 -2 -1 0 1 2 3 0.0 0.1 0.2 0.3 0.4 -3 -2 -1 0 1 2

Fig 1 .

 1 Fig 1. Distribution of the generalized least squares estimator of N = 2, 000 independent simulations of the fractional autoregressive model (14) with size n = 30 and unknown parameter a = 0.3. The Hurst index is assumed to be equal to 0.7.

Fig 2 .

 2 Fig 2. Distribution of the generalized least squares estimator of N = 2, 000 independent simulations of the fractional autoregressive model (14) with size n = 500 and unknown parameter a = 0.3. The Hurst index is assumed to be equal to 0.7.

Fig 3 .Fig 4 .Fig 5 .Fig 6 .Fig 7 .

 34567 Fig 3. Distribution of the generalized least squares estimator of N = 2, 000 independent simulations of the fractional autoregressive model (14) with size n = 2, 000 and unknown parameter a = 0.3. The Hurst index is assumed to be equal to 0.7.

Fig 8 .

 8 Fig 8. Comparison of the logarithm of quadratic risk of the generalized least squares estimator (in black) and Whittle estimator (in green) of N = 1, 000 independent simulations of the fractional autoregressive model (14) with a = -0.9 and with different values of H. The sample size n is (from left to right and from top to bottom) n = 50, 100, 150, 200.

Fig 9 .

 9 Fig 9. Comparison of the logarithm of quadratic risk of the generalized least squares estimator (in black) and Whittle estimator (in green) of N = 1, 000 independent simulations of the fractional autoregressive model (14) with a = -0.5 and with different values of H. The sample size n is (from left to right and from top to bottom) n = 50, 100, 150, 200.

Fig 10 .

 10 Fig 10. Comparison of the logarithm of quadratic risk of the generalized least squares estimator (in black) and Whittle estimator (in green) of N = 1, 000 independent simulations of the fractional autoregressive model (14) with a = 0 and with different values of H. The sample size n is (from left to right and from top to bottom) n = 50, 100, 150, 200.

Fig 11 .

 11 Fig 11. Comparison of the logarithm of quadratic risk of the generalized least squares estimator (in black) and Whittle estimator (in green) of N = 1, 000 independent simulations of the fractional autoregressive model (14) with a = 0.5 and with different values of H. The sample size n is (from left to right and from top to bottom) n = 50, 100, 150, 200.

  1≤i,j≤n . Remark 1. The covariance matrix Ω n,H depends only on the Hurst parameter H. As a symmetric real matrix, Ω n,H can be diagonalized by an orthogonal matrix. More explicitly, there exists a real orthogonal matrix P n,H such that D n,H = P T

	1/2 n,H :

n,H Ω n,H P n,H is a diagonal matrix. Since Ω n,H is positive definite (see Proposition 1 hereafter), we have Ω

Acknowledgments

The author acknowledges the support of the ANR project "Efficient inference for large and highfrequency data" (ANR-21-CE40-0021). Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40