Patrick Cattiaux

Jose Rafael León

Angie Pineda Centeno

Clémentine Prieur
email: clementine.prieur@imag.fr

A A Pineda Centeno

An overlook on statistical inference issues for stochastic damping Hamiltonian systems under the fluctuation-dissipation condition

Keywords: hypoelliptic diffusion, variance estimation, fluctuation-dissipation models AMS Subject Classification: Primary 62M05, secondary 60H10, 60F05, 35H10

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The stochastic damping hamiltonian system is defined as the following system of SDE

dX t = Y t dt dY t = σ(X t , Y t ) dW t -(c(X t , Y t )Y t + ∇V (X t ))dt.
(1.1)

We denote the vector solution as Z t = (X t , Y t ) ∈ R 2d . The matrix c is called the damping force and V the potential, the symmetric matrix σ(x, y) is the diffusion term and W is a standard Brownian motion defined in R d .

In the present paper, we will study some inference issues in the following particular case:

dX t = Y t dt dY t = (2β -1 ) 1/2 σ(X t )dW t -(σ 2 (X t )Y t + ∇V (X t ))dt. (1.2)
Observe that the diffusion term in (1.2), σ(x) = (2β -1 ) 1/2 σ(x), depends only on the x coordinate and on an unknown parameter β. Moreover, the damping force has the form c(x, y) = β 2 σ2 (x). Since Einstein, this last relationship between the damping force and the diffusion term is known as the fluctuation-dissipation relation and has numerous applications. For instance, it appears as a tool for the simulation of molecular dynamics (see [START_REF] Horenko | Likelihood-Based Estimation of Multidimensional Langevin Models and its Application to Biomolecular Dynamics[END_REF] and references therein and also the recent book by Lelièvre et al. [START_REF] Lelièvre | Free energy computations: a mathematical perspective[END_REF]Section 2.2.3]). System (1.2) also appears as limit of the Ehrenfest nuclei dynamics and is called Langevin dynamics (see [START_REF] Szepessy | Langevin molecular dynamics derived from Ehrenfest dynamics[END_REF]). Under additional assumptions (see Section 2 for details), we can prove that the Markov process solution of (1.2) has an invariant measure whose density is proportional to the Boltzmann distribution exp(-βH(x, y)) where H(x, y) = 1 2 |y| 2 + V (x). In this work, by following the spirit of our recent papers [START_REF] Cattiaux | Estimation for Stochastic Damping Hamiltonian Systems under Partial Observation. I. Invariant density[END_REF][START_REF] Cattiaux | Estimation for Stochastic Damping Hamiltonian Systems under Partial Observation. II. Drift term[END_REF][START_REF] Cattiaux | Recursive Estimation for Stochastic Damping Hamiltonian Systems[END_REF][START_REF] Cattiaux | Estimation for Stochastic Damping Hamiltonian Systems under Partial Observation. III. Diffusion term[END_REF], we aim to propose convergent estimators of the various unknown quantities involved in the formulation (1.2) of fluctuation-dissipation stochastic damping hamiltonian systems. The fluctuationdissipation relation allows to simplify the assumptions ensuring the ergodicity of the process solution of (1.2). These assumptions are stated in Section 2. In Section 3, we detail two classes of examples. We first provide an explicit construction of the Lyapunov function in the case of hypoelliptic diffusions on R 2d under a drift condition. We then consider the class of periodic potentials defined on the torus. In Section 4, we state our results concerning the inference of such models. Section 5 is devoted to numerical experiments.

Hypotheses and notation

In this section, we precise the context of our study. We will state some hypotheses implying the existence and uniqueness of the solution of (1.2). These hypothesis also imply the existence and unicity of an invariant probability measure with smooth density with respect to the Lebesgue measure. Finally, an additional assumption is provided to ensure the exponential ergodicity of the underlying process, which is crucial for establishing our estimation results. Before stating the hypothesis let us introduce some notation. We will denote by C ∞ (A, B) the infinitely differentiable functions f : A → B where A and B are finite dimensional normed spaces if B = R we denote simply C ∞ (A). The symbol M(R m , R m ) denotes the normed space of square matrices. We assume that our system satisfies the following hypotheses.

• H 0 the diffusion matrix σ is symmetric, C ∞ (R d , M(R d , R d ))
, bounded as well as its first and second partial derivatives and uniformly elliptic, i.e. ∀x ∈ R d , σ(x) ≥ σ 0 Id (in the sense of quadratic forms) for a positive constant σ 0 > 0; • H 1 the potential V is lower bounded and C ∞ (R d ).

Remark 2.1 Let σ * denote the transposed matrix of σ. Actually the law of the process solution of (1.2) depends on σσ * (which is the second order term of the infinitesimal generator). If this symmetric matrix is smooth, it is well known that one can find a smooth symmetric square root of it, which is the choice we make for σ, hence σσ * = σ 2 . In the following, we use the notation σ(x) = (σ ij (x)) 1≤i,j≤d .

Remark 2.2 Given the particular form of the damping term, we have that c(x, y) is a symmetric and uniformly elliptic matrix.

Let Ω = C(R + , R 2d ) denote the space of continuous functions from R + to R 2d , equipped with the usual compact convergence topology. According to Wu [13, Lemma 1.1], Hypotheses H 0 , H 1 and Remark 2.2 entail that for every initial state z = (x, y) ∈ R 2d , the stochastic differential equation (1.2) admits a unique weak solution P z on Ω, and that this solution is non-explosive. Let (P t (z, dz ) , t > 0) be the semigroup of transition probability kernels of the Markov process ((Z t ) t>0 , (P z ) z∈R 2d ) solution of (1.2). Let us prove that P t is strong Feller for each t > 0. Let C 1,2 (R 2d ) be the space of all functions f (x, y) such that ∂ k x f , k = 0, 1 and ∂ l y f , l = 1, 2 are continuous on R 2d . The infinitesimal generator L of (1.2) is given by: for any

f ∈ C 1,2 (R 2d ), Lf (x, y) = 1 β d i, j=1 (σ 2 ) ij (x)∂ yi ∂ yj + y • ∇ x -(σ 2 (x)y + ∇V (x)) • ∇ y f (x, y) . (2.3)
Let us define the differential operators A k , k = 0, . . . , d as

A 0 = d j=1 y j • ∂ xj - d j=1 ∂ xj V (x)∂ yj - d j=1 y j d l=1 (σ 2 ) jl (x)∂ yl A k = d j=1 σ kj (x)∂ yj k = 1, . . . , d.
We can decompose L in the Hörmander form, i.e. L = 1

β d k=1 A 2 k +A 0 . Let [A, B
] denote the Lie bracket between the operators A and B. As σ 2 is uniformly elliptic, it is possible to prove that (see [START_REF] Lelièvre | Free energy computations: a mathematical perspective[END_REF] pages 91-92 for more details):

Span(A 0 , .., A d , [A 0 , A 1 ], .., [A 0 , A d ]) = Span(∂ x1 , .., ∂ xd , ∂ y1 , .., ∂ yd ).
Then, according to Hörmander's theorem on hypoellipticity, we deduce that P t is strong Feller. Moreover, if a stationary probability distribution µ having a positive density with respect to the Lebesgue measure exists, then the irreducibility property is verified, the stationary probability distribution is unique, and pathwise ergodicity holds for any initial consition z ∈ R 2d (see [START_REF] Lelièvre | Free energy computations: a mathematical perspective[END_REF]Theorem 2.7] and references therein). Let us now consider the probability measure dµ(x, y) = p β s (x, y)dxdy with p β s (x, y) = C(β)e -βH(x,y) > 0 .

(2.4)

It is possible to prove that µ is a stationary probability distribution. Indeed, for any smooth test function ϕ:

R 2d

Lϕ(x, y)p β s (x, y)dxdy

= 1 β R 2d (∇ y p β s ∇ x ϕ -∇ x p β s ∇ y ϕ)dxdy + 1 β R 2d div y (p β s σ 2 (x)∇ y ϕ)dxdy = 0.
Thus L * (µ) = 0 and we conclude that the probability measure µ is an invariant measure. We refer to [START_REF] Lelièvre | Free energy computations: a mathematical perspective[END_REF] page 90 for more details on the above calculation.

In the following, we reinforce Hypotheses H 0 and H 1 by assuming that the stationary solution of (1.2) is exponentially β-mixing (see Hypothesis H 2 below).

• H 2 there exist two constants C > 0, λ > 0 such that for all t > 0,

R 2d sup f measurable, |f |≤1 P t f (z) - R 2d f dµ µ(dz) ≤ Ce -λt .
(2.5)

Assumption H 2 will allow the use of adequate covariance inequalities for the proof of inference results stated in this paper (see [START_REF] Doukhan | Mixing. Properties and examples[END_REF] and references therein).

In the following we present two classes of examples for which the exponential ergodicity is satisfied, by an explicit construction of the Lyapunov function.

Two important classes of examples

In this section, we consider two classes of examples for which exponential ergodicity is proved.

Hypoelliptic diffusion on R 2d under a drift condition

Let us consider System (1.2) with Assumptions H 0 and H 1 . We know from Section 2 that the process Z solution of (1.2) has a unique invariant probability measure described by (2.4). The infinitesimal generator L of System (1.2) is given by (2.3). In the sequel we assume that the potential satisfies condition (0.5) of [START_REF] Wu | Large and moderate deviations and exponential convergence for stochastic damping Hamiltonian systems[END_REF]:

< ∇V (x), x x >→ ∞ a x → ∞ with • the usual euclidean norm on R d . (3.6)
From [START_REF] Wu | Large and moderate deviations and exponential convergence for stochastic damping Hamiltonian systems[END_REF]Theorem 4.1], we could then directly conclude that exponential convergence holds. However, from a pedagogical point of view, we think it is worth detailing the different steps leading to the construction of an appropriate Lyapunov test function. The construction of the Lyapunov function we detail below follows the main lines of the proof of Theorem 4.1 in [START_REF] Wu | Large and moderate deviations and exponential convergence for stochastic damping Hamiltonian systems[END_REF], but is slightly different. For sake of simplicity in the notation, we assume without loss of generality that β = 2 in the sequel. We first introduce a function G : R d → R d defined as:

G(x) = Φ( x ) x x for x = 0 ∈ R d , G(0) = 0 ,
with Φ : R + → R + a non decreasing smooth function equal to zero on a small neighborhood of 0 and equal to 1 for x ≥ 1. By construction we have G ∈ C ∞ b (R d , R d ) and lim x →∞ ∂ xj G i (x) HS = 0, where • HS denotes the Hilbert-Schmidt norm of operators. We can now propose the folllowing form for the Lyapunov function:

Ψ(x, y) = e F (x,y)-inf R d F ≥ 1 , where F (x, y) = a( 1 2 ||y|| 2 + V (x)) + b < G(x), y > .
The constants a and b will be calibrated later on. Following the proof of [13, Theorem 3.1], we can prove that:

- LΨ Ψ = -LF - 1 2 ||σ(x)∇ y F || 2 .
October 12, 2016

Statistics: A Journal of Theoretical and Applied Statistics 11˙10˙16Fluctuationdisipation

Then, as

L a( 1 2 ||y|| 2 = a 1 2 d i=1 (σ 2 ) ii (x)-< σ 2 (x)y, y > -< y, ∇V (x) > L (aV (x)) = a < ∇V (x), y > L (b < G(x), y >) = b d j=1 d i=1 (∂ xj G i (x))y i y j -b < y, σ 2 (x)G(x) > -b < ∇V (x), G(x) > we get L(F (x, y)) = a 1 2 d i=1 (σ 2 ) ii (x)-< σ 2 (x)y, y > +b   d j=1 d i=1 ∂ xj G i (x)y i y j -< y, σ 2 (x)G(x) > -< ∇V (x), G(x) >   . Moreover 1 2 σ(x)∇ y F 2 = 1 2 σ(x) (ay + bG(x)) 2 .
Hence

- LΨ Ψ = -a 1 2 d i=1 (σ 2 ) ii (x)-< σ 2 (x)y, y > -b   d j=1 d i=1 ∂ xj G i (x)y i y j -< y, σ 2 (x)G(x)) > -< ∇V (x), G(x) >   - 1 2 σ(x) (ay + bG(x)) 2 .
The aim is now to apply [6, Theorem 5.2c] (see also [START_REF] Wu | Large and moderate deviations and exponential convergence for stochastic damping Hamiltonian systems[END_REF]Theorem 2.4]). We already know from Section 2 that we have the strong Feller property and that there exists a unique invariant probability measure whose density p β s (z) is positive on R d . By construction we have Ψ ≥ 1, and from Hypotheses H 0 and H 1 we know that Ψ is smooth. Thus, to prove Hypothesis H 2 , it remains to prove the existence of some compact subset K ⊂ R 2d and of some constants ε , C > 0 such that

- LΨ Ψ ≥ ε1 K c -C1 K . (3.7)
Proof of Inequality (3.7): From Hypothesis H 0 we know that there exists two positive constants σ 0 and σ Let us first write the following bound: for all

1 such that ∀x ∈ R d , σ 0 Id ≤ σ(x) ≤ σ 1 Id
x ∈ R d , σ(x)(ay + bG(x)) 2 ≤ σ 2 1 (2a 2 y 2 + 2b 2 G(x) 2 ) = σ 2 1 (2a 2 y 2 + 2b 2 ) . (3.8)
We then split the proof in two parts.

October 12, 2016

Statistics: A Journal of Theoretical and Applied Statistics 11˙10˙16Fluctuationdisipation

Step 1: x ∈ R d , y ≤ L . Let us first prove that there exists some positive constant ε such that:

sup x∈R d aσ 2 (x) -b∂ xj G i (x) ≥ (σ 2 0 a 2 -ε )Id for all y ∈ R d .
The above inequality is to be understood in the sense of quadratic forms. We have:

< (aσ 2 (x) -b∂ xj G i (x))y, y >≥ σ 2 0 a y 2 -b < ∂ xj G i (x)y, y > .
We have:

lim x →∞ ∂ xj G i (x) HS = 0 . (3.9)
We can thus choose b such that b sup

x∈R d ∂ xj G i (x) HS ≤ ε .
We then get:

sup x∈R d < (aσ 2 (x) -b∂ xj G i (x))y, y >≥ (aσ 2 0 -ε ) y 2 . (3.10)
Then, from the definition of G and from Assumption (3.6) we get the existence of some positive constants M , M such that:

sup x∈R d < ∇V (x), G(x) >≥ -M (3.11)
and sup

x∈R d b < y, σ 2 (x)G(x) > ≥ -M y . (3.12)
Moreover we have: for all

x ∈ R d , -a 1 2 d i=1 (σ 2 ) ii (x) ≥ -a d 2 σ 2 1 , a < σ 2 (x)y, y >≥ aσ 2 0 y 2 . (3.13)
Then, from (3.8), (3.10), (3.17), (3.12) and (3.13), we get, for all x, y ∈ R d :

- LΨ Ψ ≥ (aσ 2 0 -a 2 σ 2 1 -ε + M b) y 2 -M y - adσ 2 1 2 + b 2 σ 2 1 . (3.14)
Let us choose a such that aσ 2 0 -a 2 σ 2 1 > 0. We say that a satisfies Condition C 1 . Then, we choose b such that ε satisfies aσ 2 0 -a 2 σ 2 1 -ε > 0. We say that b satisfies Condition C 2 . We then get: for all L > 0, for all x ∈ R d and for all y ≤ L ,

- LΨ Ψ ≥ -M y - adσ 2 1 2 + b 2 σ 2 1 ≥ -C(L ) with C(L ) = M L+ adσ 2 1 2 + b 2 σ 2 1 > 0 .
October 12, 2016

Statistics: A Journal of Theoretical and Applied Statistics 11˙10˙16Fluctuationdisipation

Step 2: x ∈ R d , y > L 0 . Let ε > 0. As a satisfies C 1 , we deduce from (3.14) that there exists L 0 > 0 such that for all x ∈ R d and for all y > L 0 :

- LΨ Ψ ≥ ε . (3.15)
Step 3: x > L 0 , y ≤ L 0 . From (3.8), (3.10), (3.12) and (3.13), we get, for all x, y ∈ R d :

- LΨ Ψ ≥ (aσ 2 0 -a 2 σ 2 1 ) y 2 + b < ∇V (x), G(x) > -M y - adσ 2 1 2 + b 2 σ 2 1 . (3.16)
Then, from the definition of G and from Assumption (3.6) we get:

< ∇V (x), G(x) >→ ∞ as x → ∞ . (3.17)
Thus there exists L 0 > 0 such that for all y ≤ L 0 and for all x > L 0 we get (3.15). The compact subset is thus defined as K = {z = (x, y) , x ≤ L 0 and y ≤ L 0 }.

Conclusion:

Under Assumptions H 0 , H 1 and (3.6), it is possible to construct explicitly a Lyapunov function satisfying the minoration (3.7): we take Ψ(x, y) = e F (x,y)-inf

R d F ≥ 1 , with F (x, y) = a( 1 2 ||y|| 2 + V (x)) + b < G(x)
, y > and with a and b satisfying Conditions C 1 and C 2 .

Particle in a periodic potential

We now consider z = (x, y) ∈ T d × R d , where T d denotes the d-dimensional torus. Under Assumptions H 0 , H 1 and if V ≥ 1, it is possible to construct explicitly a Lyapunov function (see [START_REF] Mattingly | Geometric ergodicity of some hypo-elliptic diffusions for particle motions[END_REF]Section 3]) satisfying the lower bound (3.7): we take Ψ(x, y) = y 2 2 + V (x).

Main inference results

The aim of this section is to consider nonparametric inference issues for System (1.2). In the following, we only consider the framework with complete observations, that is the framework where both coordinates X and Y are observed. However, most of our results can be extended to the case where only the first coordinate X is observed.

In [START_REF] Cattiaux | Estimation for Stochastic Damping Hamiltonian Systems under Partial Observation. I. Invariant density[END_REF], a central limit theorem has been proved for the estimation of the stationary density p β s (x, y) described by (2.4) in the case where the diffusion term σ is constant. However, as the main ingredient for the proof is the exponential ergodicity, which is here ensured by Assumptions H 0 , H 1 and H 2 , Corollary 4.1 below is straightforward.

Let K be some C 2 function with compact support A such that A K(x, y)dxdy = 1. We may also assume, without loss of generality that A is a bounded ball. Moreover, we assume that there exists m ∈ N * such that for all non constant polynomial P (x, y) with degree less than or equal to m, P (u, v)K(u, v)dudv = 0. We then define

pβ s (x, y) := 1 nb d 1,n b d 2,n ln+n i=ln+1 K x -X ihn b 1,n , y -Y ihn b 2,n . ( 4
(i) b 1,n , b 2,n and h n → 0, (ii) n b d 1,n b d 2,n → +∞, (iii) b1,n b2,n h 2 n → 0 , (iv) ∃ m ∈ N * such that n b d 1,n b d 2,n max(b 1,n , b 2,n ) 2(m+1) → 0, (v) l n h n → +∞.
Then, starting from any initial point z 0 = (x 0 , y 0 ),

nb d 1,n b d 2,n pβ s (x, y) -p β s (x, y) D -----→ n→+∞ N 0, p β s (x, y) K 2 (s, t)dsdt I .
Proof. The exponential ergodicity is a consequence of Hypotheses H 0 , H 1 and H 2 . Theorem 2.10 and Corollary 2.12 in [START_REF] Cattiaux | Estimation for Stochastic Damping Hamiltonian Systems under Partial Observation. I. Invariant density[END_REF] still hold true if σ depends on x (see Theorem 2.1 in [START_REF] Konakov | Explicit parametrix and local limit theorems for some degenerate diffusion processes[END_REF]). Thus the result is a straightforward corollary of Theorem 3.3 and Theorem 6.2 in [START_REF] Cattiaux | Estimation for Stochastic Damping Hamiltonian Systems under Partial Observation. I. Invariant density[END_REF].

Remark 4.2 Further in the paper, we only need the convergence in probability of pβ s (x, y) to apply Slutsky's Theorem. From covariance inequalities stated in [START_REF] Viennet | Inequalities for absolutely regular sequences: application to density estimation[END_REF] for β-mixing sequences, we know that

nh n b d 1,n b d 2,n pβ s (x, y) -p β s (x, y) = O P (1) (4.19)
with O P (1) means stochastically bounded, and as soon as (i) -(ii), (iv) -(v) in Corollary 4.1 are satisfied.

Let us now focus our attention on the diffusion term σ(x). We know from Theorem 4.2 in [START_REF] Cattiaux | Estimation for Stochastic Damping Hamiltonian Systems under Partial Observation. III. Diffusion term[END_REF] that under Hypotheses H 0 , H 1 and H 2 , it is possible to prove a central limit theorem for the estimation of E µ σ 2 (X 0 , Y 0 ) where µ(dz) = p β s (z)dz is the unique stationary density for System (1.2) described by (2.4). Indeed, let g(x, y) = -(σ 2 (x)y + ∇V (x)). We know from the expression of p β s (z) that E µ (|g(Z 0 )| r ) < ∞ for all r > 0. It yields the first part of Assumption H 4 in [START_REF] Cattiaux | Estimation for Stochastic Damping Hamiltonian Systems under Partial Observation. III. Diffusion term[END_REF]. The second part is a straightforward consequence of our Hypothesis H 2 .

The aim is now to propose a methodology to get a pointwise estimate of the diffusion term. For that purpose, we first consider the drift estimation. Recall that g(x, y) = -[σ 2 (x)y + ∇V (x)]. We define the Nadaraya-Watson estimator :

H n (x, y) = 1 (n -1)b d 1,n b d 2,n ln+n-1 i=ln+1 K( x -X ihn b 1,n , y -Y ih b 2,n ) Y (i+1)hn -Y ihn h n . ( 4
(ii') n h n b d 1,n b d 2,n → +∞, (iii) ∃ p > 1 such that h n (b 1,n b 2,n ) d(2-1 p ) -----→ n→+∞ 0, (iv) ∃ m ∈ N * such that n b d 1,n b d 2,n max(b 1,n , b 2,n ) 2(m+1) → 0, (v) l n h n → +∞
and:

• H 3 g belongs to the domain of the infinitesimal generator L, in all L γ (µ), for 1 ≤ γ < +∞.

Then, starting from any initial condition z 0 = (x 0 , y 0 ), we get : From this formula, we can prove that Hypothesis H 3 is satisfied as soon as the potential V has a polynomial growth as x tends to infinity.

E H n pβ s (x, y) -g(x, y) 2 = O 1 nb d 1,n b d 2,n h n . ( 4
Proof. In the following, we consider the stationary case, that is we start from the stationary distribution: (X 0 , Y 0 ) ∼ p β s (z)dz with p β s (z) described by (2.4), and l n ≡ 0:

H n (x, y) = 1 (n -1)b d 1,n b d 2,n n-1 i=1 K( x -X ihn b 1,n , y -Y ihn b 2,n ) Y (i+1)hn -Y ihn h n .
The case where we start from any initial point z 0 = (x 0 , y 0 ) can then be deduced from the stationary case, following the proof in [2, Theorem 5.1].

Step 1: Let us first prove that E[H n (x, y)] → g(x, y)p β s (x, y) as n tends to infinity. (4.22)

We define the auxiliary sequence (n -1)b d 1,n b d 2,n Ĥn (x, y) as:

n-1 i=1 K( x -X ihn b 1,n , y -Y ihn b 2,n ) σ(X ihn )(W (i+1)hn -W ihn ) + (i+1)hn ihn g(X t , Y t )dt h n . Then (n -1)h n b d 1,n b d 2,n H n (x, y) -Ĥn (x, y) is equal to: n-1 i=1 K( x -X ihn b 1,n , y -Y ihn b 2,n ) (i+1)hn ihn [σ(X t ) -σ(X ihn )]dW t . (4.23)
Let • denote the euclidean norm in R d . Then, applying the Itô isometry to (4.23), we bound

h 2 n (n -1) 2 b 2d 1,n b 2d 2,n E ||H n (x, y) -Ĥn (x, y)|| 2 by: n-1 i=1 E K 2 ( x -X ihn b 1,n , y -Y ihn b 2,n ) (i+1)hn ihn Tra((σ(X t ) -σ(X ihn )) 2 )dt . (4.24)
Thus, as the first derivative of σ is uniformly bounded, we deduce from (4.24) and from a Taylor formula at order one that:

E ||H n (x, y) -Ĥn (x, y)|| 2 = O h n (n -1)b d 1,n b d 2,n . (4.25)
We then define (n -

1)b d 1,n b d 2,n Hn (x, y) as n-1 i=1 K( x -X ihn b 1,n , y -Y ihn b 2,n ) σ(X ihn )(W (i+1)hn -W ihn ) h n + g(X ihn , Y ihn ) . Then (n -1)h n b d 1,n b d 2,n Ĥn (x, y) -Hn (x, y) is equal to: n-1 i=1 K( x -X ihn b 1,n , y -Y ihn b 2,n ) (i+1)hn ihn (g(X t , Y t )dt -g(X ihn , Y ihn )) dt .
Let P t denote the semigroup associated to the Markov process Z t . Using the super index for denoting the coordinates of a vector, we get:

E Ĥk n (x, y) -Hk n (x, y) = 1 b d 1,n b d 2,n K( x -u b 1,n , y -v b 2,n ) 1 h n hn 0 P t (g k )(u, v) -g k (u, v) dt p β s (u, v)dudv = K(u 1 , v 1 )du 1 dv 1 × hn 0 t h n (P t -I)(g k ) t (x -b 1,n u 1 , y -b 2,n v 1 ) dt p β s (x -b 1,n u 1 , y -b 2,n v 1 ) = o (h n ) . ( 4
.26) Indeed, from the regularity assumptions H 0 and H 1 , each coordinate of g is in the domain of L, and then the dominated convergence Theorem allows to take the limit n → +∞ in the integral. Moreover, we have: Step 2: This step is devoted to the study of E Ĥn (x, y) -Hn (x, y)

E[ Hn (x, y)] = K(u 1 , v 1 )g(x -b 1,n u 1 , y -b 2,n v 1 )p β s (x -b 1,n u 1 , y -b 2,n v 1 )

2

. We define:

B n (i) = 1 h n (i+1)hn ihn (g(X t , Y t )dt -g(X ihn , Y ihn )) dt .
We have:

Ĥn (x, y) -Hn (x, y) = 1 (n -1)b d 1,n b d 2,n n-1 i=1 K( x -X ihn b 1,n , y -Y ihn b 2,n )B n (i).
The above equality can be written coordinatewise as:

Ĥk n (x, y) -Hk n (x, y) = 1 (n -1)b d 1,n b d 2,n n-1 i=1 K( x -X ihn b 1,n , y -Y ihn b 2,n )B k n (i) .
We have:

Var (n -1)b d 1,n b d 2,n ( Ĥk n (x, y) -Hk n (x, y)) = n-1 i=1 Var(K( x -X ihn b 1,n , y -Y ihn b 2,n )B k n (i)) + 1≤i =l≤n-1 Cov(K( x -X ihn b 1,n , y -Y ihn b 2,n )B k n (i), K( x -X lhn b 1,n , y -Y lhn b 2,n )B k n (l)) .
To bound the first term we need to bound the second moment and the square of the first moment. As the procedure is detailed in [START_REF] Cattiaux | Estimation for Stochastic Damping Hamiltonian Systems under Partial Observation. II. Drift term[END_REF] we only detail the computations of the bound for the second moment. By stationarity we have:

n-1 i=1 E (K( x -X ihn b 1,n , y -Y ihn b 2,n )B k n (i)) 2 = (n -1) E[(K( x -X 0 b 1,n , y -Y 0 b 2,n )B k n (0))] 2 ≤ (n -1) 1 h 2 n hn 0 t 2 K 2 ( x -u b 1,n , y -v b 2,n ) P t (g k )(u, v) -g k (u, v) t 2 p β s (u, v)dudvdt = O (n -1)h n b d/p 1,n b d/p 2,n . (4.28)
This latter bound has been obtained by applying Hölder's inequality with conjugate exponents p and q. Then we used that

Pt(g k )(u,v)-g k (u,v) t is bounded in L 2q (p β s (u, v)dudv) by assumption.
To bound the covariance terms, we use the covariance inequality in [5, Section 1.2.2, p. 9, Theorem 3, Item (1)], which can be applied here as the process (Z t ) is β-mixing (see (2.5)). In the following, C denotes some positive constant which may vary from line to line. Let r, p ∈ N * . We get: October 12, 2016 Statistics: A Journal of Theoretical and Applied Statistics 11˙10˙16Fluctuationdisipation

i =l Cov(K( x -X ihn b 1,n , y -Y ihn b 2,n )B k n (i), K( x -X lhn b 1,n , y -Y lhn b 2,n )B k n (l)) ≤ Cn n-2 j=1 min b d 1,n b d 2,n 1-1 r e -jhn λ r , Var K x -X 0 b 1,n , y -Y 0 b 2,n B k n (0) ≤ Cn n-2 j=1 min b d 1,n b d 2,n 1-1 r e -jhn λ r , b d 1,n b d 2,n 1 p h n .
If we choose r, p such that 1 r + 1 p = 1, then the right hand term in the above inequality

is O n (b 1,n b 2,n ) d p .
Summing the variance and the covariance terms we obtain:

Var( Ĥk n (x, y) -Hk n (x, y)) = O n -1 (b 1,n b 2,n ) d p -2d . (4.29)
Using Cauchy-Schwarz Inequality and (4.28), it is possible to prove that

E Ĥk n (x, y) -E Hk n (x, y) 2 = O h n b d p -2d 1,n b d p -2d 2,n . (4.30)
Thus, we deduce from (4.29) and (4.30),

E Ĥk n (x, y) -Hk n (x, y) 2 = O h n (b 1,n b 2,n ) d p -2d . (4.31)
Step 3: In this last step, we want to prove that E(H n (x, y) -g(x, y)p β s (x, y)) 2 tends to zero as n tends to infinity. We write:

E(H k n (x, y) -g k (x, y)p β s (x, y)) 2 ≤ 3 J k 1n + J k 2n + J k 3n (4.32) with J k 1n = E H k n (x, y) -Ĥk n (x, y) 2 , J k 2n = E Ĥk n (x, y) -Hk n (x, y) 2 and J k 3n = E Hk n (x, y) -g k (x, y)p β s (x, y) 2 .
From (4.25) and (4.31) we get:

J k 1n + J k 2n = O h n (n -1)b d 1,n b d 2,n + h n (b 1,n b 2,n ) d p -2d . ( 4

.33)

Let us now consider the third term, J k 3n , in the decomposition (4.32). The stationarity yields:

E Hk n (x, y) = K(u 1 , v 1 )g k (x -b 1,n u 1 , y -b 2,n v 1 )p β s (x -b 1,n u 1 , y -b 2,n v 1 )du 1 dv 1 .
The asymptotic behavior of this term can be classically handled with Hypothesis (iv) (see e.g., [2, page 14]). It yields:

(n -1)b d 1,n b d 2,n h n J k 3n = (n -1)b d 1,n b d 2,n h n (E[ Hk n (x, y) -g k (x, y)p β s (x, y)]) 2 = o(1) . (4.34)
From (4.33) and (4.34) we get for each k:

E(H k n (x, y) -g k (x, y)p β s (x, y)) 2 = O 1 (n -1)b d 1,n b d 2,n h n + h n (b 1,n b 2,n ) d p -2d
. (4.35)

As (i) -(ii ), (iv) -(v)
E H n pβ s (x, y) -g(x, y) 2 = O 1 nb d 1,n b d 2,n h n + h n (b 1,n b 2,n ) d p -2d
.

Then, choosing p > 1 such that (iii) in Theorem 4.3 is satisfied, we get (4.21) which achieves the proof of Theorem 4.3.

Corollary 4.5 Let (e i ) i=1,...,d denote the canonical basis of R d . Under Assumptions of Theorem 4.3 we have :

-< ĝn (x, e i ) -ĝn (x, 0), e j > Pz → σ 2 ij (x).
Proof. It is a straightforward consequence of Theorem 4.

Let us now consider an estimation procedure for -β ∇V (x). We know from (2.4) that p Then, starting from any initial condition z 0 = (x 0 , y 0 ), we get :

nb (d+2) 1,n b d 2,n ∇ x pβ s pβ s (x, y) + β∇V (x) D -----→ n→+∞ N (0, L(x, y)) ,
where L = (l i,j (x, y)) 1≤i,j≤d is a diagonal square matrix with

l jj (x, y) = 1 p β s (x, y) (∂ xj K(u, v)) 2 dudv ,
where p β s (x, y) is described by (2.4).

Proof. In the following, we consider once more the stationary case, that is we start from the stationary distribution: (X 0 , Y 0 ) ∼ p β s (z)dz with p β s described by (2.4) and l n ≡ 0. The extension to the non-stationary case is classical and is left to the reader. We have :

A n (x, y) := ( ∇ x pβ s pβ s (x, y) + β∇V (x)) = ( ∇ x pβ s pβ s (x, y) - ∇ x p β s pβ s (x, y)) + ( ∇ x p β s pβ s (x, y) + β∇V (x)) = 1 pβ s (x, y) (∇ x pβ s (x, y) -∇ x p β s (x, y)) - ∇ x p β s pβ s p β s (x, y)(p β s (x, y) -p β s (x, y)).
Then, as in the proof of Theorem 3.3 in [START_REF] Cattiaux | Estimation for Stochastic Damping Hamiltonian Systems under Partial Observation. I. Invariant density[END_REF], we can prove that the second term in the right hand side of the last equality above is

O Pz ( nb d 1,n b d 1,n -1/2
). Recaling that D denotes the convergence in distributions of probability measures, we have

D lim n→∞ nb (d+2) 1,n b d 2,n A n (x, y) = D lim n→∞ nb (d+2) 1,n b d 2,n 1 pβ s (x, y) (∇ x pβ s (x, y) -∇ x p β s (x, y)) = 1 p β s (x, y) D lim n→∞ nb (d+2) 1,n b d 2,n (∇ x pβ s (x, y) -∇ x p β s (x, y)) ,
where the last line was obtained from Slutsky's Theorem.

In the following we sketch the study of the convergence in distribution of

R n := nb (d+2)
1,n b d 2,n (∇ x ps (x, y) -∇ x p β s (x, y)). Let ∂x l denote the partial derivative with respect to the x l -variable. We get from the covariance inequality in [5, Section 1.2.2, p. 9, Theorem 3, Item (1)] :

Cov(∂ xj p(x, y), ∂ xl p(x, y)) = O( p β s (x, y) nb (d+2) 1,n b d 2,n ∂ xj K(u, v)∂ xl K(u, v)dudv).
Hence and as a consequence that the kernel K has bounded support we have (nb

(d+2) 1,n b d 2,n )Cov(∂ xj p(x, y), ∂ xl p(x, y)) → δ jl p β s (x, y) (∂ xj K(u, v)) 2 dudv,
where δ ij stand for the Kronecker symbol. The random sequence R n is a sum of a triangular arrays of β-mixing random vectors in R d . An easy modification via the Cramer-Wald device allows to extend [START_REF] Cattiaux | Estimation for Stochastic Damping Hamiltonian Systems under Partial Observation. I. Invariant density[END_REF]Theorem 4.3] to random vectors. Then, defining D(x, y) = (d ij (x, y)) as a diagonal matrix and if the sequences h n , b 1,n and b 2,n satisfy Hypotheses (i), (ii) , (iii) and (iv), we get

R n D → N (0, D(x, y)) where d jj (x, y) = p β s (x, y) (∂ xj K(u, v)) 2 dudv.
It thus yields :

nb (d+2) 1,n b d 2,n A n (x, y) D → N (0, L(x, y))
where L is also diagonal and l jj (x, y) = 1 p β s (x, y)

(∂ xj K(u, v)) 2 dudv.

Numerical experiments

This section is devoted to numerical experiments on inference issues for System (1.2).

For the simulations, we chose l n ≡ 0 in Formula (4.18).

A constant diffusion term

We first consider the particular case where V (x) = 6x 6 -12x 4 + 6x 2 + 1 7 , s(x) ≡ 2 and β = 2. This particular case enters the framework of Section 3, Item 1. We thus consider the following system:

dX t = Y t dt dY t = 2dW t -(4Y t + 36X 5 t -48X 3 t + 12X t )dt.
(5.36)

We will first implement the nonparametric estimation of the invariant density

p s (x, y) = C exp -y 2 -12x 6 + 24x 4 -12x 2 - 2 7
and of the drift term g(x, y) = -4y -36x 5 + 48x 3 -12x .

In the following, we make use of the explicit Euler-Verlet-Midpoint scheme (see e.g., [8, (2.50) p. 94]) in which the deterministic Verlet step has been corrected using a generalized Metropolis-Hastings strategy (see [START_REF] Lelièvre | Free energy computations: a mathematical perspective[END_REF]Algorithm 2.11,p. 95]) to simulate an approximated discrete sampling ( X i , Y i ) i∈N of (X t , Y t ) t∈R+ .

Estimation of the invariant density

We estimate the invariant density p β s using the kernel estimator defined by (4.18) in Section 4, and by replacing (X i , Y i ) i∈N by ( X i , Y i ) i∈N . The invariant density is estimated on a grid (z l ) l=1,...,L = (x l , y l ) l=1,...,L of size L = 55 × 39 of

min{ X i -b 1,n } n i=1 , max{ X i + b 1,n } n i=1 × min{ Y i -b 2,n } n i=1 , max{ Y i + b 2,n } n i=1 .
We use the Epanechnikov kernel defined by K(u, v) = K e (u) × K e (v) with K e (w) = each of size n,

1 M M m=1 1 L L l=1 p(m) s (z l ) -p s (z l ) 2 1 L L l=1 (p s (z l )) 2
, and we obtain: We use a regularised version of the estimator considered in Theorem 4.3: gn (x, y) = H n (x, y) ps (x, y) + ε (5.37)

Estimation of the drift term

where ps has been defined in (4.18), and H n (x, y) in (4.20), and where ε > 0 is a small term whose aim is the regularisation of the estimate. We then replace, as previously,

(X i , Y i ) i∈N by ( X i , Y i ) i∈N .
The estimation of the drift in Figure 3 We also provide in Figure 4 a plot of the absolute difference between the drift and its estimate. Figure 4 emphasizes the limit of the approach, as we can see some boundary effects in the estimation. thus have: ∇V (x) = 36x 5 -48x 3 + 12x. We now consider the estimation of -β∇V (x) = -2∇V (x) as in Theorem 4.6.

On Figure 5, we present in the left hand side a discretization of ∇ x p s (x, y), and in the right hand side its estimate ∇ x ps (x, y). On Figure 6, we present in the left hand side a discretization of -2V (x) for y 0 = -0.265, and in the right hand side its estimate ∇x ps(x,y) ps(x,y) . We choose here a gaussian kernel. The plot of the absolute difference for the estimation of the derivative of the potential in Figure 7 also shows some boundary effects. The correction of boundary effects could be the next aim in this framework.

Finally, we have drawn in Figure 8 below gn (x, y) -gn (x, 0) for two different values of y, to illustrate the result of Corollary 4.5.

We conclude that the quality of gn (x, y) -gn (x, 0) to estimate β 2 depends on the value of y. It would be interesting to investigate further the asymptotic properties of this estimator, in the more general case where the diffusion term σ(x, y) depends on the position x.

A diffusion term that depends on the position

In this section, we consider the following system:

   dX t = Y t dt dY t = exp - 1 X 2 t + 1 dW t -exp - 2 X 2 t + 1 Y t + 36X 5 t -48X 3 t + 12X t dt.
(5.38) We will first implement the nonparametric estimation of the invariant density

p s (x, y) = C exp -y 2 -12x 6 + 24x 4 -12x 2 - 2 7 ,
which is unchanged with respect to Section 5.1, and of the drift term g(x, y) = -exp -1 x 2 + 1 y -36x 5 + 48x 3 -12x .

The choice of the kernel and of the schemes for simulation are similar to the ones in Section 5.1 and are thus not recalled here.

Estimation of the invariant density For a sample size of n = 10 6 we obtain the following graphics: the discretized theoretical bivariate invariant density and its estimated version are drawn in Figure 9. In Figure 10 is plotted the absolute difference between the invariant density and its estimate.

  in the sense of quadratic forms. Remark that the right hand side of Inequality (3.7) is composed with two terms, the first one being concerned with what happens on the non bounded domain K c and the second one with what happens on the compact set K.

  October 12, 2016 Statistics: A Journal of Theoretical and Applied Statistics 11˙10˙16Fluctuationdisipation

. 21 )

 21 Remark 4.4 From the fluctuation-dissipation relation, we have an analytical formula for the stationary density p β s .

Theorem 4 . 6

 46 β s (x, y) = C(β)e -β(|y| 2 +V (x)) , thus ∇xp β s p β s (x, y) = -β∇V (x). Theorem 4.6 below stated a central limit theorem for ∇x pβ Assume Hypotheses H 0 , H 1 and H 2 are fulfilled. Assume that the bandwidths b 1,n , b 2,n and the discretization step h n satisfy Hypotheses (i), (ii), (iii), (iv) and (v) of Corollary 4.1. Assume that the kernel K is twice continuously differentiable on its bounded support.

  Journal of Theoretical and Applied Statistics 11˙10˙16Fluctuationdisipation

3 4 ( 1 -

 31 w 2 )1 |w|≤1 , and we consider h n = n -γ , b 1,n = n -α1 and b 2,n = n -α2 , with α 1 = α 2 = 0.2 and γ = 0.15. The relative mean integrated squared error below is computed on M = 30 samples, October 12, 2016 Statistics: A Journal of Theoretical and Applied Statistics 11˙10˙16Fluctuationdisipation

6

 6 Relative mean integrated squared error 0.2225 0,03841 0,0071 0,0012 As expected, the relative error decreases as n increases. For a sample size of n = 10 6 we obtain the following graphics: the discretized theoretical bivariate invariant density and its estimated version are drawn in Figure1. The estimation of the marginal invariant density for the position (resp. velocity) is drawn in the left (resp. right) hand part of Figure2.

Figure 1 .

 1 Figure 1. Discretized theoretical (left) and estimated (right) invariant bivariate density for V (x) = 6x 6 -12x 4 + 6x 2 + 1/7, σ(x) ≡ 2.

Figure 2 .

 2 Figure 2. Theoretical (plain line) and estimated (dashed line) invariant density of the position (left) and of the velocity (right) for the potential V (x) = 6x 6 -12x 4 + 6x 2 + 1/7.

  Journal of Theoretical and Applied Statistics 11˙10˙16Fluctuationdisipation

  is drawn on a discretized grid of (x, y) ∈ [-1.15, 1.15] × [-1.7, 1.7]. Note that [-1.15, 1.15] × [-1.7, 1.7] on which the estimated invariant density of the position and the one of the velocity are both positive.

Figure 3 .

 3 Figure 3. Discretized theoretical (left) and estimated (right) drift on [-1.15, 1.15] × [-1.7, 1.7].

Figure 4 .

 4 Figure 4. Absolute value of the difference between the drift and its estimate at left on the restricted domain [-0.85, 0.85] × [-1.7, 1.7], and a same graphic with a different orientation at right

Figure 5 .

 5 Figure 5. Discretization of ∇xps(x, y) (left), and its estimate (right).

Figure 6 .

 6 Figure 6. Discretization of the -2V (x) (plain line) and its estimate (dashed line) for the fixed value y = -0.265.

Figure 7 .

 7 Figure 7. Absolute value of the difference between xps and its estimate.

Figure 8 .

 8 Figure 8. The dashed line corresponds to gn(x, y) -gn(x, 0) for y = 1 (left) and for y = 0.085 (right), the red line corresponds to the corresponding value β 2 × y = 4y.

Figure 9 .

 9 Figure 9. Discretized theoretical (left) and estimated (right) invariant bivariate density for V(x) = 6x 6 -12x 4 + 6x 2 + 1/7, σ(x) = exp(-1 x 2 +1).

Figure 10 .

 10 Figure 10. Absolute value of the difference between the theoretical invariant density and its estimate.

Figure 12

 12 Figure 12 plots the absolute difference between the drift and its estimate.

Figure 11 .

 11 Figure 11. Discretized theoretical (left) and estimated (right) drift on [-1, 1] × [-1.7, 1.7].

Figure 12 .

 12 Figure 12. Absolute value of the difference between the drift and its estimate at left on the restricted domain [-0.85, 0.85] × [-1.7, 1.7], and the same graphic with a different orientation at right.

  .18) Corollary 4.1 Assume Hypotheses H 0 , H 1 and H 2 are fulfilled. Assume that the bandwidths b 1,n , b 2,n and the discretization step h n satisfy the following assumptions :

  Assume that Hypotheses H 0 , H 1 and H 2 are satisfied. Moreover we assume that (i) b 1,n , b 2,n and h n → 0,

	.20)
	We get Theorem 4.3 stated below.
	Theorem 4.3

  in Theorem 4.3 imply (i) -(ii), (iv) -(v) in Corollary 4.1, we get the bound stated in (4.19). Thus, from Slutsky's Theorem and from (4.35) we get :

Acknowledgments. This work has been partially supported by the project ECOS NORD/FONACIT (France-Venezuela) V12M01 Equations différentielles stochastiques et analyse de sensibilité : applications environnementales, by the 16-MATH-03 SIDRE project of the MATH-AmSud Program and by the Inria International Chairs program.