Jean-Marc Azaïs

François Bachoc

Agnès Lagnoux

Thi Mong

Ngoc Nguyen

Semi-parametric estimation of the variogram scale parameter of a Gaussian process with stationary increments

Keywords: quadratic variations, scale covariance parameter, asymptotic normality, moment method, aggregation of estimators

We consider the semi-parametric estimation of the scale parameter of the variogram of a onedimensional Gaussian process with known smoothness. We suggest an estimator based on quadratic variations and on the moment method. We provide asymptotic approximations of the mean and variance of this estimator, together with asymptotic normality results, for a large class of Gaussian processes. We allow for general mean functions and study the aggregation of several estimators based on various variation sequences. In extensive simulation studies, we show that the asymptotic results accurately depict the finite-sample situations already for small to moderate sample sizes. We also compare various variation sequences and highlight the efficiency of the aggregation procedure.

Introduction

General context and state of the art Gaussian process models are widely used in statistics. For instance, they enable to interpolate observations by Kriging, notaby in computer experiment designs to build a metamodel [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF][START_REF] Stein | Interpolation of Spatial Data: Some Theory for Kriging[END_REF]. A second type of application of Gaussian processes is the analysis of local characteristics of images [START_REF] Richard | Anisotropy of hölder gaussian random fields: characterization, estimation, and application to image textures[END_REF] and one dimensional signals (e.g. in finance, see [START_REF] Wu | Gaussian process volatility model[END_REF][START_REF] Han | Financial time series volatility analysis using gaussian process state-space models[END_REF] and the references therein). A central problem with Gaussian processes is the estimation of the covariance function or the variogram. In this paper, we consider a real-valued Gaussian process (X(t)) t∈R with stationary increments. Its semi-variogram is well-defined and given by

V (h) := 1 2 E (X(t + h) -X(t)) 2 . (1)
Ideally, one aims at knowing perfectly the function V or at least estimate it precisely, either in a parametric setting or in a nonparametric setting. The parametric approach consists in assuming that the mean function of the Gaussian process (the drift) is a linear combination of known functions (often polynomials) and that the semi-variogram V belongs to a parametric family of semi-variograms {V θ , θ ∈ Θ ⊂ R p } for a given p in N * . Furthermore, in most practical cases, the semi-variogram is assumed to stem from a stationary autocovariance function k defined by k(h) = Cov(X(t), X(t + h)). In that case, the process is supposed to be stationary, and V can be rewritten in terms of the process autocovariance function k: V (h) = k(0) -k(h). Moreover, a parametric set of stationary covariance functions is considered of the form {k θ , θ ∈ Θ} with Θ ⊂ R p . In such a setting, several estimation procedures for θ have been introduced and studied in the literature. Usually in practice, most of the software packages (like, e.g. DiceKriging [START_REF] Roustant | Dicekriging, diceoptim: Two r packages for the analysis of computer experiments by kriging-based metamodeling and optimization[END_REF]) use the maximum likelihood estimation method (MLE) to estimate θ (see [START_REF] Stein | Interpolation of Spatial Data: Some Theory for Kriging[END_REF][START_REF] Santner | The design and analysis of computer experiments[END_REF][START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] for more details on MLE). Unfortunately, MLE is known to be computationally expensive and intractable for large data sets. In addition, it may diverge in some complicated situations (see Section 5.3). This has motivated the search for alternative estimation methods with a good balance between computational complexity and statistical efficiency. Among these methods, we can mention low rank approximation [START_REF] Stein | Limitations on low rank approximations for covariance matrices of spatial data[END_REF], sparse approximation [START_REF] Hensman | Gaussian processes for big data[END_REF], covariance tapering [START_REF] Furrer | Covariance tapering for interpolation of large spatial datasets[END_REF][START_REF] Kaufman | Covariance tapering for likelihood-based estimation in large spatial data sets[END_REF], Gaussian Markov random fields approximation [START_REF] Datta | Hierarchical nearest-neighbor Gaussian process models for large geostatistical datasets[END_REF][START_REF] Rue | Gaussian Markov random fields, Theory and applications[END_REF], submodel aggregation [START_REF] Cao | Generalized product of experts for automatic and principled fusion of Gaussian process predictions[END_REF][START_REF] Deisenroth | Distributed Gaussian processes[END_REF][START_REF] Hinton | Training products of experts by minimizing contrastive divergence[END_REF][START_REF] Rullière | Nested Kriging predictions for datasets with a large number of observations[END_REF][START_REF] Tresp | A Bayesian committee machine[END_REF][START_REF] Van Stein | Optimally weighted cluster Kriging for big data regression[END_REF] and composite likelihood [START_REF] Bachoc | Fixed-domain asymptotic properties of composite likelihood estimators for Gaussian processes[END_REF].

Framework and motivation The approaches discussed above are parametric. In this paper, we consider a more general semi-parametric context. The Gaussian process X is only assumed to have stationary increments and no parametric assumption is made on the semi-variogram V in [START_REF] Adler | Uniform quadratic variation for Gaussian processes[END_REF]. Assume, to simplify, that the semi-variogram is a C ∞ function outside 0. This is the case for most of the models even if the sample paths are not regular, see the examples in Section 2.1. Let D be the order of differentiability in quadratic mean of (X(t)) t∈R . This is equivalent to the fact that V is 2D differentiable and not 2D + 2 differentiable. Let us assume that the 2D'th derivative of V has the following expansion at the origin:

V (2D) (h) = V (2D) (0) + C(-1) D |h| s + r(h), (2)
where C 0, 0 < s < 2, and the remainder function r satisfies some hypothesis detailed further (see Section 2.1) and is a o(|h| s ) as h → 0. Note that, since s < 2, V is indeed not (2D + 2) differentiable. The quantity s is the smoothness parameter and we call C the scale parameter. In this paper, we assume D and s to be known and we focus on the theoretical study of the semi-parametric estimation of C defined in [START_REF] Azaïs | Level sets and extrema of random processes and fields[END_REF] in dimension one. Remark that this makes it possible to test whether the Gaussian process stems from a white noise or not. Notice that we also perform some additional simulations in higher dimensions. The value of C may lead to significantly different behaviors of the process X as one can see in Figure 1 which represents several realizations of a Gaussian process with exponential covariance function for different values of C (V (h) = 1 -exp(-C|h|), which satisfies [START_REF] Azaïs | Level sets and extrema of random processes and fields[END_REF] with D = 0). More concretely, for instance when D = 0, C provides the first order approximation of E (X(t + h) -X(t))

2 when h is small. Moreover, when D = 0, C = +∞ traduces independence, i.e. the process X reduces to a white noise whereas C = 0 corresponds to a constant process X.

As a motivating example, consider the following case where the estimation of C is beneficial. Assume that we observe a signal S depending on a vector of initial parameters x given by a computer code described by the application:

S : E → R [0,1] (3)
x → (S x (t), t ∈ [0, 1]),

where E stands for the initial parameter space. In order to interpret the output curve t → S x (t) for a given parameter x, it is useful to consider that this curve is the realization of a Gaussian process, with scale parameter C x . It is then insightful to study the black-box x → C x , for instance by means of a sensitivity analysis, or in the aim of finding which inputs x lead to C x < +∞, that is to output curves with a dependence structure. A necessary step to such a study is the estimation of the value of C x , given a discretized version of the curve t → S x (t). More generally, estimating C enables to assess if an observed signal is composed of independent components (C = +∞) or not (C < ∞), and to quantify the level of dependence. We refer to the real data sets studied in Section 5.3 for further discussion.

State of the art on variogram estimation Nonparametric estimation of the semi-variogram function is a difficult task since the resulting estimator must necessarily lead to a "valid" variogram (conditional negative definiteness property) [13, p. 93]. This requirement usually leads to complicated and computationally involved nonparametric estimators of the variogram [START_REF] Hall | On the nonparametric estimation of covariance functions[END_REF][START_REF] Hall | Properties of nonparametric estimators of autocovariance for stationary random fields[END_REF] that may need huge data sets to be meaningful. A simpler estimator based on the moment method has been proposed in [START_REF] Matheron | Traité de géostatistique appliquée, Tome I[END_REF][START_REF] Cressie | Robust estimation of the variogram: I[END_REF] but does not always conduce to a valid variogram. A classical approach to tackle this problem has been proposed in the geostatistics literature [START_REF] David | Geostatistical ore reserve estimation[END_REF][START_REF] Journel | Mining geostatistics[END_REF][START_REF] Clark | Practical geostatistics[END_REF] and consists in fitting a parametric model of valid semi-variograms to a pointwise nonparametric semi-variogram estimator by minimizing a given distance between the nonparametric estimator and the semi-variograms at a finite number of lags. The reader is refered to [START_REF] Cressie | Statistics for spatial data[END_REF]Chapter 2] for further details on semi-variogram model fitting and to [START_REF] Lahiri | On asymptotic distribution and asymptotic efficiency of least squares estimators of spatial variogram parameters[END_REF] for least-squares methods.

State of the art on quadratic variations

In order to remedy the drawbacks of the MLE previously mentioned, we focus on an alternative estimation method using quadratic variations based on the observations of the process X at a triangular array of points (t j ) j , where j = 1, ..., n. This method is also an alternative to the existing methods mentioned in the previous paragraph. Quadratic variations have been first introduced by Levy in [START_REF] Lévy | Le mouvement brownien plan[END_REF] to quantify the oscillations of the Brownian motion. Then a first result on the quadratic variation of a Gaussian non-differentiable process is due to Baxter (see e.g. [START_REF] Baxter | A strong limit theorem for Gaussian processes[END_REF], [START_REF] Grenander | Abstract inference[END_REF]Chap. 5] and [START_REF] Gladyšev | A new limit theorem for stochastic processes with Gaussian increments[END_REF]) that ensures (under some conditions) the almost sure convergence (as n tends to infinity) of

n j=1 (X(t j ) -X(t j-1 ) 2 , (4)
where t j = j/n, for j = 1, ..., n (by convention t 0 = 0 and X 0 = 0). A generalization of the previous quadratic variations has been introduced in Guyon and Léon [START_REF] Guyon | Convergence en loi des H-variations d'un processus gaussien stationnaire sur R[END_REF]: for a given real function H, the H-variation is given by

V H,n := n j=1 H X(t j ) -X(t j-1 ) (Var(X(t j ) -X(t j-1 )) 1/2 .
(5)

In [START_REF] Guyon | Convergence en loi des H-variations d'un processus gaussien stationnaire sur R[END_REF], technical conditions are assumed and a smoothness parameter 0 < s < 2, similar to the one in (2) when D = 0, is considered. Then the most unexpected result of [START_REF] Guyon | Convergence en loi des H-variations d'un processus gaussien stationnaire sur R[END_REF] is that (V H,n /n) n has a limiting normal distribution with convergence rate n 1/2 when 0 < s < 3/2 whereas the limiting distribution is non normal and the convergence rate is reduced to n 2-s when 3/2 < s < 2. Moreover, for statistical purposes, it has been proved by Coeurjolly that quadratic variations are optimal (details and precisions can be found in [START_REF] Coeurjolly | Estimating the parameters of a fractional Brownian motion by discrete variations of its sample paths[END_REF]). In [START_REF] Istas | Quadratic variations and estimation of the local Hölder index of a Gaussian process[END_REF], Istas and Lang generalized the results on quadratic variations. They allowed for observation points of the form t j = jδ n for j = 1, . . . , n, with δ n depending on n and tending to 0 as n goes to infinity. They studied the generalized quadratic variations defined by:

V a,n := n-1 i=1 k a k X(i + kδ n ) 2 , (6)
where the sequence a = (a k ) k has a finite support and some vanishing moments. Then they built estimators of the smoothness parameter and the scale parameter C and showed that these estimators are almost surely consistent and asymptotically normal. In the more recent work of Lang and Roueff [START_REF] Lang | Semi-parametric estimation of the Hölder exponent of a stationary Gaussian process with minimax rates[END_REF], the authors generalized the results of Istas and Lang [START_REF] Istas | Quadratic variations and estimation of the local Hölder index of a Gaussian process[END_REF] and Kent and Wood [START_REF] Kent | Estimating the fractal dimension of a locally self-similar Gaussian process by using increments[END_REF] on an increment-based estimator in a semi-parametric framework with different sets of hypothesis. Another generalization for non-stationary Gaussian processes and quadratic variations along curves is done in [START_REF] Adler | Uniform quadratic variation for Gaussian processes[END_REF]. See also the studies of [START_REF] Perrin | Quadratic variation for Gaussian processes and application to time deformation[END_REF] and [START_REF] Coeurjolly | Estimating the parameters of a fractional Brownian motion by discrete variations of its sample paths[END_REF].

Contributions of the paper Now let us present the framework considered in our paper. We assume that the Gaussian process X has stationary increments and is observed at times t j = jδ n for j = 1, . . . , n with δ n tending to zero. Note that t j also depends on n but we omit this dependence in the notation for simplicity. We will only consider δ n = n -α with 0 < α 1 throughout the article. Two cases are then considered: α = 1 (δ n = 1/n, infill asymptotic setting [START_REF] Cressie | Statistics for spatial data[END_REF], that we call the infill situation throughout) and 0 < α < 1 (δ n → 0 and nδ n → ∞, mixed asymptotic setting [START_REF] Cressie | Statistics for spatial data[END_REF], that we call the mixed situation throughout). The paper is devoted to the estimation of the scale parameter C from one or several generalized quadratic a-variations V a,n defined in [START_REF] Bachoc | Fixed-domain asymptotic properties of composite likelihood estimators for Gaussian processes[END_REF]. Calculations show that the expectation of V a,n is a function of C so that C can be estimated by the moment method.

Our study is related to the study of Istas and Lang [START_REF] Istas | Quadratic variations and estimation of the local Hölder index of a Gaussian process[END_REF] in which they estimate both the scale parameter C and the local Hölder index (a function of D and s in (2)). Our main motivation for focusing on the case where the local Hölder index is known is, on the one hand, to provide a simpler method to implement and analyze the estimator, and on the other hand to address more advanced statistical issues, such as efficiency and aggregation of several estimators of C. In addition, our results hold under milder technical conditions than in [START_REF] Istas | Quadratic variations and estimation of the local Hölder index of a Gaussian process[END_REF], and in particular apply to most semi-variogram models commonly used in practice. In particular, we also show that a necessary condition in [START_REF] Istas | Quadratic variations and estimation of the local Hölder index of a Gaussian process[END_REF], namely the fact that the quantity in ( 9) is non-zero when the variation used has a large enough order, in fact always holds. Thus, our study has a larger scope of application, in terms of necessary technical conditions, than that in [START_REF] Istas | Quadratic variations and estimation of the local Hölder index of a Gaussian process[END_REF].

We establish asymptotic approximations of the expectation and the variance and a central limit theorem for the quadratic variations under consideration and for the estimators deduced from them. In particular, given a finite number of sequences a, we prove a joint central limit theorem (see Corollary 3.8). In addition, our method does not require a parametric specification of the drift (see Section 3.4); therefore it is more robust than MLE.

For a finite discrete sequence a with zero sum, we define its order as the largest integer M such that

k a k k = 0 for = 1, . . . , M -1.
Roughly speaking k a k f (kδ n ) is an estimation of the M th derivative of the function f at zero. The order of the simplest sequence: -1, 1 is M = 1. Natural questions then arise. What is the optimal sequence a? In particular, what is the optimal order? Is it better to use the elementary sequence of order 1 (-1, 1) or the one of order 2 (-1, 2, -1)? For a given order, for example M = 1, is it better to use the elementary sequence of order 1 (-1, 1) or a more general one, for example (-1, -2, 3) or even a sequence based on discrete wavelets? Can we efficiently combine the information of several quadratic a-variations associated to several sequences? As far as we know, these questions are not addressed yet in the literature. Unfortunately, the asymptotic variance we give in Proposition 3.1 or Theorem 3.7 does not allow either to address theoretically this issue. However, by Corollary 3.8, one may gather the information of different quadratic a-variations with different orders. In order to validate such a procedure, an important Monte Carlo study is performed. The main conclusion is that gathering the information of different quadratic a-variations with different orders M produces closer results to the optimal Cramér-Rao bound computed in Section 4.2. The simulations are illustrated in Figure 4. We also illustrate numerically the convergence to the asymptotic distribution considering different models (exponential and Matérn models).

Finally, we show that our suggested quadratic variation estimator can be easily extended to the twodimensional case and we consider two real data sets in dimension two. When comparing our suggested estimator with maximum likelihood estimation, we observe a very significant computational benefit for our estimator.

Organization of the paper The paper is organized as follows. In Section 2, we detail the framework and present the assumptions on the process. In Section 3, we introduce our quadratic variation estimator and provide its asymptotic properties, together with discussion. Section 4 is devoted to the analysis of the statistical efficiency of our estimator. In Section 5, we provide the results of the Monte Carlo simulation and on the real data sets. A conclusion is provided in Section 6 together with some perspectives. All the proofs have been postponed to the Appendix.

2 General setting and assumptions

Assumptions on the process

In this paper, we consider a Gaussian process (X(t)) t∈R which is not necessarily stationary but only has stationary increments. The process is observed at times jδ n for j = 0, . . . , n with δ n going to 0 as n goes to infinity. As mentioned in the introduction, we will only consider δ n = n -α with 0 < α 1 throughout the article. Two cases are then considered: α = 1 (δ n = 1/n, infill situation) and 0 < α < 1 (δ n → 0 and nδ n → ∞, mixed situation). The semi-variogram of X is defined by

V (h) := 1 2 E (X(t + h) -X(t)) 2 .
In the sequel, we denote by (Const) a positive constant which value may change from one occurrence to another. For the moment, we assume that X is centered, the case of non-zero expectation will be considered in Section 3.4. Now, we introduce the following assumptions. The form of (H 1 ) and (H 2 ) change following whether we are in infill situation or in the particular mixed situation (δ n = n -α with 0 < α < 1).

(H 0 ) V is a C ∞ function on (0, +∞].
Infill situation:

δ n = 1/n. (H 1 )
The semi-variogram is 2D times differentiable with D 0 and there exists C > 0 and 0 < s < 2 such that for any h ∈ R, we have

V (2D) (h) = V (2D) (0) + C(-1) D |h| s + r(h), with r(h) = o(|h| s ) as |h| → 0. (7)
In (H 1 ), the integer D is the greatest integer such that V is 2D-times differentiable everywhere. We recall that, when X is assumed to be a stationary process, we have

V (h) = k(0) -k(h).
If the covariance function k belongs to a parametric set of the form {k θ ; θ ∈ Θ} with Θ ⊂ R p , then C is a deterministic function of the parameter θ.

(H 2 ) For some β < -1/2, for |h| < 1:

• when s < 3/2, r (2) (h) (Const) |h| β ;
• when s 3/2, r

(h) (Const) |h| β . (H 3 ) As h → 0, |r(h)| = o |h| s+1/2 . (3)
Mixed situation : δ n = n -α with 0 < α < 1.

We must add to (H 1 ):

|r(h)| (Const) |h| s |h| > 1.
The new expression of (H 2 ) is

• when s < 3/2, there exists β with s -2 < β < -1/2 such that, for all h ∈ R, r (2) (h) (Const) |h| β ;
• when s 3/2, there exists β with s -3 < β < -1/2 such that, for all h ∈ R,

r (3) (h) (Const) |h| β .
Here

(H 3 ) writes |r(h)| = o |h| s+(1/2α) ,
as h → 0.

Remark 2.1.

• When D > 0, the D-th derivative X (D) in quadratic mean of X is a Gaussian stationary process with autocovariance function k given by k(h) = (-1) D+1 V (2D) (h). This implies that the Hölder exponent of the paths of X (D) is s/2. Because s < 2, D is exactly the order of differentiation of the paths of X.

• If we denote H = D + s/2, H represents the local Hölder index of the process [START_REF] Ibragimov | Gaussian Random Processes[END_REF].

• Note that in the infill situation (δ n = 1/n), (H 2 ) is almost minimal. Indeed, the condition β < -1/2 does not matter since the smaller β, the weaker the condition. And for example, when s < 3/2, the second derivative of the main term is of order |h| s-2 and we only assume that β > s -2.

Examples of processes that satisfy our assumptions

We present a non exhaustive list of examples in dimension one that satisfy our hypotheses. In these examples, we provide a stationary covariance function k, and we recall that this defines V , with

V (h) = k(0) -k(h).
• The exponential model:

k(h) = exp(-C|h|) (D = 0, s = 1, C = C). For this model, (H 0 ) to (H 2 )
always hold and (H 3 ) holds when α > 1/2, that is when the observation domain does not increase too fast.

• The generalized exponential model:

k(h) = exp(-C|h| s ), s ∈ (0, 2) (D = 0, s = s, C = C).
For this model, (H 0 ) to (H 2 ) always hold and (H 3 ) holds when 1/(2α) < s. Hence, in the infill situation, we need s > 1/2 and, in the mixed situation, the observation domain needs to increase slowly enough.

• The generalized Slepian model [START_REF] Slepian | On zeros of Gaussian noise[END_REF]:

k(h) = (1 -C|h| s ) + , s ∈ (0, 1] (D = 0, s = s, C = C).
For this model, (H 0 ) to (H 3 ) hold in the infill situation and when C < 1. For (H 0 ), we remark that, in this case, V is smooth on (0, 1] and not on (0, ∞), but this is sufficient for all the results to hold.

• The Matérn model:

k(h) = 2 1-ν Γ(ν) √ 2νθh ν K ν ( √ 2νθh),
where ν > 0 is the regularity parameter of the process. The function K ν is the modified Bessel function of the second kind of order ν. See, e.g., [START_REF] Stein | Interpolation of Spatial Data: Some Theory for Kriging[END_REF] for more details on the model. In that case, D = ν and s = 2ν -2D. Here, it requires tedious computations to express the scale parameter C as a function of ν and θ. However, in Section 5.1, we derive the value of C in two settings (ν = 3/2 and ν = 5/2). For this model, (H 0 ) to (H 2 ) always hold and (H 3 ) holds when s < 2 -1/(2α).

Hence, in the infill situation we need s < 3/2, and in the mixed situation the observation domain needs to increase slowly enough.

All the previous examples are stationary (and thus have stationary increments). The following one is not stationary.

• The fractional Brownian motion (FBM) process denoted by (B s (t)) t∈R and defined by

Cov(B s (u), B s (t)) = C |u| s + |t| s -|u -t| s .
A reference on this subject is [START_REF] Cohen | Fractional fields and applications[END_REF]. This process is classically indexed by its Hurst parameter H = s/2. Here, D = 0, s = s and C = C. We call the FBM defined by C = 1 the standard FBM.

We remark that the Gaussian model, or square-exponential, defined by k(h) = σ 2 e -h 2 θ 2 , with (σ 2 , θ) ∈ (0, ∞), does not satisfy our assumptions, because it is too regular (i.e. it is C ∞ everywhere).

Detailed verification of the assumptions with the generalized exponential model We consider the generalized exponential model, where k(h) = exp(-C|h| s ) for some fixed s ∈ (0, 2). Since we have

V (h) = 1 -exp(-C|h| s ) for h ∈ R, (H 0 ) is trivially satisfied. Now we show that (H 1 ) holds for D = 0.
Indeed, V is a continuous function and we have

V (h) = 0 + C|h| s + 1 -C|h| s -exp(-C|h| s ) := V (0) + C(-1) 0 |h| s + r(h), with r(h) = 1 -C|h| s -exp(-C|h| s ). As h → 0, r(h) = O(|h| 2s ) = o(|h| s ) and thus (H 1 ) holds in the infill situation. Furthermore, as |h| → ∞, r(h) = -C|h| s + o(1) = O(|h| s
) and so (H 1 ) holds also in the mixed situation. Let us now show that (H 2 ) is also satisfied. First consider the case where s < 3/2 and let

β = 1 2 ((s -2) + (-1/2)) if 2s -2 > -1/2 (s > 3/4), 1 2 ((s -2) + (2s -2)) if 2s -2 -1/2 (s 3/4).
One can check that s -2 < β and that β < -1/2 since s -2 < -1/2. Consider the case where |h| 1. Since

r(h) = ∞ k=2 (-1) k+1 C k k! |h| sk , one has r (2) (h) = |h| 2s-2 ∞ k=2 (-1) k+1 C k k! sk(sk -1)|h| s(k-2) = |h| 2s-2 g(h)
where g is a bounded function on

[-1, 1]. Hence, |r (2) (h)| (Const)|h| β since 2s -2 β.
Then the case s < 3/2 is complete in the infill situation. Consider now the case where |h| 1. Simply, one can show that

r (2) (h) = -Cs(s -1)|h| s-2 + -Cs(s -1)|h| s-2 + C 2 s 2 |h| 2s-2 exp(-C|h| s ). Hence |r (2) (h)| (Const)|h| s-2 (Const)|h| β since β > s -2.
The case where s 3/2 can be treated analogously. Finally, it is simple to show that (H 3 ) holds when 1/(2α) < s.

Discrete a-differences

Now, we consider a non-zero finite support sequence a of real numbers with zero sum. Let L(a) be its length. Since the starting point of the sequence plays no particular role, we will assume when possible that the first non-zero element is a 0 . Hence, the last non-zero element is a L(a)-1 . We define the order M (a) of the sequence as the first non-zero moment of the sequence a:

L(a)-1 j=0 a j j k = 0, for 0 k < M (a)
and

L(a)-1 j=0 a j j M (a) = 0.
To any sequence a, with length L(a) and any function f , we associate the discrete a-difference of f defined by

∆ a,i (f ) = L(a)-1 j=0 a j f ((i + j)δ n ), i = 1, . . . n , (8)
where n stands for n -L(a) + 1. As a matter of fact, in the case of the simple quadratic a-variation given by a 0 = 0 and a 1 = -1, the operator ∆ a is a discrete differentiation operator of order one. More generally,

L(a)-1 j=0
a j f (jδ n ) is an approximation (up to some multiplicative coefficient) of the M (a)-th derivative (when it exists) of the function f at zero. We also define ∆ a (X) as the Gaussian vector of size n with entries ∆ a,i (X) and Σ a its variance-covariance matrix.

Examples -Elementary sequences. The simplest case is the order 1 elementary sequence a (1) defined by a

(1) 0 = -1 and a

(1)

1 = 1 We have L(a (1) ) = 2, M (a (1) ) = 1.
More generally, we define the k-th order elementary sequence a (k) as the sequence with coefficients a

(k) j = (-1) k-j k j , j = 0, . . . , k. Its length is given by L(a (k) ) = k + 1.
For two sequences a and a , we define their convolution b = a * a as the sequence given by b j = k-l=j a k a l . In particular, we denote by a 2 * the convolution a * a. Notice that the first non-zero element of b is not necessarily b 0 but b L(a )-1 as mentioned in the following properties.

Properties 2.2. The following properties of convolution of sequences are direct.

(i) The support of a * a (the indices of the non-zero elements) is included in -(L(a ) -1), (L(a) -1) while its order is M (a) + M (a ). In particular, a 2 * has length 2L(a) -1, order 2M (a) and is symmetrical.

(ii) The composition of two elementary sequences gives another elementary sequence.

The main result of this section is Proposition 2.6 that is required to quantify the asymptotic behaviors of the two first moments of the quadratic a-variations defined in [START_REF] Clark | Practical geostatistics[END_REF] (see Proposition 3.1). In order to prove (9), we establish two preliminary tools (Proposition 2.4 and Lemma 2.5). In that view, we need to define the integrated fractional Brownian motion (IFBM). We start from the FBM defined in Section 2.2 which has the following non anticipative representation:

B s (u) = u -∞ f s (t, u)dW (t),
where dW (t) is a white noise defined on the whole real line and

f s (t, u) = (Const) ((u -t) + ) (s-1)/2 -((-t) + ) (s-1)/2 .
For m 0 and t 0, we define inductively the IFBM by

B (-0) s (u) = B s (u) B (-m) s (u) = u 0 B (-(m-1)) s (t)dt.

Definition 2.3 (Non degenerated property).

A process Z has the ND property if for every k > 0 and every

t 1 < t 2 < • • • < t k belonging to the domain of definition of Z, the distribution of Z(t 1 ), . . . , Z(t k ) is non degenerated.
We have the following results.

Proposition 2.4. The IFBM has the ND property.

Lemma 2.5. The variance function of the IFBM satisfies, for all m ∈ N,

Var B (-m) s (u) -B (-m) s (v) = Nm i=1 P m,i (v)h m,i (u) + P m,i (u)h m,i (v) + (-1) m 2|u -v| s+2m (s + 1) . . . (s + 2m) ,
where N m ∈ N, for i = 1, ..., N m , P m,i is a polynomial of degree less or equal to m and h m,i is some function.

Proposition 2.6. If the sequence a has order M (a) > D, then

j a 2 * j |j| 2D+s = 0 (i.e. (-1) D j a 2 * j |i| 2D+s < 0). (9)
Note that ( 9) is stated as an hypothesis in [START_REF] Istas | Quadratic variations and estimation of the local Hölder index of a Gaussian process[END_REF].

3 Quadratic a-variations

Definition

Here, we consider the discrete a-difference applied to the process X and we define the quadratic avariations by

V a,n = ∆ a (X) 2 = n i=1 (∆ a,i (X)) 2 , ( 10
)
recalling that n = n -L(a) + 1. When no confusion is possible, we will use the shorthand notation L and M for L(a) and M (a).

Main results on quadratic a-variations

The basis of our computations of variances is the identity

E[∆ a,i (X)∆ a ,i (X)] = -∆ a * a ,i-i (V ), (11)
for any sequences a and a . A second main tool is the Taylor expansion with integral remainder (see, for example, [START_REF] Journel | Mining geostatistics[END_REF]). So we introduce another notation. For a sequence a, a scale δ, an order q and a function f , we define

R(i, δ, q, f, a) = - j a j j q 1 0 (1 -η) q-1 (q -1)! f ((i + jη)δ)dη. (12)
By convention, we let

R(i, δ, 0, f, a) = -∆ a,i (f ). Note that R(-i, δ, 2q, |•| s , a * a ) = R(i, δ, 2q, |•| s , a * a)
. One of our main results is the following.

Proposition 3.1 (Moments of V a,n ). Assume that V satisfies (H 0 ) and (H 1 ). 1) If we choose a sequence a such that M > D, then E[V a,n ] = nC(-1) D δ 2D+s n R(0, 1, 2D, |•| s , a 2 * ) (1 + o(1)), (13)
as n tends to infinity. Furthermore, (-1

) D R(0, 1, 2D, |•| s , a 2 * ) is positive.
2) If V satisfies additionally (H 2 ) and if we choose a sequence a so that M > D + s/2 + 1/4, then as n tends to infinity:

Var(V a,n ) = 2nC 2 δ 4D+2s n i∈Z R 2 (i, 1, 2D, |•| s , a 2 * )(1 + o(1)) (14)
and the series above is positive and finite.

Remark 3.2. (i) Notice that ( 13) and ( 14) imply concentration in the sense that

V a,n E[V a,n ] L 2 -→ n→+∞ 1.
(ii) In practice, since the parameters D and s are known, it suffices to choose M such that M D + 1 when s < 3/2 and M D + 2 when 3/2 s < 2.

(iii) The expression of the asymptotic variance appears to be complicated. Anyway, in practice, it can be easily approximated. Some explicit examples are given in Section 5.

Following the same lines as in the proof of Proposition 3.1 and using the identities (a * a ) j = (a * a) -j and R(i,

1, 2D, |•| s , a * a ) = R(-i, 1, 2D, |•| s , a * a)
, one may easily derive the corollary below. The proof is omitted. Corollary 3.3 (Covariance of V a,n and V a ,n ). Assume that V satisfies (H 0 ), (H 1 ), and (H 2 ). Let us consider two sequences a and a so that M (a) ∧ M (a ) > D + s/2 + 1/4. Then, as n tends to infinity, one has

Cov(V a,n , V a ,n ) = 2nC 2 δ 4D+2s n i∈Z R 2 (i, 1, 2D, |•| s , a * a ) (1 + o(1)). ( 15
)
Particular case -D = 0:

(i) We choose a as the first order elementary sequence (a 0 = -1, a 1 = 1 and M = 1). As n tends to infinity, one has

E[V a,n ] = nCδ s n (2 + o(1)); Var(V a,n ) = 2nC 2 δ 2s n i∈Z (|i -1| s -2 |i| s + |i + 1| s ) 2 (1 + o(1)), s < 3/2.
(ii) General sequences. We choose two sequences a and a so that M (a) ∧ M (a ) > s/2 + 1/4. Then, as n tends to infinity, one has

E[V a,n ] = -nCδ s n   j a 2 * j |j| s   (1 + o(1)); Var(V a,n ) = 2nC 2 δ 2s n i∈Z   j a 2 * j |i + j| s   2 (1 + o(1)); Cov(V a,n , V a ,n ) = 2nC 2 δ 2s n   |j| L a * a j |j| s   2 (1 + o(1)) + nC 2 δ 2s n i∈Z *      |j| L a * a j |i + j| s   2 +   |j| L a * a j |i + j| s   2    (1 + o(1)).
Now we establish the central limit theorem.

Theorem 3.4 (Central limit theorem for V a,n ). Assume (H 0 ), (H 1 ) and (H 2 ) and M > D + s/2 + 1/4. Then V a,n is asymptotically normal in the sense that

V a,n -E[V a,n ] Var(V a,n ) D -→ n→+∞ N (0, 1). ( 16
)
Remark 3.5.

• If M = D + 1, the condition M > D + s/2 + 1/4 in Proposition 3.1 implies s < 3/2.
However, when M = D + 1 and s 3/2, it is still possible to compute the variance but the convergence is slower and the central limit theorem does not hold anymore. More precisely, we have the following.

-If s > 3/2 and M = D + 1 then, as n tends to infinity,

Var(V a,n ) = (Const) × δ 4D+2s n × n 2s-4(M -D)+2 × (1 + o(1)). (17)
-If s = 3/2 and M = D + 1 then, as n tends to infinity

Var(V a,n ) = (Const) × δ 4D+2s n × n log n × (1 + o(1)). (18)
We omit the proof. Analogous formula for the covariance of two variations can be derived similarly.

• Since the work of Guyon and León [START_REF] Guyon | Convergence en loi des H-variations d'un processus gaussien stationnaire sur R[END_REF], it is a well known fact that in the simplest case (D = 0, L = 2, M = 1) and in the infill situation (δ n = 1/n, α = 1), the central limit theorem holds true for quadratic variations if and only if s < 3/2. Hence assumption M > D + s/2 + 1/4 is minimal.

Corollary 3.6 (Joint central limit theorem). Assume that V satisfies (H 0 ), (H 1 ) and (H 2 ). Let a (1) , . . . , a (k) be k sequences with order greater than D + s/2 + 1/4. Assume also that, as n → ∞, the k × k matrix with term i, j equal to

1 nδ 4D+2s n Cov V a (i) ,n , V a (j) ,n converges to an invertible matrix Λ ∞ . Then, V a (1) ,...,a (k) ,n = (V a (1) ,n , . . . , V a (k) ,n ) is asymptotically normal in the sense that n → ∞ V a (1) ,...,a (k) ,n -E V a (1) ,...,a (k) ,n n 1/2 δ 2D+s n D -→ n→+∞ N (0, Λ ∞ ).

Estimators of C based on the quadratic a-variations

Guided by the moment method, we define

C a,n := V a,n n(-1) D δ 2D+s n R(0, 1, 2D, |•| s , a 2 * ) . ( 19
)
Then C a,n is an estimator of C which is asymptotically unbiased by Proposition 3.1. Now our aim is to establish its asymptotic behavior.

Theorem 3.7 (Central limit theorem for C a,n ). Assume (H 0 ) to (H 3 ) and that M (a) > D + s/2 + 1/4. Then C a,n is asymptotically normal. More precisely, we have

C a,n -C Var(C a,n ) D -→ n→+∞ N (0, 1), (20)
with

Var(C a,n ) = (Const)n -1 (1 + o(1)).
The following corollary is of particular interest: it will give theoretical results when one aggregates the information of different quadratic a-variations with different orders. As one can see numerically in Section 5.2, such a procedure appears to be really promising and circumvents the problem of the determination of the optimal sequence a.

Corollary 3.8. Under the assumptions of Theorem 3.7, consider k sequences a (1) , . . . , a (k) so that, for i = 1, . . . , k, M (a (i) ) > D + s/2 + 1/4. Assume furthermore that the covariance matrix of

(C a (i) ,n /Var(C a (i) ,n ) 1/2 ) i=1,...,k converges to an invertible matrix Γ ∞ as n → ∞. Then, ([C a (i) ,n - C]/Var(C a (i) ,n ) 1/2 ) i=1,...,k converges in distribution to the N (0, Γ ∞ ) distribution.

Adding a drift

In this section, we do not assume anymore that the process X is centered and we set for t 0,

f (t) = E[X(t)].
We write X the corresponding centered process: X(t) = X(t) -f (t). As it is always the case in statistical applications, we assume that f is a C ∞ function. We emphasize on the fact that our purpose is not proposing an estimation of the mean function.

Corollary 3.9. Assume the same assumptions as in Theorem 3.7, and recall that δ n = n -α for α ∈ (0, 1].

Define K α M,n = sup t∈[0,n 1-α ] |f (M ) (t)|.
and if we assume in addition that

K α M,n = o(n -1/4 δ D-M +s/2 n ), (21)
then (20) still holds for X.

Note that in the infill situation (δ n = 1/n, α = 1), K 1 M,n does not depend on n. Obviously, ( 21) is met if f is a polynomial up to an appropriate choice of the sequence a (and M ). In the infill situation, a sufficient condition for ( 21) is M > D + s/2 + 1/4 which is always true. Moreover, it is worth noticing that we only assume regularity on the M -th derivative of the drift. No parametric assumption on the model is required, unlike in the MLE procedure.

Elements of comparison with existing procedures

Quadratic variations versus MLE

In this section, we compare our methodology to the very popular MLE method. For details on the MLE procedure, the reader is referred to, e.g. [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF][START_REF] Santner | The design and analysis of computer experiments[END_REF].

Model flexibility As mentioned in the introduction, the MLE methodology is a parametric method and requires the covariance function to belong to a parametric family of the form {k θ , θ ∈ Θ}. In the procedure proposed in this paper, it is only assumed that the semi-variogram satisfies the conditions given in Section 2.1, and that D and s are known. In this latter case, the suggested variation estimator is feasible, while the MLE is not defined.

Adding a drift In order to use the MLE estimator, it is necessary to assume that the mean function of the process is a linear combination of known parametric functions:

f (t) = q i=1 β i f i (t),
with known f 1 , . . . , f q and where β 1 , . . . , β q need to be estimated. Our method is less restrictive and more robust. Indeed, we only assume the regularity of the M -th derivative of the mean function in assumption [START_REF] Gladyšev | A new limit theorem for stochastic processes with Gaussian increments[END_REF]. It does not require parametric assumptions neither on the semi-variogram nor the mean function. Furthermore, no estimation of the mean function is necessary.

Computational cost

The cost of our method is only O(n) (the method only requires the computation of a sum) while the cost of the MLE procedure is known to be O(n 3 ).

Practical issues In some real data frameworks, it may occur that the MLE estimation diverges as can be seen in Section 5.3. Such a dead end can not be possible with our procedure.

Quadratic variations versus other methods

Least-square estimators In [START_REF] Lahiri | On asymptotic distribution and asymptotic efficiency of least squares estimators of spatial variogram parameters[END_REF], the authors propose a Least Square Estimator (LSE). More precisely, given an estimator V n of the semi-variagram assumed to belong to a parametric family of semi-variograms {V θ , θ ∈ Θ}, the Ordinary Least Square Estimation (OLSE) consists in minimizing in the parameter θ the quantity

h∈I (V n (h) -V θ (h)) 2 ,
where I is a set of lags and V n is a non-parametric estimator of the semi-variogram. Several variants as the Weighted Least Square Estimation and the General Least Square Estimation have been then introduced. Then the authors of [START_REF] Lahiri | On asymptotic distribution and asymptotic efficiency of least squares estimators of spatial variogram parameters[END_REF] provide necessary and sufficient conditions for these estimators to be asymptotically efficient and they show that when the number of lags used to define the estimators is chosen to be equal to the number of variogram parameters to be estimated, the ordinary least squares estimator, the weighted least squares and the generalized least squares estimators are all asymptotically efficient. Similarly as for the MLE, least square estimators require a parametric family of variograms, while quadratic variation estimators do not. [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of Gaussian processes[END_REF][START_REF] Bachoc | Asymptotic analysis of covariance parameter estimation for gaussian processes in the misspecified case[END_REF][START_REF] Bachoc | Cross validation and maximum likelihood estimations of hyper-parameters of gaussian processes with model misspecification[END_REF][START_REF] Zhang | Kriging and cross-validation for massive spatial data[END_REF] are based on minimizing scores based on the leave one out prediction errors, with respect to covariance parameters θ, when a parametric family of semi-variograms {V θ , θ ∈ Θ} is considered. Hence, as the MLE, they require a parametric family of variograms. Furthermore, as the MLE, the computation cost is in O(n 3 ), while this cost is O(n) for quadratic variation estimators.

Cross validation estimators Cross validation estimators

Composite likelihood Maximum composite likelihood estimators follow the principle of the MLE, with the aim of reducing its computational cost [START_REF] Varin | An overview of composite likelihood methods[END_REF][START_REF] Mateu | Fitting negative spatial covariances to geothermal field temperatures in Nea Kessani (Greece)[END_REF][START_REF] Pardo-Igúzquiza | AMLE3D: A computer program for the inference of spatial covariance parameters by approximate maximum likelihood estimation[END_REF][START_REF] Stein | Approximating likelihoods for large spatial data sets[END_REF][START_REF] Vecchia | Estimation and model identification for continuous spatial processes[END_REF]. In this aim, they consist in optimizing, over θ, the sum, over i = 1, . . . , n of the conditional likelihoods of the observation X i = X(t i ), given a small number of observations which observation locations are close to t i when a parametric family of semi-variograms {V θ , θ ∈ Θ} is considered. The computation cost of an evaluation of this sum of conditional likelihood is O(n), in contrast to O(n 3 ) for the full likelihood. Nevertheless, the composite likelihood estimation requires to perform a numerical optimization, while our suggested estimator does not. Furthermore, a parametric family of variograms is required for the composite likelihood but not for our estimator. Finally, [START_REF] Bachoc | Fixed-domain asymptotic properties of composite likelihood estimators for Gaussian processes[END_REF] recently showed that the composite likelihood estimator has rate of convergence only n s when D = 0 and 0 < s < 1/2 in (2). Hence, the rate of convergence of quadratic variation estimators (n 1/2 ) is larger in this case.

Already known results on quadratic a-variations

In [START_REF] Istas | Quadratic variations and estimation of the local Hölder index of a Gaussian process[END_REF], Istas and Lang consider a Gaussian process with stationary increments in the infill case and assume 2 as in our paper. Then they establish the asymptotic behavior of V a,n under more restrictive hypothesis of regularity on V than ours (in particular on r and on δ n ). Then they propose an estimation of both the local Hölder index H = D + s/2 and the scale parameter C, based on quadratic a-variations and study their asymptotic behavior. The expression of the estimation of C is much more complex than our that simply stems from the moment method. More precisely, they consider I sequences (a j ) j=1,...,I with length p j and the vector U of length I whose coordinate j is given by V a j ,n . Noticing that the vector U/n converges to the product AZ where A is a I × p matrix derived from the sequences a with p = max j p j and Z is the vector of (C(-1) D (jδ n ) 2h ) j=1,...,p , they estimate Z by Ẑ = (A A) -1 A U and derive their estimators of h and C from Ẑ.

As explained in the introduction in Section 1, Lang and Roueff in [START_REF] Lang | Semi-parametric estimation of the Hölder exponent of a stationary Gaussian process with minimax rates[END_REF] generalize the results of Istas and Lang in [START_REF] Istas | Quadratic variations and estimation of the local Hölder index of a Gaussian process[END_REF] and [START_REF] Kent | Estimating the fractal dimension of a locally self-similar Gaussian process by using increments[END_REF]. They consider the infill situation and use quadratic a-variations to estimate both the scale parameter C and smoothness parameter s under a similar hypothesis as in 2. Furthermore, they assume three types of regularity assumptions on V : Hölder regularity of the derivatives at the origin, Besov regularity and global Hölder regularity. Nevertheless, estimating both C and s leads to a more complex estimator of C and to proofs significantly different and more complicated.

To summarize, our contributions, additionally to the existing references [START_REF] Lang | Semi-parametric estimation of the Hölder exponent of a stationary Gaussian process with minimax rates[END_REF][START_REF] Istas | Quadratic variations and estimation of the local Hölder index of a Gaussian process[END_REF][START_REF] Kent | Estimating the fractal dimension of a locally self-similar Gaussian process by using increments[END_REF], is to provide an estimation method for C which definition, implementation and asymptotic analysis are simpler. As a result, we need fewer technical assumptions. In fact, our assumptions can be easily shown to hold in many classical examples. This also enables us to study the aggregation of quadratic variation estimators from different sequences, see Section 5.2.

Efficiency of our estimation procedure

In this section, in order to decrease the asymptotic variance, we propose a procedure to combine several quadratic a-variations leading to aggregated estimators. Then our goal is to evaluate the quality of these proposed estimators. In that view, we compare their asymptotic variance with the theoretical Cramér-Rao bound in some particular cases in which this bound can be explicitly computed.

Aggregation of estimators

Now in order to improve the estimation procedure, we suggest to aggregate a finite number of estimators:

k j=1 λ j C a (j) ,n
based on k different sequences a (1) , ..., a (k) with weights λ 1 , . . . , λ k . Ideally, one should provide an adaptive statistical procedure to choose the optimal number k * of sequences, the optimal sequences and the optimal weights λ * . Such a task is beyond the scope of this paper. Nevertheless, in this section, we consider a given number k of given sequences a (1) , ..., a (k) leading to the estimators C a (1) ,n , . . . , C a (k) ,n defined by [START_REF] Deisenroth | Distributed Gaussian processes[END_REF]. Then we provide the optimal weights λ * . Using [START_REF] Lavancier | A general procedure to combine estimators[END_REF] or [START_REF] Bates | The combination of forecasts[END_REF], one can establish the following lemma. Lemma 4.1. We assume that for j = 1, ..., k, the conditions of Corollary 3.8 are met. Let R be the k × k asymptotic variance-covariance matrix of the vector of length k whose elements are given by (n 1/2 /C)C a (j) ,n , j = 1, . . . , k. Then for any λ 1 , . . . , λ k ,

(n 1/2 /C)( k j=1 λ j C a (j) ,n -C) D -→ n→+∞ N (0, λ T Rλ).
Let 1 k be the "all one" column vector of size k and define

λ * = R -1 1 k 1 T k R -1 1 k .
One has k j=1 λ * j = 1 and λ * T Rλ * λ T Rλ.

As will be shown with simulations in Section 5, the aggregated estimator considerably improves each of the original estimators C a (1) ,n , ..., C a (k) ,n . We call v a,s its normalized asymptotic variance.

Cramér-Rao bound

To validate the aggregation procedure, we want to compare the obtained asymptotic variance with the theoretical Cramér-Rao bound. In that view, we compute in the following section the Cramér-Rao bound in two particular cases.

We consider a family Y C (C ∈ R + ) of centered Gaussian processes. Let R C be the (n -1) × (n -1) variance-covariance matrix defined by

(R C ) i,j = Cov (Y C (iδ n ) -Y C ((i -1)δ n ) , Y C (iδ n ) -Y C ((i -1)δ n )) .
Assume that C → R C is twice differentiable and R C is invertible for all C ∈ R + . Then, let

I C = 1 2 Tr R -1 C ∂ ∂C R C R -1 C ∂ ∂C R C ( 22
)
be the Fisher information. The quantity 1/I C is the Cramér-Rao lower bound for estimating C based on

Y C (iδ n ) -Y C ((i -1)δ n ) i=2,...,n
(see for instance [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of Gaussian processes[END_REF][START_REF] Dahlhaus | Efficient parameter estimation for self-similar processes[END_REF]). Now we give two examples of families of processes for which we can compute the Cramér-Rao lower bound explicitly. The first example is obtained from the IFBM defined in Section 2.2.

Lemma 4.2. Let 0 < s < 2 and let X be equal to √ CB

(-D) s where B (-D) s is the IFBM. Then Y C = X (D) is a FBM whose semi-variogram V C is given by V C (h) = 1 2 E (Y C (t + h) -Y C (t)) 2 = C|h| s . ( 23
)
Hence in this case, we have 1/I C = 2C 2 /(n -1).

Proof. [START_REF] Guyon | Convergence en loi des H-variations d'un processus gaussien stationnaire sur R[END_REF] implies that ∂R C /∂C = R 1 then [START_REF] Grenander | Abstract inference[END_REF] gives the result. Now we consider a second example given by the generalized Slepian process defined in Section 2.2.

Let s 1 and Y C with stationary autocovariance function k C defined by

k C (h) = (1 -(C/2)|h| s ) + , for any h ∈ R. ( 24
)
This function is convex on R and it follows from Pólya's theorem [START_REF] Pólya | Remarks on characteristic functions[END_REF] that k C is a valid autocovariance function. We thus easily obtain the following lemma whose proof is omitted.

Lemma 4.3. Let X be the integration D times of Y C defined via [START_REF] Hall | On the nonparametric estimation of covariance functions[END_REF]. Then, in the infill situation (δ n = 1/n, α = 1) and for C < 2, the semi-variogram of Y C is given by [START_REF] Guyon | Convergence en loi des H-variations d'un processus gaussien stationnaire sur R[END_REF] and by consequence 1/I C = 2C 2 /(n -1).

Numerical results

In this section, we first study to which extent the asymptotic results of Proposition 3.1 and Theorem 3.7 are representative of the finite sample behavior of quadratic a-variations estimators. Then, we study the asymptotic variances of these estimators provided by Proposition 3.1 and that of the aggregated a-variations estimators of Section 4.1.

Simulation study of the convergence to the asymptotic distribution

We carry out a Monte Carlo study of the quadratic a-variations estimators in three different cases. In each of the three cases, we simulate N = 10, 000 realizations of a Gaussian process on [0, 1] with zero mean function and stationary autocovariance function k. In the case D = 0, we let k(h) = exp(-C|h|).

Hence (H 1 ) holds with D = 0 and s = 1. In the case D = 1, we use the Matérn 3/2 autocovariance [START_REF] Roustant | DiceOptim: Two R packages for the analysis of computer experiments by Kriging-based metamodeling and optimization[END_REF] :

k(h) = 1 + √ 3 |h| θ e - √ 3 |h| θ .
One can show, by developing k into power series, that (H 1 ) holds with D = 1, s = 1 and C = 6 √ 3/θ 3 . Finally, in the case D = 2, we use the Matérn 5/2 autocovariance function:

k(h) = 1 + √ 5 |h| θ + 5|h| 2 3θ 2 e - √ 5 |h| θ .
Also (H 1 ) holds true with D = 2, s = 1 and C = 200 √ 5/3θ 5 .

In each of the three cases, we set C = 3. For n = 50, n = 100 and n = 200, we observe each generated process at n equispaced observation points on [0, 1] and compute the quadratic a-variations estimator C a,n of Section 3.3. When D = i, i = 0, 1, 2, we choose a to be the elementary sequence of order i + 1.

In Figure 2, we display the histograms of the 10, 000 estimated values of C for the nine configurations of D and n. We also display the corresponding asymptotic Gaussian probability density functions provided by Proposition 3.1 and Theorem 3.7. We observe that there are few differences between the histograms and limit probability density functions between the cases (D = 0, 1, 2). In these three cases, the limiting Gaussian distribution is already a reasonable approximation when n = 50. This approximation then improves for n = 100 and becomes very accurate when n = 200. Naturally, we can also see the estimators' variances decrease as n increases. Finally, the figures suggest that the discrepancies between the finite sample and asymptotic distributions are slightly more pronounced with respect to the difference in mean values than to the difference in variances. As already pointed out, these discrepancies are mild in all the configurations.

Analysis of the asymptotic distributions

Now we consider the normalized asymptotic variance of C a,n obtained from ( 14) in Proposition 3.1. We consider the infill situation (δ n = 1/n, α = 1) and we let

v a,s = 2 i∈Z R 2 (i, 1, 2D, |•| s , a 2 * ) R 2 (0, 1, 2D, |•| s , a 2 * , (25)
so that (n 1/2 /C)(C a,n -C) converges to a N (0, v a,s ) distribution as n → ∞, where v a,s already defined in Section 4.1 does not depend on C (nor on n).

First, we consider the case D = 0 and we plot v a,s as a function of s for various sequences a in Figure 3. The considered sequences are the following:

• the elementary sequence of order 1: a (1) given by (-1,1);

• the elementary sequence of order 2: a (2) given by (1,-2,1);

• the elementary sequence of order 3: a (3) given by (-1, 3, -3, 1);

• the elementary sequence of order 4, a (4) given by (1,-4, 6,-4,1);

• a sequence of order 1 and with length 3: a (5) given by (-1,-2,3);

• a Daubechies wavelet sequence [START_REF] Daubechies | Orthonormal bases of compactly supported wavelets[END_REF] with M = 2 as in [START_REF] Istas | Quadratic variations and estimation of the local Hölder index of a Gaussian process[END_REF]: a (6) given by (-0.1830127,-0.3169873,1.1830127,-0.6830127);

• a second Daubechies wavelet sequence with M = 3: a (7) given by (0.0498175,0.12083221,-0.19093442,-0.650365,1.14111692,-0.47046721). From Figure 3, we can draw several conclusions. First, the results of Section 4.2 suggest that 2 is a plausible lower bound for v a,s . We shall call the value 2 the Cramér-Rao lower bound. Indeed, we observe numerically that v a,s 2 for all the s and a considered here. Then we observe that, for any value of s, there is one of the v a,s which is close to 2 (below 2.5). This suggests that quadratic variations can be approximately as efficient as maximum likelihood, for appropriate choices of the sequence a. We observe that, for s = 1, the elementary sequence of order 1 (a 0 = -1, a 1 = 1) satisfies v a,s = 2. This is natural since for s = 1, this quadratic a-variations estimator coincides with the maximum likelihood estimator, when the observations stem from the standard Brownian motion. Except from this case s = 1, we could not find other quadratic a-variations estimators reaching exactly the Cramér-Rao lower bound 2 for other values of s.

Second, we observe that the normalized asymptotic variance v a,s blows up for the two sequences a satisfying M = 1 when s reaches 1.5. This comes from Remark 3.5: the variance of the quadratic avariations estimators with M = 1 is of order larger than 1/n when s 1.5. Consequently, we plot v a,s for 0.1 s 1.4 for these two sequences. For the other sequences satisfying M 2, we plot v a,s for 0.1 s 1.9.

Third, it is difficult to extract clear conclusions about the choice of the sequence: for s smaller than, say, 1.2 the two sequences with order M = 1 have the smallest asymptotic variance. Similarly, the elementary sequence of order 2 has a smaller normalized variance than that of order 3 for all values of s. Also, the Daubechies sequence of order 2 has a smaller normalized variance than that of order 3 for all values of s. Hence, a conclusion of the study in Figure 3 is the following. When there is a sequence of a certain order for which the corresponding estimator reaches the rate 1/n for the variance, there is usually no benefit in using a sequence of larger order. Finally, the Daubechies sequences appear to yield smaller asymptotic variances than the elementary sequences (the orders being equal). The sequence of order 1 given by (a 0 , a 1 , a 2 ) = (-1, -2, 3) can yield a smaller or larger asymptotic variance than the elementary sequence of order 1, depending on the value of s. For two sequences of the same order M , it seems nevertheless challenging to explain why one of the two provides a smaller asymptotic variance. Now, we consider aggregated estimators, as presented in Section 4.1. A clear motivation for considering aggregation is that, in Figure 3, the smallest asymptotic variance v a,s corresponds to different sequences a, depending on the values of s. In Figure 4 left, we consider the case D = 0 and we use four sequences: a (1) , a (5) a (2) and a (6) . We plot their corresponding asymptotic variances v a (i) ,s as a function of s, for 0.1 s 1.4 as well as the variance of their aggregation. It is then clear that aggregation drastically improves each of the four original estimators. The asymptotic variance of the aggregated estimator is very close to the Cramér-Rao lower bound 2 for all the values of s. In Figure 4 right, we perform the same analysis but with sequences of order larger than 1. The four considered sequences are now a (6) , a (2) a (3) and a (4) . The value of s varies from 0.1 to 1.9 Again, the aggregation is clearly the best.

Eventually, Figures 5 and6 explore the case D = 1. Conclusions are similar.

Real data examples

In this section, we consider real data of spatially distributed processes in dimension two. In this setting, we extend the estimation procedure based on the quadratic a-variations that is then compared to the MLE procedure.

A moderate size data set

We compare two methods of estimation of the autocovariance function of a separable Gaussian model on a real data set of atomic force spectroscopy1 . The data consist of observations taken on a grid of step 1/15 on [0, 1] 2 , so they consist of 256 points of the form X(i/15, j/15) i = 0, . . . , 15, j = 0, . . . , 15.

The first method is maximum likelihood estimation in a Kriging model, obtained from the function km or the R toolbox DiceKriging [START_REF] Roustant | DiceOptim: Two R packages for the analysis of computer experiments by Kriging-based metamodeling and optimization[END_REF]. For this method, the mean and autocovariance functions are assumed 4 but for D = 1. On the left, from top to bottom in the legend, the sequences are the elementary sequence of order 2, the Daubechies sequence of order 2 and the elementary sequence of order 3. On the right, from top to bottom in the legend, the sequences are the elementary sequences of orders 3 and 4 and the Daubechies sequence of order 3. 0.8770 and θ2 = 0.6547 and takes about 22 minutes on a personal computer. The second method provides θ1 = 0.607 and θ2 = 0.107 and takes about 0.05 seconds on a personal computer. Hence, our suggested quadratic variation estimator provides a very significant computational benefit. Both estimators conclude that the spatial correlation is more important along the x-axis that along the y-axis, which is graphically confirmed in Figure 7. As in Section 5.3.1, the maximum likelihood estimator provides less correlation than the quadratic variation estimator. Finally, if the field of deformation is considered with no preliminary subsampling, its size is 400 × 600. In this case, the MLE can not be directly implemented while the quadratic variation estimator can be.

Conclusion

We have provided an in-depth analysis of the estimation of the scale parameter of a one-dimensional Gaussian process by quadratic variations. Indeed, the knowledge of this scale parameter is essential when studying a Gaussian process, as it enables to quantify its dependence structure, or to test independence. We have addressed a semi-parametric setting, where no parametric family of variograms needs to be assumed to contain the unknown variogram. We have suggested an estimator, based on previous references, which numerical implementation is straightforward. Our theoretical analysis follows the principles of previous references, but is significantly simpler and holds under mild and simple to check technical assumptions. Based on this theoretical analysis, we have been able to tackle more advanced statistical topics, such as the aggregation of estimators based on different sequences, in the aim of improving the statistical efficiency. Our analysis paves the way for further research topics. For instance, it would be interesting to estimate the variances and covariances of a set of quadratic variation estimators, in order to estimate the optimal aggregation of them.

A Proofs

A.1 Proof of the results of Section 2.3 Proof of Proposition 2.4. By the stochastic Fubini theorem,

B (-m) s (u 1 ) = u1 0 du 2 • • • um 0 du m+1 um+1 -∞ dW (t)f s (t, u m+1 ) = u1 -∞ dW (t) u1 t du 2 • • • um t du m+1 f s (t, u m+1 ) =: u1 0 g m,s (u 1 , t)dW (t).
The positiveness of f s (t, u) for u > 0 implies that of g m,s (t, u). As a consequence, for 0 < t 1 < . . . t k , B (-m) s (t k ) includes a non-zero component: so that the lemma holds with the convention (s + 1) . . . (s + 0) = 1. Thus we prove it by induction on m and assume that it holds for m ∈ N. We have, with

t k t k-1 g m,s (u, t)dW (t), which is independent of (B (-m) s (t 1 ), . . . , B (-m) s (t k-1 )) implying that B (-m) s (t k ) is not
K (-r) (u, v) = E B (-r) s (u)B (-r) s (v) , for r ∈ N, K (-m) (u, v) = 1 2 Var B (-m) s (u) -B (-m) s (0) + Var B (-m) s (v) -B (-m) s (0) -Var B (-m) s (u) -B (-m) s (v) =ψ(u) + ψ(v) - 1 2 Nm i=1 P m,i (v)h m,i (u) - 1 2 Nm i=1 P m,i (u)h m,i (v) - 1 2 (-1) m 2|u -v| s+2m (s + 1) . . . (s + 2m) ,
where ψ is some function. Since we have K

(-(m+1)) (u, v) = u 0 v 0 K (-m) (x, y)dxdy, K (-(m+1)) (u, v) = Nm+1 i=1 P m+1,i (v) h m+1,i (u) + Nm+1 i=1 P m+1,i (u) h m+1,i (v) + (-1) m+1 1 (s + 1) . . . (s + 2m) v 0 u 0 |x -y| s+2m dx dy, (27)
where N m+1 ∈ N, where for i = 1, ..., N m+1 , P m+1,i is a polynomial of degree less or equal to m + 1 and h m+1,i is some function. For v u, we have

v 0 u 0 |y -x| s+2m dx dy = v 0 y 0 (y -x) s+2m dx + u y (x -y) s+2m dx dy = v 0 y s+2m+1 2m + 1 + (u -y) s+2m+1 2m + 1 dy = v s+2m+2 (2m + 1)(2m + 2) - (u -v) s+2m+2 (2m + 1)(2m + 2) + u s+2m+2 (2m + 1)(2m + 2)
.

By symmetry, we obtain, for u, v ∈ N,

u 0 v 0 |x -y| s+2m dx dy = u s+2m+2 (2m + 1)(2m + 2) + v s+2m+2 (2m + 1)(2m + 2) - |u -v| s+2m+2 (2m + 1)(2m + 2) . (28)
Hence, from the relation We conclude using the ND property of the IFBM stated in Proposition 2.4.

Var B (-(m+1)) s (u) -B (-(m+1)) s (v) = K (-(m+1)) (v, v) + K (-(m+1)) (u, u) -2K (-(m+1)) (v, u), (27)

A.2 Preliminary results

Lemma A.1. Let Z = (X, Y ) be a centered Gaussian vector of dimension 2 then Cov X 2 , Y 2 = 2Cov 2 (X, Y ) .

Proof of Lemma A.1. This Lemma is a consequence of the so called Mehler formula [START_REF] Azaïs | Level sets and extrema of random processes and fields[END_REF]. Its proof is immediate using the cumulant method. 2) Using Lemma A.1, [START_REF] Journel | Mining geostatistics[END_REF] with q = 2D, the fact that D M , and the vanishing moments of the sequence a 2 * , we obtain Var(V a,n ) =2 n i,i =1

Cov 2 (∆ a,i (X), ∆ a,i (X)) = 2 Using [START_REF] Kaufman | Covariance tapering for likelihood-based estimation in large spatial data sets[END_REF] , with q = 2M , |•| s instead of V (2D) and δ n = 1, and using the vanishing moments of the sequence a 2 * , we get, for i large enough so that i and i + j always have the same sign in the sum below,

which is the general term of a convergent series.

(ii) Now we show that the term C n is negligible compared to A n . This will imply in turn that B n is negligible compared to A n , from the Cauchy-Schwarz inequality. We have to give bounds to the series with general term R 2 (i, δ n , 2D, r, a 2 * ) with R(i, δ n , 2D, r, a 2 * ) =j a 2 * j j 2D 1 0

(1 -η) 2D-1 (2D -1)! r ((i + jη)δ n ) dη.

For fixed i, the assumptions (7) on r in (H 1 ) are sufficient to build a dominated convergence argument to prove that R 2 (i, δ n , 2D, r, a 2 * ) = o(δ 2s n ) which leads to the required result. So we concentrate our attention on indices i such that |i| > 2L. Now we use (H 2 ) and the notation d = 2 if s < 3/2 and d = 3 if s 3/2. Consider β as in (H 2 ) and remark that, in the mixed situation, we have s -d < β < -1/2. In the infill situation, it is assumed that for |h| < 1, |r (d) |(h) (Const)|h| β with β < -1/2. Since h is restricted to [-1, 1], we may also consider without loss of generality that β has been chosen such that s -d < β < -1/2. Using [START_REF] Kaufman | Covariance tapering for likelihood-based estimation in large spatial data sets[END_REF] as in the proof of item 1), if 2D + d 2M , one gets R(i, δ n , 2D, r, a 2 * ) = - Proof of Theorem 3.4. By a diagonalization argument, V a,n can be written as

V a,n = n i=1 λ i Z 2 i ,
where λ 1 , . . . , λ n are the non-zero eigenvalues of variance-covariance matrix Σ a of ∆ a (X) and the Z i are independent and identically distributed standard Gaussian variables. Hence,

V a,n -E(V a,n ) Var(V a,n ) = n i=1 λ i n r=1 λ 2 r (Z 2 i -1). (35)
In such a situation, Lemma 2 in [START_REF] Istas | Quadratic variations and estimation of the local Hölder index of a Gaussian process[END_REF] implies that the Lindeberg condition is a sufficient condition required to prove the central limit theorem and is equivalent to Proof of Corollary 3.6. To prove the asymptotic joint normality it is sufficient to prove the asymptotic normality of any non-zero linear combination

LC(γ) = k j=1 γ j V a (j) ,n ,
where γ j ∈ R for j = 1, . . . , k. We have again the representation

LC(γ) = n i=1 λ i Z 2 i ,
where the λ i 's are now the non-zero eigenvalues of the variance-covariance matrix σ = k j=1 γ j Σ a (j) ,n , and the Z i 's are as before. The Lindeberg condition has the same expression. On one hand, as n goes to infinity,

1 nδ 4D+2s n n i=1 λ i → γ Λ ∞ γ

Figure 1 :

 1 Figure 1: Several realizations of Gaussian processes with exponential covariance function with parameter C = 0.1 (left), 10 (middle), and 100 (right).

Figure 2 :

 2 Figure 2: Comparison of the finite sample distribution of C a,n (histograms) with the asymptotic Gaussian distribution provided by Proposition 3.1 and Theorem 3.4 (probability density function in blue line). The vertical red line denotes the true value of C = 3. From left to right, n = 50, 100, 200. From top to bottom, D = 0, 1, 2.

Figure 3 :

 3 Figure3: Case D = 0. Plot of the normalized asymptotic variance v a,s of the quadratic a-variations estimator, as a function of s, for various sequences a. The legend shows the values a 0 , ..., a l of these sequences (rounded to two digits). From top to bottom in the legend, the sequences are the elementary sequence of order 1, the sequence (-1, -2, 3) which has order 1, the d Daubechies sequences of order 2 and 3 and the elementary sequences of orders 2 and 3. The horizontal line corresponds to the Cramér-Rao lower bound 2.

Figure 4 :

 4 Figure 4: Case D = 0. Plot of the normalized asymptotic variance v a,s of the quadratic a-variations estimator, as a function of s, for various sequences a and for their aggregation. On the left, including the order one elementary sequence, on the right without. The horizontal line corresponds to the Cramér-Rao lower bound 2.

  18,-0.32,1.18,-0.68 0.05,0.12,-0.19,-0.65,1.14,-0.47

Figure 5 :

 5 Figure 5: Same setting as in Figure 3 but for D = 1. From top to bottom in the legend, the sequences are the elementary sequences of order 2, 3 and 4 and the Daubechies sequences of order 2 and 3.

Figure 6 :

 6 Figure6: Same setting as in Figure4but for D = 1. On the left, from top to bottom in the legend, the sequences are the elementary sequence of order 2, the Daubechies sequence of order 2 and the elementary sequence of order 3. On the right, from top to bottom in the legend, the sequences are the elementary sequences of orders 3 and 4 and the Daubechies sequence of order 3.

  collinear to this set of variables. By induction, this implies in turn that B (-m) s (t 1 ), . . . , B (-m) s (t k ) are not collinear. Proof of Lemma 2.5. For m = 0, we have Var B (-0) s (u) -B (-0) s (v) = 2|u -v| s

  , and[START_REF] Hinton | Training products of experts by minimizing contrastive divergence[END_REF], we conclude the proof of the lemma.Proof of Proposition 2.6. Using Lemma 2.5 (with m = D) and the vanishing moments of a of order less or equal than D, we havek,l a k a l |k -l| 2D+s = (Const)(-1) D k,l a k a l Var(B (-D) s (k) -B (-D) s (l)) = (Const)(-1) D+1 Var k a k B (-D) s (k) .

Lemma A. 2 .

 2 Assume that V satisfies (H 0 ), (H 1 ) and (H 2 ). One has, whenM > D + s + 1/4, max i=1,...,n   i =1,...,n |Σ a (i, i )|   = o Var(V a,n ) 1/2 .

  |i|)R 2 (i, δ n , 2D, V (2D) , a 2 * ) |i|) C(-1) D δ s n R(i, 1, 2D, |•| s , a 2 * ) + R(i, δ n , 2D, r, a 2 * ) 2 . =:A n + B n + C n ,where B n comes from the double product.(i) We show that A n converges. Indeed,A n = 2C 2 δ 4D n n -1 i=-n +1 (n -|i|)δ 2s n R 2 (i, 1, 2D, |•| s , a 2 * ) = 2C 2 n δ 4D+2s n i∈Z f n (i), with f n (i) := n -|i| n R 2 (i, 1, 2D, |•| s , a 2 * )1 |i| n -1 . Since f n (i) ↑ R 2 (i, 1, 2D, |•| s , a2* ) for fixed i and n going to infinity, it suffices to study the series i∈Z R 2 (i, 1, 2D, |•| s , a 2 * ).

R(i, 1 , 2 * j j 2M 1 0( 1 -

 1211 2D, |•| s , a 2 * ) = R(i, 1, 2M, g, a 2 * ) =j a η) 2M -1 (2M -1)! g((i + jη))dη,where g is the 2(M -D)-th derivative of |•| s (defined on R \ {0}). For i sufficiently large, g(i + jη) is bounded by (Const)|i| s-2(M -D) so thatR 2 (i, 1, 2D, |•| s , a2* ) is bounded by (Const)i 2(s-2(M -D)) ,

j a 2 * j j 2D+d δ d n 1 0( 1

 211 -η) 2D+d-1 (2D + d -1)! r (d) ((i + jη)δ n ) dη.The condition |i| > 2L ensures that the integral is always convergent. Then we haveR 2 (i, δ n , 2D, r, a 2 * ) (Const)δ 2d+2β n i 2β .(34)Since β < -1/2 , the series in i converges and the contribution to C of the indices i such that |i| > 2L is bounded by (Const)δ 4D+2d+2β n which is negligible compared to δ 4D+2s n since d + β > s.

  max i=1,...,n |λ i | = o Var(V a,n ) .From Lemma A.2, one has max follows using the following classical linear algebra result (see for instance[START_REF] Luenberger | Introduction to Dynamic Systems: Theory, Models, and Applications[END_REF] Ch. 6

Personal communication from C. Gales and J. M. Senard.

Acknowledgments This work has been partially supported by the French National Research Agency (ANR) through project PEPITO (no ANR-14-CE23-0011). The authors are grateful to Laurent Risser, for providing them the image registration data of Section 5.3.2.

to be E(X(i/15, j/15)) = µ and Cov(X(i/15, j/15), X(i /15, j /15)) = σ 2 e -θ1|i-i |/15 e -θ2|j-j |/15 .

(

The parameters µ, σ 2 , θ 1 , θ 2 are estimated by maximum likelihood.

The second method assumes the same autocovariance model [START_REF] Han | Financial time series volatility analysis using gaussian process state-space models[END_REF] and consists in the following steps.

(1) Estimate σ 2 the sum of square

with m = 15 i,j=0 X(i/15, j/15).

(2) For each column j of [X(i/15, j/15)] i,j=0,...,15 , the vector of 16 observations obey our model with s = 1 and C 1 = σ 2 θ 1 . Hence, we can estimate C 1 by Ĉ1,j with the estimator [START_REF] Deisenroth | Distributed Gaussian processes[END_REF], with the elementary sequence of order 1. We thus obtain an estimate Ĉ1 by averaging the Ĉ1,j for j = 0, . . . , 15.

(3) We perform the same analysis row by row to obtain an estimate Ĉ2 .

(

The first method, based on maximum likelihood, provides infinite values for θ 1 and θ 2 , so that it considers the 256 observed values as completely spatially independent. On the other hand, the second method provides the values θ1 = 14.72 and θ2 = 15.73. This corresponds to a correlation of approximately 1/e ≈ 0.36 between direct neighbors on the grid. Hence, the second method, based on our suggested quadratic variation estimator, is able to detect a weak correlation (that can be checked graphically), but not the maximum likelihood estimator.

A large size data set

The second data set consists in a two-dimensional field of deformation amplitude, corresponding to the registration of two real images. The deformation field is obtained from the software presented in [START_REF] Risser | ADNI: Simultaneous multiscale registration using large deformation diffeomorphic metric mapping[END_REF]. Figure 7 displays the two images to be registered and the deformation field. After a subsampling of the field of deformation amplitude, the data consist of observations taken on a rectangular grid of steps 1/56 and 1/59 on [0, 1] 2 , so they consist of 3420 points of the form X(i/56, j/59) i = 0, . . . , 56, j = 0, . . . , [START_REF] Vecchia | Estimation and model identification for continuous spatial processes[END_REF].

With these data, we consider the same autocovariance model as in Section 5.3.1. We estimate the parameters θ 1 and θ 2 from the same two methods as in Section 5.3.1. The first method provides θ1 = Proof of Lemma A.2. Using the stationary increments of the process, one has

Recall that

We have seen in the proof of Proposition 3.1 (( 33) and ( 34)) that for i sufficiently large

Thus the sum in ( 29) is bounded by

On the other hand, we have proved also in the proof of Proposition 3.1 that

giving the result. Thus, one has to check that

) which is true by the assumptions made. We skip the details.

A.3 Proof of the main results

Proof of Proposition 3.1. 1) By definition of V a,n in [START_REF] Clark | Practical geostatistics[END_REF] and identity [START_REF] Coeurjolly | Estimating the parameters of a fractional Brownian motion by discrete variations of its sample paths[END_REF], we get

Recall that n = n -L + 1 is the size of the vector ∆ a (X). In all the proof, j is assumed to vary from -L + 1 to L -1. We use a Taylor expansion of V ((i + j)δ n ) at (iδ n ) and of order q 2D:

Note that this expression is "telescopic" in the sense that if q < q 2D,

Combining [START_REF] Journel | Mining geostatistics[END_REF] (with i = 0 and q = 2D), the vanishing moments of the sequence a 2 * and (H 1 ) yields:

The first term is non-zero by (9) in Proposition 2.6 and a dominated convergence argument together with (H 1 ) shows that the last term is o(δ 2D+s n ) giving [START_REF] Cressie | Statistics for spatial data[END_REF].

where stands for the transpose. On the other hand, by the triangular inequality for the operator norm (which is the maximum of the |λ i |'s), one gets max i=1,...,n

In the proof of Theorem 3.4, we have established that Σ a (j) ,n op = o(n 1/2 δ 2D+s n ) leading to the result.

A.4 Proof of the remaining results in Section 3

Proof of Theorem 3.7. We use the definition of C a,n and the following decomposition:

Following the proof of Proposition 3.1, the second term is proportional to

which is negligible compared to (Const)n 1/2+αs-α(s+1/2α) by (H 3 ) and thus goes to 0 as n goes to infinity. Then Slutsky's lemma and Theorem 3.4 lead straightforwardly to the required result.

Proof of Corollary 3.9. Obviously, one has

Using the triangular inequality A + B 2 -A 2 B 2 + 2 A B , it suffices to have ∆ a (f ) 2 = o(Var(V a,n (X) 1/2 ) = o(n 1/2 δ 2D+s n ) to deduce the central limit theorem for X from that for X. By application of the Taylor-Lagrange formula, one gets