F Comte

C Lacour

Y Rozenholc

Adaptive estimation of the dynamics of a discrete time stochastic volatility model

Keywords: J.E.L. Classification number: C13-C14-C22 Keywords. Adaptive Estimation, Autoregression, Deconvolution, Heteroscedastic, Hidden Markov Model, Nonparametric Projection Estimator

This paper is concerned with the discrete time stochastic volatility model Yi = exp(Xi/2)ηi, Xi+1 = b(Xi) + σ(Xi)ξi+1, where only (Yi) is observed. The model is re-written as a particular hidden model: Zi = Xi + εi, Xi+1 = b(Xi) + σ(Xi)ξi+1, where (ξi) and (εi) are independent sequences of i.i.d. noise. Moreover, the sequences (Xi) and (εi) are independent and the distribution of ε is known. Then, our aim is to estimate the functions b and σ 2 when only observations Z1, . . . , Zn are available. We propose to estimate bf and (b 2 + σ 2 )f and study the integrated mean square error of projection estimators of these functions on automatically selected projection spaces. By ratio strategy, estimators of b and σ 2 are then deduced. The mean square risk of the resulting estimators are studied and their rates are discussed. Lastly, simulation experiments are provided: constants in the penalty functions defining the estimators are calibrated and the quality of the estimators is checked on several examples.

Introduction

In this paper, we consider the following model:

(1)

Y i = exp(X i /2)η i , X i+1 = b(X i ) + σ(X i )ξ i+1 ,
where (η i ) and (ξ i ) are two independent sequences of independent and identically distributed (i.i.d.) random variables (noise processes). Only Y 1 , . . . , Y n are observed, while the process of interest is the unobserved volatility V i = exp(X i /2), and in particular the functions driving its dynamics, b(.) and σ(.). We will describe an estimation method leading to nonparametric estimates of these functions.

For identifiability of the model, the density of η must be known (e.g. N (0, 1)). This model is often called a discrete time stochastic volatility process. Examples of such representation for economic or financial processes can be found in [START_REF] Shephard | Stochastic volatility. Selected readings. Advanced Texts in Econometrics[END_REF] or Shephard (2006). (1) Université Paris Descartes, MAP5, UMR CNRS 8145. fabienne.comte@parisdescartes.fr, yves.rozenholc@parisdescartes.fr. (2) Laboratoire de Probabilités et Statistique, Université Paris Sud-Orsay. claire.lacour@math.u-psud.fr.

1.1. Comparison with continuous time stochastic volatility models. Let us emphasize that our model is written directly in discrete time with fixed (to 1) step of observation.

This makes its structure different from stochastic volatility models when first considered in continuous time and in a second step discretely sampled. It would be conceivable to study a system [START_REF] Bennett | Errors-in-variables in joint population pharmacokinetic/pharmacodynamic modeling[END_REF] d log(S t ) = √ U t dW t dU t = m(U t )dt + (U t )dB t ,

where (W t , B t ) is a two dimensional (standard) Brownian motion and U is a positive diffusion process. In that case, define Y i = log(S i+1 /S i ). Then, Y i has the same distribution as V i η i where V i = ( i+1 i U s ds) 1/2 and η i are i.i.d. N (0, 1) random variables. Consequently, the first equation of [START_REF] Bennett | Errors-in-variables in joint population pharmacokinetic/pharmacodynamic modeling[END_REF] can lead exactly to the first equation of [START_REF] Baraud | Adaptive estimation in an autoregressive and a geometrical β-mixing regression framework[END_REF] with specific Gaussian distribution for η. The tools for estimating the common density of the V i 's are thus common to both models. But the second equations of models ( 1) and ( 2), even if they represent the same idea of a time dynamics, do not coincide.

Indeed, many statistical tools have been developed for the estimation of µ and , see Renò (2008), Bandi and Renò (2008), [START_REF] Comte | Nonparametric estimation for a stochastic volatility model[END_REF] in the nonparametric setting, or Gloter (2008) in the parametric context. But all these authors consider the high frequency data context, a set of assumptions which are often meaningful, but specific to observations with small step ∆ = ∆ n that tends to zero when n grows to infinity, and long interval of time: T n = n∆ n tends to infinity with n. Most strategies in this context require a huge sample of observations, which may not be available. Moreover, the discrete time approximations of the variables crucially exploit the small step assumption. Therefore, this does not correspond to the discrete time setting with fixed step of observation. Very few papers use another context to provide estimators of the continuous time model. Following ideas developed in Hansen and Scheinkman (1987) and later in [START_REF] Hansen | Spectral methods for identifying scalar diffusions[END_REF], Gobet et al. (2004) propose a strategy for low frequency data (i.e. fixed sample step ∆ of observation). They obtain nonparametric estimators in this context, but their method is specific to the underlying continuous time model. The methodology is similar to the operator eigenvalues approximation studied in [START_REF] Carrasco | Linear Inverse Problems in Structural Econometrics[END_REF] but can not be applied to our model.

Consequently, Model [START_REF] Baraud | Adaptive estimation in an autoregressive and a geometrical β-mixing regression framework[END_REF], and in particular the autoregressive equation in [START_REF] Baraud | Adaptive estimation in an autoregressive and a geometrical β-mixing regression framework[END_REF], is not an approximation of the continuous time model [START_REF] Bennett | Errors-in-variables in joint population pharmacokinetic/pharmacodynamic modeling[END_REF] (it would require a small sample step for an approximation to hold) and does not share any underlying continuous time structure. 1.2. Logarithmic transformation of (1) and related models. When the model is described in discrete time only, estimation strategies have been proposed for the estimation of the stationary density of (X i ) only, see van Es et al. (2005), [START_REF] Comte | Adaptive density deconvolution with dependent inputs[END_REF]. In this case, the model is rewritten as

(3) Z i = X i + ε i X i+1 = b(X i ) + σ(X i )ξ i+1
where ε i = ln(η 2 i ) -E(ln(η 2 i )) and Z i = ln(Y 2 i ) -E(ln(η 2 i )). In our setting, E(ln(η 2 i )) is known. Clearly, the logarithmic transformation of Y 2 i implies that the sign of Y i can not be recovered. When we are interested in volatility dynamics, the sign of the innovation in Y does not matter since we have the additional assumption of independent volatility and price innovation, i.e. (η) and (ξ) are independent. All our results require this independence assumption, and thus, this context does not authorize the so-called leverage effect (see also Remark 3.1).

It also happens that in some contexts, model ( 3) is considered directly. Such a non linear autoregressive model observed with additive noise is also called an errors-in-variables model and is commonly considered in economic and biological applications (see Schenmach (2007), [START_REF] Hong | A simple estimator for nonlinear error in variable models[END_REF] or Benett and Wakefield (2001)). Statistical methods used to study these models are mainly related to parametric or semiparametric specifications of the model, see [START_REF] Chanda | Large sample analysis of autoregressive moving average models with errors-in-variables[END_REF], [START_REF] Comte | Semi-parametric estimation in the nonlinear errors-invariables autoregressive model[END_REF].

Lastly, Model (3) belongs to the general class of hidden Markov models (HMM). These models constitute a very famous class of discrete time processes with applications in various areas (see [START_REF] Cappé | Inference in hidden Markov models[END_REF]). Here our model is simpler in the sense that our noise is additive, but in standard HMMs it is assumed that the joint density of (X i , Z i ) has a parametric form.

Statistical bibliography for model (3).

To our knowledge, the question of estimating b and σ 2 in Model (3) on the basis of observations Z 1 , . . . , Z n has not been studied yet.

Only the following regressive model

Z i = X i + ε i , Y i = b(X i ) + ξ i , in which (Y i
) and (Z i ) for i = 1, . . . , n + 1 are observed, has received attention. In that case, two processes are observed, all sequences (X i ), (ξ i ), (ε i ) can be supposed independent identically distributed (i.i.d.) and independent from each other, and (Y i ) is homoscedastic (σ(x) ≡ 1). In this context, Fan and Truong (1990), and Comte and Taupin (2007) study the problem of the estimation of b. See also [START_REF] Fan | Nonparametric function estimation involving errors-in-variables[END_REF], [START_REF] Fan | Multivariate regression estimation with errors-in-variables: asymptotic normality for mixing processes[END_REF], [START_REF] Ioannides | Nonparametric regression with errors in variables and applications[END_REF], [START_REF] Koo | B-spline estimation of regression functions with errors in variable[END_REF]. Most authors propose estimators of b based on the ratio of two estimators, namely an estimator of bf divided by an estimator of f , where f denotes the common density of the (i.i.d. in their context) (X i ).

Several papers develop estimation methods for f , see [START_REF] Fan | On the optimal rates of convergence for nonparametric deconvolution problems[END_REF], Pensky and Vidakovic (1999), [START_REF] Comte | Penalized contrast estimator for density deconvolution[END_REF], [START_REF] Carrasco | Linear Inverse Problems in Structural Econometrics[END_REF], and the optimality of the rates are studied in [START_REF] Fan | On the optimal rates of convergence for nonparametric deconvolution problems[END_REF], [START_REF] Butucea | Deconvolution of supersmooth densities with smooth noise[END_REF] and Butucea and Tsybakov (2007).

1.4. Statistical strategy. The quotient strategy is also adopted in our more general setting. More precisely, we assume that the process (X i ) is stationary, with stationary density denoted by f , and we estimate b (resp. b 2 + σ 2 ) as a ratio of an estimator of bf (resp. (b 2 +σ 2 )f ) divided by an estimator of f . We estimate f with the adaptive estimator proposed by [START_REF] Comte | Penalized contrast estimator for density deconvolution[END_REF]. To estimate the numerators, we adopt the same type of strategy as for the estimation of f , namely: first, we study a projection estimator, then we propose a model selection criterion to choose the best projection space as possible, in term of mean integrated risk bound.

In the setting of (3), regarding the identifiability of the model, it must be assumed that the distribution of ε, f ε (or equivalently of η) is fully known. For instance, the process η is often modeled as a standard Gaussian i.i.d. sequence, and then ε i has the distribution of ln(N (0, 1) 2 ) + ln(2) + C where C is the Euler constant. [START_REF] Van Es | Nonparametric volatility density estimation for discrete time models[END_REF] specifically study this case in terms of density estimation, however more general distributions can also be considered (see [START_REF] Comte | Penalized contrast estimator for density deconvolution[END_REF]Comte et al. ( , 2007))).

In this respect, this allows to consider general classes of noise density f ε and also various classes of regularities for the functions to estimate (bf , (b 2 + σ 2 )f , f ). Then, when all functions belong to fixed (but user-unknown) regularity spaces, our risk bounds provide rates of convergence of the estimators.

We want here to emphasize the following point. After the pioneering work of [START_REF] Fan | On the optimal rates of convergence for nonparametric deconvolution problems[END_REF], people mainly remembered that the estimation of a twice differentiable density measured with an additive gaussian noise ε had a deplorable logarithmic rate. This is true, but often, the functions to recover are much more regular than only twice differentiable. If the signal has the same regularity as the noise, then the rate can become polynomial again (see Example 3.1 below). This was found by [START_REF] Carrasco | Linear Inverse Problems in Structural Econometrics[END_REF] with their specific method and also by [START_REF] Pensky | Adaptive wavelet estimator for nonparametric density deconvolution[END_REF], [START_REF] Comte | Penalized contrast estimator for density deconvolution[END_REF], or [START_REF] Butucea | Sharp optimality for density deconvolution with dominating bias[END_REF] with wavelet, projection or kernel estimators respectively.

The plan of the paper is the following: we first give the notations, the assumptions and describe projection spaces in Section 2. Next, Section 3 explains the estimation strategy for b and gives bounds of the integrated mean square risk of the estimators. Section 4 develops the same study for the estimation of σ 2 . Simulation experiments are conducted in Section 5 in order to illustrate the method. Lastly, proofs are gathered in Sections 6-7-8 and an appendix (section 9) describes auxiliary tools.

General setting and assumptions

2.1. The principle. Let us assume that the sequence (X i ) is stationary and let us denote by f the common density of the X i 's. The principle of the estimation methods relies in all cases on a "Nadaraya-Watson-strategy" in the sense that b or b 2 + σ 2 are estimated as ratio of an estimator of = bf (respectively ϑ = (b 2 + σ 2 )f ) and an estimator of f . In all cases, we use the adaptive estimator of f described in [START_REF] Comte | Penalized contrast estimator for density deconvolution[END_REF] or [START_REF] Comte | Adaptive density deconvolution with dependent inputs[END_REF] which study independent and β-mixing contexts.

2.2. Notations and Assumptions. Subsequently we denote by u * the Fourier transform of the function u defined by u * (t) = e itx u(x)dx, and by u , u ∞ , u ∞,K , < u, v >, u * v the quantities

u 2 = u 2 (x)dx, u ∞ = sup x∈R |u(x)|, u ∞,K = sup x∈K |u(x)|, < u, v >= u(x)v(x)dx with zz = |z| 2 and u * v(x) = u(t)v(x -t)dt.
Moreover, we recall that for any integrable and square-integrable functions u, u 1 , u 2 , (4)

(u * ) * (x) = 2πu(-x) and u 1 , u 2 = (2π) -1 u * 1 , u * 2 .
We consider the autoregressive model [START_REF] Birgé | Minimum contrast estimators on sieves: Exponential bounds and rates of convergence[END_REF]. The assumptions are the following: A1 (a) The ε i 's are i.i.d. centered (E(ε 1 ) = 0) random variables with finite variance,

E(ε 2 1 ) = s 2 ε . The density of ε 1 , f ε , belongs to L 2 (R), and for all x ∈ R, f * ε (x) = 0. (b)
The ξ i 's are i.i.d. centered with unit variance (E(ξ 2 1 ) = 1). A2 The X i 's are stationary and absolutely regular. A3 The sequences (ξ i ) i∈N and (ε i ) i∈N are independent. It follows from A3 and the data generating process of the X i that the sequences (X i ) i∈N and (ε i ) i∈N are independent.

The Z i 's are observed but the X i 's are not, the stationary density f of the X i 's is unknown but the density f ε of the ε i 's is known.

Standard assumptions on b, σ and the ξ i 's ensure that the sequence (X i ) i∈Z is stationary with stationary density denoted by f . This sequence is also absolutely regular, with βmixing coefficients denoted by β(k), see Doukhan (1994) or [START_REF] Comte | Adaptive estimation of mean and volatility functions in (auto-) regressive models[END_REF] for precise sets of conditions. We shall consider that the mixing is at least arithmetical with rate θ, i.e. that there exists θ > 0 such that [START_REF] Butucea | Deconvolution of supersmooth densities with smooth noise[END_REF] ∀k ∈ N,

β(k) ≤ (1 + k) -(1+θ) ,
or, more often, geometrical, i.e. ∃θ > 0, ∀k ∈ N, β(k) ≤ e -θk . The definition of the β-mixing coefficients and related properties are recalled in Section 9.

As we develop a L 2 -strategy, we need the target functions to be square-integrable.

A4 The function to estimate ( = bf , ϑ = (b 2 + σ 2 )f , or f ) is square-integrable.

In the following, we also assume that f ε is such that

A5 For all t in R, A 0 (t 2 + 1) -γ/2 exp{-µ|t| δ } ≤ |f * ε (t)| ≤ A 0 (t 2 + 1) -γ/2 exp{-µ|t| δ }, with γ > 1/2 if δ = 0.
Under Assumption A5, when δ = 0, the errors are usually called "ordinary smooth" errors, and "super smooth" errors when δ > 0, µ > 0. The standard examples are the following : Gaussian or Cauchy distributions are super smooth of order (γ = 0, µ = 1/2, δ = 2) and (γ = 0, µ = 1, δ = 1) respectively, and the Laplace (symmetric exponential) distribution is ordinary smooth (δ = 0) of order γ = 2. When ε = ln(η 2 ) -E(ln(η 2 )) with η ∼ N (0, 1) as in van Es et al. (2005), then E(ln(η 2 )) =ln(2) -C, where C is the Euler Constant, and Var(ln(η 2 )) = π 2 /2 and ε is super-smooth with γ = 0, µ = π/2 and δ = 1:

f * ln(η 2 1 ) (x) = 2 ix √ π Γ(1 + ix), |f * ln(η 2 1 ) (x)| ∼ x→+∞ 2/ee -π|x|
S m = Span{ϕ m,j , j ∈ Z} = {f ∈ L 2 (R), supp(f * ) ⊂ [-mπ, mπ]}.
Moreover, (S m ) m∈Mn , with M n = {1, . . . , m n }, denotes the collection of linear spaces.

In practice, we should consider the truncated spaces S

(n) m = Span{ϕ m,j , j ∈ Z, |j| ≤ K n },
where K n is an integer depending on n, and the associated estimators under the additional assumption: x 2 ψ 2 (x)dx < A ψ < ∞, where ψ = bf, (b 2 + σ 2 )f or f is the function to estimate. This is done in [START_REF] Comte | Penalized contrast estimator for density deconvolution[END_REF] and does not change the main part of the study. For the sake of simplicity, we write the sums over Z in the theoretical part of the present study. For t belonging to a space S m of the collection (S m ) m∈Mn , let

γ n (t) = 1 n n i=1 ( t 2 -2Z i+1 u * t (Z i )), u t (x) = 1 2π t * (-x) f * ε (x)
.

The following sequence of equalities, relying on the Fourier equalities (4), explains the choice of the contrast γ n :

E(Z 2 u * t (Z 1 )) = E(b(X 1 )u * t (Z 1 )) = u * t * f ε (-.), bf = 1 2π t * f * ε (-.) f * ε (-.), (bf ) * = 1 2π t * , (bf ) * = t, bf = E(b(X 1 )t(X 1 )) = t(x)b(x)f (x)dx = t, , (7) using that
(8) E(σ(X 1 )ξ 2 u * t (X 1 + ε 1 )) = E(ξ 2 )E(σ(X 1 )u * t (X 1 + ε 1 )) = 0.
Therefore, we find that

E(γ n (t)) = t 2 -2 , t = t -2 - 2
is minimal when t = . Thus, we define

(9) ˆ m = arg min t∈Sm γ n (t)
As γ n is minimized over S m only, ˆ m is in fact an unbiased estimator of m , which in turn is expected to be near of . Indeed, the estimator can also be written

(10) ˆ m = j∈Z âm,j ( )ϕ m,j , with âm,j ( ) = 1 n n i=1 Z i+1 u * ϕ m,j (Z i ),
an clearly, [START_REF] Cappé | Inference in hidden Markov models[END_REF] implies that E(â m,j ( )) = a m,j ( ). Now, the decomposition of the Mean Integrated Squared Error (MISE) will show that m plays here the role of a bandwidth parameter. Thus, we have to explain how to select an adequate value of m. To this end, we define ˆ m, by setting m = arg min

m∈Mn γ n ( ˆ m ) + pen(m) ,
where the penalty function is given by pen(m) = κE(Z 2 2 )Ψ(m) where ( 11)

Ψ(m) =          ∆(m) n if 0 ≤ δ < 1/3 m [(3δ-1)/2]∧δ ∆(m) n if δ ≥ 1/3, and ∆(m) = 1 2π πm -πm dx |f * ε (x)| 2 ,
where x ∧ y := inf(x, y).

The idea behind this criterion is that γ n ( ˆ m ) estimates the squared bias part of the MISE, and pen(m) has the order of the variance term. A kind of cross-validation method is thus performed here.

In practice

E(Z 2 2
) is unknown and is replaced by its empirical version, (1/n) n i=1 Z 2 i . Then, the resulting penalty function, pen, becomes random. We note that γ n ( ˆ m ) = -j∈Z [â m,j ( )] 2 , which explains (13) below.

Remark 3.1. We can see here why we can not omit the independence assumption between (ε i ) and (ξ i ), and consider the so-called leverage effect. Indeed, assume that we have

ξ i = ρε i + 1 -ρ 2 ε ⊥ i , for an i.i.d. centered noise (ε ⊥ i ), independent of (ε i ).
Then, because of the dynamics of the autoregressive X i sequence, if the noises (ε i ) and (ξ i ) are not independent, then the sequences (X i ) and (ε i ) are not independent either. Then, as for a given i, X i and ε i are no longer independent in model (3), we can not estimate f (the basic convolution link is lost). To keep X i and ε i independent for a given i (and then f can still be estimated, see Comte et al. (2007[START_REF] Comte | Adaptive density deconvolution with dependent inputs[END_REF])), we have to write the dynamics of X i as follows: X i+1 = b(X i ) + σ(X i )ξ i . But then, Equality [START_REF] Cappé | Inference in hidden Markov models[END_REF] which justifies the definition of the contrast γ n (t) used for the estimation of , is no longer true. Indeed, instead of (8), we have

E(σ(X 1 )ξ 1 u * t (X 1 + ε 1 )) = ρE(σ(X 1 )ε 1 u * t (X 1 + ε 1 )
), and this last term is not necessarily zero. As a conclusion, the independence assumption between (ε i ) and (ξ i ) is necessary.

3.1.2. Second step: the estimators of f . The second stage of the estimation procedure is to estimate f . In fact, [START_REF] Comte | Penalized contrast estimator for density deconvolution[END_REF] explain how to estimate f in an adaptive way and the mixing context is studied in [START_REF] Comte | Adaptive density deconvolution with dependent inputs[END_REF]. The estimator of f on S m is defined by

(12) fm = j∈Z âm,j (f )ϕ m,j with âm,j (f ) = 1 n n i=1 u * ϕ m,j (Z i ).
Then we define f m,

(13) m = arg min m∈Mn    - j∈Z [â m,j (f )] 2 + p en(m)    ,
where the penalty function is given by p en(m) = κΨ(m) with Ψ(m) given by [START_REF] Comte | Adaptive estimation of the stationary density of discrete and continuous time mixing processes[END_REF]. For the properties of f m we refer to [START_REF] Comte | Penalized contrast estimator for density deconvolution[END_REF]. Up to the multiplicative constants, the control of the mean square risk of the estimator is the same as the one obtained for here.

3.1.3. Last step: the estimator of b. We estimate b on a compact set B only and the following additional assumption is required:

A6 (a) ∀x ∈ B, f 0 ≤ f (x) ≤ f 1 for two positive constants f 0 and f 1 . (b) b is bounded on B.
Then we can define:

(14) b = b m, m = ˆ m f m if ˆ m/ f m ≤ a n , b = b m, m = 0 else
, where a n is a sequence to be specified later.

3.2.

Risk bound for ˆ m and ˆ m.

3.2.1.

Risk bound for ˆ m . We define the following empirical centered process

ν n (t) = 1 n n k=1 (Z k+1 u * t (Z k ) -t, ),
and with ( 6) and ( 10), we note that the following equalities hold

-ˆ m 2 = -m 2 + m -ˆ m 2 = -m 2 + j∈Z (a m,j ( ) -âm,j ( )) 2 = -m 2 + j∈Z ν 2 n (ϕ m,j ). Therefore E -ˆ m 2 ≤ -m 2 + j∈Z Var[ν n (ϕ m,j )].
The termm 2 is a deterministic integrated squared bias term. The part corresponding to j∈Z Var[ν n (ϕ m,j )] is the variance term, which has to be bounded. Then, the risk bound on the estimate ˆ m is as follows: Proposition 3.1. Consider the estimator ˆ m of defined by ( 9) where = bf with b and f as in Model [START_REF] Birgé | Minimum contrast estimators on sieves: Exponential bounds and rates of convergence[END_REF]. Then under Assumptions A1-A4, if E(b 4 (X 1 )) < +∞ and θ > 1 for arithmetical mixing (see ( 5)), we have

(15) E( -ˆ m 2 ) ≤ -m 2 + 2E(Z 2 2 ) ∆(m) n + 8K m n .
Note that the last term m/n in ( 15) is always smaller that ∆(m)/n and might have been omitted, up to a less precise constant before ∆(m)/n.

Rate of the estimator. When f *

ε satisfies A5, then the order of the variance term is bounded by: C∆(m)/n ≤ C m 2γ+1-δ exp(2µ(πm) δ )/n. For the squared bias term, we have

-m 2 = (2π) -1 |x|≥πm |f * (x)| 2 dx.
To evaluate this term, regularity conditions must be considered for . We shall assume that belongs to the space:

(16) S s,a,r (A) = {u : +∞ -∞ |u * (x)| 2 (x 2 + 1) s exp{2a|x| r }dx ≤ A}, δ = 0 δ > 0 f ε ordinary smooth f ε supersmooth r = 0 Sobolev(s) π m = O(n 1/(2s+2γ+1) ) rate = O(n -2s/(2s+2γ+1) ) π m = [ln(n)/(2µ + 1)] 1/δ rate = O((ln(n)) -2s/δ ) r > 0 C ∞ π m = [ln(n)/2a] 1/r rate = O ln(n) (2γ+1)/r n m solution of m2s+2γ+1-r exp{2µ(π m) δ + 2aπ r mr } = O(n)
Table 1. Choice of m and corresponding rates under A1-A5 and if belongs to S s,a,r (A) defined by [START_REF] Comte | Semi-parametric estimation in the nonlinear errors-invariables autoregressive model[END_REF].

for nonnegative constants s, a, r and A > 0. When r = 0, this corresponds to Sobolev spaces of order s, which are classically considered. But we can also think that may be as regular as f ε can be. When r > 0, a > 0, this corresponds to analytic functions, which are called "super-smooth" functions. When belongs to a space S s,a,r (A) defined by ( 16), then the order of the squared bias is

-m 2 = (2π) -1 |x|≥πm |f * (x)| 2 dx ≤ Cm -2s exp(-2a(πm) r ).
In this setting, we obtain from [START_REF] Comte | Finite sample penalization in adaptive density deconvolution[END_REF] the risk bound terms:

Cm -2s exp(-2a(πm) r ) + C m 2γ+1-δ exp(2µ(πm) δ )/n. Example 3.1. Consider the example of a Gaussian η, η ∼ N (0, 1), i.e. ε = ln(η 2 ) - E(ln(η 2 )
) is super-smooth with γ = 0, µ = π/2 and δ = 1. It is true that if is in a Sobolev i.e. belongs to S s,a,r (A) with r = a = 0, then the best attainable rate is of order [log(n)] -2s . But if is as regular as f ε and belongs to S s,a,r (A) with s = 0, r = 1, a > 0, then the optimal choice is π m = log(n)/(π + 2a) and the resulting rate is polynomial, of order n -2a/(π+2a) . This is summarized in Proposition 3.2 and Table 1.

Proposition 3.2. Assume that the Assumptions of Proposition 3.1 hold. In addition, assume that belongs to a space S s,a,r (A) defined by [START_REF] Comte | Semi-parametric estimation in the nonlinear errors-invariables autoregressive model[END_REF] and that Assumption A5 is fulfilled. Then the estimate ˆ m with m as in Table 1, has the rates given in Table 1 in terms of its mean square integrated risk E( ˆ m -2 ).

When r > 0, δ > 0 the value of m is not explicitly given. It is obtained as the solution of the equation m2s+2γ+1-r exp{2µ(π m) δ + 2aπ r mr } = O(n). For explicit general formulae of the rates in these cases, we refer to [START_REF] Lacour | Rates of convergence for nonparametric deconvolution[END_REF].

We just recall the particular case δ = r, for positive other parameters: this gives the (almost) polynomial rate

n -µ/(a+µ) (log(n)) -τ with τ = [2µδ + (δ -2γ -1)a]/[(a + µ)δ].

3.3.

Risk bound of ˆ m. These rates enhance the interest of building an estimator for which the choice of the relevant model m is automatically performed. This is done with ˆ m, and we can prove the following result: Theorem 3.1. Assume that Assumptions A1-A4 hold, that E(b 8 (X 1 )), E(σ 8 (X 1 )) and E(ξ 8 1 ) are finite and that E(ε 6 1 ) < +∞. Assume moreover that the process X is geometrically β-mixing, (or arithmetically β-mixing with θ > 14) and that the collection M n is such that for all m ∈ M n , pen(m) ≤ 1, then

E( ˆ m -2 ) ≤ C inf m∈Mn -m 2 + pen(m) + C n .
The proof of Theorem 3.1 is sketched in Section 6. Theorem 3.1 shows that the estimator automatically selects the optimal m when δ ≤ 1/3: indeed, in that case, the penalty has exactly the same order as the variance (namely ∆(m)/n). When δ > 1/3, a compromise is still performed, but the penalty is slightly greater than the variance. In an asymptotic setting, this implies a loss in the rate of convergence of the estimator, but this loss can be shown to be negligible with respect to the rates. For discussions on this point, see [START_REF] Comte | Penalized contrast estimator for density deconvolution[END_REF]. Theorem 3.2. Assume that Assumptions A1-A4 hold. Assume that the process X is geometrically β-mixing, (or arithmetically β-mixing with θ > 3 in ( 5)) and that the collection

M n is such that for all m ∈ M n , p en(m) ≤ 1, then E( f m -f 2 ) ≤ C inf m∈Mn f -f m 2 + p en(m) + C n
where f m denotes the orthogonal projection of f on S m .

Then it is common (see e.g. [START_REF] Lacour | Nonparametric estimation of the stationary density and the transition of a Markov chain[END_REF] or Comte and Taupin (2007)) to obtain a bound on the MISE computed on a compact set, denoted here by B. In the following result, we denote by h 2 B := B h 2 (x)dx, for any square integrable function h on B. Theorem 3.3. Assume that the assumptions of Theorems 3.1 and 3.2 and A5-A6 hold, that f belongs to a space S s,a,r (A) with s > 1/2 if r = 0 and to a space S s ,a ,r (A ), that ln(ln(n)) ≤ m n ≤ (n/ ln(n)) 1/(2γ+1) for f m. Then, for n great enough, we have

E( b -b 2 B ) ≤ C 1 E( f m -f 2 ) + C 2 E( ˆ m -2 ) + C 3 n
where a n = n ω with ω > 1/2 and C 1 , C 2 , C 3 are constants.

The MISE bound in Theorem 3.3 means that the rate for estimating b will be the worst between the rates of estimation of f and bf , which are optimal or near-optimal.

The proof of this result is omitted. The reader is referred to [START_REF] Lacour | Nonparametric estimation of the stationary density and the transition of a Markov chain[END_REF] or Comte and Taupin (2007).

Estimation of σ 2

4.1. Steps of the estimation. We now aim at estimating σ 2 and we also follow the strategy described in Section 2.1.

First step. We set ϑ = (b 2 + σ 2 )f and we first estimate ϑ. To this end, we consider the following contrast:

(17) γn (t) = t 2 - 2 n n k=1 Z 2 k+1 -s 2 ε u * t (Z k ).
As f ε is assumed to be known, then so is the variance s 2 ε . Suggestions on how to estimate this quantity can be drawn from Butucea and Matias (2004), but this has a cost and they obtain logarithmic rates of convergence for both the variance estimator and the induced plugged in estimator of the density. Then we define [START_REF] Delyon | Limit theorem for mixing processes[END_REF] θm = arg min where p en(m) is a penalty function given by: p en(m) = κE((Z 2 2s 2 ε ) 2 )Ψ(m), with Ψ(m) given by [START_REF] Comte | Adaptive estimation of the stationary density of discrete and continuous time mixing processes[END_REF]. Again, the expectation

E[(Z 2 2 -s 2 ε ) 2 ] is replaced by its empirical version in practice, namely (1/n) n i=1 (Z 2 i -s 2 ε ) 2 .
Second step. As previously, we use as an estimator of f , the estimator f m defined by ( 12)-( 13). Its risk is controlled by Theorem 3.2.

Third step. We obtain, by defining, similarly to ( 14),

b 2 + σ 2 = θ m f m if θ m/ f m ≤

4.2.

Risk bounds for θm and θ m. It is not difficult to check that E(γ n (t)) = t 2 -2 ϑ, t which justifies the choice of γn given in [START_REF] Comte | Adaptive estimation in a nonparametric regression model with errors-in-variables[END_REF]. We can also easily obtain the decomposition γn (t)γn (s

) = t -ϑ 2 -s -ϑ 2 -2ν n (t -s) where νn (t) = 1 n n k=1 (Z 2 k+1 -s 2 ε )u * t (Z k ) -t, ϑ .
As for b previously, we can write that

θm -ϑ 2 = ϑ m -ϑ 2 + j∈Z ν2 n (ϕ m,j ).
With the same tools as for the study of , using a relevant decomposition of the empirical process νn , we prove (see Section 7) that:

Proposition 4.1. Consider the estimator θm of ϑ defined by [START_REF] Delyon | Limit theorem for mixing processes[END_REF] where ϑ = (b 2 + σ 2 )f with b, σ and f as in Model (3). Then under Assumptions A1-A4, and if ξ 2 , ε 1 , b 2 (X 1 ) and σ 2 (X 1 ) admit moments of order 4, then

E( ϑ -θm 2 2 ) ≤ ϑ -ϑ m 2 + 4E[(Z 2 2 -s 2 ε ) 2 ] ∆(m) n + K m n where K = 16 2 k≥0 (k + 1)β(k)E((b 2 (X 1 ) + σ 2 (X 1 )) 4 ) if k kβ(k) < +∞.
The empirical processes involved in the decomposition of νn are of the same type as the processes studied for the estimation of . Therefore, we give the risk bound for θ m but we omit the proof.

Theorem 4.1. Assume that Assumptions A1-A4 hold, that E(b p (X 1 )), E(σ p (X 1 )) and E(ξ p 1 ) are finite for a p ≥ 16 and that E(ε 12 1 ) < +∞. Assume that the process X is geometrically β-mixing and that the collection M n is such that for all m ∈ M n , p en(m) ≤ 1, then

E( θ m -ϑ 2 ) ≤ C inf m∈Mn ϑ -ϑ m 2 + p en(m) + C n .

Simulation results

5.1. The models. We start from the first equation of the multiplicative model, Y i = exp(X i /2)η i ,, i = 1, . . . , n. Then we compute:

(19) Z i = ln(Y 2 i ) -E(ln(η 2 1
)), i = 1, . . . , n and the Z i 's follow the equation

Z i = X i + ε i , where ε i is centered with variance s 2 ε . We have ε i = ln(η 2 i ) -E(ln(η 2 1
)). We consider three types of noises: the case ε i Laplace with variance s 2 ε,L , the case ε i Gaussian N (0, σ 2 ε,G ), and the case called hereafter "logχ 2 (1)" where η i in N (0, 1); in that case, E(ln(η 2 1 )) = -Cln(2) where C is the Euler constant, and s 2 ε,χ = π 2 /2. In practice, we generate directly the ε i 's. To compute the Z i 's, we need also to generate the X i 's.

We want to experiment different models, thus we have to choose functions b, σ. In all cases, the ξ i 's are i.i.d. N (0, 1) random variables. Then we consider the following functions, defining three models (20)

   b 0 1 (x) = 0.25x σ 0 1 (x) = 3, b 0 2 (x) = 0.25 sin(2πx + π/3) σ 0 2 (x) = ( √ 2/7)(0.31 + 0.7 exp(-5x 2 )), b 0 3 (x) = -0.25(x + 2 * exp(-16x 2 )) σ 0 3 (x) = 0.2 + 0.4 exp(-2x 2
). The first model is a simple linear model, the two other ones are arbitrary chosen to lead to different types of curves.

We compute recursively, up to initial conditions, a sample U 1 , . . . , U n with U i+1 = b 0 (U i ) + σ 0 (U i )ξ i+1 , for (b 0 , σ 0 ) taken as one of the (b 0 j , σ 0 j ), j = 1, 2, 3. Note that we drop out the first 100 observations to have a stationary sample.

At last, we take (X i , b, σ) = (U i , b 0 , σ 0 ) when the ε i 's are Laplace or Gaussian, and

X i = U i /s ε , b(x) = b 0 (s ε x)/s ε , σ(x) = σ 0 (s ε x)/s ε when the ε i 's are log χ 2 (1)
. The factor sε is a scaling factor computed to have the same signal to noise ratio (i.e. the ratio Var(X)/Var(ε)) for the three different noises ε. Therefore we take s ε,L = s ε,G = 1 and thus sε = √ 2/π for Model 1, and s ε,L = s ε,G = 0.05 and thus sε = 0.05 √ 2/π for Models 2 and 3.

Description of the estimation algorithm.

Let us describe the way the procedure works. Consider data (Y i ) 1≤i≤n for a given noise η in Model [START_REF] Baraud | Adaptive estimation in an autoregressive and a geometrical β-mixing regression framework[END_REF]. Then E(ln(η 2 1 )) is known and we compute Z i as given by [START_REF] Van Es | Nonparametric volatility density estimation for discrete time models[END_REF]. For instance if η 1 ∼ N (0, 1), we take

Z i = ln(Y 2 i ) + C + ln(2)
where C is the Euler constant.

The first step is to estimate . Given the Z i , we compute

âm,j ( ) = 1 n n i=1 Z i+1 u * ϕ m,j (Z i ) = 1 2πn n i=1 Z i+1 e -iuZ i ϕ m,j (u) f * ε (u)
du for j = -N, . . . , N , and N great enough. Note that the a m,j ( ) are real. Indeed, let φZ (x) = (1/n) n k=1 Z k+1 e ixZ k . Then with an elementary change of variable, we obtain

âm,j ( ) = 1 2π π -π φZ (-mu)e iju f * ε (mu) √ mdu = 1 π π 0 Re φZ (-mu)e iju f * ε (mu) √ mdu,
where Re(z) denotes the real part of the complex number z. Then m is selected by minimizing over a grid of possible m's the quantity:

γ n ( ˆ m ) + pen(m) = - j [â m,j ( )] 2 + pen(m) with pen(m) = 1 2n 1 n n i=1 Z 2 i 1 + ln 2.5 (πm) (πm)(1 + s 2 ε ) + C(m, s ε ) ∆(m), ∆(m) = 2 πm 0 du |f * ε (u)| 2 where (21) f * ε (x) = 1 1 + 1 2 s 2 ε,L x 2 (Laplace), f * ε (x) = exp(- 1 2 s 2 ε,G x 2 ) (Gaussian)
or, for the log-χ 2 (1):

(22) f * ε (x) = 1 √ π 2 i(x+C+ln(2)) Γ(1 + i(x + C + ln(2))) and C(m, s ε ) = 0 if ε 1 is Laplace, C(m, s ε ) = 2 3 (πm) 2 s 2 ε,G if ε 1 is Gaussian, C(m, s ε ) = πms ε,χ if ε 1 ∼ ln(χ 2 (1)).
Let us mention that the Fourier transforms f * ε above correspond to the densities The penalties are calibrated so that we recover standard orders m/n when s ε → 0. Indeed, setting s ε → 0 amounts to set the convolution noise to zero and is a limit case where X i would be observed.

• f ε (x) = exp(- √ 2|x|/s ε,L )/( √ 2s ε,L ) for the centered Laplace with variance s 2 ε,L , • f ε (x) = exp(-x 2 /(2s 2 ε,G ))/( √ 2πs ε,G ) for the N (0, s 2 ε,G ), • f ε (x) = K(x -ln(2) -C)
Moreover, the penalties correspond to the variance order ∆(m)/n, with no loss for Laplace errors, with a loss of order m for the log χ 2 (1) case (where δ = 1) and a loss of order m 2 in the Gaussian case (δ = 2), see formula [START_REF] Comte | Adaptive estimation of the stationary density of discrete and continuous time mixing processes[END_REF]. The common term ln 2.5 (πm) is an additional adjustment for small m's, which has no asymptotic weight compared to the other terms.

For the estimation of f , the same procedure applies with

• the function φZ replaced by f * Z (x) = (1/n) n k=1 e ixZ k , • the factor (1/n) n i= Z 2
i in the penalty replaced by 1. For the estimation of ϑ, the same procedure applies with

• the function φZ replaced by

φZ 2 (x) = 1 n n k=1 (Z 2 k+1 -s 2 ε )e ixZ k ,
• the factor (1/n) n i= Z 2 i in the penalty replaced by (1/n) n i=1 (Z 2 is 2 ε ) 2 . In the three cases, we select the three values of m among 125 values ranging from 10 ln(n)/(πn) to 10 and of course, they are not the same in general.

Simulation results. Let us present the results of our simulation experiments.

First, we performed a Monte Carlo study which is reported in Table 1. For each simulation, we compute the average squared error (ASE) at 101 grid points of our adaptive estimator and average these ASEs over 100 replications. It must be noticed that we give here the results for the functions b 0 j , σ 0 j , so that the ASE's can be compared in function of the distribution of the noise ε. We can see that, as already noticed in the density deconvolution setting, there is little difference between Laplace, logχ 2 (1) and Gaussian ε i 's, in spite of the difference between the theoretical rates. Moreover, it appears clearly that increasing the sample size leads to noticeable improvements of the results. Figures 2 to 6 illustrate our results in the three cases of couples (b j , σ 2 j ) given in (20) with rescaling for logχ 2 (1) errors ε. The true curves are bold dotted and are given together with 20 estimated curves. This set of twenty estimated curves shows, in each case, that the confidence band around the true function has very small width. Moreover, in all cases, and for all the noise distributions, the estimation are very good. With the financial setting in mind, we took samples with large sizes n = 1000 or n = 5000 which are clearly good sizes for nonparametric estimation.

We can see that b and b 2 + σ 2 are well estimated by the ratio strategy. The extraction of σ 2 sometimes suffers from scale problems (if σ 2 is much smaller that b 2 or if both are very small). In particular, we plot σ 2 with the same vertical scale as b 2 + σ 2 to take this into account.

We can therefore conclude that the method we propose is very precise to estimate b and b 2 + σ 2 , and can give an interesting idea of the shape of σ

n (t) + ν (2) n (t) + ν (3) n (t) with (23) ν (1) n (t) = 1 n n k=1 ε k+1 u * t (Z k ), ν (2) n (t) = 1 n n k=1 ξ k+1 σ(X k )u * t (Z k ), ( 24
) ν (3) n (t) = 1 n n k=1 (b(X k )u * t (Z k ) -t, ).
Here the terms ν

n and ν

(2)

n can be kept together and benefit from the uncorrelatedness of the variables involved in the sums. The term ν we find Var[ν n (ϕ m,j )] ≤ 2Var ν (1) n (ϕ m,j ) + ν (2) n (ϕ m,j ) + 2Var ν (3) n (ϕ m,j ) . The first variance involves uncorrelated and centered terms and leads to

Var 1 n n i=1 (ε i+1 + σ(X i )ξ i+1 )u * ϕ m,j (Z i ) = 1 n 2 n i=1 E[(s 2 ε + σ 2 (X i ))|u * ϕ m,j (Z i )| 2 ] so that j∈Z Var ν (1) n (ϕ m,j ) + ν (2) n (ϕ m,j ) = 1 n j∈Z E[(s 2 ε +σ 2 (X 1 ))|u * ϕ m,j (Z 1 )| 2 ] = (s 2 ε + E(σ 2 (X 1 )))∆(m) n .
We use here the following useful property of our basis (resulting from a Parseval's formula):

∀x ∈ R, j |u * ϕ m,j (x)| 2 = ∆(m),
where ∆(m) is defined by [START_REF] Comte | Adaptive estimation of the stationary density of discrete and continuous time mixing processes[END_REF] and the u * ϕ m,j (x) are just rewritten as Fourier coefficients. For the second term, we use the standard tools specific to the β-mixing context (namely [START_REF] Viennet | Inequalities for absolutely regular sequences: application to density estimation[END_REF] covariance Inequality) and we can easily prove the following Lemma: Lemma 6.1. Under Assumptions A1-A3, and if E(b 4 (X 1 )) < ∞ and k kβ(k) < +∞, then j∈Z Var ν (3) n (ϕ m,j ) ≤ E(b 2 (X 1 ))

∆(m) n + 4Km n ,
where

K = 2 k≥0 (k + 1)β(k)E(b 4 (X 1 )).
The result of Proposition 3.1 follows. 2 Proof of Lemma 6.1.

Var ν (3) n (ϕ m,j ) =

1 n 2 n k=1 Var b(X k )u * ϕ m,j (Z k ) + 1 n 2 1≤k =l≤n cov(b(X k )u * ϕ m,j (Z k ), b(X l )u * ϕ m,j (Z l )). F. COMTE, C. LACOUR, Y. ROZENHOLC Then E(b(X k )u * ϕ m,j (Z k )) = 1 2π E(b(X 1 )e ixX 1 )ϕ * m,j (-x)dx = bf, ϕ m,j = E(b(X 1 )ϕ m,j (X 1 ))
and note that, for k = l,

E(b(X k )u * ϕ m,j (Z k )b(X l ) ū * ϕ m,j (Z l )) = 1 4π 2 E(b(X k )b(X l )e ixX k -iyX l )ϕ * m,j (-x)ϕ * m,j (y)dxdy = E[b(X k )b(X l )ϕ m,j (X k )ϕ m,j (X l )]. Therefore, Var 1 n n k=1 b(X k )u * ϕ m,j (Z k ) ≤ 1 n Var(b(X 1 )u * ϕ m,j (Z 1 )) + Var 1 n n k=1 b(X k )ϕ m,j (X k ) .
The last term requires a covariance inequality for mixing variables [START_REF] Delyon | Limit theorem for mixing processes[END_REF], Viennet (1997), Theorem 9.1 in the appendix) and uses the fact that the X i 's are β-mixing with coefficients β(k).

j∈Z

Var 1 n n i=1 b(X i )ϕ m,j (X i ) ≤ j∈Z 4 n β(x)b 2 (x)|ϕ m,j (x)| 2 f (x)dx ≤ 4m n β(x)b 2 (x)f (x)dx
where β is a nonnegative function such that E(β p (X)) ≤ p k≥0 (k + 1) p-1 β(k) and by using that

j |ϕ m,j | 2 (.) ∞ = m. Therefore if E(b 4 (X 1 )) < ∞ and θ > 1, then j∈Z Var 1 n n i=1 b(X i )ϕ m,j (X i ) ≤ 4m 2 k≥0 (k + 1)β(k)E(b 4 (X 1 )) n . Moreover j∈Z Var b(X 1 )u * ϕ m,j (Z 1 ) ≤ E b 2 (X 1 ) j∈Z (u * ϕ m,j (Z 1 )) 2 so that j∈Z 1 n Var b(X 1 )u * ϕ m,j (Z 1 ) ≤ E(b 2 (X 1 ))∆(m) n ,
which gives the result.

(t) -γ n (s) = t -2 -s -2 -2ν n (t -s), and
2ν n ( ˆ m -m ) ≤ 1 4 ˆ m -m 2 + 4 sup t∈B m, m (0,1) ν 2 n (t).
Thus, we obtain, as

ˆ m -m 2 ≤ 2 ˆ m -2 + 2 m -2 that (25) 1 2 ˆ m -2 ≤ 3 2 m -2 + pen(m) + 4 sup t∈B m, m (0,1) ν 2 n (t) -pen( m).
Then we need to find a function p(m, m ) such that

(26) E sup t∈B m, m(0,1) ν 2 n (t) -p(m, m) + ≤ C n
which in turn will fix the penalty function through the requirement: ∀m, m ∈ M n , [START_REF] Hong | A simple estimator for nonlinear error in variable models[END_REF] 4p(m, m ) ≤ pen(m) + pen(m ).

Gathering ( 25), ( 26) and ( 27) will lead to, ∀m ∈ M n ,

1 2 E ˆ m -2 ≤ 3 2 m -2 + 2pen(m) +
4C n which is the result. Now, if ν n is split into several terms, deduced from the first decomposition given by ( 23)-( 24), say ν n (t) = 3

i=1 p i j=1 ν (i,j)
n (t) where p i ≤ 3, then, up to some multiplicative constants, inequality (26) will follow from inequalities

E sup t∈B m, m (0,1) [ν (i,j) n (t)] 2 -p i,j (m, m) + ≤ C i,j n ,
with C = 9 i,j C i,j and p(m, m ) = 9 i,j p i,j (m, m ). The study of the ν

(i,j) n (t) is explained below.
First we split ν

n in two parts, so that both expressions involve independent variables, conditionally to (X): ν

(1) n = ν (1,odd) n + ν (1,even) n where ν (1,even) n (t) = 1 n 1≤2k≤n ε 2k+1 u * t (Z 2k ), ν (1,odd) n (t)
= 1 n 1≤2k+1≤n ε 2k+2 u * t (Z 2k+1 ).
Now, we shall study ν

(1,even) n only since both terms lead to the same type of result. As Talagrand's Inequality requires the random variables involved to be bounded, we have an additional step that allows to obtain the result under a moment condition on the ε i 's: ν

(1,even) n = ν (1,1) n + ν (1,2) n + ν (1,3) n with ν (1,1) n (t) = 1 n 1≤2k≤n ε 2k+1 1 I |ε 2k+1 |≤n 1/4 u * t (Z 2k ) -E X (ε 2k+1 1 I |ε 2k+1 |≤n 1/4 u * t (Z 2k )) , ν (1,2) n (t)
= 1 n 1≤2k≤n E(ε 2k+1 1 I |ε 2k+1 |≤n 1/4 ) [t(X 2k ) -E(t(X 2k ))]
and

ν (1,3) n (t) = 1 n 1≤2k≤n ε 2k+1 1 I |ε 2k+1 |>n 1/4 u * t (Z 2k ) -E(ε 2k+1 1 I |ε 2k+1 |>n 1/4 u * t (Z 2k )) ,
where E X denotes the conditional expectation given (X k ) 1≤k≤n+1 . It is worth noticing that ν

(1,2) n vanishes if the ε's are symmetric and ν

(1,3) n (t) is negligible under adequate moment conditions on the ε's.

The following lemmas are proved below:

Lemma 6.2. E sup t∈B m, m (0,1) [ν (1,1) n (t)] 2 -p 1,1 (m, m) + ≤ C n ,
where p 1,1 (m, m ) = KE(ε 2 1 )Ψ(m ∨ m ) where K is a numerical constant and Ψ(m) is defined by [START_REF] Comte | Adaptive estimation of the stationary density of discrete and continuous time mixing processes[END_REF]. Lemma 6.3. If the process (X k ) is geometrically β-mixing (or arithmetically with θ > 3), then

E sup t∈B m, m (0,1) [ν (1,2) n (t)] 2 -KE(|ε 1 |) k β(k) m + m n + ≤ C n . Lemma 6.4. If E(ε 6 1 ) < +∞, and m n is the largest value of m such that ∆(m n )/n ≤ 1, then, E sup t∈B m, m (0,1)
[ν (1,3)

n (t)] 2 ≤ E sup t∈Bm n (0,1) (ν (1,3) n (t)) 2 ≤ 2E(ε 6 1 ) n .
For the study of ν

n (t), a result is given, whose proof is detailed in Section 8:

Lemma 6.5. Let τ n (t) = ν (2)
n (t) = (1/n) n k=1 ξ k+1 σ(X k )u * t (Z k ). Under the assumptions of Theorem 3.1, E sup t∈B m, m(0,1) [τ n (t)] 2 -p τ (m, m) + ≤ C n ,
where

p τ (m, m ) = κE(σ 2 (X 1 ))Ψ(m ∨ m ).
For ν

n (t) we write ν

n (t) = ν (3,1) n (t) + ν (3,2) n (t) with ν (3,1) n (t) = 1 n n k=1 [b(X k )u * t (Z k ) -b(X k )t(X k ))], ν (3,2) n (t) = 1 n n k=1 [b(X k )t(X k ) -t, ], where b(X k )t(X k ) = E (X) [b(X k )u * t (Z k )]. For ν (3)
we can apply Talagrand's Inequality conditionally to (X), for ν

, we can use approximation techniques. More precisely, using the same techniques as previously, we get Lemma 6.6. If E(b 8 (X 1 )) < +∞, and (X i ) i∈N is arithmetically β-mixing with θ > 14, then

E sup t∈B m, m (0,1) [ν (3,1) n (t)] 2 -p 3,1 (m, m) + ≤ C n , where p 3,1 (m, m ) = KE(b 2 (X 1 ))Ψ(m ∨ m ), and
E sup t∈B m, m (0,1) [ν (3,2) n (t)] 2 -K (E(b 4 (X 1 )) k (k + 1)β(k)) 1/2 m + m n + ≤ C n ,
where κ and κ are numerical constants.

The proof of the result concerning ν

(3,1) n follows the same line as the proof of Lemma 6.5, which is detailed in section 8 and is therefore omitted here. For ν (3,2) n , the bound can be obtained directly by applying Talagrand's inequality (see Theorem 9.2) to this process, if b is bounded. As this is not assumed, we write

b = b 1 + b 2 with b 1 (x) = b(x)1 I |b(x)|≤n 1/4 and b 2 (x) = b(x)1 I |b(x)|>n 1/4
. This allows to split the process in two parts and consequently to obtain the result under E(|b(X 1 )| 8 ) < +∞ and m n ≤ √ n, where m n is the largest over the m ∈ M n (a condition which is fulfilled in our problem).

Proof of Lemma 6.2.

We apply Lemma 9.2 to process ν

(1,1) n (t) conditionally to the sequence (X k ) 1≤k≤n . Given the X i 's, the variables (Z 2k , ε 2k+1 ) k≥1 are independent and we have, for m = m ∨ m ,

E X sup t∈B m,m (0,1) (ν (1,1)
n (t)) 2 ≤ j∈Z E X     1 n 1≤2k≤n ε 2k+1 1 I |ε 2k+1 |≤n 1/4 u * ϕ m ,j (Z 2k ) -E X (ε 2k+1 1 I |ε 2k+1 |≤n 1/4 u * ϕ m ,j (Z 2k ))   2   ≤ j∈Z 1 n 2 1≤2k≤n Var X ε 2k+1 1 I |ε 2k+1 |≤n 1/4 u * ϕ m ,j (Z 2k ) ≤ j∈Z 1 n 2 1≤2k≤n E X ε 2k+1 1 I |ε 2k+1 |≤n 1/4 u * ϕ m ,j (Z 2k ) 2 = 1 n 2 1≤2k≤n E X   ε 2 2k+1 1 I |ε 2k+1 |≤n 1/4 j∈Z (u * ϕ m ,j (Z 2k )) 2   ≤ E(ε 2 1 ) ∆(m ) n := H 2 , as j |u * ϕ m,j (.)| 2 ∞ ≤ ∆(m) by using Parseval's formula. Next sup x,y |y1 |y|≤n 1/4 u * t (x)| ≤ n 1/4 u * t (.) ∞ ≤ n 1/4 ∆(m * ) := M 1 .
Lastly, following the same method as in [START_REF] Comte | Penalized contrast estimator for density deconvolution[END_REF], Lemma 4, we get sup t∈B m,m (0,1)

1 n k Var X (ε 2k+1 1 I |ε 2k+1 |≤n 1/4 u * t (Z 2k ))
≤ sup t∈B m,m (0,1)

1 n k E(ε 2 2k+1 )E X [(u * t (Z 2k )) 2 ] ≤ E(ε 2 1 ) ∆ 2 (m )/(2π), where (28)
∆ 2 (m) = m 2 ϕ * (x)ϕ * (y) f * ε (mx)f * ε (my) f * ε (m(x -y)) 2 dxdy.
Then the usual bounds for ∆ 2 hold, namely, ∆

2 (m ) ≤ ∆(m ) if δ > 1 and if δ ≤ 1, ∆ 2 (m )/2π ≤ κ∆(m )/(m ) (1-δ)/2 . This gives v = c∆(m * )(m ) -(1-δ) + /2 .
Given that the orders are the same as in Comte et al. for v and H 2 and inserting the slight difference on M 1 , it can easily be checked that the conclusion still holds and therefore the result of Lemma 6.2 follows. 2

Proof of Lemma 6.3. The result given in Lemma 6.3 is a standard result of density estimation for mixing variables. We refer the reader to [START_REF] Tribouley | L p adaptive density estimation in a β mixing framework[END_REF]

E sup

t∈Bm n (0,1)

(ν (1,3) n (t)) 2 ≤ j∈Z Var 1 n k e 2k+1 u * ϕ mn,j (Z 2k ) = 1 n 2 j∈Z   k Var e 2k+1 u * ϕ mn,j (Z 2k ) + k =l cov(e 2k+1 u * ϕ mn,j (Z 2k ), e 2l+1 u * ϕ mn,j (Z 2l ))   ≤ E(e 2 3 )∆(m n ) n + 1 n 2 j∈Z E(e 3 ) 2 k =l cov(ϕ mn,j (X 2k ), ϕ mn,j (X 2l )) ≤ E(e 2 3 )∆(m n ) n + 1 n 2 j∈Z E(e 3 ) 2 Var k ϕ mn,j (X 2k ) ≤ E(e 2 3 ) ∆(m n ) n + m n k β k n ≤ 2∆(m n ) n E ε 2 1 1 I |ε 1 |≥n 1/4 ≤ 2E ε 6 1 /n .2 7. Proof of Proposition 4.1.
We bound the expectations of the empirical processes involved in order to obtain the bound of j∈Z E(ν n (ϕ m,j )), using the decomposition νn = 4 i=1 ν(i) n with ν( 1)

n (t) = (1/n) n k=1 (b 2 (X k ) + σ 2 (X k ))u * t (Z k ) -t, ϑ , ν(2) n (t) = (1/n) n k=1 (ξ 2 k+1 -1)σ 2 (X k )u * t (Z k ), ν (3)
n (t) = 1 n n k=1 (ε 2 k+1 -s 2 ε )u * t (Z k ), and ν(4) n (t) = (2/n) n k=1 [ε k+1 ξ k+1 σ(X k )+b(X k )ε k+1 +σ(X k )b(X k )ξ k+1 ]u * t (Z k ). But it is clear that ν(1)
n is the same process as ν

(3)

n with b(X k ) replaced by (b 2 + σ 2 )(X k ), that ν(2) n is of the same type as ν (2) n with σ(X k ) replaced by σ 2 (X k ) and ξ k+1 by ξ 2 k+1 -1. Next, ν (3)
n corresponds to ν

(1)

n with ε k+1 replaced by ε 2 k+1 -s 2 ε . Lastly j∈Z E[(ν (4) n (ϕ m,j )) 2 ] = 4 n j∈Z E[(ε 2 ξ 2 σ(X 1 ) + b(X 1 )ε 2 + (σb)(X 1 )ξ 2 ) 2 (u * ϕ m,j (Z k )) 2 ] ≤ 4∆(m) n [s 2 ε E(σ 2 (X 1 )) + E(b 2 (X 1 )(s 2 ε + σ 2 (X 1 )))].
The last step is to gather the terms.2 8. Proof of Lemma 6.5

If t = t 1 + t 2 with t 1 in S m and t 2 in S m , then t is such that t * has its support included in [-π max(m, m ), π max(m, m )] and therefore t belongs to S m where m = max(m, m ). We recall that B m,m (0, 1) = {t ∈ S m / t = 1}. Denote by

H 2 τ (m, m ) = (n -1 n i=1 ξ 2 i+1 σ 2 (X i ))∆(m )/n, ( 29
)
and let σ 2 τ = E(ξ 2

2 )E(σ 2 (X 1 )) = E(σ 2 (X 1 )). We have

H 2 τ (m, m ) = (n -1 n i=1 ξ 2 i+1 σ 2 (X i ) -σ 2 τ )∆(m )/n + σ 2 τ ∆(m )/n,
which is bounded by H τ,1 (m, m ) + H τ,2 (m, m ) where

H τ,1 (m, m ) = (n -1 n i=1 ξ 2 i+1 σ 2 (X i ) -σ 2 τ )1 I {|n -1 n i=1 ξ 2 i+1 σ 2 (X i )-σ 2 τ |>σ 2 τ /2}
∆(m ) n and H τ,2 (m, m ) = 3σ 2 τ ∆(m )/(2n). We infer that τ n (t) = τ

(1)

n (t) + τ (2)
n (t) with

τ (1) n (t) = 1 n n k=1 ξ k+1 σ(X k ) [u * t (Z k ) -t(X k )] , τ (2) n (t) = 1 n n k=1 ξ k+1 σ(X k )t(X k ),
and

E[ sup t∈B m,m (0,1) |τ n (t)| 2 -p τ (m, m )] + ≤ 2E[ sup t∈B m,m (0,1)
|τ (1)

n (t)| 2 -2(1 + 2 2 (m, m ))H 2 τ (m, m )] + +2E[ sup t∈B m,m (0,1) |τ (2) n (t)| 2 -6p 2 (m, m )] + +E[4(1 + 2 2 (m, m ))H 2 τ (m, m ) + 12p 2 (m, m ) -p τ (m, m )] + , where (m, m ) = 1 if 0 ≤ δ < 1/3 and (m, m ) = C(f ε )(π(m ∨ m ))
p 2 (m, m ) = E(σ 2 (X 1 )) m n .
Clearly, p 2 (m, m ) is negligible with respect to p τ (m, m ) (compare their orders in m ), so that for simplicity we consider that (12p 2 (m, m )p τ (m, m )/2) + ≤ C/n.

E[4(1 + 2 2 (m, m ))H 2 τ (m, m ) + 12p 2 (m, m ) -p τ (m, m )] + (31) ≤ 4(1 + 2 2 (m, m ))E|H τ,1 (m, m )| + E[4(1 + 2 2 (m, m ))H τ,2 (m, m ) -p τ (m, m )/2] + .
Since we only consider values of m such that the penalty are bounded by some constant K, we obtain that for some p ≥ 2, E|H τ,1 (m, m )| is bounded by

CE | 1 n n i=1 ξ 2 i+1 σ 2 (X i ) -σ 2 τ |1 I {n -1 | n i=1 (ξ 2 i+1 σ 2 (X i )-σ 2 τ )|>σ 2 τ /2} ≤ C2 p-1 E |n -1 n i=1 ξ 2 i+1 σ 2 (X i ) -σ 2 τ | p /σ 2(p-1) τ .
Moreover, we shall see below that (m, m ) is constant (if δ = 0 or 0 < δ < 1/3) or at most of order (ln(n)) δ (if δ > 1/3). According to Rosenthal's inequality (see Rosenthal (1970)) generalized to the mixing case (see Doukhan (1994) and Inequality [START_REF] Pensky | Adaptive wavelet estimator for nonparametric density deconvolution[END_REF] recalled in Lemma 9.1), we find that,

E|n -1 n i=1 ξ 2 i+1 σ 2 (X i ) -σ 2 τ | p ≤ C (p, ξ, σ(X)) n 1-p + n -p/2 .
Now, Assumption A1(i)-A5 implies that γ > 1/2, therefore |M n | ≤ √ n if δ = 0 and has logarithmic order if δ > 0 and thus, choosing p = 3 leads to m ∈Mn E|(1 + 2 2 (m, m ))H τ,1 (m, m )| ≤ C(ξ, σ(X))/n, where C(ξ, σ(X)) is a constant depending on the moments of ξ 1 and σ(X 1 ). In particular this requires that ξ admit a moment of order 8.

The last term of the inequality (31) vanishes as soon as

p τ (m, m ) = 8(1 + 2 2 (m, m ))H τ,2 (m, m ) = 12(1 + 2 2 (m, m ))E(σ 2 (X 1 ))∆(m )/n.
For this choice of p τ (m, m ), we obtain that

E[ sup t∈B m, m (0,1) |τ n (t)| 2 -p τ (m, m)] + ≤ 2 m ∈Mn E[ sup t∈B m,m (0,1) (n -1 n i=1 ξ i σ(X i )(u * t (Z i ) -t(X i ))) 2 -2(1 + 2 2 (m, m ))H 2 τ (m, m )] + +2 m ∈Mn E[ sup t∈B m,m (0,1) |τ (2) n (t)| 2 -6p 2 (m, m )] + + C n .
Then we apply the following Lemma.

Lemma 8.1. Under the assumptions on the model, if E|ξ 1 | 8 < ∞ and E(σ 8 (X 1 )), then for some given > 0:

m ∈Mn E   sup t∈B m,m (0,1) 1 n n i=1 ξ i+1 σ(X i )(u * t (Z i ) -t(X i )) 2 -2(1 + 2 2 )H 2 τ (m, m )   + ≤ K 1    m ∈Mn σ 2 τ λ 2 Γ 2 (m ) n exp -K 2 2 ∆(m ) λ 2 Γ 2 (m ) + 1 + ln 4 (n) √ n 1 n    ,
where λ 2 is a constant, Γ 2 (m) is defined by

(32) Γ 2 (m) = (m) 2γ+min[(1/2-δ/2),(1-δ)] exp{2µ(πm) δ }
and K 1 and K 2 are constants depending on the moments of ξ and σ(X).

Moreover, it also follows from [START_REF] Baraud | Adaptive estimation in an autoregressive and a geometrical β-mixing regression framework[END_REF] and [START_REF] Comte | Adaptive estimation of mean and volatility functions in (auto-) regressive models[END_REF], that the process τ [START_REF] Bennett | Errors-in-variables in joint population pharmacokinetic/pharmacodynamic modeling[END_REF] n is a standard process of the auto-regressive context and satisfies, for p 2 (m, m ) defined by [START_REF] Koo | B-spline estimation of regression functions with errors in variable[END_REF],

2 m ∈Mn E[ sup t∈B m,m (0,1) |τ (2) n (t)| 2 -6p 2 (m, m )] + ≤ c n .
We denote by

A(m) = K 1 σ 2 τ n λ 2 Γ 2 (m) exp -K 2 2 ∆(m) λ 2 Γ 2 (m) = K 1 σ 2 τ λ 2 Γ 2 (m) n exp -κ 2 2 m (1/2-δ/2) + .
The We work conditionally to the (ξ i , X i )'s and E X,ξ and P X,ξ denote the conditional expectations and probability for fixed ξ 1 , . . . , ξ n , ξ n+1 , X 1 , . . . , X n .

We apply Lemma 9.2 with f t (ξ i , X i , Z i ) = ξ i+1 σ(X i )u * t (Z i ), conditionally to the ξ i 's and X i 's to the random variables (ξ 2 , X 1 , Z 1 ), . . . , (ξ n+1 , X n , Z n ) which are independent but non identically distributed since the ξ i 's and the X i s are fixed constants.

Straightforward calculations give that for H τ (m, m ) defined in [START_REF] Klein | Concentration around the mean for maxima of empirical processes[END_REF] we have

E 2 X,ξ [ sup t∈B m,m (0,1) n -1 n l=1 ξ l+1 σ(X l )(u * t (Z l ) -t(X l ))] ≤ H 2 τ (m, m ). Let P (l) j,k (m) = E X,ξ [u * ϕ m,j (Z l )u * ϕ m,k (-Z l )]. Write sup t∈B m,m (0,1) 1 n n l=1 Var X,ξ (ξ l+1 σ(X l )u * t (Z l )) ≤ 1 n n l=1 ξ 2 l+1 σ 2 (X l )   j,k∈Z |P (l) j,k (m )| 2   1/2 .
We argue as in [START_REF] Comte | Penalized contrast estimator for density deconvolution[END_REF]. Let recall that ∆ 2 (m) is defined by [START_REF] Ioannides | Nonparametric regression with errors in variables and applications[END_REF]. We have ∆ 2 (m) ≤ λ 2 2 Γ 2 2 (m), with Γ 2 defined by (32) and λ 2 = λ 2 (γ, A 0 , δ, µ, f ε ). Now, write P (l) j,k as

P (l) j,k (m) = m 4π 2 e -ixj-iyk e im(x-y)X l ϕ * (-x)ϕ * (-y) f * ε (mx)f * ε (my) f * ε (m(x -y))dxdy.
By applying Parseval's formula we get that j,k |P

(l) j,k (m)| 2 equals ∆ 2 (m). We now write that sup t∈B m,m (0,1) (n -1 n i=1 Var X,ξ (ξ i+1 σ(X i )u * t (Z i ))) ≤ (n -1 n i=1 ξ 2 i+1 σ 2 (X i ))λ 2 Γ 2 (m ), and thus we take v τ (m, m ) = (n -1 n i=1 ξ 2 i+1 σ 2 (X i ))λ 2 Γ 2 (m ). Lastly, since sup t∈B m,m (0,1) f t ∞ ≤ 2 max 1≤i≤n |ξ i+1 σ(X i )| ∆(m ), we take M 1,τ (m, m ) = 2 max 1≤i≤n |ξ i+1 σ(X i )| ∆(m ).
By applying Lemma 9.2, we get for some constants κ

1 , κ 2 , κ 3 E X,ξ [ sup t∈B m,m (0,1) ν 2 n,1 (t) -2(1 + 2 2 )H 2 τ ] + ≤ K 1 λ 2 Γ 2 (m ) n 2 ( n i=1 ξ 2 i+1 σ 2 (X i )) exp -K 2 2 ∆(m ) λ 2 Γ 2 (m * ) + ∆(m ) n 2 ( max 1≤i≤n ξ 2 i+1 σ 2 (X i )) exp    -K 3 C( 2 ) n i=1 ξ 2 i+1 σ 2 (X i ) max i |ξ i+1 σ(X i )|     
To relax the conditioning, it suffices to integrate with respect to the law of the (ξ i+1 , X i )'s the above expression. The first term in the bound simply becomes:

σ 2 τ λ 2 Γ 2 (m ) exp[-κ 2 ∆(m )/(λ 2 Γ 2 (m ))]/n.
The second term is bounded by

(33) ∆(m ) n 2 E   (max |ξ i+1 σ(X i )| 2 ) exp   -κ 3 C( 2 ) n i=1 ξ 2 i+1 2 (X i ) max 1≤i≤n |ξ i+1 σ(X i )|     .
Since we only consider integers m such that the penalty term is bounded, we have ∆(m)/n ≤ K and the sum of the above terms for m ∈ M n and |M n | ≤ √ n is less than

K √ n E   max 1≤i≤n ξ 2 i+1 σ 2 (X i ) exp   -κ 3 C( 2 ) n i=1 ξ 2 i+1 σ 2 (X i ) max 1≤i≤n |ξ i+1 σ(X i )|     .
We need to study when such a term is less than c/n for some constant c. We bound max i |ξ i+1 σ(X i )| by m ξ,σ on the set {max i |ξ i+1 σ(X i )| ≤ m ξ,σ } and the exponential by 1 on the set {max i |ξ i+1 σ(X i )| > m ξ,σ } and by denoting µ = κ 3 C( 2 ), this yields √ n/[2 √ 2 ln(n)] and for any n ≥ 3, and for C 1 and C 2 some constants depending on the moments of ξ and σ(X), we find that

E max 1≤i≤n ξ 2 i+1 σ 2 (X i ) exp -µ n i=1 ξ 2 i+1 σ 2 (X i ) max 1≤i≤n ξ 2 i+1 σ 2 (X i ) ≤ m 2 ξ,σ E   exp(-µ n i=1 ξ 2 i+1 σ 2 (X i ) m ξ,σ )   +E max 1≤i≤n ξ 2 i+1 σ 2 (X i )1 I {max 1≤i≤n |ξ i+1 σ(X i )|>m ξ,σ } ≤ m 2 ξ,σ E exp(-µ nσ 2 τ /(2m 2 ξ,σ )) + P | 1 n n i=1 ξ 2 i+1 σ 2 (X i ) -σ 2 τ | ≥ σ 2 τ /2 +m -
E    max 1≤i≤n ξ 2 i+1 σ 2 (X i ) exp   -κ 3 C( 2 ) n i=1 ξ 2 i+1 σ 2 (X i ) max 1≤i≤n ξ 2 i+1 σ 2 (X i )      ≤ C 1 √ n + C 2 ln 4 (n) √ n 1 √ n .
Then the sum over M n with cardinality less than √ n of the terms in ( 33) is bounded by C(1 + ln(n) 4 / √ n)/n for some constant C, by using again that ∆(m )/n is bounded. 2

Appendix

As a reminder, some definitions and properties related to β-mixing sequences are given in this section. Let (Ω, A, P) be a probability space. Let Y be a random variable with values in a Banach space (B, • B ), and let M be a σ-algebra of A. (1960). Let X = (X i ) i≥1 be a strictly stationary sequence of real-valued random variables. For any k ≥ 0, the coefficients β X,1 (k) are defined by β X,1 (k) = β(σ(X 1 ), σ(X 1+k )), Let M i = σ(X k , 1 ≤ k ≤ i). The coefficients β X,∞ (k) are defined by β X,∞ (k) = sup i≥1,l≥1 sup {β(M i , σ(X i 1 , . . . , X i l )), i

+ k ≤ i 1 < • • • < i l } ,
In the paper, we do not distinguish between the two types of mixing and denote the coefficients of the process X by β(k) or β X (k). It is implicit that when only covariance inequality are involved, then the milder mixing β X,1 (k) is required, and we shall assume that stronger β X,∞ (k) mixing coefficients are used in the general case. Now, a Rosenthal-type inequality for mixing variables can be deduced from Doukhan (1994), Theorem 2 p.26 and the following result holds: We also use [START_REF] Delyon | Limit theorem for mixing processes[END_REF] covariance Inequality, successfully exploited by Viennet (1997) for partial sums of strictly stationary processes. Lastly, we recall the version of the Talagrand inequality that is required in the paper. Mention must be made that it is valid for independent but non necessarily identically distributed random variables, which is useful here when we work conditionally to one or two of the sequences.

3. Estimation of b 3 . 1 . 3 . 1 . 1 .

 31311 The steps of the estimation. First step: the estimators of = bf . The orthogonal projection of = bf on S m , m , is given by[START_REF] Butucea | Sharp optimality for density deconvolution with dominating bias[END_REF] m = j∈Z a m,j ( )ϕ m,j with a m,j ( ) = R ϕ m,j (x) (x)dx = ϕ m,j , .

3. 4 .

 4 Risk bounds for b.[START_REF] Comte | Adaptive density deconvolution with dependent inputs[END_REF] prove that f m satisfies the same inequality as ˆ m. We recall their result.

  t∈Sm γn (t) and m = arg min m∈Mn γn ( θn ) + p en(m)

  ȃn and 0 otherwise, and ȃn is a sequence to be specified in the same way as a n for the estimation of b. Clearly, b 2 + σ 2 is an estimator of b 2 + σ 2 . For the study of steps 2 and 3, see Section 3.4.Fourth step. The estimator of σ 2 must be built by settingσ2 = b 2 + σ 2 -( b) 2 .Clearly, as σ2 -σ 2 2 ≤ 2 b 2 + σ 2 -(b 2 + σ 2 ) 2 + 2 b + b 2 b -b 2 ,the risk of the final estimator is the sum of the risks of the estimators of b 2 + σ 2 and b, provided that b is bounded (by M b ) and b is bounded (by 2M b ) with probability near of one. The latter step is studied from an empirical point of view only.

Figure 1 .

 1 Figure1. ASE for the estimation of b 0 j and (σ 0 j ) 2 , j = 1, 2, 3 given in[START_REF] Fan | On the optimal rates of convergence for nonparametric deconvolution problems[END_REF], for 100 replications of the estimation procedure.

Figure 2 .Figure 3 .

 23 Figure 2. True curves: b 1 (left), b 2 1 +σ 2 1 (center), σ 2 1 (right) in dotted-bold and 20 estimated curves, n = 1000, ε Laplace.

  Figures 2 to 6 illustrate our results in the three cases of couples (b j , σ 2 j ) given in (20) with rescaling for logχ 2 (1) errors ε. The true curves are bold dotted and are given together with 20 estimated curves. This set of twenty estimated curves shows, in each case, that the confidence band around the true function has very small width. Moreover, in all cases, and for all the noise distributions, the estimation are very good. With the financial setting in mind, we took samples with large sizes n = 1000 or n = 5000 which are clearly good sizes for nonparametric estimation.We can see that b and b 2 + σ 2 are well estimated by the ratio strategy. The extraction of σ 2 sometimes suffers from scale problems (if σ 2 is much smaller that b 2 or if both are very small). In particular, we plot σ 2 with the same vertical scale as b 2 + σ 2 to take this into account.We can therefore conclude that the method we propose is very precise to estimate b and b 2 + σ 2 , and can give an interesting idea of the shape of σ 2 .

Figure 4 .Figure 5 .

 45 Figure 4. True curves: b 2 (left), b 2 2 +σ 2 2 (center), σ 2 2 (right) in dotted-bold and 20 estimated curves, n = 1000, ε Laplace.

Figure 6 .

 6 Figure 6. True curves: b 3 (left), b 2 3 +σ 2 3 (center), σ 2 3 (right) in dotted-bold and 20 estimated curves, n = 5000, Gaussian ε.

2 6. 2 .

 22 Proof of Theorem 3.1. The proof could be sketched as follows. Let us define for m, m ∈ M n , B m (0, 1) = {t ∈ S m , t = 1} and B m,m (0, 1) = {t ∈ S m + S m , t = 1}. Under the definition of m, ∀m ∈ M n , γ n ( ˆ m) + pen( m) ≤ γ n ( m ) + pen(m). For all functions s and t, γ n

  or Comte and Merlevède (2002), p.217.2 Proof of Lemma 6.4. Let e k = ε k 1 I |ε k |>n 1/4 .

  σ ) + C (p, ξ, σ(X))m 2 ξ,σ [n 1-p + n -p/2 ] + nE(|ξ 1 | r+2 )E(|σ(X 1 )| r+2 )m -r ξ,σ . Since E|ξ 1 | 8 < ∞, we take p = 3, c = 4 in Lemma 9.1, r = 4, m ξ,σ = σ τ C( 2 )κ 3

  Let P Y |M be a conditional distribution of Y given M, and let P Y be the distribution of Y . Let B(B) be the borel σ-algebra on (B, • B ). Define now β(M, σ(Y )) = E sup A∈B(X ) |P Y |M (A) -P Y (A)| The coefficient β(M, σ(Y )) is the usual mixing coefficient, introduced by Volkonskiȋ and Rozanov

Lemma 9 . 1 .≤

 91 Let (Y k ) 1≤k≤n be a sequence of centered and stationary β-mixing variables with coefficients β(k), admitting moments of order r + 1 and r ≥ 2, then if∃c ∈ 2N, c ≥ r, such that k≥1 (k + 1) c-2 β(k) 1 c+1 < +∞, F. COMTE, C. LACOUR, Y. ROZENHOLC we have the bound C(r) n 1-r [E|Y 1 | r+1 ] r/(r+1) + n -r/2 [(E|Y 1 | 3 )] r/3 .

Theorem 9 . 1 .

 91 [START_REF] Delyon | Limit theorem for mixing processes[END_REF],[START_REF] Viennet | Inequalities for absolutely regular sequences: application to density estimation[END_REF]) Let P be the distribution of Z 0 on a probability space X , f dP = E P (f ) for any function f P -integrable. For r ≥ 2, let L(r, β, P ) be the set of functions b Z :X → R + such that b Z = l≥0 (l + 1) r-2 b l,Z with 0 ≤ b l,Z ≤ 1 and E P (b l,Z ) ≤ β Z (l)We define B r as B r = l≥0 (l + 1) r-2 β Z (l). Then for 1 ≤ p < ∞ and any function b Z in L(2, β, P ), E P (b p Z ) ≤ pB p+1 , as soon as B p+1 < ∞. The following result holds for a strictly stationary absolutely regular sequence, (Z i ) i∈Z , with β-mixing coefficients(β Z (k)) k≥0 : if B 2 < +∞,there exists b Z ∈ L(2, β, ∞) such that for any positive integer n and any measurable function f ∈ L 2 (P ), we have Var n i=1 f (Z i ) ≤ 4nE P (b Z f 2 ) = 4n b Z (x)f 2 (x)dP (x).

Lemma 9 . 2 . 1 ,

 921 Let Y 1 , . . . , Y n be independent random variables, letν n,Y (f ) = (1/n) n i=1 [f (Y i )-E(f (Y i ))] and let F be a countable class of uniformly bounded measurable functions. Then for ξ 2 > 0E sup f ∈F |ν n,Y (f )| 2 -2(1 + 2 2 with C( 2 ) = √ 1 + 2 -1, K 1 = 1/6,andsup f ∈F f ∞ ≤ M 1 , E sup f ∈F |ν n,Y (f )| ≤ H, sup f ∈F 1 n n k=1 Var(f (Y k )) ≤ v.This result follows from the concentration inequality given in[START_REF] Klein | Concentration around the mean for maxima of empirical processes[END_REF] and arguments in Birgé and Massart (1998) (see the proof of their Corollary 2 page 354).

  /2 .

	2.3. The projection spaces. As projection estimators are used in all cases, we hereby
	provide a description of the projection spaces. Let us define
	ϕ(x) =	sin(πx) πx	and ϕ m,j (x) =	√	mϕ(mx -j),
	where m can be replaced by 2 m . The key point is that ϕ * (x) = 1 I [-π,π] (x). It is well
	known (see Meyer (1990), p.22) that {ϕ m,j } j∈Z is an orthonormal basis of the space of square integrable functions having a Fourier transform with compact support included
	into [-πm, πm]. Such a space is denoted by S m .

  study of A(m ) is standard in deconvolution (see[START_REF] Comte | Penalized contrast estimator for density deconvolution[END_REF]) and leads to choose 2 (m, m ) as a constant if δ ≤ 1/3 and of orderm δ-(1/2-δ/2) + if δ > 1/3, to ensure that m ∈Mn A(m ) is less than C/n.With p τ (m, m ) given in Lemma 6.5, by gathering all terms we find the result. 2 Proof of Lemma 8.1.

  ] + nE(|ξ 2 σ(X 1 )| r+2

	Again by applying Rosenthal's inequality (see Lemma 9.1), we obtain that
	E max 1≤i≤n	ξ 2 i+1 σ 2 (X i ) exp -µ	n i=1 ξ 2 i+1 σ 2 (X i ) max 1≤i≤n ξ 2 i+1 σ 2 (X i )
	≤ m 2 ξ,σ e -µ	√ nστ /( √	2m ξ,σ ) + m 2 ξ,σ	C(p, ξ, σ(X)) n p	[n + n p/2
	r ξ,σ E( max 1≤i≤n	|ξ i+1 σ(X i )| r+2 )
	≤ m 2 ξ,σ e -µ	√ nστ /( √	2m ξ,σ ) + m 2 ξ,σ 2 p σ -2p τ	E	n 1	i=1 n	ξ 2 i+1 σ 2 (X i ) -σ 2 τ	p
	+m -r ξ,σ E( max 1≤i≤n	|ξ i+1 σ(X i )| r+2 ).


A P Rinaldi
email: rinaldi@bo.ingv.it

M Todesco
email: todesco@bo.ingv.it

M Bonafede

Hydrothermal instability and ground displacement at the Campi Flegrei caldera

Ground deformation is commonly observed at active volcanoes, where it represents a reliable sign of unrest and a potential precursor of eruptive activity. The source of deformation, however, is not always unequivocally constrained. Magma ascent and differentiation are generally involved, but hydrothermal fluids may play a role, due to thermal expansion and pore

pressure acting on rocks. The identification of mechanisms driving ground displacement bears important consequences for hazard evaluation. The aim of this work is to evaluate mechanical effects associated with pressurization and heating of hydrothermal fluids. We first simulate the heat and fluid flow driven by the arrival of magmatic fluids from greater depth. Then, we calculate the rock deformation arising from simulated pressure and temperature changes within a shallow hydrothermal system. We employ a mathematical model, based on the linear theory of thermo-poro-elasticity and on a system of distributed equivalent forces. Results show that stronger degassing of a magmatic source may cause several centimeters of uplift.

Introduction

This research work stems from the recent evolution of the Campi Flegrei (CF) caldera (Italy): although the last eruption took place in 1538 AD, secondary activity has always been present and, in particular, changes of geochemical and geophysical parameters highlighted periods of unrest during the last 50 years. Ground deformation (or bradyseism) characterized these unrests, and more than 3 meters of uplift were recorded at the end of the two major episodes, in 1969-1971 and 1982-1984. A slow subsidence begun in 1985 and, since then, only minor uplifts (a few cm each) occurred in the caldera [START_REF] De Natale | A mechanical fluid-dynamical model for ground movements at Campi Flegrei caldera[END_REF]. Each uplift phase has been accompanied by seismic activity and gravity changes, and followed by variation of discharge rate and gas composition in fumaroles. Minor crises were similar to the larger ones in terms of deformation pattern, compositional change and seismicity.

Both large and small unrests have been carefully studied for their actual and potential consequences on the densely inhabited surrounding region. Early models ascribed the observed deformation to a pressure or volume change inside a magma chamber [START_REF] Bonafede | Displacement and stress fields produced by a centre of dilatation and by a pressure source in a viscoelastic half-space: application to the study of ground deformation and seismic activity at Campi Flegrei, Italy[END_REF][START_REF] Bianchi | Modeling of surface deformation in volcanic areas: The 1970-1972 and 1982-1984 crises of Campi Flegrei, Italy[END_REF][START_REF] De Natale | Geophysical and geochemical modelling of the 1982-1984 unrest phenomena at Campi Flegrei caldera (southern Italy)[END_REF]. [START_REF] Casertano | Hydrodynamics and geodynamics in the Phlegrean fields area of Italy[END_REF] made a pioneer study considering the importance of fluids in bradyseismic events. More recently, other authors recognized the effects of heating and expansion of hydrothermal fluids during an unrest [START_REF] Bonafede | Axi-symmetric deformation of a thermo-poro-elastic half-space: inflation of a magma chamber[END_REF][START_REF] Bonafede | Hot fluid migration: an efficient source of ground deformation: application to the 1982-1985 crisis at Campi Flegrei-Italy[END_REF][START_REF] De Natale | Geophysical and geochemical modelling of the 1982-1984 unrest phenomena at Campi Flegrei caldera (southern Italy)[END_REF][START_REF] De Natale | A mechanical fluid-dynamical model for ground movements at Campi Flegrei caldera[END_REF]Gaeta et al., 1998;[START_REF] Hurwitz | Hydrothermal fluid flow and deformation in large calderas: Inferences from numerical simulations[END_REF][START_REF] Hutnak | Numerical models of caaldera deformation: Effects of multiphase and multicomponent hydrothermal fluid flow[END_REF][START_REF] Orsi | Mechanical and thermofluid behaviour during unrest at the Campi Flegri caldera (Italy)[END_REF]. The occurrence of seismicity in the CF caldera (Troise et al., 1997) and gravity changes [START_REF] Bonafede | Modelling gravity variations consistent with ground deformation in the Campi Flegrei caldera (Italy)[END_REF] were also modeled. These models The use of a more sophisticated model of hydrothermal circulation showed that several aspects of the complex unrest dynamics (including deformation) are related to the intensity of magmatic degassing [START_REF] Chiodini | Magma degassing as a trigger of bradyseismic events: the case of Phlegrean fields (Italy)[END_REF]Todesco et al., 2003aTodesco et al., ,b, 2004;;[START_REF] Todesco | Modelling hydrothermal fluid circulation and gravity signals at the Phlegrean Fields caldera[END_REF].

The involvement of both magma and hydrothermal fluids in the recent evolution of the CF caldera was confirmed by further studies relating deformation and gravity data to the density of the source of deformation [START_REF] Gottsmann | Unrest at Campi Flegrei: A contribution to the magmatic versus hydrothermal debate from inverse and finite element modelling[END_REF][START_REF] Amoruso | Simultaneous inversion of deformation and gravity changes in a horizontally layered half-space: Evidences for magma intrusion during 1982-1984 unrest at Campi Flegrei caldera (Italy)[END_REF]Bonafede and Ferrari, 2008).

In this work, we want to quantify the amount of deformation that can be caused by hydrothermal fluids during a generic unrest period at Campi Flegrei.

Following a well established approach [START_REF] Todesco | Multiphase fluid circulation and ground deformation: a new perspective on[END_REF](Todesco et al., , 2004;;[START_REF] Hurwitz | Hydrothermal fluid flow and deformation in large calderas: Inferences from numerical simulations[END_REF], we first simulate the evolution of the hydrothermal system; then we compute the deformation arising from the simulated changes in pore pressure and temperature.

Our simulation of hydrothermal circulation is based on a conceptual model that was developed for Campi Flegrei in previous papers [START_REF] Chiodini | Magma degassing as a trigger of bradyseismic events: the case of Phlegrean fields (Italy)[END_REF]Todesco et al., 2003aTodesco et al., ,b, 2004;;[START_REF] Todesco | Modelling hydrothermal fluid circulation and gravity signals at the Phlegrean Fields caldera[END_REF]): the hydrothermal system is fed by fluids of magmatic origin, and unrest events correspond to periods of increased magmatic degassing.

Page 4 of 28 A c c e p t e d M a n u s c r i p t

Pore pressure and temperature changes, arising from a given unrest event, are considered here as axially symmetric distributed sources of deformation, and described in terms of a system of equivalent forces. The resulting displacement of the ground surface is then computed analytically, according to the linear theory of thermo-poro-elasticity. Results show that the simulated unrest period can generate centimeters of vertical displacement at the free surface of an elastic homogeneous half-space.

The ground deformation model

The mathematical model presented here is based on the linear theory of poro-elasticity, which describes the elastic deformation of a porous medium taking into account the flow of a hot fluid that propagates through the pores.

Here we follow the formulation proposed by [START_REF] Rice | Some basic stress diffusion solutions for fluid-saturated elastic porous media with compressible constituents[END_REF] and, given the high temperatures of the fluid, we also take into account the thermoelastic response of the rock as described by [START_REF] Mctigue | Thermoelastic response of fluid-saturated porous rock[END_REF]. The deformation e ij due to a given change in pore pressure and temperature is:

e ij = 1 2µ σ ij - ν 1 + ν σ kk δ ij + ∆p 3H δ ij + α s 3 ∆T δ ij (1)
where σ ij is the stress tensor; σ kk is the trace of the stress tensor; µ is the shear modulus; ν is the Poisson's ratio in free-drainage conditions; 1/H = 1/K -1/K s is the Biot's constant (K is the isothermal, drained bulk modulus and K s is the bulk modulus of the solid constituent); α s is the volumetric thermal expansion coefficient for the solid matrix; δ ij is the Kronecker delta.

∆p and ∆T are the pore pressure and temperature changes, respectively.

In order to evaluate the displacement field arising from pressure and temperature changes at a given point, we follow the approach proposed by Aki We first consider a porous medium in an initial configuration of vanishing stress and strain. We identify a volume source, with dimensions dx, dy, dz, and remove it from its surrounding, without affecting the porous matrix.

We then increase pore pressure and temperature in the source volume by supplying heat and fluid, at constant (vanishing) stress. At this stage (A), the source remains unstressed (τ

(A) ij = 0
), but its volume is isotropically strained ("Aki's stress-free strain") by an amount:

e (A) ij = 1 3 ∆θδ ij
where ∆θ is the relative volume change, given by the trace of the strain tensor (1)

∆θ = ∆V V 0 = e (A) kk = ∆p H + α s ∆T
In order to place the source back in its original position in the porous medium (stage B), we need to restore its original volume V 0 : we apply a stress field τ

τ ij = τ (A) ij + τ (B) ij = τ (B) ij , e ij = e (A) ij + e (B) ij = 0
Now the source is back in its place, but a traction discontinuity exists over its surface: the matrix is still unstressed and unstrained while the source is subject to the artificial "stress glut" τ ij applied to restore its initial volume.

The difference between the outer and the inner values of traction is then

Page 6 of 28

A c c e p t e d M a n u s c r i p t -τ ij n j . Removal of this traction discontinuity will cause the source to expand again but, this time, the surrounding matrix acts against the expansion and prevents the source from reaching the stress-free volume. The displacement field due to the applied pressure and temperature changes at the source can be expressed in terms of this traction discontinuity [START_REF] Aki | Quantitative Seismology -Theory and Methods[END_REF]:

u i (x) = V K∆θ(x )G ik,k (x, x ) dx dy dz (2)
where K is the bulk modulus, x is the observation point, x is the source position and G ik is the Green's Tensor.

In order to evaluate the divergence of the Green's Tensor we use a system of distributed equivalent forces. The isotropic point source can be mathematically described as a system of three equal orthogonal dipoles placed at the source point (Wang, 2000). The displacement field G ik arising from a single force in a homogeneous half-space with a traction-free boundary was provided by [START_REF] Mindlin | Force at a point in the interior of a semi-infinite solid[END_REF].

In this work, we first calculate the displacement due to a single dipole set along the x, ŷ or ẑ axes, and then we sum the single components of each field in order to have the total displacement at the observation point. The vertical component of the displacement for this case is:

u tot z = F h (1 -2ν) 8πµ(1 -ν) (z -z ) R 3 1 + -(3 -4ν) z + z R 3 2 + 2z R 3 2 - 6z(z + z ) 2 R 5 2 (3)
where

F
= (x -x ) 2 + (y -y ) 2 + (z -z ) 2
) and R 2 is the distance between the mirror-source and

the observation point (R 2 2 = (x -x ) 2 + (y -y ) 2 + (z + z ) 2 ).
Since we are using a free surface calculation, the horizontal component of the displacement will be different than the vertical:

u tot r = F h (1 -2ν) 8πµ(1 -ν) 1 R 3 1 + -(3 -4ν) 1 R 3 2 - 6z(z + z ) R 5 2 (r -r ) (4)
where r 2 = x 2 + y 2 and r 2 = x 2 + y 2 represent the radial distance of the observation point and of the source from the ẑ axes, respectively.

The fluid flow model

Hydrothermal circulation is simulated with the multi-phase and multicomponent TOUGH2 model [START_REF] Pruess | TOUGH2 user's guide, version 2.0[END_REF], which describes the coupled flow of heat and fluids through the porous medium. In the present application, the considered fluid components are water and carbon dioxide.

The computational domain, two-dimensional and axisymmetric, is 10 km wide and 1.5 km deep (Fig. 1 shows the axial region only, since most of changes happen near the symmetry axis).

Bottom and side boundaries are impervious and adiabatic. Atmospheric conditions are fixed along the upper boundary, which is open to heat and fluid flows. The properties of the porous medium are listed in Table 1.

The shallow hydrothermal circulation is driven by the injection of a hot (ca. collected at the CF caldera, and corresponds to a CO 2 /H 2 O molar ratio of 0.17 [START_REF] Chiodini | Magma degassing as a trigger of bradyseismic events: the case of Phlegrean fields (Italy)[END_REF]. As described in previous works [START_REF] Chiodini | Magma degassing as a trigger of bradyseismic events: the case of Phlegrean fields (Italy)[END_REF]Todesco et al., 2003aTodesco et al., ,b, 2004)), the prolonged activity of the fluid source generates, at the steady state, a wide two-phase plume, with a shallow and a deep single-phase gas regions region (Fig. 1a). High temperature characterizes the entire plume (Fig. 1b), which is also slightly pressurized with respect to the hydrostatic gradient.

Starting from these steady-state initial conditions, the simulation presented here describes an initial unrest phase, provided by a sudden increase of fluid flow, followed by a longer quiet period. The unrest phase lasts 20 months, during which both the input of magmatic fluids and the carbon dioxide content increase (6000 tons/day of CO 2 and 6100 tons/day of H 2 O, CO 2 /H 2 O=0.4) at same, constant input temperature (ca. 623 K). During the following quiet phase, flow rate and composition return to the initial values.

Inlet conditions corresponding to the different phases are listed in Table 2. A sequence of similar unrest and quiet periods was shown to be consistent with observed changes in gas composition and gravity at the CF caldera during the last 30 years [START_REF] Chiodini | Magma degassing as a trigger of bradyseismic events: the case of Phlegrean fields (Italy)[END_REF][START_REF] Todesco | Modelling hydrothermal fluid circulation and gravity signals at the Phlegrean Fields caldera[END_REF].

The imposition of a larger discharge rate during the unrest is accompanied by a pressure build-up, which gradually spreads out as the fluids propagate

Page 9 of 28 A c c e p t e d M a n u s c r i p t through the system. At the end of the unrest, pressure changes range from +5 MPa, near the fluid inlet, to +0.1 MPa near the ground surface (Fig. 2a).

The unrest is also associated with temperature changes, even if the specific enthalpy of injected fluids does not changes. Temperature changes, up to tens of Kelvin, are mostly confined along the edges of the two-phase plume (Fig. 2b). Most of these temperature changes are related to the lateral spreading of the two-phase region during the unrest. As the hot, gas-rich fluids replace the colder liquid water, temperature may increase up to 40 K (Fig. 3a). Temperature changes within the two-phase plume may occur as well, as a consequence of water phase changes: pressure build-up near the inlet can cause the condensation of water vapor (Fig. 3a) and the associated release of latent heat increases the temperature in the region above the fluid source (+5 K). At shallower depths, some fraction of liquid water in the twophase region evaporates (Fig. 3a), subtracting latent heat and causing minor cooling (ca. -1.5 K).

When the unrest is over, pore pressure drops and the two-phase plume begins to shrink gradually. Decompression enhances water evaporation and associated cooling. Fluids entered during the unrest still propagate upwards, affecting shallower portions of the domain during the quiet (Fig. 3b). At the end of the simulation, pressure and temperature anomalies persist: localized overpressure (up to +0.1 MPa) is still present at shallow depth (< 200 m), while most of the domain undergoes decompression (up to -0.3 MPa), particularly near the plume border (Fig. 2c). Temperature anomalies are mostly confined along the edges of the two-phase plume (Fig. 2d). A slow temperature decline occurs at depth, where heating (+20 K) still persists

Page 10 of 28

A c c e p t e d M a n u s c r i p t at the end of the simulation. Moderate heating (a few K) takes place at shallower depths, as hot fluids rise toward the surface (Fig. 3). The heated region is not uniform along the plume due to the competing effects of fluid displacement and phase changes during the simulation (Fig. 2d). It is interesting to note that pressure and temperature changes affect different regions of the domain, which do not necessarily reflect the position of the inlet of magmatic fluids (i.e. the location of the "source" that drives the unrest, Fig. 1).

Resulting displacement

To compute the vertical ground displacement arising from the simulated evolution, each element of the computational domain is considered as a potential source of deformation. Pressure and temperature changes at different times are calculated with respect to the initial pressure and temperature values in each element. The corresponding vertical displacement can be computed from the contribution of each element (el) of the grid (5):

u el z = - (1 + ν) 12π(1 -ν) ∆θ el V el (z -z el ) R 3 1 + -(3 -4ν) z + z el R 3 2 + 2z R 3 2 - 6z(z + z el ) 2 R 5 2 (5)
where ∆θ el = ∆p el /H + α s ∆T el represent the dilatation (compression) of a single element of the grid. As in previous papers (Todesco et al., 2004;[START_REF] Hurwitz | Hydrothermal fluid flow and deformation in large calderas: Inferences from numerical simulations[END_REF], at this time we do not account for the effects of rock deformation back onto fluid propagation (full coupling). Rock deformation is expected to affect fluid circulation by altering the value and distribution of permeability. Nevertheless, major effects are expected only when pore pressure approaches lithostatic values (Zencher et al., 2006).

The temporal evolution of ground displacement, calculated at the surface, on the symmetry axis, is shown in Figure 4a. The values for thermo-poroelastic parameters are listed in Table 3. Uplift begins as soon as the injection rate is increased, at the beginning of the unrest, and reaches the maximum value (9 cm) at the end of the crisis. The quiet period is characterized by a slow subsidence that reflects the lower inflow of magmatic fluids. The minimum ground elevation (-1 cm) is attained at the end of the simulation. This trend of deformation arises from the temporal evolution of pressure and temperature anomalies (Fig. 4b). The faster uplift phase is due to the combined effects of increasing average pressure and temperature. When the unrest is over, the average pressure quickly declines, but the average temperature keeps increasing until hot fluids begin to discharge at the surface. The effects of fast decompression are therefore mitigated by thermal expansion, and the resulting subsidence rate is slower than the uplift rate. After ca. 60 months

Page 12 of 28 A c c e p t e d M a n u s c r i p t of simulation, the average pore pressure drops below its initial value and later on declines at a smaller rate. As a result, the ground level falls below its initial elevation, and then tends to stabilize at the end of the simulation.The pressure drop is associated with an overall decline of fluid density, mainly due to the injection of larger amounts of carbon dioxide during the unrest, and to the evaporation of liquid water associated with decompression.On a larger time scale, the system tends to return to its initial configuration and therefore to restore the initial ground level, which correspond to the steady activity of the fluid source during quiet times.

The radial distribution of ground uplift, calculated at different times, is shown in Figure 5. At the end of unrest, the pattern of deformation is similar to a Mogi-type source (Mogi, 1958): the maximum uplift is right above the source and deformation is mostly confined within 4 km from the symmetry axis (Fig. 5a). A different pattern characterizes the subsidence. With time, changes along the border of the plume become shallower (Fig. 2), and some spikes in vertical displacement develop at 500-700 m from the symmetry axis.

At the end of simulation (Fig. 5b), minor subsidence affects the axial region, and vertical displacement vanishes within the first 2 km. For K values above 15 GPa, the initial ground elevation is not restored within the simulation time. Lower bulk moduli lead to larger vertical displacement during the unrest, and drive ground subsidence below the initial elevation.

Differences in ground displacement scale quite linearly with increasing bulk modulus during the uplift: 5 GPa changes in bulk modulus produce about 1,5 cm difference in ground displacement (Fig. 6a). This linearity is conserved to the end of simulation, although following a different pattern. The shear modulus of the rocks also has a substantial effect on ground deformation.

Low rock rigidity leads to larger vertical displacement (up to 18 cm) and faster subsidence during the quiet. In this case, displacements do not scale linearly with increased rigidity (Fig. 6b).

Bulk modulus K s changes, within the range 20-40 GPa, and thermal expansion coefficient changes, within the range 10 -7 -10 -5 K -1 , have minor effect on the maximum calculated ground displacements.

Conclusion

In this work we present a model of ground deformation based on the linear theory of thermo-poro-elasticity. The model was applied to compute the vertical ground displacement associated with heating and pressurization of hydrothermal fluids during an unrest event at the CF caldera. The evolution of the hydrothermal system was simulated with the TOUGH2 numerical code [START_REF] Pruess | TOUGH2 user's guide, version 2.0[END_REF]. Then, the corresponding ground deformation was calculated under the assumption of homogeneous mechanical properties of the porous medium. Our approach does not require the definition of boundary

Page 14 of 28

A c c e p t e d M a n u s c r i p t conditions for the displacement along the periphery of the computational domain. This is particularly relevant in a small system such as ours, where the presence of a shallow bottom boundary certainly affects the mechanical response of the free surface.

Results show that increasing the fluid injection rate in our system (by a factor 3.5) leads to a maximum uplift of the order of 10 cm, which may double depending on the choice of mechanical properties of the porous rock.

A slow subsidence takes place during the quiet phase, when hot fluids reach shallower depths and heating mitigates the effects of pressure drop, leading to a slower evolution.

The evolution of the two-phase region during the quiet is associated with a complex pattern of ground deformation along the surface, with peak values characterizing the edges of the plume. A complex radial pattern of deformation was already emphasized by Todesco et al. (2004) A direct comparison with previous results [START_REF] Todesco | Multiphase fluid circulation and ground deformation: a new perspective on[END_REF](Todesco et al., , 2004) ) is not possible, due to a different choice of rock properties and inlet conditions.

Nevertheless, results obtained by [START_REF] Todesco | Multiphase fluid circulation and ground deformation: a new perspective on[END_REF], with a fivefold increase of the flow rate at the inlet, indicate a maximum uplift (15 cm) of the same order of magnitude.

The calculated displacements reflects the flow rate at the inlet, and the permeability of the system, which directly controls the pressure build-up and Initial conditions also play a role: the deformation depends on pressure and temperature changes with respect to the initial conditions. The same input of magmatic fluids entering a colder region would cause much larger vertical displacement. If the system becomes heated by subsequent unrest events, as in the case of CF, ground deformation is expected to become progressively less pronounced over time, even if the input of magmatic fluids is unchanged. The remarkable difference between the large and the small unrest events at the CF caldera could be due to different sources of deformation acting at different times, but it could also derive from significant changes in rock properties and system conditions after the first major episodes.

With the present choice of initial conditions, feeding rate, and rock properties, based on a sound conceptual model that was developed for CF over time, the calculated displacement is consistent with the deformation observed during the recent, minor uplift events at CF. Our results therefore support the importance of hydrothermal fluids in these recent episodes. The strong influence of system properties and conditions on the calculated deformation pattern should be carefully considered when trying to infer properties of the deformation source from the observed pattern of ground displacement.

Page 16 of 28 Specific Heat 1000 J kg -1 K -1 .

  describe the maximum uplift (1.8 m) observed in 1982-84, and authors sug-Page 3 of 28 A c c e p t e d M a n u s c r i p t gest that the following subsidence could be due to lateral diffusion of fluids. Using a 1D model, Gaeta et al. (2003) explain minor uplifts as a consequence of hydrothermal fluid circulation. Although these models provide important conceptual insights, they are based on very simplified descriptions of the fluid dynamics and of the embedding medium.

  and h are the intensity and arm of each dipole; (x, y, z) is the observation point and (x , y , z ) is the source point; R 1 is the distance Page 7 of 28 A c c e p t e d M a n u s c r i p t between the deformation source and the observation point (R 2 1

  623 K) mixture of water and carbon dioxide. This mixture represents the Page 8 of 28 A c c e p t e d M a n u s c r i p t magmatic component that enters the domain near the symmetry axis (Fig. 1). Initial conditions are obtained by simulating a long-lasting (i.e.: thousands of years) injection of magmatic fluids. The prescribed flow rate at the inlet (1000 tons/day of CO 2 and 2400 tons/day of H 2 O) reflects data

  V el and z el are the volume and depth of the element. R 1 and R 2 are the distance of the observation point from the element and from the image-element. The total displacement is calculated by summing the contributions of all grid block. Contrary to previous works, our approach for the computation of ground deformation does not require the specification of arbitrary boundary conditions over the boundary of the computational domain, since Mindlin's solutions provide vanishing stress over the free surface and vanishing displacements and stress at remote distance.

  The temporal evolution of vertical displacement does not depend only on the dynamics of the hydrothermal fluids. The mechanical properties of the porous rock control both the magnitude and the temporal evolution of ground deformation. This is a relevant aspect, as mechanical properties may change significantly with rock type and physical condition, and site-specific values are often poorly constrained. The effects of different values of bulk and shear moduli are shown in Figure 6. Increasing the bulk modulus reduces the Page 13 of 28 A c c e p t e d M a n u s c r i p t uplift during the unrest, and leads to a slower subsidence during the quiet.

  heat propagation. Stronger inputs of magmatic fluids would cause larger uplifts. Similarly, different choices of rock permeability (or any change in permeability during the unrest) are going to affect system conditions and its response to unrest periods. Higher permeability favors fluid propagation, and Page 15 of 28 A c c e p t e d M a n u s c r i p t hence heating of the system but, at the same time, reduces pressure build up. On the contrary, low permeability hinders fluid propagation and favors the pressure increase, which however would remain confined within a smaller region. The overall outcome in terms of ground displacement is difficult to estimate and further research should focus on this aspect.

  the Phlegrean Fields (Italy). In: Proceedings of the 28 th Workshop on Geothermal Research Engineering, Stanford, CA, USA. Todesco, M., J. Rutqvist, G. Chiodini, K. Pruess, and C. M. Oldenburg (2004), Modeling of recent volcanic episodes at Phlegrean Fields (Italy): geochemical variations and ground deformation, Geothermics, 33, 531-547. Troise, C., G. De Natale, and F. Pingue (1997), A model for earthquake generation during unrest episodes at Campi Flegrei and Rabaul calderas, Geophys. Res. Lett., 24(13), 1575-1578. Wang, H. F. (2000), Theory of Linear Poroelasticity with Applications to Geomechanics and Hydrology, 276 pp., Princeton University Press, Princeton and Oxford. Zencher, F., M. Bonafede, and R. Stefansson (2006), Near-lithostatic pore pressure at seismogenic depths: a thermoporoelastic model, Geophys. J.Int., 166, 1318Int., 166, -1334. .

Figure 1 :

 1 Figure 1: Initial conditions: a) volumetric gas fraction and b) temperature (K).The star indicates the inlet of magmatic fluids, which has a radius of 150 m. The computational domain was discretized into 2580 elements, with radial dimensions ranging from 25 to 3196 m and thickness from 5 to 25 m.

Figure 2 :

 2 Figure 2: Simulated changes with respect to the initial conditions at different times: a,c) pore pressure (MPa) and b,d)temperature (K) changes at the end of unrest and at the end of simulation, respectively. The figures refer to axial region. The star indicates the inlet of magmatic fluids.

Figure 3 :

 3 Figure 3: Simulated changes with respect to the initial conditions: volumetric gas fraction at the end of unrest (a) and at the end of simulation (b). Arrows describe the pattern of flow for the liquid (black) and gas phase (magenta). Flow velocity ranges from 0 to 1.25 g s -1

Figure 4 :Figure 5 :

 45 Figure 4: (a) Temporal evolution of vertical ground displacement at the symmetry axis.(b) Temporal variation of the average temperature (solid line) and average pore pressure (dashed line). Values are computed over the entire domain and normalized with respect to the maximum value. Shaded area in both figures highlights the unrest period. Rock mechanical properties: µ=2 GPa; K=5 GPa; K s =30 GPa; α=10 -5 K -1

Figure 6 :

 6 Figure 6: Temporal evolution of vertical ground displacement at the simmetry axis for different values of the rigidity (a) and bulk modulus (b). The other parameters are as listed in Table3. Shaded are in both figures highlights the unrest period.

  Figure 1

  Figure 2

Table 1 :

 1 Rock properties considered in the simulation of the hydrothermal system. These properties remain constant during the simulation.

		CO 2 H 2 O Molar Ratio
	Quiet	1000 2400	0.17
	Unrest 6000 6100	0.40

Table 2 :

 2 Flow rate (ton/day) and CO 2 /H 2 O molar ratio at the inlet during the unrest and the following quiet. The enthalpy of injected fluids corresponds to a temperature of ca. 623 K and does not change during the simulation.

	µ	2 GPa
	K	5 GPa
	K s	30 GPa
	α s 10 -5 K -1

Table 3 :

 3 Rock mechanical properties considered in the computation of the vertical ground displacement. µ is the shear modulus, K is the drained bulk modulus, K s is the bulk modulus of the solid constituent and α s is the volumetric thermal expansion coefficient for the solid matrix.

	Page 22 of 28

Acknowledgements

We wish to thank Andrew Newman, whose careful review greatly improved the manuscript.