domain: stat math

text:
F Comte

NONPARAMETRIC MOMENT METHOD FOR MCKEAN-VLASOV STOCHASTIC DIFFERENTIAL EQUATIONS

Keywords: Invariant density, interaction function, Hermite basis, McKean-Vlasov stochastic dierential equation, moment method, projection estimators, nonparametric drift estimation

We study the nonparametric estimation of both the potential and the interaction terms of a McKean-Vlasov stochastic dierential equation (SDE) in stationary regime from a continuous observation on a time interval [0, T ], with asymptotic framework T → +∞. The problem is quite dierent from the case of usual diusions with no interaction term and the observation of only one sample path is not enough to estimate both functions. We consider the observation of four i.i.d. sample paths. The observation of two sample paths could be enough at the cost of much more computations. Estimators of the potential and the interaction functions are built using a combination of a moment method and a projection method on sieves. We dene a specic risk tted to this estimation problem and obtain a bound for it. A nonparametric estimator of the invariant density also is proposed. The method is implemented on simulated data for several examples of McKean-Vlasov SDEs and a model selection procedure is experimented.

Introduction

McKean-Vlasov processes are described by a stochastic dierential equation where the coecients depend on the distribution of the solution itself. They are also called nonlinear or self-stabilizing processes. These processes have been extensively studied since their rst description by [START_REF] Mckean | A class of Markov processes associated with nonlinear parabolic equation[END_REF], [START_REF] Mckean | A class of Markov processes associated with nonlinear parabolic equation[END_REF], and appear in various applications: for the modelling of granular media in statistical physics [START_REF] Benedetto | A kinetic equation for granular media[END_REF], [START_REF] Benedetto | A kinetic equation for granular media[END_REF]), in neurosciences (see e.g. [START_REF] Baladron | Mean-eld description and propagation of chaos in networks of Hodgkin-huxley and FitzHugh-Nagumo neurons[END_REF], [START_REF] Baladron | Mean-eld description and propagation of chaos in networks of Hodgkin-huxley and FitzHugh-Nagumo neurons[END_REF], [START_REF] Dawson | Critical Dynamics and Fluctuations for a Mean-Field Model of Cooperative Behavior[END_REF], [START_REF] Dawson | Critical Dynamics and Fluctuations for a Mean-Field Model of Cooperative Behavior[END_REF]), for population dynamics and ecology [START_REF] Carrillo | The derivation of swarming models: Mean-eld limit and Wasserstein distances[END_REF], [START_REF] Carrillo | The derivation of swarming models: Mean-eld limit and Wasserstein distances[END_REF], Molginer and Edelstein-Keshet (1999), [START_REF] Molginer | A non-local model for a swarm[END_REF]), for epidemics dynamics [START_REF] Ball | Stochastic SIR in Structured Populations[END_REF], [START_REF] Ball | Stochastic SIR in Structured Populations[END_REF], [START_REF] Forien | Household epidemic models and McKean-Vlasov Poisson driven SDEs[END_REF], [START_REF] Forien | Household epidemic models and McKean-Vlasov Poisson driven SDEs[END_REF]) and in nance (see e.g. [START_REF] Giesecke | Inference for large nancial systems[END_REF], [START_REF] Giesecke | Inference for large nancial systems[END_REF] and the references therein). A large number of contributions is devoted to probabilistic properties of these models. Existence and uniqueness of solutions under several dierent sets of assumptions can be found in e.g., Snitzman (1991), [START_REF] Sznitman | Topics in propagation of chaos[END_REF], [START_REF] Méléard | Asymptotic behaviour of some interacting particle systems; McKean-Vlasov and Boltzmann models[END_REF], [START_REF] Méléard | Asymptotic behaviour of some interacting particle systems; McKean-Vlasov and Boltzmann models[END_REF], [START_REF] Funaki | A certain class of diusion processes associated with nonlinear parabolic equations[END_REF], [START_REF] Funaki | A certain class of diusion processes associated with nonlinear parabolic equations[END_REF], [START_REF] Gärtner | On the Mc-Kean-Vlasov limit for interacting diusions[END_REF], [START_REF] Gärtner | On the Mc-Kean-Vlasov limit for interacting diusions[END_REF], Benachour (1) : Université de Paris, MAP5, UMR 8145 CNRS,F-75006, FRANCE, email: fabienne.comte@u-paris.fr, valentine.genon-catalot@mi.parisdescartes.fr (2) Université Paris Cité, LPSM UMR 8001 CNRS,F-75006, FRANCE, catherine.laredo@inrae.fr. et al. (1998a), [START_REF] Benachour | Nonlinear self-stabilizing processes -I Existence, invariant probability, propagation of chaos[END_REF], [START_REF] Benachour | Nonlinear self-stabilizing processes -II Convergence to invariant probability[END_REF], [START_REF] Benachour | Nonlinear self-stabilizing processes -II Convergence to invariant probability[END_REF], Hermann et al. (2008), [START_REF] Herrmann | Large deviations and a Kramers'type low for selfstabilizing diusions[END_REF], [START_REF] Kolokoltsov | Non linear Markov processes and kinetic equations 182[END_REF], [START_REF] Kolokoltsov | Non linear Markov processes and kinetic equations 182[END_REF]. More references are given below.

We are interested in a one dimensional McKean-Vlasov process of the form [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF] dX(t) = -[b(X(t)) + φ µ t (X(t))]dt + σdW (t), X(0) = η ∼ ν, where b(.), φ(.) are deterministic functions, denotes the convolution product, µ t (.) is the law of X(t), W is a Brownian motion and η is a random variable independent of W . Such process can be obtained as the limit of a system of N interacting particles as N tends to innity (propagation of chaos). More precisely, consider the N -dimensional process given by the stochastic dierential system:

dX N i (t) = -b(X N i (t))dt - 1 N N j=1 φ(X N i (t) -X N j (t))dt + σdW i (t), (2)
= -b(X N i (t))dt -φ µ N t (X N i (t))dt + σdW i (t), X N i (0) = X i 0 , i = 1, . . . , N,
where

µ N t = N -1 N j=1 δ X N j (t)
is the empirical measure associated with (X N j (t), j = 1, . . . , N ), (W i , i = 1 . . . , N ) are N independent Brownian motions, X i 0 , i = 1, . . . , N are independent and identically distributed (i.i.d.) random variables, independent of (W j , j = 1, . . . , N ). The function φ describes the interaction between the N particles. For any xed integer p, the process (X N i (t), i = 1, . . . , p) converges in distribution as N tends to innity to p i.i.d. McKean-Vlasov processes (X (i) (t)), given by

dX (i) (t) = -[b(X (i) (t)) + φ µ t (X (i) (t))]dt + σdW i (t), X (i) (0) = X i 0 , i = 1, . . . , p.
Thus, the process (X (i) (t)) provides a good approximation of the behaviour of the i-th particle of the system (see e.g., [START_REF] Bolley | Uniform convergence to equilibrium for granular media[END_REF], [START_REF] Bolley | Uniform convergence to equilibrium for granular media[END_REF], Eberle et al. (2019), [START_REF] Eberle | Quantitative Harris type theorems for diusions and Mc-Kean-Vlasov processes[END_REF], [START_REF] Fournier | On the rate of convergence in Wasserstein distance of the empirical measure[END_REF], [START_REF] Fournier | On the rate of convergence in Wasserstein distance of the empirical measure[END_REF], [START_REF] Malrieu | Convergence to equilibrium for granular media equations and their Euler schemes[END_REF], [START_REF] Malrieu | Convergence to equilibrium for granular media equations and their Euler schemes[END_REF], Cattiaux et al. (2008), [START_REF] Cattiaux | Probabilistic approach for granular media equations in the non-uniformly convex case[END_REF]). Moreover, existence and uniqueness of invariant distributions for [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF] are studied in Benachour et al. (1998a,b), [START_REF] Benachour | Nonlinear self-stabilizing processes -I Existence, invariant probability, propagation of chaos[END_REF], [START_REF] Benachour | Nonlinear self-stabilizing processes -II Convergence to invariant probability[END_REF], [START_REF] Malrieu | Convergence to equilibrium for granular media equations and their Euler schemes[END_REF], [START_REF] Malrieu | Convergence to equilibrium for granular media equations and their Euler schemes[END_REF], Butkovski (2014), [START_REF] Butkovsky | On ergodic properties of nonlinear Markov chains and Mc-Kean-Vlasov equations[END_REF], Cattiaux et al. (2008), [START_REF] Cattiaux | Probabilistic approach for granular media equations in the non-uniformly convex case[END_REF], Eberle et al. (2019), [START_REF] Eberle | Quantitative Harris type theorems for diusions and Mc-Kean-Vlasov processes[END_REF], Hermann and Tugaut (2010), [START_REF] Herrmann | Non uniqueness of stationary measures for self-stabilizing diusions[END_REF]. Statistical inference for interacting particle systems and their mean-eld limits has started more recently. On the one hand, several papers deal with inference based on the observation on a xed time interval of an interacting particle system with asymptotic properties as the number N of particles tends to innity (see e.g. parametric inference in Kasonga, 1990, [40], Giesecke et al. (2020), [START_REF] Giesecke | Inference for large nancial systems[END_REF], Chen (2021), [START_REF] Chen | Maximum likelihood estimation of potential energy in interacting particle systems from single-trajectory data[END_REF], Li et al. (2021), [START_REF] Li | On the identiability of interaction functions in systems of interacting particles[END_REF], Sharrock et al. (2023), [START_REF] Sharrock | Online parameter estimation for the McKean-Vlasov stochastic dierential equation Stochastic Process[END_REF], Lu et al. (2022), [START_REF] Lu | Learning Interaction Kernels in Stochastic Systems of Interacting Particles from Multiple Trajectories[END_REF], Amorino et al. ( 2023), [START_REF] Amorino | Parameter estimation of discretely observed interacting particle systems[END_REF], Della Maestra and Homann (2022b), [START_REF] Della Maestra | The LAN property for McKean-Vlasov models in a mean-eld regime[END_REF], Pavlotis and Zanoni (2022), [START_REF] Pavliotis | Eigenfunction martingale estimators for interacting particle systems and their mean eld limit[END_REF], semiparametric inference in Belomestny et al. (2023), [START_REF] Belomestny | Semiparametric estimation of McKean-Vlasov SDEs[END_REF], and nonparametric inference in Della Maestra and Homann (2022a), [START_REF] Della Maestra | Nonparametric estimation for interacting particle systems: McKean-Vlasov models[END_REF], Comte and Genon-Catalot (2023), [START_REF] Comte | Nonparametric adaptive estimation for interacting particle systems[END_REF]).

On the other hand, several papers consider inference based on the observation of the mean-eld limit process (1) on a time interval [0, T ], that is inference for the limiting process of one typical particle for large N . In Genon-Catalot and Larédo (2021a), [START_REF] Genon-Catalot | Probabilistic properties and parametric inference of small variance nonlinear self-stabilizing stochastic dierential equations[END_REF] and (2021b), [START_REF] Genon-Catalot | Parametric inference for small variance and long time horizon McKean-Vlasov diusion models[END_REF], estimation of unknown parameters in the potential term b and the interaction term φ are studied under the asymptotic framework σ → 0. More recently, statistical inference has been investigated under the assumption of stationarity for (1) when the process is observed in the stationary regime with asymptotic framework T → +∞. In Genon-Catalot and Larédo (2023a), [START_REF] Genon-Catalot | Inference for ergodic McKean-Vlasov stochastic dierential equations with polynomial interactions[END_REF], the case of b = 0 and φ equal to an odd polynomial is studied and the estimation of the coecients of the polynomial is treated by a pseudo-likelihood approach. In Genon-Catalot and Larédo (2023b), [START_REF] Genon-Catalot | Parametric inference for ergodic McKean-Vlasov stochastic differential equations[END_REF], general functions b and φ depending on unknown parameters are considered and a dierent pseudo-likelihood function is proposed. Pavlotis and Zanoni (2023), [START_REF] Pavliotis | A method of moments estimator for interacting particle systems and their mean eld[END_REF], use a moment method to estimate the coecients of polynomial functions b and φ.

In here, our aim is the nonparametric estimation of the functions b(.), φ(.) from the continuous observation of a process distributed as (1) on a time interval [0, T ] when the process is in stationary regime and T tends to innity. When φ ≡ 0 (usual diusion process), nonparametric inference for the function b, based on the observation of one sample path on a time interval [0, T ], when the process is in stationary regime and T is large, is an extensively developped subject. We can quote Homann (1999), [START_REF] Homann | Adaptive estimation in diusion processes[END_REF], Comte et al. (2007), [START_REF] Comte | Penalized nonparametric mean square estimation of the coecients of diusion processes[END_REF], see also the book Kessler et al. (2012), [START_REF] Kessler | Statistical methods for stochastic dierential equations[END_REF], Dalalyan andReiss (2006, 2007), [START_REF] Dalalyan | Asymptotic statistical equivalence for scalar ergodic diusions[END_REF], [START_REF] Dalalyan | Asymptotic statistical equivalence for ergodic diusions: the multidimensional case[END_REF], Schmisser (2013), [START_REF] Schmisser | Penalized nonparametric drift estimation for a multidimensional diusion process[END_REF], [START_REF] Strauch | Sharp adaptive drift estimation for ergodic diusions: the multivariate case[END_REF][START_REF] Strauch | Exact adaptive pointwise drift estimation for multidimensional ergodic diusions[END_REF], [START_REF] Strauch | Sharp adaptive drift estimation for ergodic diusions: the multivariate case[END_REF], [START_REF] Strauch | Exact adaptive pointwise drift estimation for multidimensional ergodic diusions[END_REF], Nickl and Ray (2020), [START_REF] Nickl | Nonparametric statistical inference for drift vector elds of multi-dimensional diusions[END_REF], Comte and Genon-Catalot (2021), [START_REF] Comte | Drift estimation on non compact support for diusion models[END_REF].

In the case of the McKean-Vlasov process (1) with φ = 0 observed in the stationary regime, the joint nonparametric estimation of (b, φ) with general functions (b, φ) has never been investigated up to our knowledge. It diers completely from the case of a usual diusion. Indeed, although there are clear assumptions (see below) for existence and uniqueness of an invariant distribution, this distribution is not explicitely given as a function of b and φ. It is solution of an implicit equation. The convolution term depending on the unknown marginal distribution of the process introduces new diculties for the nonparametric estimation of b and φ linked with identiability problems. For these reasons, the observation of only one sample path is not enough to infer both b and φ. The novelty, in what follows, is that we assume that the observation is composed of four independent copies (X ( ) (t), t ∈ [0, T ], = 1, 2, 3, 4) in stationary regime and our asymptotic framework is T → +∞. Actually, only two copies could be sucient but this would induce cumbersome and tedious additional computations. Our estimation method is a moment method inspired by the work in parametric setting of [START_REF] Pavliotis | A method of moments estimator for interacting particle systems and their mean eld[END_REF], [START_REF] Pavliotis | A method of moments estimator for interacting particle systems and their mean eld[END_REF].

In Section 2, we present our set of assumptions ensuring that equation (1) admits a unique strong solution and that the model has a unique invariant distribution ρ(x)dx. We give useful properties of the density ρ and of the process (1) in the stationary regime. Proposition 3 states a result on the innitesimal generator of (1) often used in the proofs. In Section 3, the principle of the estimation method is described. It is based on properties of the innitesimal generator of (1), and it combines a moment method with a projection method on nite dimensional subspaces of L 2 (R, dx) := L 2 (dx). Nonparametric estimators of b and φ based on the observation of four i.i.d. sample paths (X (i) (t), t ∈ [0, T ]) distributed as [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF] in the stationary regime are built. Dening an adequate risk for the estimators is an additional diculty since b, φ are not square integrable in our framework. We propose a denition of the risk well tted to the problem which is as follows. For a couple (u, v) of functions, we set, when dened

(3) (u, v) 2 V := (u(x) + v ρ(x)) 2 ρ(x)dx.
The risk of an estimator ( b, φ) is dened by ( 4)

E( ( b -b, φ -φ) 2 V ).
Theorem 1 gives the risk bound of our projection estimators on a xed space. Section 4 is devoted to the nonparametric estimation of the invariant density which exhibits a parametric rate. In Section 5, we present numerical simulation results on several models for the estimation of the invariant density and of b and φ. In the last case, we experiment a model selection procedure. Section 6 contains concluding remarks and Section 7 is devoted to proofs.

Assumptions and preliminaries.

In the sequel, the notation means ≤ up to a constant. Our set of assumptions is essentially the same as the one given in Genon-Catalot and Larédo (2023b), [START_REF] Genon-Catalot | Parametric inference for ergodic McKean-Vlasov stochastic differential equations[END_REF]. Let us set

B(x) = x 0 b(y)dy, Φ(x) = x 0 φ(y)dy.
The function b is the potential term, the function φ is the interaction function. We assume that b ≡ 0 and that these functions satisfy:

• [H1] The function Φ is even. The functions x → B(x) and x → Φ(x) are C 2 and convex, one of the two being strictly convex: there exist constants K and λ such that

∀x, B (x) ≥ K ≥ 0, Φ (x) ≥ λ ≥ 0, K + λ > 0.
• [H2] The functions B (x) = b(x), Φ (x) = φ(x) are locally Lipschitz with polynomial growth, i.e. there exist c > 0,

∈ N * = {1, 2, 3, . . .} such that ∀x, y ∈ R, |b(x) -b(y)| + |φ(x) -φ(y)| ≤ c|x -y|(1 + |x| + |y| ).
• [H3] The functions b(x), b (x), φ(x), φ (x) have polynomial growth: there exists a constant c such that

|b(x)| + |b (x)| + |φ(x)| + |φ (x)| ≤ c(1 + |x| ).
According to [START_REF] Malrieu | Convergence to equilibrium for granular media equations and their Euler schemes[END_REF], [START_REF] Malrieu | Convergence to equilibrium for granular media equations and their Euler schemes[END_REF], Benachour et al. (1998a), [START_REF] Benachour | Nonlinear self-stabilizing processes -I Existence, invariant probability, propagation of chaos[END_REF], Cattiaux et al. (2008), [START_REF] Cattiaux | Probabilistic approach for granular media equations in the non-uniformly convex case[END_REF], Herrmann et al. (2008), [START_REF] Herrmann | Large deviations and a Kramers'type low for selfstabilizing diusions[END_REF], under assumptions [H1]-[H3], equation (1) admits a unique solution ((X(t), u t (x)dx)), t ≥ 0), where u t (x)dx = L(X(t)) is the marginal distribution of X(t). Model (1) admits a unique invariant distribution ρ such that R x 2 ρ(x)dx < +∞. And if the initial variable X(0) of (1) follows the invariant distribution ρ, then, for all t, L(X(t)) = L(X(0)). Thus, we assume

• [H4] The initial variable X(0) follows the invariant distribution ρ. Under [H1]-[H3], the invariant distribution ρ has density ρ(x) given as the solution of the implicit equation ( 5)

ρ(x) = 1 M exp - 2 σ 2 x 0 [b(u) + φ ρ(u)]du , where M = R exp{-2 σ 2 y 0 [b(u) + φ ρ(u)]du]}dy. As (Φ ρ) = φ ρ, x 0 φ ρ(u)du = Φ ρ(x) - Φ ρ(0), (6) ρ(x) ∝ exp [- 2 σ 2 (B(x) + Φ ρ(x))].
Note that if b ≡ 0, there is a one-parameter family of invariant distributions (see e.g. Genon-Catalot and Larédo (2023a), [START_REF] Comte | Nonparametric adaptive estimation for interacting particle systems[END_REF]). The invariant distribution is not unique unless its expectation is specied. We exclude this case here. Contrary to classical stochastic dierential equations, the invariant distribution is not explicit. Nevertheless, we can prove:

Proposition 1. Under [H1]-[H3], (7)
ρ(x) exp [- (K + λ) σ 2 x + b(0) + φ ρ(0) (K + λ) 2 ],
Moreover, for 0 ≤ z ≤ x and for x ≤ z ≤ 0,

ρ(x) ρ(z) ≤ exp [C(x -z)] where C = -(2/σ 2 )(b(0) + φ ρ(0)). (8)
Note that ( 7) is proved in Genon-Catalot and Larédo (2023b), [START_REF] Genon-Catalot | Parametric inference for ergodic McKean-Vlasov stochastic differential equations[END_REF]. Proposition 1 has the obvious consequence that the invariant distribution has moments of any orders and [START_REF] Benachour | Nonlinear self-stabilizing processes -II Convergence to invariant probability[END_REF] implies that, for all k ∈ R and all ε > 0,

(9) R exp (kx)(ρ(x)) ε dx < +∞.
Note that, by [H3], φ ρ has polynomial growth as ρ has moments of any order. Under [H1]-[H4], the initial variable η follows the invariant distribution ρ(x)dx, which implies that the distribution L(X(t)) = µ t (dx) satises ∀t ≥ 0, µ t (dx) = ρ(x)dx.

Therefore, equation (1) becomes: [START_REF] Butkovsky | On ergodic properties of nonlinear Markov chains and Mc-Kean-Vlasov equations[END_REF]

dX(t) = -(b(X(t)) + φ ρ(X(t)))dt + σdW (t), X(0) = η ∼ ρ(x)dx.
As noted in Genon-Catalot and Larédo (2023b), the following holds: Proposition 2. Assume [H1]-[H4] and consider the stochastic dierential equation

dY (t) = -(b(Y (t)) + φ ρ(Y (t)))dt + σdW (t).
Then (Y (t)) is a positive recurrent diusion with stationary density given by [START_REF] Belomestny | Semiparametric estimation of McKean-Vlasov SDEs[END_REF]

. If Y (0) ∼ ρ(x)dx, it is ergodic. Moreover, -If Y (0) = X(0), (Y (t)) ≡ (X(t)). -If Y (0) = X(0) = η ∼ ρ(x)dx, then X(t) = Y (t) for all t ≥ 0.
The result simply derives from the uniqueness of solutions but it has important consequences. In the stationary regime, the process (X(t)) given by ( 10) is identical to the classical SDE (Y (t)) in the stationary regime and is ergodic.

The innitesimal generator of the SDE (10) is given by: (11)

Lg = σ 2 2 g -(b + φ ρ)g = σ 2 2ρ (g ρ) .
The operator L acts on L 2 (ρ(x)dx) := L 2 (ρ) and is dened on the domain D, (

) D = {g ∈ L 2 (ρ), g abs. continuous, Lg ∈ L 2 (ρ), lim |x|→∞ g (x)ρ(x) = 0}. 12
(see e.g. Genon-Catalot et al. (2000), [START_REF] Genon-Catalot | Stochastic volatility models as hidden Markov models and statistical applications[END_REF]). Below, we use that, for all function Ψ ∈ D, LΨ(x)ρ(x)dx = 0. This relation obviously holds as, by ( 11)- [START_REF] Cattiaux | Probabilistic approach for granular media equations in the non-uniformly convex case[END_REF],

LΨ(x)ρ(x)dx = σ 2 2ρ (Ψ ρ) ρ = σ 2 2 (Ψ ρ) = [Ψ ρ] +∞ -∞ = 0.
In order to apply our estimation method, we need the following lemma.

Lemma 1. Consider a function ψ such that ψ ∈ L 2 (dx), ψ is absolutely continuous and bounded and ψ has polynomial growth. Then, x 0 ψ(t)dt and x 0 ψ ρ(t)dt both belong to D. Note that, if a function h has exponential growth, i.e. there exists a constant a > 0 such that, for all u ∈ R, |h(u)| ≤ ae au , then by Proposition 1 (see [START_REF] Benachour | Nonlinear self-stabilizing processes -II Convergence to invariant probability[END_REF] and ( 9)), h ∈ L 2 (ρ). Now, we can state a key property for the statistical strategy. Proposition 3. Assume [H1]-[H4] and that h has exponential growth. Set

h c = h- R h(x)ρ(x)dx
and dene the function Γ h (x) by Γ h (x) = x 0 (Γ h ) (y)dy where [START_REF] Chen | Maximum likelihood estimation of potential energy in interacting particle systems from single-trajectory data[END_REF] Then, Γ h ∈ D, LΓ h = h c , and (Γ h ) ∈ L 2 (ρ). Moreover, for (X(t)) satisfying [START_REF] Butkovsky | On ergodic properties of nonlinear Markov chains and Mc-Kean-Vlasov equations[END_REF], [START_REF] Comte | Drift estimation on non compact support for diusion models[END_REF]

(Γ h ) (y) = 2 σ 2 ρ(y) y -∞ h c (u)ρ(u)du = - 2 σ 2 ρ(y) +∞ y h c (u)ρ(u)du
) T 0 h c (X(s))ds = σ T 0 (Γ h ) (X(s))dW (s) + Γ h (X(0)) -Γ h (X(T )), (15)
Var ( 1

T T 0 h c (X(s))ds) ≤ 2σ 2 T [(Γ h ) (x)] 2 ρ(x)dx + 8 T 2 Γ 2 h (x)ρ(x)dx.
Note that relation [START_REF] Chen | Maximum likelihood estimation of potential energy in interacting particle systems from single-trajectory data[END_REF] holds since h c (x)ρ(x)dx = 0.

3. Estimation method for the potential and the interaction functions.

Notation. For h a function, we set h 2 = h 2 (x)dx, h 2 ρ = h 2 (x)ρ(x)dx. For x a vector in R d , x R d denotes the Euclidean norm of the vector. For a matrix M , M ⊥ denotes its transpose. For M a square matrix, Tr(M ) is the trace of M , and M op is the square-root of the largest eigenvalue of M M ⊥ . If M is symmetric, M op is simply equal to sup i {|λ i |} where the λ i are the eigenvalues of M . The symbol means ≤ up to a multiplicative xed constant.

3.1. Principle of the method. Following Pavliotis and Zanoni (2023), [START_REF] Pavliotis | A method of moments estimator for interacting particle systems and their mean eld[END_REF], we notice that, for a function ψ satisfying the assumptions of Lemma 1, it follows from the denition of ρ and (11)- [START_REF] Cattiaux | Probabilistic approach for granular media equations in the non-uniformly convex case[END_REF], that for Ψ(x) = x 0 ψ(y)dy, as LΨ(x)ρ(x)dx = 0, we have ( 16)

ψ(x)b(x)ρ(x)dx + ψ(x)φ ρ(x)ρ(x)dx = 1 2 ψ (x)σ 2 ρ(x)dx.
To estimate b and φ, we must take into account that the functions b, φ and φ ρ do not belong to L 2 (dx). By [H3], these three functions have polynomial growth and therefore, as ρ has moments of any order, belong to L 2 (ρ). So we proceed as follows.

Consider two orthonormal bases of L 2 (dx), (ϕ j ) j≥0 , (θ ) ≥0 composed of functions derivable, bounded, with derivatives having polynomial growth, thus satisfying the assumptions of Lemma 1. We denote for τ = ϕ and τ = θ:

(17) L τ (q) := sup x∈R q-1 k=0 τ 2 k (x),
and assume that L ϕ (m) < +∞, L θ (p) < +∞ for all xed m ≥ 1, p ≥ 1. The quantity L τ (q) only depends on the space generated by (τ 0 , . . . , τ q-1 ). Indeed,

L τ (q) = sup h∈span(τ 0 ,...,τ q-1 ) h 2 ∞ / h 2 .
It is nite for classical examples of bases. Generally, for all q, L τ (q) ≤ c τ q α with c τ , α constants linked with the basis (see Section 3.5). As ρ is bounded, for all j, ϕ j belongs to L 2 (ρ). For m ≥ 1, the functions ϕ 0 , . . . , ϕ m-1 are linearly independent and generate a m-dimensional space

S m = Vect(ϕ 0 , . . . , ϕ m-1 )
which is a subspace of L 2 (ρ). Analogously, we dene the p-dimensional subspace of L 2 (ρ) Σ p = Vect(θ 0 , . . . , θ p-1 ).

Recall that, by denition (3), for a couple (u, v) of functions, we have set, when dened,

(u, v) 2 V = u + v ρ 2 ρ . Note that (u, v) 2 V is nite for (u, v) = (b, φ) and for (u, v) ∈ S m × Σ p . Indeed for (u, v) ∈ S m × Σ p , u 2 ρ ≤ ρ ∞ u 2 and by the Young inequality, v ρ ≤ ρ 1 v = v , so that v ρ belongs to L 2 (dx) ⊂ L 2 (ρ). Note also that (u, v) → (u, v) V is only a semi- norm on L 2 (ρ) × L 2 (ρ). Indeed, (u, v) V = 0 is equivalent to u(x) + v ρ(x) = 0 for all x ∈ R.
But this does not always implies u ≡ 0, v ≡ 0. Nevertheless, we denote (u 1 , v 1 ), (u 2 , v 2 ) V the scalar product associated to the semi-norm . V . Now, we dene the V -orthogonal projection of (b, φ) on S m × Σ p : [START_REF] Comte | Non compact estimation of the conditional density from direct or noisy data[END_REF]

(b V m , φ V p ) = arg min (u,v)∈Sm×Σp (b -u, φ -v) 2 V . Setting b V m = m-1 k=0 b k ϕ k and φ V p = p-1
=0 c θ , we obtain [START_REF] Dalalyan | Asymptotic statistical equivalence for scalar ergodic diusions[END_REF] V m,p (b, c) m,p = Z m,p using the notation

(b, c) m,p = (b 0 , . . . , b m-1 , c 0 , . . . , c p-1 ) ⊥ = b m c p ∈ M m+p,1 (R)
with b m = (b 0 , . . . , b m-1 ) ⊥ and c p = (c 0 , . . . , c p-1 ) ⊥ , and where V m,p is the

(m + p) × (m + p) matrix, element of M m+p,m+p (R), ( 20
) V m,p =       V 1 = ϕ j ϕ k ρ 0≤j,k≤m-1 V 1,2 = θ ρ ϕ k ρ 0 ≤ k ≤ m -1 0 ≤ ≤ p -1 V 2,1 = θ r ρ ϕ j ρ 0 ≤ r ≤ p -1, 0 ≤ j ≤ m -1 V 2 = θ ρ θ r ρ ρ 0≤ ,r≤p-1       , where V 2,1 = (V 1,2 ) ⊥ , (21)
Z m,p =     (b + φ ρ)ϕ k ρ 0≤k≤m-1 (b + φ ρ)(θ ρ)ρ 0≤ ≤p-1     ∈ M m+p,1 (R).
As the functions ψ = ϕ 0 , . . . , ϕ m-1 , θ 0 ρ, . . . , θ p-1 ρ satisfy the assumptions of Lemma 1, we may apply formula [START_REF] Comte | Nonparametric adaptive estimation for interacting particle systems[END_REF] to these functions and we get, using that

(θ ρ) = θ ρ, (22)
Z m,p = σ 2 2     ϕ k ρ 0≤k≤m-1 θ ρ ρ 0≤ ≤p-1     ∈ M m+p,1 (R).
3.2. Theoretical constraints and identiability. Below, vectors are denoted using bold letters and coordinates or functions are denoted with usual letters (not bold). For (x 0 , . . . , x m-1 , y 0 , . . . , y p-1 ) ∈ R m+p , we denote [START_REF] Della Maestra | The LAN property for McKean-Vlasov models in a mean-eld regime[END_REF] (x, y) m,p = (x 0 , . . . , x m-1 , y 0 , . . . , y p-1 ) ⊥ = x m y p with x m = (x 0 , . . . , x m-1 ) ⊥ , y p = (y 0 , . . . , y p-1 ) ⊥ . The vector (b, c) m,p in ( 19) is uniquely dened if and only if the matrix V m,p is invertible. We emphasize that the matrix V m,p is symmetric and nonnegative. Indeed,

(x, y) ⊥ m,p V m,p (x, y) m,p =   m-1 j=0 x j ϕ j (x) + p-1 =0 y θ ρ(x)   2 ρ(x)dx ≥ 0.
Note that

(x, y) ⊥ m,p V m,p (x, y) m,p = (u, v) 2 V for u = m-1 j=0 x j ϕ j , v = p-1 =0 y θ l .
As, by [START_REF] Benachour | Nonlinear self-stabilizing processes -I Existence, invariant probability, propagation of chaos[END_REF], the support of ρ is R, V m,p is invertible if:

(24)    ∀x ∈ R, m-1 j=0 x j ϕ j (x) + p-1 =0 y θ ρ(x) = 0    ⇒ (x, y) m,p = (0, 0) m,p
Condition ( 24) is an identiability constraint linked with the choice of the bases. In Section 3.5, we provide an example of bases for which [START_REF] Eberle | Quantitative Harris type theorems for diusions and Mc-Kean-Vlasov processes[END_REF] is fullled. Nevertheless, [START_REF] Eberle | Quantitative Harris type theorems for diusions and Mc-Kean-Vlasov processes[END_REF] is no enough for our purpose and we need to reinforce the identiability constraint and set the following assumption:

• [H5] There exists a constant f 0 > 0 xed and not depending on m, p, such that

∀x ∈ R m+p , x R m+p = 1, x ⊥ V m,p x ≥ f 0 > 0.
Assumption [H5] ensures that, for all (m, p), the eigenvalues of V m,p are all larger than f 0 . As a consequence, under [H5], there exists a positive constant s 0 independent of m, p such that

V -1 m,p op ≤ s 0 .
Clearly s 0 = 1/f 0 suits.

b m (x) = m-1 j=0 b j ϕ j (x), φ p (x) = p-1 =0 c θ (x), x ∈ R.
We assume that we have at our disposal four independent paths of (10), (X (i) (t)) t∈[0,T ] , for i = 1, 2, 3, 4, i.e.

dX (i) (t) = -(b(X (i) (t)) + φ ρ(X (i) t)))dt + σdW i (t), X (i) (0) = η i ∼ ρ(x)dx,
where (W i ) are independent Brownian motions, (η i ) are independent random variables, independent of the Brownian motions (W i ). For m, p ≥ 0, we dene an empirical version of the theoretical matrix V m,p and of the vector vector Z m,p , denoted by V m,p and Z m,p . The matrix V m,p is built using the two sample paths (X (i) (t)) t∈[0,T ] , for i = 1, 2 and the vector Z m,p is built using the two sample paths (X (i) (t)) t∈[0,T ] , for i = 3, 4. In this way, the estimators V m,p and Z m,p are independent. We could use only two sample paths to estimate both, but this independence simplies considerably the study. Let us now dene V m,p and Z m,p . For this, we introduce the following notation. For a function h, let us set

h(x) = 1 T T 0 h(x -X (2) (t))dt, h(x) = 1 T T 0 h(x -X (4) (t))dt.
Note that Eh(x) = Eh(x) = h ρ(x) and h ρ(X (1) (s)) can be estimated by h(X (1)

(s)) = 1 T T 0 h(X (1) (s) -X (2) (t)
)dt, using two trajectories. Although h(x), h(x) are random and depend respectively on X (2) and X (4) , we omit this dependence in the notation for the sake of simplicity. We set

V m,p = V 1 V 1,2 ( V 1,2 ) ⊥ V 2 with V 1 = 1 T T 0 ϕ j (X (1) (s))ϕ k (X (1) (s))ds 0≤j,k≤m-1 V 1,2 = 1 T T 0 θ (X (1) (s)) ϕ k (X (1) (s))ds 0≤k≤m-1,0≤ ≤p-1 V 2 = 1 T T 0 θ (X (1) (s)) θ r (X (1) (s))ds 0≤ ,r≤p-1
.

The matrix V m,p is built using the sample paths

(X (i) (t)) t∈[0,T ] , for i = 1, 2 as for = 0, . . . , p-1, θ (X (1) (s)) = 1 T T 0 θ (X (1) (s) -X (2) (t))dt.
Next, we set

Z m,p = σ 2 2      1 T T 0 ϕ k (X (3) (s))ds 0≤k≤m-1 1 T T 0 θ (X (3) (s))ds 0≤ ≤p-1     
, which depends on the sample paths (X (i) (t)) t∈[0,T ] , for i = 3, 4, as, for = 0, . . . , p -1,

θ (X (3) (s)) = 1 T T 0 θ (X (3) (s) -X (4) (t)))dt.
Moments strategies often bring unbiased estimators. Here, due to the use of the functions h(x), h(x) which are unbiased estimators of h ρ, it is not exactly the case, and we can prove:

Proposition 4. Under Assumptions [H1]-[H4], it holds that E( Z m,p ) = Z m,p , E( V m,p ) = V m,p +     0 m×m 0 m×p 0 p×m O 1 T 1 I p×p     .,
where 1 I p×p denote the p × p matrix with all coecients equal to 1.

Thus Z m,p is unbiased and V m,p is asymptotically unbiaised when T grows to innity.

Provided that V m,p is invertible, using notation [START_REF] Della Maestra | The LAN property for McKean-Vlasov models in a mean-eld regime[END_REF], relation [START_REF] Dalalyan | Asymptotic statistical equivalence for scalar ergodic diusions[END_REF] suggests to dene the estimator of (b, c) m,p by

(b, c) m,p = V -1 m,p Z m,p = b m c p .
Although under [H5], the theoretical matrix V m,p is invertible with minimum eigenvalue away from 0, this does not guaranteee the invertibility of V m,p . Therefore, to get a proper denition of our estimator, we propose the standard cuto strategy leading to replace V -1 m,p by ( 26)

V -1 m,p := V -1 m,p 1 ∆m,p with ∆ m,p = { V -1 m,p op ≤ 2s 0 }.
Then we set ( 27)

(b, c) m,p = V -1 m,p Z m,p = b m c p .
where

b m = ( b 1 , . . . , b m-1 ) ⊥ , c p = ( c 1 , . . . , c p-1 ) ⊥ .
The estimator of (b, φ) is thus given by

b m (x) = m-1 j=0 b j ϕ j (x), φ p (x) = p-1 =0 c θ (x).
By construction, V -1 m,p and Z m,p are independent. The risk of the estimator ( b m , φ p ) is dened by: E( ( b m -b, φ p -φ) 2

V ), (see ( 3), ( 4)).

3.4. Study of estimator risk. In the sequel, we use the following notations.

Let (b, φ)(x) = (b(x), φ(x)) and set (b, φ) m,p (x) = (b V m (x), φ V p (x))
where b V m and φ V p given by [START_REF] Comte | Non compact estimation of the conditional density from direct or noisy data[END_REF]. Analogously, let

(b, φ) m,p (x) = ( b m (x), φ p (x)),
and recall that (see ( 3)

) (b, φ) 2 V = (b + φ ρ 2 ρ .
For (x, y) m,p =

x m y p a vector of R m × R p , dene the square norm associated with V m,p dened by [START_REF] Dalalyan | Asymptotic statistical equivalence for ergodic diusions: the multidimensional case[END_REF], by

(x, y) m,p 2
Vm,p := (x, y) ⊥ m,p V m,p (x, y) m,p =   m-1 j=0 x j ϕ j (x) + p-1 =0 y θ ρ(x)   2 ρ(x)dx ≤ 2 ρ ∞     m-1 j=0 x j ϕ j (x)   2 dx + p-1 =0 y θ ρ(x) 2 dx   .
According to Young's inequality, we have

p-1 =0 y θ ρ 2 ≤ ρ 2 1 p-1 =0 y θ 2 = p-1 =0 y θ 2 .
Now using that the functions (ϕ j ) 0≤j≤m-1 and (θ ) 0≤ ≤p-1 are orthonormal in L 2 , we get

(x, y) m,p (28)
Vm,p ≤ 2 ρ ∞ (x, y) m,p 2 R m+p . 2
Now, we can link the V -norm of a couple of functions to the . Vm,p -norm of a vector of

R m ×R p . Indeed, consider two functions ω(x) = m-1 j=0 ω j ϕ j (x) ∈ S m and ζ(x) = p-1 =0 ζ θ (x) ∈ Σ p with coecients respectively ω m = (ω j , 0 ≤ j ≤ m -1) ⊥ and ζ p = (ζ , 0 ≤ ≤ p -1) ⊥ . We have (ω, ζ) 2 V = (ω, ζ) ⊥ m,p V m,p (ω, ζ) m,p = (ω, ζ) m,p 2
Vm,p .

Thus, using (28), we have

(ω, ζ) 2 V ≤ 2 ρ ∞ (ω, ζ) 2 where (ω, ζ) 2 = m-1 j=0 ω 2 j + p-1 =0 θ 2 = (ω, ζ) m,p 2 R m+p .
The following theorem gives a bound for the risk of the estimator

(b, φ) m,p (x) = ( b m (x), φ p (x))
with xed m, p.

Theorem 1. Assume that Assumptions [H1]-[H5] hold and that

L ϕ (m)+L θ (p) ≤ T , with L ϕ (m)
and L θ (p) dened in [START_REF] Comte | Regression function estimation as a partly inverse problem[END_REF]. Then

E (b, φ) m,p -(b, φ) 2 V ≤ inf (u,v)∈Sm×Σp (b -u, φ -v) 2 V + C L ϕ (m) + L θ (p) T ,
where C is a constant depending on b, φ, ρ, s 0 .

Let us comment this result. The risk bound shows a decomposition into a bias term

inf (u,v)∈Sm×Σp (b -u, φ -v) 2
V and a variance term (L ϕ (m) + L θ (p))/T (see [START_REF] Comte | Regression function estimation as a partly inverse problem[END_REF]). This result justies the choice of the V -norm to dene the risk. The bias term decreases with (m, p) and the variance term increases with m, p. This decomposition may be used to realize the square bias-variance compromise provided that the rate of the bias term can be assessed on specic regularity spaces.

We need the following Propositions 5, 6 and Lemma 2 to prove Theorem 1.

Proposition 5. Under [H1]-[H4], it holds

E[ Z m,p -Z m,p 2 R m+p ] L ϕ (m) + L θ (p) T . Proposition 6. Under [H1]-[H4], E Tr( V m,p -V m,p ) 2 L ϕ (m) T
.

Therefore, we have

E V m,p -V m,p 2 op ≤ E Tr( V m,p -V m,p ) 2 L ϕ (m) T .
Note that this bound does not depend on p.

Lemma 2. Under [H1]-[H5]

,

E V -1 m,p -V -1 m,p 2 op 1 ∆m,p ≤ (40 + 1/4)s 4 0 E V m,p -V m,p 2 op
Consequently,

E V -1 m,p -V -1 m,p 2 op 1 ∆m,p L ϕ (m) T .
Remark 1. Following ideas in Genon-Catalot and Larédo (2023a,b), [START_REF] Genon-Catalot | Inference for ergodic McKean-Vlasov stochastic dierential equations with polynomial interactions[END_REF], [START_REF] Genon-Catalot | Parametric inference for ergodic McKean-Vlasov stochastic differential equations[END_REF], we might have proposed (1) (s)ds instead of the term in V 2 , with ρ [T ] computed with the path X (2) as in Section 4, formula [START_REF] Genon-Catalot | Parametric inference for small variance and long time horizon McKean-Vlasov diusion models[END_REF]. It could be used in the denition of Z m,p also, but would yield additional bias, compared to our strategy. In other words h(x) is an unbiased estimator of h ρ(x) which is not the case of h ρ [T ] (x).

T -1 T 0 ϕ ρ [T ] (X (1) (s)ϕ r ρ [T ] (X
3.5. Identiability and risk bound when using the Hermite basis. Let us discuss the respective orders of the bias and variance terms in the bound obtained in Theorem 1.

The variance term order depends on the bases. Let us consider the Hermite basis which is dened as follows.The Hermite polynomial and the Hermite function of order j are given, for j ≥ 0, by: (29)

H j (x) = (-1) j e x 2 d j dx j (e -x 2 ), h j (x) = c j H j (x)e -x 2 /2 , c j = 2 j j! √ π -1/2
The sequence (h j , j ≥ 0) is an orthonormal basis of L 2 (dx). Hermite functions with odd (resp. even) index are odd (resp. even) and h *

j (x) = √ 2πi j h j (x)
, where h * j is the Fourier transform of h j . Moreover (see Abramowitz and Stegun (1964, 22.14.17), [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF], [START_REF] Szegö | Orthogonal polynomials[END_REF], [START_REF] Szegö | Orthogonal polynomials[END_REF] p.242, Indritz (1961), [START_REF] Indritz | An inequality for Hermite polynomials[END_REF]), h j ∞ ≤ Φ 0 , Φ 0 1/π 1/4 0.7511, so that the basis (ϕ j = h j ) is constituted of bounded functions which satisfy the assumptions of Lemma 1.

Moreover, the identiability condition ( 24) is fullled by the Hermite basis. Indeed, replacing in [START_REF] Eberle | Quantitative Harris type theorems for diusions and Mc-Kean-Vlasov processes[END_REF] ϕ j by h j and θ by h 2 -1 , and applying Fourier transform yields [START_REF] Genon-Catalot | Probabilistic properties and parametric inference of small variance nonlinear self-stabilizing stochastic dierential equations[END_REF]

∀x ∈ R, m-1 j=0 y j i j h j (x) + p-1 =0 y +m i 2 +1 h 2 +1 (x)ρ * (x) = 0.
Multiplying by e x 2 /2 gives P (x) + Q(x)ρ * (x) = 0, ∀x ∈ R, where P and Q are polynomials. Let us moreover assume that V (.) and W (.) are C ∞ and let f (n) (x) the n-th derivative of f . Inequality [START_REF] Benachour | Nonlinear self-stabilizing processes -II Convergence to invariant probability[END_REF] yields that ρ belongs to the Schwarz class

S(R) = {f ∈ C ∞ (R, R), s.t.∀m ≥ 1, ∀n ≥ 0, f (n) (x) = O(x -m )
at ∞}. Therefore ρ also belongs to S(R). Letting x → ∞ in (30) yields that P (x) → 0, which implies that the term with highest degree of P (x), i.e. of y m-1 i m-1 H m-1 (x), tends to 0. Therefore y m-1 = 0. Iterating the procedure yields successively that y m 1 , y m-2 , . . . , y 0 are all equal to 0. Equation ( 30) reduces to the equation:

∀x ∈ R, p-1 =0 y +m i 2 +1 h 2 +1 (x) = 0.
The functions h j (.) being orthogonal, this implies that y +m = 0 for = 0, . . . , p -1. Lastly, it is proved in Lemma 1 of Comte and Lacour (2023), [START_REF] Comte | Non compact estimation of the conditional density from direct or noisy data[END_REF]

that L ϕ (m) ≤ C √ m.
As a consequence, if we use Hermite functions (h 0 , h 1 , . . . , h m-1 ) to span S m and (h 1 , h 3 , . . . , h 2p-1 ) to span Σ p (recall that φ is odd), the variance term has order

L ϕ (m) + L θ (p) T ≤ C √ m + √ p T .
Note that the assumptions of Lemma 1 imply that we cannot use compactly supported bases whose functions are not derivable on R.

Let us now discuss the square bias term. If the functions b and φ were square integrable on R, we would upper bound the square bias by

ρ ∞ inf u∈Sm b -u 2 + inf v∈Σp φ -v 2 = ρ ∞   j≥m b, ϕ j 2 + ≥p φ, θ 2 +1 2  
For b and φ in Sobolev-Hermite spaces with regularity s 1 and s 2 (see e.g. Comte and Genon-Catalot (2020), section 3.4, [START_REF] Comte | Regression function estimation as a partly inverse problem[END_REF]), the resulting order would be m -s 1 + p -s 2 and the resulting bias rate would be T -2s /(2s +1) where s = min(s 1 , s 2 ). This is the optimal rate for estimating a function with regularity s . Thus, our bias term is meaningful. However, the order of our specic bias term taking into account that b, φ, φ ρ are in L 2 (ρ) and not in L 2 (dx) would require the denition of specic regularity spaces to assess the rate of

inf u∈Sm b -u 2 ρ + inf v∈Σp (φ -v) ρ 2 ρ .

Estimation of the stationary density

As ρ ∈ L 2 (dx), we consider an orthonormal basis of L 2 (dx), still denoted (ϕ j ) j≥0 and write ρ(x) = j≥0 a j ϕ j (x) with a j = ϕ j , ρ = ϕ j (x) ρ(x)dx. For an integer D, we consider the estimator [START_REF] Genon-Catalot | Parametric inference for small variance and long time horizon McKean-Vlasov diusion models[END_REF]

ρ D = D-1 j=0 a j ϕ j , a j = 1 T T 0 ϕ j (X(s))ds.
We can prove the following result Proposition 7. Under Assumptions [H1]-[H4],

E ρ D -ρ 2 ≤ ρ -ρ D 2 + C T ,
where ρ D = D-1 j=0 a j ϕ j , C is a positive constant (not depending on D, T ) and . is the L 2 -norm.

It follows from Proposition 7 that the variance of the estimator does not depend on D. Therefore, taking D as large as possible will make the bias term negligible under weak regularity conditions, and ρ can be estimated with (parametric) rate 1/T . In Comte and Merlevède (2005), the same result is obtained under specic conditions, mainly the so-called Castellana and Leadbetter condition which is dicult to check in general. The interesting point here is that the parametric rate is obtained without such specic condition in our framework.

Remark 2. In most cases, as mentioned above, the invariant density has no explicit expression. However, in the special case where b is odd and φ is linear, ρ is explicit. We follow Hermann and Tugaut (2010) to explain this case. Let us assume that b is odd and φ(x) = cx. Then, φ ρ(x) = c(x -m) where m is the expectation of ρ. The invariant density ρ is given by

ρ(x) = [K(m)] -1 exp [-2σ -2 (B(x) + c 2 (x -m) 2 )]
where [K(m)] -1 is the norming constant and m is solution of

x exp [-2σ -2 (B(x) + c 2 (x -m) 2 ]dx = m exp [-2σ -2 (B(x) + c 2 (x -m) 2 ]dx = mK(m)
As B is even, this equation has an obvious solution m = 0. As the invariant distribution is unique, we nd

ρ(x) = exp [-2σ -2 (B(x) + c 2 x 2 )][ exp [-2σ -2 (B(y) + c 2 y 2 )]dy] -1 .
Moreover, equation [START_REF] Butkovsky | On ergodic properties of nonlinear Markov chains and Mc-Kean-Vlasov equations[END_REF] is equal to

dX(t) = -[b(X(t)) + cX(t)]dt + σdW (t), X(0) ∼ ρ(x)dx.
Models 1,2,3 presented below in Section 5 correspond to such cases and we use them to evaluate the performances of the invariant density estimator.

Simulation experiments

We implemented the method by simulating N = 200 particules following equation ( 2) with Euler-type discretization scheme with step ∆ = 0.1. The initial value is taken equal to 0 and we exclude the 50 rst values of the process to reach approximately the stationary regime. Two values of T presented below are T = 200 and T = 2000, involving n = 2000 and n = 20000 generated X N i (j∆), j = 1, . . . , n and i = 1, . . . , N . In all cases, we took σ = 0.25. The estimation uses the sample paths X N i (j∆), j = 1, . . . , n for i = 1, . . . , 4. First, we present 3 models for which we can compute the true density ρ so that we can illustrate the performance of the density estimator described in Section 4 (see Remark 2). With b = 0.5, c = 0.5, we consider (1

) b(x) = bx, φ(x) = cx, ρ(x) is a N (0, σ 2 /[2(b + c)]), ( 2
) b(x) = b tanh(x), φ(x) = cx, ρ(x) = c 1 exp(-cx 2 /σ 2 )/(cosh(x)) 2b/σ 2 , ( 3
) b(x) = bx 3 , φ(x) = cx, ρ(x) = c 2 exp -(bx 4 /2 + cx 2 )/σ 2 .
The constants c 1 and c 2 are computed numerically so that the function ρ integrates to one. We compute estimators corresponding to the Hermite basis (ϕ j = h j ) 0≤j≤D-1 (see [START_REF] Genon-Catalot | Stochastic volatility models as hidden Markov models and statistical applications[END_REF]). In all cases, the true is in bold red.

choice of the parameters b, c seems important as it inuences the simulation results. Note that if D is too large (D = T , e.g.), it implies numerical problems. We can see on Figure 1 that the method works very well and is improved when T grows. The three densities, though dierent, have very similar forms. Figure 2 presents the estimation of b and φ for Example 2. The detail about our protocol of estimation is detailed below. Two other models are illustrated hereafter where b is not odd and/or φ is not linear:

(4

) b(x) = 2(x -1) 3 + 2x, φ(x) = 4 tanh(x), ( 5
) b(x) = 2x, φ(x) = 2(x + x 3 ).
To compute the estimators of b, we use the Hermite basis (ϕ j ) 0≤j≤m-1 for m = 1, . . . , 6 (see Section 3.5). To compute the estimators of φ, we need to take into account that φ is odd. So we consider the basis (θ ) 0≤ ≤p-1 = (h 2 +1 ) 0≤ ≤p-1 as Hermite functions of odd index are odd functions. We consider dimensions m ∈ {1, . . . , m max = 6} and p ∈ {1, . . . , p max = 3}. Larger proposals systematically involve estimators with an obvious variance eect. The fact that we consider small dimensions replaces the cuto of V -1 m,p . For examples (1) to (3), we could have simulated by Euler scheme 4 independent sample paths of dX(t) = -(b(X(t)) + cX(t))dt + σdB(t) instead of the particles system. We found that the estimation results were of the same type, so we kept the particles system which is suitable for all models.

Though we have no theoretical result concerning an adaptive choice for the dimensions m, p, we propose a criterion following the standard method to select data-driven dimensions. The selection of (m, p) is done by choosing the couple which minimizes

Crit(m, p) = -(b, c) ⊥ m,p Z m,p + κ × √ m + 4 √ p T ,
where (b, c)

⊥ m,p Z m,p = (b, c) m,p 2
Vm,p

= Z ⊥ m,p V -1
m,p Z m,p . The criterion is inspired by an estimation of the empirical bias, as usually performed for least-squares contrast minimization, and a second term which has the order of the variance of the estimator. The constant κ is roughly calibrated and chosen here as κ = 5.

Figure 2 shows 20 estimators computed for 20 independent systems of particules of particles, in the case of Examples 2, for T = 2000. Figure 3 illustrates, for example 4, the improvement which occurs when increasing T from 200 to 2000: it is less striking than in the density case, but still existing. Figure 4 shows, for example 5, on the rst line all the proposals which are considered for b and φ for all couples (m, p) of the list. We can see that the selection of a relevant estimator is crucial. As Hermite bases are parsimonious, small values of m and p perform well. The choice m = 1 is generally not good, and the most selected couple for (m, p) is (2, 1).

To conclude, the performance of the density estimator in nice, but the problem of estimating b and φ is obviously more dicult. In all cases, the choice of the functions b, φ is important because too small or too large ranges of values of the process can generate numerical diculties.

6. Concluding remarks 6.1. Summary and remaining questions. In this paper, we consider the nonparametric estimation of the functions b, φ of model [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF] based on the continuous observation on a time interval [0, T ] of four i.i.d. sample paths in stationary regime. The asymptotic framework is T → +∞. Due to the presence of the interaction term in equation ( 1) which contains a convolution with the unknown marginal law of the process, the problem is completely dierent from the case of a usual diusion (φ ≡ 0) and the observation of only one sample path is not enough to estimate both b and φ. Actually, two sample paths could be enough but induces much more tedious computations. The estimation method is inspired by a paper in parametric setting of [START_REF] Pavliotis | A method of moments estimator for interacting particle systems and their mean eld[END_REF]. To build the estimators, we use a combination of a moment method and a projection method. Dening a specic risk for the estimators tted to the problem, we obtain a bound for the risk of our projection estimators involving as usual a variance term and a square bias term. The order of the variance term is precised for the Hermite basis. The rate of the square bias term which is very specic to the problem and is not simple to evaluate. We sketch a discussion on this point. Note that no benchmark is available for comparison. We also provide a nonparametric estimator of the invariant density of the process which has the parametric rate without any condition such as Castellana and Leadbetter's one. The method is implemented on simulated data using an empirical model selection criterion, and performs reasonably well.

Remaining questions are worth of further research. It would be of interest to propose a theoretical background for assessing the rate of the bias term and for the selection procedure. Moreover, the estimation procedure is based on the estimation of the vector Z m,p suggested by formula [START_REF] Della Maestra | Nonparametric estimation for interacting particle systems: McKean-Vlasov models[END_REF], using the moment relation [START_REF] Comte | Nonparametric adaptive estimation for interacting particle systems[END_REF]. But the vector Z m,p is also given by ( 21) and could have been estimated by

Z m,p = 1 T T 0 ϕ k (X (3) (s))dX (3) (s) ⊥ 0≤k≤m-1 , 1 T T 0 θ (X (3) (s))dX (3) (s) ⊥ 0≤ ≤p-1 ⊥
keeping the same matrix V m,p . The other estimator V -1 m,p Z m,p requires another study which is worth being done. 6.2. Estimation from the observation of the particle systems. Suppose that the observation is the particles system (2) for i = 1, . . . , 4, t ∈ [0, T ] and X N i (0) ∼ ρ(x)dx with N, T large. Then, we can apply the same estimation method to obtain estimators of b and φ. This is illustrated in Section 5 as the easiest way to simulate four independent sample paths X (i) (t), i = 1, . . . , 4 is to simulate the particles system (2) for large N . Of course, the theoretical study would have to take into account this point by adding a bias term due to the convergence of the particle system to the limit processes. 6.3. Discrete observations. In practice, only discrete time observations are available. Therefore, it is worth of interest to study the same estimation problem based on discrete observations (X (i) (j∆), j = 1 . . . , n, i = 1, . . . , 4) with ∆ = ∆ n tending to 0 and n∆ n tending to innity (high frequency framework). One can consider the discretized versions of our estimators and this is actually what is done in the simulation section. The extension of the estimation results would be to the price of some more computations and some additional bias terms under conditions on n, ∆ n . The usual condition for diusions, n∆ 2 n = o(1), would probably be the same.

7. Proofs 7.1. Proof of Proposition 1. We have b (x) ≥ K, φ (x) ≥ λ. Thus, (φ ρ) (x) = φ (x - y)ρ(y)dy ≥ λ. Therefore, for x ≥ 0, b(x) ≥ Kx + b(0), φ ρ(x) ≥ λx + φ ρ(0).

This implies, for

0 ≤ z ≤ x, x z [b(u) + φ ρ(u)]du ≥ (K + λ)(x 2 -z 2 )/2 + (b(0) + φ ρ(0))(x -z).
Thus,

ρ(x) ρ(z) ≤ exp [- (K + λ) σ 2 (x 2 -z 2 ) + C(x -z)], C = -(2/σ 2 )(b(0) + φ ρ(0)).
Analogously, for x ≤ z ≤ 0,

ρ(x) ρ(z) ≤ exp [- (K + λ) σ 2 (z 2 -x 2 ) + C(x -z)].
Therefore, for all x,

ρ(x) exp [- (K + λ) σ 2 x 2 + Cx] ∝ exp - (K + λ) σ 2 x - Cσ 2 2(K + λ) 2 , Moreover, for 0 ≤ z ≤ x or x ≤ z ≤ 0, ρ(x) ρ(z) ≤ exp [C(x -z)]. 2 7.2. Proof of Lemma 1. Set Ψ(x) = x 0 ψ(t)dt.
We prove that Ψ and . 0 ψ ρ both belong to D dened by [START_REF] Cattiaux | Probabilistic approach for granular media equations in the non-uniformly convex case[END_REF]. We have

Ψ 2 (x)ρ(x)dx = +∞ 0 ... + 0 -∞ ..., where +∞ 0 Ψ 2 (x)ρ(x)dx ≤ +∞ 0 x x 0 ψ 2 (t)dtρ(x)dx = 0<t<x xψ 2 (t)ρ(x)dxdt = +∞ 0 ψ 2 (t)( +∞ t xρ(x)dx)dt ≤ +∞ 0 ψ 2 (t)dt +∞ 0 xρ(x)dx < +∞.
We proceed analogously for the integral on (-∞, 0) and thus Ψ ∈ L 2 (ρ).

To prove that x 0 ψ ρ(t)dt belongs to L 2 (ρ), it is therefore enough to prove that ψ ρ belongs to L 2 (R). By the Young inequality,

(ψ ρ) 2 (x)dx ≤ ( ρ(x)dx) 2 ψ 2 (x)dx < +∞.
Next, we must prove that ψ(x)ρ(x) and ψ ρ(x)ρ(x) tend to 0 as x tends to +∞ and -∞. This holds as ψ is bounded by assumption. And ψ ρ is also bounded as

|ψ ρ(x)| ≤ ψ 2 (t)dt ρ 2 (y)dy 1/2 ≤ ρ ∞ ψ 2 (t)dt 1/2
. Now, we prove that LΨ and L( . 0 (ψ ρ) belong to L 2 (ρ). We have

LΨ = σ 2 2 ψ -(b + φ ρ)ψ.
As ψ has polynomial growth, [ψ (x)] 2 ρ(x)dx < +∞ by Proposition 1. By our assumptions, b = V and φ ρ have polynomial growth. So, it is enough to check that (1 + x )ψ(x) belongs to L 2 (ρ) which holds since ψ is bounded and ρ has moments of any order. We have L( .

0 (ψ ρ) = σ 2 2 ψ ρ -(b + φ ρ)ψ ρ.
We use that ψ ρ is bounded and b + φ ρ has polynomial growth to deduce that (b + φ ρ)ψ ρ belongs to L 2 (ρ). Next, as ψ has polynomial growth, ψ ρ has also polynomial growth. This implies that ψ ρ belongs to L 2 (ρ). 2 7.3. Useful Lemma. We state here a useful Lemma. Lemma 3. Assume [H1]-[H4] and consider a function h such that there exists a positive constant a such that ∀u ∈ R, |h(u)| ≤ ae au . The following integrals are nite and bounded by constants depending on ρ and h only:

(1) I -= 0 -∞ y -∞ h(u)ρ(u)du 2 dy ρ(y) , I + = +∞ 0 +∞ u h(u)ρ(u)du 2 dy ρ(y) .
(

) J -= x<y<0 |x|ρ(x) ρ 2 (y) y -∞ h(u)ρ(u)du 2 dxdy, J + = 0<y<x xρ(x) ρ 2 (y) +∞ y h(u)ρ(u)du 2 dxdy. (3) 0 -∞ ( x -∞ ρ 2 (y)dy) dx ρ(x) , +∞ 0 ( +∞ x ρ 2 (y)dy) dx ρ(x) . ( 2
) +∞ 0 xρ(x)dx x 0 du ρ 2 (u) +∞ u ρ 2 (y)dy, 0 -∞ |x|ρ(x)dx 0 x du ρ 2 (u) u -∞ ρ 2 (y)dy. ( 4
) +∞ 0 xρ(x)dx x 0 du ρ 2 (u) +∞ u ρ(y)dy 2 , 0 -∞ |x|ρ(x)dx 0 x du ρ 2 (u) u -∞ 5

ρ(y)dy

Proof of Lemma 3. We start with the rst item. It is enough to look at I -, the other one being analogous. We have:

I -= 2 u<v<y<0 dy ρ(y) h(u)ρ(u)h(v)ρ(v)dudv.
On the integration set,

ρ(u)ρ(v) ρ(y) = [ρ(u)ρ(v)] 1/2 [ ρ(u) ρ(y) ] 1/2 [ ρ(v) ρ(y) ] 1/2 ≤ [ρ(u)ρ(v)] 1/2 e (C/2)(u+v)-Cy .
Thus,

I -≤ 2 u<v<y<0 h(u)h(v)[ρ(u)ρ(v)] 1/2 e (C/2)(u+v)-Cy dydudv = 2 u<v<0 h(u)h(v)[ρ(u)ρ(v)] 1/2 e (C/2)(u+v) ( 0 v e -Cy dy)dudv = 2 C u<v<0 [ρ(u)ρ(v)] 1/2 e (C/2)u (e (C/2)v -e -(C/2)v )h(u)h(v)dudv ≤ 2 C u<0 [ρ(u)] 1/2 |h(u)|e (C/2)u du v<0 |(e (C/2)v -e -(C/2)v )h(v)|[ρ(v)] 1/2 dv < +∞
by Proposition 1 and the assumption on h.

We prove the second item. We have:

J -= 2 x<y<0,u<v<y |x| ρ(x)ρ(u)ρ(v) ρ 2 (y) h(u)h(v)dxdydudv.
We can write:

ρ(x)ρ(u)ρ(v) ρ 2 (y) = [ρ(x)] 1/3 [ ρ(x) ρ(y) ] 2/3 [ ρ(u) ρ(y) ] 2/3 [ ρ(v) ρ(y) ] 2/3 [ρ(u)ρ(v)] 1/3 .
On the integration set, we have:

ρ(x)ρ(u)ρ(v) ρ 2 (y) ≤ [ρ(x)ρ(u)ρ(v)] 1/3 e (2/
3)C(x-y) e (2/3)C(u-y) e (2/3)C(v-y) .

Thus,

J -≤ x<0,u<v |x||h(u)h(v)|[ρ(x)ρ(u)ρ(v)] 1/3 e (2C/3)(x+u+v) ( +∞ sup(x,v) e -2Cy dy)dxdudv < +∞,
by Proposition 1 and the assumption on h.

We treat [START_REF] Baladron | Mean-eld description and propagation of chaos in networks of Hodgkin-huxley and FitzHugh-Nagumo neurons[END_REF]. By Proposition 1, and Inequality (8), we get straightforwardly that For (4), we have:

+∞ 0 xρ(x)dx x 0 du 1 ρ 2 (u) +∞ u ρ 2 (y)dy = x>0,0<u<x,0<u<y xρ(x) ρ 2 (y) ρ 2 (u) dxdudy.
We can write:

ρ(x) ρ 2 (y) ρ 2 (u) = [ρ(x)] 1/3 [ ρ(x) ρ(u) ] 2/3 [ ρ(y) ρ(u) ] 4/3 [ρ(y)] 2/3
Therefore, using (8) again,

x>0,0<u<x,0<u<y xρ(x) ρ 2 (y) ρ 2 (u) dxdudy ≤ x>0,0<u<x,0<u<y x[ρ(x)] 1/3 [ρ(y)] 2/3 exp [(2/3)C(x -u)] exp [(4/3)C(y -u)]dxdudy x>0,y>0 x[ρ(x)] 1/3 [ρ(y)] 2/3 exp [(2/3)Cx + (4/3)Cy] x∧y 0 exp (-2u)dudxdy < +∞.
For ( 5),

+∞ 0 xρ(x)dx x 0 du ρ 2 (u) +∞ u ρ(y)dy 2 = 1 0<u<x 1 u<y,u<y xρ(x)ρ(y)ρ(y ) ρ 2 (u) dxdydy du.
As above, on the integration set,

ρ(x) ρ(y)ρ(y ) ρ 2 (u) = [ρ(x)] 1/2 [ ρ(x) ρ(u) ] 1/2 [ ρ 3/4 (y)ρ 3/4 (y ) ρ 3/2 (u) ][ρ(y)ρ(y )] 1/4
≤ [ρ(x)] 1/2 ρ(y)] 1/4 ρ(y )] 1/4 e (C/2)x+(3C/4)(y+y ) e -(5C/4)u , and conclude that the rst integral of ( 5) is nite. We proceed analogously for the second one. 11) and ( 12), we obtain, taking into account that we want lim |y|→∞ g (y)ρ(y) = 0,

g (y) = 2 σ 2 ρ(y) y -∞ h c (u)ρ(u)du = - 2 σ 2 ρ(y) +∞ y h c (u)ρ(u)du,
where the equality above comes from the fact that +∞ -∞ h c (u)du = 0. We choose g(x) =

x 0 g (u)du which gives Lg = h c . We must now prove (32) I = [g (y)] 2 ρ(y)dy < +∞, and J = g 2 (x)ρ(x)dx < +∞.

We can split I into I = 4 σ 4 (I + + I -) with

I -= 0 -∞ dy ρ(y) y -∞ h c (u)ρ(u)du 2 , I + = +∞ 0 dy ρ(y) +∞ y h c (u)ρ(u)du 2 .
It is enough to look at the rst integral, the other one being analogous. We have, using

h c = h -hρ, I -≤ 2 0 -∞ dy ρ(y) y -∞ h(u)ρ(u)du 2 + 2 0 -∞ dy ρ(y) y -∞ ρ(u)du 2 h(u)ρ(u)du 2 .
These integrals are nite by Lemma 3 (with h(u) = 1 for the second one). Now, we also split J = J + + J -, with

J + = +∞ 0 [g(x)] 2 ρ(x)dx, J -= 0 -∞ [g(x)] 2 ρ(x)dx,
and only treat J -. We have

J -≤ 0 -∞ |x|ρ(x) 0 x [g (y)] 2 dydx = x<y<0 |x|[g (y)] 2 ρ(x)dxdy ≤ 8 σ 4 x<y<0 |x| dy ρ 2 (y) y -∞ h(u)ρ(u)du 2 ρ(x)dxdy + 8 σ 4 x<y<0 |x| dy ρ 2 (y) y -∞ ρ(u)du 2 h(z)ρ(z)dz 2 := J (1) -+ J (2) -× h(z)ρ(z)dz 2
These integrals are nite by Lemma 3. Thus, ( 32) is satised and Γ h = g belongs to D and satises LΓ h = h c . Then ( 14) is obtained thanks to the Ito formula, and ( 15) is deduced using that X(t) is stationary. 2 7.5. Proof of Proposition 4. Obviously, E 1 T T 0 ϕ k (X (3) (s))ds = ϕ k (x)ρ(x)dx. For the term 1 T T 0 θ (X (3) (s))ds, the expectation is computed in two steps by using the independence of X (3) and X (4) .

E 1 T T 0 θ (X (3) (s))ds = E(θ (x))ρ(x)dx = θ ρ(x)ρ(x)dx.

Gathering both terms yields that

E( Z m,p ) = Z m,p . Concerning V m,p , it is clear that E 1 T T 0 ϕ j (X (1) (s))ϕ k (X (1) (s))ds = ϕ j (x)ϕ k (x)ρ(x)dx
and with the same two-step computation as previously, we have

E 1 T T 0 θ (X (1) (s)) ϕ k (X (1) (s))ds = E( θ (x))ϕ k (x)ρ(x)dx = θ ρ(x)ϕ k (x)dx.

Now we prove that

E 1 T T 0 θ (X (1) (s)) θr (X (1) (s))ds = θ ρ(x)θ r ρ(x)ρ(x)dx + O( 1
T
).

First, we split:

1 T T 0 θ (X (1) (s)) θr (X (1) (s))ds-θ ρ(x)θ r ρ(x)ρ(x)dx = A 1 ( , r)+A 2 ( , r)+A 3 ( , r)+A 4 ( , r),
where

A 1 ( , r) = T -1 T 0
( θ (X (1) (s)) -θ ρ(X (1) (s)))( θr (X (1) (s)) -θ r ρ(X (1) (s))ds.

A 2 ( , r) = T -1 T 0 ( θ (X (1) (s)) -θ ρ(X (1) (s)))θ r ρ(X (1) (s)ds, A 3 ( , r) = A 2 (r, ), A 4 ( , r) = T -1 T 0 θ ρ(X (1) (s))θ r ρ(X (1) (s))ds -θ ρ(x)θ r ρ(x)ρ(x)dx,
The term A 4 ( , r) is centered as E(θ ρ(X (1) (s))θ r ρ(X (1)

(s)) = θ ρ(x)θ r ρ(x)ρ(x)dx.
By the two-step computation already used, the terms A 2 ( , r), A 3 (r, ) are centered too by the independence of X (1) and X (2) , as

E X (1) (s) ( θ (X (1) (s)) -θ ρ(X (1) (s)) = 0.
For the term A 1 ( , r), we have:

EA 1 ( , r) = ρ(u)duE ( θ (u) -θ ρ(u))( θr (u) -θ r ρ(u)) ;
We intend to use Proposition 3. Set θ ,u (x) = θ (u -x) and note that θ c ,u (x) = θ ,u (x)θ ,u (x)ρ(x)dx = θ (u -x) -θ ρ(u). Now with Γ θ ,u such that LΓ θ ,u = θ c ,u , it holds that

1 T T 0 θ c ,u (X (2) (t))dt = θ (u) -θ ρ(u) = σ T T 0 Γ θ ,u (X (2) (s))dW 2 (s) + 1 T Γ θ ,u (X (2) (0)) -Γ θ ,u (X (2) (T )) . (33)
Consequently,

E ( θ (u) -θ ρ(u))( θr (u) -θ r ρ(u)) = T (1)
,r (u) + T

,r (u) + T

,r (u) + T

,r (u)

with

T (1) ,r (u) = σ 2 T 2 E T 0 Γ θ ,u (X (2) (s))Γ θr,u (X (2) (s))ds = σ 2 T Γ θ ,u (y)Γ θr,u (y)ρ(y)dy T (2) ,r (u) = 1 T 2 E Γ θ ,u (X (2) (0)) -Γ θ ,u (X (2) (T )) Γ θr,u (X (2) (0)) -Γ θr,u (X (2) (T )) T (3) ,r (u) = σ T 2 E T 0 Γ θ ,u (X (2) (s)dW 2 (s) Γ θr,u (X (2) (0)) -Γ θr,u (X (2) (T )) T (4) ,r (u) = T (3) r, (u)
By the Cauchy-Schwarz inequality, we get:

T (2) ,r (u) ≤ 2 T 2 Γ 2 θ ,u (y)ρ(y)dy Γ 2 θr,u (y)ρ(y)dy 1/2 , T (3)
,r (u) ≤ 2σ T 3/2 (Γ θ ,u (y)) 2 ρ(y)dy Γ 2 θr,u (y)ρ(y)dy

1/2
.

Next, we must integrate all these terms w.r.t. ρ(u)du.

Using the Cauchy-Schwarz inequality for the second and the third inequalities leads to:

ρ(u)T (1)
,r (u)du =

σ 2 T Γ θ ,u (y)Γ θr,u (y)ρ(y)ρ(u)dudy, ρ(u)T (2)
,r (u)du ≤

σ 4 T 2 Γ 2 θ ,u (y)ρ(y)ρ(u)dudy Γ 2 θr,u (y)ρ(y)ρ(u)dudy 1/2 ρ(u)T (3) ,r (u)du ≤ 2σ T 3/2 [Γ θ ,u (y)] 2 ρ(y)ρ(u)dudy Γ 2 θr,u (y)ρ(y)ρ(u)dudy 1/2 .
Computations analogous to the ones of Lemma 3 allow to prove that the integrals above are nite. We only treat the rst one. Recall that

Γ θ ,u (y) = 2 σ 2 ρ(y) y -∞ θ c ,u (z)ρ(z)dz = - 2 σ 2 ρ(y) +∞ y θ c ,u (z)ρ(z)dz.
Therefore, we split T

,r (u) = which is nite by Lemma 3.

Gathering the four terms, we get that

E 1 T ( θ (X (1) (s) -θ ρ(X (1) (s))( θr (X (1) (s) -θ r ρ(X (1) (s))ds = σ 2 T Γ θ ,u (y)Γ θr,u (y)ρ(y)dy ρ(u)du + 1 T 3/2 O(1).
Finally,

E 1 T T 0 θ (X (1) (s)) θr (X (1) (s))ds = θ ρ(x) θ r ρ(x) ρ(x)dx + σ 2 T Γ θ ,u (y)Γ θr,u (y)ρ(y)dy ρ(u)du + 1 T 3/2 O(1).
This states a more precise result and ends the proof. 2 7.6. Proof of Theorem 1. We have

(b, φ) m,p -(b, φ) 2 V = (b, φ) m,p -(b, φ) 2 V 1 ∆m,p + (b, φ) 2 V 1 ∆ c m,p = (b, φ) m,p -(b, φ) m,p 2 V + (b, φ) m,p -(b, φ) 2 V 1 ∆m,p + (b, φ) 2 V 1 ∆ c m,p
where we use Pythagoras theorem: (b, φ) m,p -(b, φ) is orthogonal to S m × Σ p w.r.t the V -scalar product and thus V -orthogonal to (b, φ) m,p -(b, φ) m,p . Therefore

(b, φ) m,p -(b, φ) 2 V ≤ (b, φ) m,p -(b, φ) m,p 2 V 1 ∆m,p + (b, φ) m,p -(b, φ) 2 V + (b, φ) 2 V 1 ∆ c m,p
The middle rhs term is the squared bias term. Let us study the two others. We write, on

∆ m,p , (b, φ) m,p -(b, φ) m,p 2 V = (b, c) m,p -(b, c) m,p 2
Vm,p = V -1 m,p Z m,p -V -1 m,p Z m,p 2
Vm,p

= ( V -1 m,p -V -1 m,p )( Z m,p -Z m,p ) + V -1 m,p ( Z m,p -Z m,p ) + ( V -1 m,p -V -1 m,p )Z m,p 2 Vm,p ≤ 3 ( V -1 m,p -V -1 m,p )( Z m,p -Z m,p ) 2 Vm,p + 3 V -1 m,p ( Z m,p -Z m,p ) 2 Vm,p + 3 ( V -1 m,p -V -1 m,p )Z m,p 2 Vm,p
We set

T 1 := ( V -1 m,p -V -1 m,p )( Z m,p -Z m,p ) 2 Vm,p 1 ∆m,p T 2 := V -1 m,p ( Z m,p -Z m,p ) 2 Vm,p 1 ∆m,p T 3 := ( V -1 m,p -V -1 m,p )Z m,p 2
Vm,p 1 ∆m,p T 4 := (b, φ) 2 V 1 ∆ c m,p
We have

E( (b, φ) m,p -(b, φ) 2 V ) ≤ (b, φ) m,p -(b, φ) 2 V + 3E(T 1 + T 2 + T 3 ) + E(T 4 ).
The proof of Theorem 1 is structured in a study of the terms, E(T i ) for i = 1, . . . , 4 and relies on the results stated in Propositions 5 , 6 and Lemma 2.

• Study of T 4 = (b, φ) 2 V 1 ∆ c m,
E(T 4 ) = E( (b, φ) 2 V 1 ∆ c m,p ) ≤ (4 + 1/4)s 2 0 (b, φ) 2 V E V m,p -V m,p 2 op .
The bound for E(T 4 ) follows then from Proposition 6,

E(T 4 ) s 2 0 (b, φ) 2 V L ϕ (m) T Proof of Lemma 4. Let (35)
Ω m,p = V -1/2 m,p V m,p V -1/2 m,p -Id m+p op < 1 2 .
We write that

P(∆ c m,p ) = P(∆ c m,p ∩ Ω m,p ) + P(∆ c m,p ∩ Ω c m,p ).
We have

P(Ω c m,p ) ≤ 4E V -1/2 m,p V m,p V -1/2 m,p -Id m+p 2 op ≤ 4 V -1 m,p 2 op E V m,p -V m,p 2 op .
Therefore, using [H5],

P(∆ c m,p ∩ Ω c m,p ) ≤ P(Ω c m,p ) ≤ 4s 2 0 E V m,p -V m,p 2 op .
Next, on ∆ c m,p , due to [H5], it holds that V -1 m,p -V -1 m,p op > s 0 and

P(∆ c m,p ∩ Ω m,p ) ≤ E 1 { V -1 m,p -V -1 m,p op>s0} 1 Ωm,p ≤ 1 s 2 0 E V -1 m,p -V -1 m,p 2 op 1 Ωm,p .
Using (47), we get

P(∆ c m,p ∩ Ω m,p ) ≤ s 2 0 4 E V m,p -V m,p 2 op .
Adding both bounds gives [START_REF] Giesecke | Inference for large nancial systems[END_REF] and concludes the proof of Lemma 4. 2

• Study of

T 3 = ( V -1 m,p -V -1 m,p )Z m,p 2
Vm,p 1 ∆m,p . Using [START_REF] Gärtner | On the Mc-Kean-Vlasov limit for interacting diusions[END_REF] and Lemma 2 yields

E(T 3 ) ≤ 2 ρ ∞ E ( V -1 m,p -V -1 m,p ) 2 op 1 ∆m,p Z m,p 2 R m,p ≤ 80.5 s 4 0 ρ ∞ E( V m,p -V m,p 2 op ) Z m,p 2 R m+p . (36)
Now, we can prove the following result:

Lemma 5. We have

Z m,p 2
R m+p ≤ 2 ρ ∞ b + φ ρ 2 ρ = 2 ρ ∞ (b, φ) 2 V .
Thus, by [START_REF] Herrmann | Non uniqueness of stationary measures for self-stabilizing diusions[END_REF], Proposition 6 and Lemma 5, we obtain,

E(T 3 ) s 4 0 L ϕ (m) T .
Proof of Lemma 5. First, we note that, using [START_REF] Comte | Nonparametric adaptive estimation for interacting particle systems[END_REF] we get

Z m,p 2
R m+p = m-1 k=0 ϕ k (x)[b(x) + φ ρ(x)]ρ(x)dx 2 + p-1 =0 θ ρ(x)[b(x) + φ ρ(x)]ρ(x)dx 2 .
Using again the relation

(u v) w = u (v -w) with v -(x) = v(-x) yields θ ρ(x)[(b(x) + φ ρ(x)]ρ(x)dx = θ (x)(b + φ ρ)ρ ρ -)(x)dx.
Therefore, by the projection argument and the Young inequality

( u v ≤ u 1 v ) as ρ 1 = 1, we get Z m,p 2
R m+p ≤ (b + φ ρ) 2 ρ 2 + [(b + φ ρ)ρ ρ -] 2 ≤ (b + φ ρ) 2 ρ 2 + ρ - 2 1 (b + φ ρ) 2 ρ 2 ≤ 2 ρ ∞ (b + φ ρ) 2 ρ = 2 ρ ∞ (b, φ) 2 V . 2
• Study of E(T 2 ). Using [H5] and ( 28) yields

E(T 2 ) = E V -1 m,p ( Z m,p -Z m,p )1 ∆m,p 2 Vm,p ≤ V -1 m,p op E Z m,p -Z m,p 2 R m+p 1 ∆m,p ≤ s 0 E Z m,p -Z m,p 2 R m+p .
Applying Proposition 5 yields that ( 37)

E(T 2 ) s 0 L ϕ (m) + L θ (p) T .
• Study of

T 1 = ( V -1 m,p -V -1 m,p )( Z m,p -Z m,p ) 2
Vm,p 1 ∆m,p . By (28),

T 1 ≤ 2 ρ ∞ ( V -1 m,p -V -1 m,p 2 op 1 ∆m,p ) Z m,p -Z m,p 2
R m+p .
As the two factors of the right-hand-side are independent

E(T 1 ) ≤ 2 ρ ∞ E V -1 m,p -V -1 m,p 2 op 1 ∆m,p E Z m,p -Z m,p 2 R m+p .
Using Lemma 2 we get

E V -1 m,p -V -1 m,p 2 op 1 ∆m,p ≤ (40 + 1/4)s 4 0 E V m,p -V m,p 2 op 1 ∆m,p , E(T 1 ) ≤ 2(40 + 1/4)s 4 0 ρ ∞ E V m,p -V m,p 2 op 1 ∆m,p E Z m,p -Z m,p 2 R m+p .
By applying Proposition 5 and 6 and using L ϕ (m) ≤ T , we obtain

E(T 1 ) s 4 0 L ϕ (m) + L θ (p) T .
Gathering the ve bounds ends the proof of Theorem 1.2 7.7. Proof of Proposition 5. We have

E[ Z m,p -Z m,p 2 R m+p ] = σ 4 4 E m-1 k=0 1 T T 0 ϕ k (X (3) (s))ds -ϕ k (x)ρ(x)dx 2 + p-1 =0 1 T T 0 θ k (X (3) (s))ds -θ k ρ(x) ρ(x)dx 2 . (38)
We prove two lemmas for each term of the r.h.s. Lemma 6. Assume that ϕ k is a dierentiable and bounded basis. Let g k := Γ ϕ k be dened by

Lg k = [ϕ k -ϕ k (x)ρ(x)dx] (see Proposition 3). Then, ( 39
) E 1 T T 0 ϕ k (X (3) (s))ds -ϕ k (x)ρ(x)dx 2 ≤ 2σ 2 T [g k (x)] 2 ρ(x)dx + 8 T 2 g 2 k (x)ρ(x)dx.
Moreover,

m-1 k=0 E 1 T T 0 ϕ k (X (3) (s))ds -ϕ k (x)ρ(x)dx 2 L ϕ (m) T . Lemma 7. Recall θ (x) = 1 T T 0 θ (x -X (4) (t))dt.
Let G x, be dened by

LG x, (y) = θ (x -y) -θ ρ(x) and h be dened by Lh = θ ρ -θ ρ(x)ρ(x)dx. Then,

E 1 T T 0 θ (X (3) (s))ds -θ ρ(x)ρ(x)dx 2 ≤ 4σ 2 T G x,k (y) 2 ρ(x)ρ(y)dxdy + 8 T 2 (G x, (y)) 2 ρ(x)ρ(y)dxdy + 4σ 2 T [h (x)] 2 ρ(x)dx + 8 T 2 [h (x)] 2 ρ(x)dx.
Moreover,

p-1 =0 E 1 T T 0 θ (X (3) (s))ds -θ ρ(x)ρ(x)dx 2 L θ (p) T .
From the two Lemmas and formula [START_REF] Hutton | Interchanging the order of dierentiation and stochastic integration[END_REF], we get the result of Proposition 5. 2

Proof of Lemma 6. By denition of L and g k := Γ ϕ k and Proposition 3, we have

1 T T 0 ϕ k (X (3) (s))ds -ϕ k (x)ρ(x)dx) = σ T T 0 g k (X (3) (s))dW 3 (s) + 1 T [g k (X (3) (0)) -g k (X (3) (T ))],
and consequently by [START_REF] Comte | Penalized nonparametric mean square estimation of the coecients of diusion processes[END_REF],

E 1 T T 0 ϕ k (X (3) (s))ds -ϕ k ρ 2 ≤ 2σ 2 T [g k (x)] 2 ρ(x)dx + 8 T 2 g 2 k (x)ρ(x)dx.
This yields inequality [START_REF] Indritz | An inequality for Hermite polynomials[END_REF]. Now we have

g k (x) = 2 σ 2 ρ(x) x -∞ ϕ k (y) -ϕ k (z)ρ(z)dz ρ(y)dy = - 2 σ 2 ρ(x) +∞ x ϕ k (y) -ϕ k (z)ρ(z)dz ρ(y)dy.
We split the term The second term is easier and analogous.

[g k (x)] 2 ρ(x)dx = +∞ 0 [g k (x)] 2 ρ(x)dx + 0 -∞ [g k (x)]
The term I -is treated analogously using the second formula for g k . Therefore,

σ 2 T m-1 k=0 [g k (x)] 2 ρ(x)dx L ϕ (m) T .
Now, with g k (x) = x 0 g k (u)du, we look at g 2 k (x)ρ(x)dx and split it again into +∞ 0 . . . + 0 -∞ . . .. We use that for x ≥ 0,

[g k (x)] 2 = x 0 g k (u)du 2 ≤ x x 0 [g k (u)] 2 du.
We can write +∞ 0

g 2 k (x)ρ(x)dx ≤ 2 σ 4 +∞ 0 xρ(x)dx x 0 du 1 ρ(u) +∞ u ϕ k (y)ρ(y)dy 2 + 8 σ 4 +∞ 0 xρ(x)dx x 0 du 1 ρ(u) +∞ u ρ(y)dy ϕ k (z)ρ(z)dz) 2 .
Then, with the same integration by part for ϕ k ρ as previously, we get Therefore, by Lemma 3, all integrals being nite, we obtain the result. 2 Proof of Lemma 7. We have θ (X (3) (s))-θ ρ(x)ρ(x)dx = θ (X (3) (s)-θ ρ(X (3) (s))+θ ρ(X (3) (s))-θ ρ(x)ρ(x)dx.

Then, for all x,

θ (x) -θ ρ(x) = 1 T T 0 θ (x -X (4) (t))dt -θ ρ(x) = σ T T 0 G x, (X (4) (t))dW 4 (t) + 1 T [G x, (X (4) (0)) -G x, (X (4) (T ))].
And,

1 T T 0 θ ρ(X (3) (s))ds-θ ρ(x)ρ(x)dx = σ T T 0 h (X (3) (s)))dW 3 (s)+ 1 T [h (X (3) (0)))-h (X (3) (T )))].
Thus,

1 T T 0 θ (X (3) (s))ds -θ ρ(x)ρ(x)dx = 1 T T 0 ds σ T T 0 G X (3) (s)), (X (4) (t))dW 3 (t) + 1 T T 0 ds 1 T [G X (3) (s)), (X (4) (0)) -G X (3) (s)), (X (4) (T ))] + σ T T 0 h (X (3) (s)))dW 3 (s) + 1 T [h (X (3) (0))) -h (X (3) (T )))]
It follows that

E 1 T T 0 θ (X (3) (s))ds -θ ρ(x)ρ(x)dx 2 1 T T 0 ds σ 2 T 2 EE X (3) (s) T 0 dt[G X (3) (s)), (X (4) (t))] 2 + 1 T T 0 ds 1 T 2 EE X (3) (s) [G X (3) (s)), (X (4) (0)) -G X (3) (s)), (X (4) (T ))] 2 + σ 2 T 2 T 0 E[h (X (3) (s))] 2 ds + 1 T 2 E[h (X (3) (0))) -h (X (3) (T )))] 2 .
This yields

E 1 T T 0 θ (X (3) (s))ds -θ ρ(x)ρ(x)dx 2 1 T [G x, (y)] 2 ρ(x)ρ(y)dxdy + 1 T 2 [G x, (y)] 2 ρ(x)ρ(y)dxdy + 1 T [h (x)] 2 ρ(x)dx + 1 T 2 [h (x)] 2 ρ(x)dx
Next, we have to study p-1 =0 for each of the above terms following the lines of Lemma 6. By formula [START_REF] Chen | Maximum likelihood estimation of potential energy in interacting particle systems from single-trajectory data[END_REF], we have for y ≥ 0,

G x, (y) = 2 σ 2 ρ(y) +∞ y θ (x -z) -θ ρ(x) ρ(z)dz.
By integration by part, we obtain, for y ≥ 0 Following the proof of Lemma 3, we prove that all integrals are nite. The same holds for the case y ≤ 0. The order of the term is L θ (p).

G x, (y) = 2 σ 2 ρ(y) θ (x -y)ρ(y) + +∞ y θ (x -z)ρ (z)dz -θ ρ (x) +∞ y ρ(z)dz . Thus p-1 =0 [G x, (y)] 2 = 4 σ 4 p-1 =0 θ (x -y) + 1 ρ(y) +∞ y θ (x -z)ρ (z)dz - θ ρ (x) ρ(y) +∞ y ρ(z)dz 2 ≤ 12 σ 4 L θ (p) + 1 ρ 2 (y) +∞ y [ρ (z)] 2 dz + (ρ ) 2 +∞ y ρ(z)dz ρ 2 (y)
The terms [G x, (y)] 2 ρ(x)ρ(y)dxdy are treated in a similar way with an additional integration due to the bound G 2

x, (y) ≤ |x| [0,x] [G x, (u)] 2 du. This yields the same order L ϕ (p) for the sum over . Now, for the terms in h , we note that θ ρ -θ ρρ = θ ρ -θ ρ ρ.

By Formula (13), we get, for x ≥ 0,

h (x) = 2 σ 2 ρ(x) +∞ x θ ρ (u) -θ ρ (z)ρ(z)dz ρ(u)du so that p-1 =0 [h (x)] 2 ≤ 4 σ 4 ρ 2 (x) +∞ x θ ρ (u) -θ ρ (z)ρ(z)dz ρ(u)du 2 . Now, note that for x ≥ 0 +∞ x θ ρ (u)ρ(u)du = θ (v) 1 u≥x ρ (u -v)ρ(u)du dv which yields p-1 =0 +∞ x θ ρ (u)ρ(u)du 2 ≤ 1 u≥x ρ (u -v)ρ(u)du 2 dv.
This implies that

p-1 k=0 ∞ 0 [h (x)] 2 ρ(x)dx ≤ 8 σ 4 +∞ 0 1 u≥x ρ (u -v)ρ(u)du 2 dv dx ρ(x) + ρ (u -v)ρ(u)du 2 dv +∞ 0 +∞ x ρ(u)du 2 dx ρ(x) .
As a consequence, this term is bounded by C(ρ)/σ 4 where C(ρ) is a constant depending on ρ only.

We proceed as done several times previously with p-1 =0 ∞ 0 h 2 k (x)ρ(x)dx, which yields the same order. 2 7.8. Proof of Proposition 6. Recall that the matrices V m,p and V m,p are symmetric nonnegative. To obtain Proposition 6, we prove that

E Tr( V m,p -V m,p ) 2 L ϕ (m) T .

Let us recall that for a function

h, h(x) = 1 T T 0 h(x -X (2) (s))ds.
We have Tr( V m,p -V m,p ) 2 = T 11 + 2T 12 + T 22 where

T 11 = Tr( V 1 -V 1 ) 2 ) = m-1 j,k=0 1 T T 0 ϕ j (X (1) (s))ϕ k (X (1) (s))ds -ϕ j ϕ k ρ 2 T 12 = Tr( V 1,2 -V 1,2 )( V 1,2 -V 1,2 ) ⊥ = m-1 k=0 p-1 =0 1 T T 0 θ (X (1) (s)) ϕ k (X (1) (s))ds -θ ρϕ k ρ 2 T 22 = Tr( V 2 -V 2 ) 2 ) = p-1 ,r=0 1 T T 0 θ (X (1) (s)) θ r (X (1) (s))ds -θ ρ θ r ρ ρ 2 , Lemma 8. We have E(T 11 ) Lϕ(m) T + Lϕ(m) T 2
Proof of Lemma 8. Let g jk (x) = Γ ϕ j ϕ k where Γ h is dened in Proposition 3. We have:

E 1 T T 0 ϕ j (X (1) (s))ϕ k (X (1) (s))ds -ϕ j ϕ k ρ 2 ≤ 2σ 2 T [g jk (x)] 2 ρ(x)dx)+ 8 T 2 g 2 jk (x)ρ(x)dx.
For each integral of the rhs, we write = 0 -∞ + +∞ 0

and use the appropriate expression of g jk in each integral. Thus, We proceed analogously for the second term. This yields

4 σ 4 +∞ 0 [g jk (x)] 2 ρ(x)dx = +∞ 0 dx ρ(x) +∞ x [ϕ j (y)ϕ k (y) -ϕ j (z)ϕ k (zρ(z)dz]ρ(y)dy 2 ) ≤ 2 +∞ 0 dx ρ(x) +∞ x [ϕ j (y)ϕ k (y)ρ(y)dy 2 + +∞ x ρ(y)dy 2 ϕ j (z)ϕ k (z)ρ(z)dz]
0≤j,k≤m-1 +∞ 0 [g jk (x)] 2 ρ(x)dx ≤ c 1 L ϕ (m), with σ 4 4 c 1 = 2 × +∞ 0 dx ρ(x) +∞ x ρ 2 (y)dy + 2 +∞ 0 dx ρ(x) +∞ x ρ(y)dy 2 ρ 2 (x)dx.
The integral on (-∞, 0) is treated analogously. This yields 0≤j,k≤m-1

[g jk (x)] 2 ρ(x)dx ≤ C 1 L ϕ (m). Next, +∞ 0 ρ(x)g 2 jk (x)dx ≤ +∞ 0 ρ(x)x x 0 [g jk (y)] 2 dy = 0<y<x xρ(x)[g jk (y)] 2 dxdy.
We can compute as above

ρ 2 (y) m-1 j=0 m-1 k=0 [g jk (y)] 2 ≤ 2L ϕ (m) +∞ y ρ 2 (z)dz + 2 +∞ y ρ(z)dz 2 ρ 2 (u)du.
Thus,

σ 4 4 0≤j,k≤m-1 +∞ 0 ρ(x)g 2 jk (x)dx ≤ 2L ϕ (m) 0<y<x,y<<z x ρ(x)ρ 2 (z) ρ 2 (y) dxdydz + 4L ϕ (m) 0<y<x,y<z<z x ρ(x)ρ(z)ρ(z ) ρ 2 (y) dxdydzdz ρ 2 (u)du,
where all the integrals are nite. Analogously, 0

-∞ ρ(x)g 2 jk (x)dx L ϕ (m). This concludes the proof. 2 Lemma 9. ET 12 Lϕ(m) T + Lϕ(m) T 2 + 1 T + 1 T 2 . Proof of Lemma 9. 1 T T 0 θ (X (1) (s)) ϕ k (X (1) (s))ds -θ ρ(x)ϕ k (x)ρ(x)dx = T 1 (k, ) + T 2 (k, ) T 1 (k, ) = 1 T T 0 [θ (X (1) (s)) -θ ρ(X (1) (s))] ϕ k (X (1) (s))ds T 2 (k, ) = 1 T T 0 θ ρ(X (1) (s)) ϕ k (X (1) (s))ds -θ ρ(x)ϕ k (x)ρ(x)dx.
The second term is easier. Set h k = Γ ϕ k θ ρ so that

T 2 (k, ) = σ T T 0 h k (X (1) (s))dW 1 (s)) + 1 T [h k (X (1) (0)) -h k (X (1) (T ))].
Therefore,

E[T 2 (k, )] 2 1 T [h k (x)] 2 ρ(x)dx + 1 T 2 h 2 k (x)ρ(x)dx.
We must compute, for x > 0, (and then for x < 0)

σ 2 4 ρ 2 (x) p-1 =0 m-1 k=0 [h k (x)] 2 ≤ 2 p-1 =0 m-1 k=0 +∞ x θ ρ(y)ϕ k (y)ρ(y)dy 2 +2 p-1 =0 m-1 k=0 +∞ x ρ(y)dy 2 θ ρ(y)ϕ k (y)ρ(y)dy 2 ≤ 2 p-1 =0 +∞ x [θ ρ(y)] 2 ρ 2 (y)dy + 2 [θ ρ(y)] 2 ρ 2 (y)dy( +∞ x ρ(y)dy) 2 . Now, p-1 =0 +∞ x [θ ρ(y)] 2 ρ 2 (y)dy = p-1 =0 +∞ x θ (z)ρ(y -z)dz 2 ρ 2 (y)dy ≤ +∞ x ρ 2 (y -z)dzρ 2 (y)dy = +∞ x ρ 2 (y)dy ρ 2 (u)du. And p-1 =0 [θ ρ(y)] 2 ρ 2 (y)dy ≤ ρ 2 (y)dy 2 .
Thus,

p-1 =0 m-1 k=0 +∞ 0 [h k (x)] 2 ρ(x)dx +∞ 0 dx ρ(x) +∞ x ρ 2 (y)dy ρ 2 (u)du + ρ 2 (y)dy 2 +∞ 0 dx ρ(x) +∞ x ρ(y)dy 2 .
We proceed analogously for +∞ 0 .

The other term

[h jk (x)] 2 ρ(x)dx is treated analogously. This means that E( k, T 2 2 (k, )) 1 T + 1 T 2 .
Let us now study T 1 (k, ). This term contains the two trajectories X (1) and X (2) . For each x, set θ (x -z) = θ ,x (z) and G x, = Γ θ ,x so that

θ (x) -θ ρ(x) = σ T T 0 G x, (X (2) (t))dW 2 (t) + 1 T [G x, (X (2) (0)) -G x, (X (2) (T ))].
Recall that

G x, (y) = 2 σ 2 ρ(y) +∞ y [θ (x -z) -θ ρ(x)]ρ(z)dz = - 2 σ 2 ρ(y) y -∞ [θ (x -z) -θ ρ(x)]ρ(z)dz Thus, T 1 (k, ) = 1 T T 0 ϕ k (X (1) (s))ds σ T T 0 G X (1) (s)), (X (2) (t))dW 2 (t) + 1 T T 0 ϕ k (X (1) (s))ds 1 T [G X (1) (s)), (X (2) (0)) -G X (1) (s)), (X (2) (T ))] := S 1 (k, ) + S 2 (k, ).
We have:

S 1 (k, ) = 1 T T 0 1 T T 0 ϕ k (X (1) (s))G X (1) (s)), (X (2) (t))ds dW 2 (t)
Using the independance of the two trajectories, we get

ES 2 1 (k, ) = 1 T 4 T 0 dtE T 0 ϕ k (X (1) (s))G X (1) (s)), (X (2) (t))ds 2 = 1 T ϕ 2 k (x)[G x, (y)] 2 ρ(x)ρ(y)dxdy. Moreover, ES 2 2 (k, ) ≤ 2 T 2 ϕ 2 k (x)[G x, (y
)] 2 ρ(x)ρ(y)dxdy. We have to compute k, of the two terms. We have

p-1 =0 m-1 k=0 ϕ 2 k (x)[G x, (y)] 2 ρ(x)ρ(y)dxdy ≤ L ϕ (m) p-1 =0 [G x, (y)] 2 ρ(x)ρ(y)dxdy L ϕ (m) p-1 =0 ρ(x)dx[ +∞ 0 [G x, (y)] 2 ρ(y)dy] + L ϕ (m) p-1 =0 ρ(x)dx[ 0 -∞ [G x, (y)] 2 ρ(y)dy]
We only treat the term containing +∞ 0 . We have:

p-1 =0 ρ(x)dx[ +∞ 0 [G x, (y)] 2 ρ(y)dy] p-1 =0 ρ(x)dx[ +∞ 0 dy ρ(y) +∞ y θ (x -z)ρ(z)dz 2 ] + p-1 =0 ρ(x)dx[ +∞ 0 dy ρ(y) θ (x -z))ρ(z)dz 2 +∞ y ρ(z)dz 2 p-1 =0 ρ(x)dx[ +∞ 0 dy ρ(y) θ (u)1 [y,+∞) (x -u)ρ(x -u)du 2 ] + p-1 =0 ρ(x)dx +∞ 0 dy ρ(y) +∞ y ρ(z)dz 2 θ (u)ρ(x -u)du 2 ρ(x)ρ(x -u) ρ(x -u) ρ(y) 1 x-u>y>0 )dxdydu + ρ(x)ρ 2 (x -u)dxdu +∞ 0 dy ρ(y) +∞ y ρ(z)dz 2
Using the bound on ρ and a previous study, we nd that the second term above is nite. For the rst term, we have

ρ(x)ρ(x -u) ρ(x -u) ρ(y) 1 x-u>y>0 )dxdydu ≤ ρ(x)ρ(x -u)e C(x-u-y)) 1 x-u>y>0 )dxdydu = ρ(x)ρ(x -u)e C(x-u) x-u 0 e (-Cy) dy1 x-u>0 )dxdydu = ρ(x)ρ(x -u)e C(x-u) 1 x-u>0 )dxdu -ρ(x)ρ(x -u)e -C(x-u) 1 x-u>0 )dxdu = due -Cu +∞ u ρ(x)ρ(x -u)e Cx dx -due Cu +∞ u ρ(x)ρ(x -u)e -Cx dx
The two terms of the sum are analogous. We look at the rst one. For this, we compute, using Proposition 1 (where b is given):

ρ(x)ρ(x -u) exp - K + λ σ 2 [(x + b) 2 + (x -u + b) 2 ] exp - K + λ σ 2 (u 2 /2) exp - 2(K + λ) σ 2 [(x + b -(u/2))) 2 ]
Thus,

+∞ u ρ(x)ρ(x -u)e Cx exp - K + λ σ 2 (u 2 /2) e Cx exp - 2(K + λ) σ 2 [(x + b -(u/2))) 2 ] dx exp - K + λ σ 2 (u 2 /2) e C((u/2)-b) e Cz exp - 2(K + λ)z 2 σ 2 dz exp - K + λ σ 2 (u 2 /2) e C((u/2)-b) .
Finally,

due -Cu +∞ u ρ(x)ρ(x -u)e Cx dx du exp - K + λ σ 2 (u 2 /2) e -C((u/2)) < +∞.
This means that 1) (s))θ r ρ(X (1) (s))ds -θ ρ θ r ρ ρ.

k,l ES 2 1 (k, ) L ϕ (m) T . 0 θ ρ(X (
We only treat the term A 1 ( , r) as the other ones are easier and analogous terms have been already treated previously.

Using that

[(1/T ) T 0 Φ(u)du] 2 ≤ (1/T ) T 0 Φ 2 (u)du, we have E r, A 2 1 ( , r) ≤ E    r, (θ (u) -θ ρ(u))(θ r (u) -θ r ρ(u)) 2 ρ(u)du    ≤ E p-1 =0 (θ (u) -θ ρ(u)) 2 p-1 r=0 (θ r (u) -θ r ρ(u)) 2 ρ(u)du = E    p-1 =0 (θ (u) -θ ρ(u)) 2 2 ρ(u)du    := A 1 (40)
We rst work on the sums and write each as e.g. p-1 =0 ϕ , Ξ2 for some complicated and random Ξ, in order to upper bound the term by Ξ 2 .

Recall that θ (u) -θ ρ(u) = σ T T 0 Γ θ ,u (X (2) (s))dW 2 (s) + 1 T (Γ θ ,u (X (2) (0)) -Γ θ ,u (X (2) (T ))),

Therefore, using the above equality, we get: ρ(X (2) (s)) dW 2 (s).

More simply, θ ρ(u)

T 0 1 (X (2) (s)<0)
ρ(X (2) (s))

X (2) (s) -∞ ρ(v)dv dW 2 (s) 2 = ρ 2 (z)dz T 0 1 (X (2) (s)<0)
ρ(X (2) (s))

X (2) (s) -∞ ρ(v)dv dW 2 (s) 2
We proceed in the same way for the term containing 1 (X (2) (s)<0) using the second expression of Γ θ ,u (y). So that by gathering all terms, we have obtained We only look at the rst term:

ρ(u)B 2 11 (u) du = ρ(u)ρ 2 (v)ρ 2 (v )1 v<0,v <0 M 2 T (v)M 2 T (v )dudvdv = ρ 2 (v)ρ 2 (v )1 v<0,v <0 M 2 T (v)M 2 T (v )dvdv .
We have to take the expectation of the above term. Hence, we deal with E(M 2 T (v)M 2 T (v ). Using the Cauchy-Schwarz and the Burkholder-Davis-Gundy inequalities yields

E(M 2 T (v)M 2 T (v )) ≤ (E(M 4 T (v))E(M 4 T (v ))) 1/2 (E( M (v) 2 T )E( M (v ) 2 T )) 1/2 .
We have (see [START_REF] Mckean | A class of Markov processes associated with nonlinear parabolic equation[END_REF]):

E( M (v) 2 T ) = E T 0 1 (v<X (2) (s)<0) ρ 2 (X (2) (s) ds 2 ≤ T T 0 E 1 (v<X (2) (s)<0)
ρ 4 (X (2) (s) ds = T 2 1 (v<x<0) ρ 3 (x) dx Thus, we obtain:

ρ(u)B 2 11 (u) du T 2 ρ 2 (v)ρ 2 (v )1 v<0,v <0 1 (v<x<0) ρ 3 (x) 1 (v <x <0) ρ 3 (x ) dxdx 1/2 dvdv = T 2 ρ 2 (v)1 v<0 1 (v<x<0) ρ 3 (x) dx 1/2 dv 2 .
We split ρ 2 (v) = ρ 1/2 (v)ρ 3/2 (v) and get: Dealing analogously with the other terms, we nally nd (see [START_REF] Kasonga | Maximum likelihood theory for large interacting systems[END_REF] and ( 43)):

EA 1 1 T 2 .
The proof of Lemma 10 is complete. 2 7.9. Proof of Lemma 2. We split the expectation along

Ω m,p = V -1/2 m,p V m,p V -1/2 m,p -Id m+p op < 1 2 .
(see [START_REF] Herrmann | Large deviations and a Kramers'type low for selfstabilizing diusions[END_REF]) and its complement. On the complement, we have, by the denition of ∆ m,p ((26)), Thus, gathering the last two inequalities, we get, using [H5], For the other term, we have with [H5], We have A -1 B op = B op < 1/2 on Ω m,p , so that

E V -1 m,p -V -1 m,p
E V -1 m,p -V -1 m,p 2 op 1 ∆m,p 1 Ωm,p ≤ s 2 0 E V 1/2 m,p V -1 m,p V 1/2 m,p -Id m+p
A = A + B = V 1/2 m,p V -1 m,p V 1/2
m,p and Therefore we obtain (47) By [START_REF] Malrieu | Convergence to equilibrium for granular media equations and their Euler schemes[END_REF] and (47), we get Lemma 2. Both terms are handled similarly, so we consider the rst one only. where C(ρ) is a nite constant coming from Lemma 3.

A -1 -A -1 op = V 1/2 m,p V -1 m,p V 1/2 m,p -Id m+p op ≤ B op A -1 2 op 1 -A -1 B op = V -1/2 m,p V m,p V -1/2 m,p -Id m+p op 1 -V -1/2 m,p V -1 m,p V
E V -1 m,p -V -1 m,p 2
For T 2 , we write T 2 = m-1 j=0 ( 0 -∞ g 2 j ρ + +∞ 0 g 2 j ρ) and consider only the last sum, both terms being again similar. For x ≥ 0, we write that This term is nite by Lemma 3. So T 2 is bounded by say C (ρ)/σ 4 . This yields

E( ρ m -ρ m 2 ) ≤ C(ρ) σ 4 1 T (1 + 1 T ).
Plugging this in (50) ends the proof of Proposition 7. 2

3. 3 .

 3 Denition of the estimators. To dene a nonparametric estimator of b, φ, we estimate the functions b V m , φ V p , i.e. we build estimators b 0 , . . . , b m-1 and c 0 , . . . , c p-1 of the coecients b 0 , . . . , b m-1 and c 0 , . . . , c p-1 and set[START_REF] Forien | Household epidemic models and McKean-Vlasov Poisson driven SDEs[END_REF]

Figure 1 .

 1 Figure 1. Density estimation for Examples 1 to 3. 20 estimated curves in dotted green, rst line T = 200 (D = 70), last line T = 2000 (D = 225). In all cases, the true is in bold red.

Figure 2 .Figure 3 .

 23 Figure 2. Example 2, T = 2000. 20 estimated b (left) and φ right, all selected by penalization. The true functions are in bold red.

Figure 4 .

 4 Figure 4. Example 5, T = 200. First line: the 6x3 proposals for b (left) and the 3x6 proposals for φ (right). Second line: b 2 (right) and φ 1 (right) selected by penalization. In all cases, the true is in bold red.

2 7. 4 .

 24 Proof of Proposition 3. Recall that h having exponential growth, belongs to L 2 (ρ) as noted before the proposition. Let us nd Γ h = g such that g ∈ D and Lg = h c . Using (

  As θ is bounded, θ c ,u (.) is bounded too, so y<0 . . . y<0 dy ρ(y) z<y,z <y ρ(z)ρ(z )dzdz ∝

2

 2

2 2 ≤ 2 ≤

 222 2 ρ(x)dx := 4 σ 4 (I + + I -) y) -ϕ k (z)ρ(z)dz ρ(y)dy )dy ϕ k (z)ρ(z)dz) Now we write by integration by part, +∞ x ϕ k (y)ρ(y)dy = -ϕ k (x)ρ(x) -+∞ x ϕ k (y)ρ (y)dy as ϕ j is bounded and lim x→±∞ ρ(x) = 0.As a consequence for all x > 0, it holds y)ρ (y)1 [x,+∞] (y)dy 2L ϕ (m)ρ 2 (x) + 2 +∞ 0 [ρ (y)] 2 1 [x,+∞] (y)dy by using the usual projection argument. Thus 2L ϕ (m) + 2 y)] 2 dywhere, thanks to some adaptation of. Lemma 3, noting that ρ (y) = ρ(y) × k(y) where k(y) has polynomial growth, we can prove

σ 4 (

 4 [ρ (z)] 2 dz)

2 .

 2

  , (y)] 2 ρ(y)ρ(x)dxdy ≤ 12 σ 4 L θ (p) +

2 . 2 ≤

 22 (y)ϕ k (y)ρ(y)dy

T 0 1 (X ( 2 ) 1 (

 121 (s)<0) ρ(X (2) (s)) 1 (u-z≤X (2) (s)) θ (z)ρ(u -z)dz dW 2 (s) v≤X (2) (s)<0)

B 1 ( 0 Γ 2 B 0 1 2 = ρ 2 2 B 1 ( 2 B 14 (u) = ρ 2 ( 2 .

 1020222121422 θ ,u (X(2) (s))dW 2 (s) 11 (u) + B 12 (u) + B 13 (u) + B 14 (u) withB 11 (u) = ρ 2 (u -z)dz T (u-z≤X (2) (s)<0) ρ(X (2) (s)) dW 2 (s) (v)1 v<0 M 2 T (v)dv B 12 (u) = ρ 2 (z)dz 13 (u) = ρ 2 (u -z)dz T 0 0≤X (2) (s)<u-z) ρ(X (2) (s)) dW 2 (s)Therefore, we have ρ(u)B 2 1 (u)du 4 4 i=1 ρ(u)B 2 1i (u) du.

2 ≤ T 2 ρ 1 / 2

 212 (v)1 v<0 exp [-(3/2)Cv]

2 op ) 1 2 op ≤ 4 V - 1 m,p 2 op

 212412 ∆m,p 1 Ω c m,p ≤ 10s 2 0 P(Ω c m,p ) Moreover P(Ω c m,p ) ≤ 4E V -1/2 m,p V m,p V -1/2 m,p -Id m+p E V m,p -V m,p2op .

2 op 1 - 1 / 2 m 2 op 1 -

 211221 ∆m,p 1 Ωm,p .Let A = Id m+p and B = V -1/2 m,p V m,p V ,p -Id m+p . Then we use the following theorem:Theorem.[START_REF] Stewart | Matrix perturbation theory[END_REF])) Let A, B be (m × m) matrices. If A is invertible and A -1 B op < 1, then à := A + B is invertible and it holds Ã-1 -A -1 op ≤ B op A -1 A -1 B op .

-1/ 2 m,p -Id m+p op ≤ 1 2 V- 1 m,p op 2 V

 2212 -1/2 m,p V m,p V -1/2 m,p -Id m+p op ≤ V m,p -V m,p op .

op 1 ∆m,p 1 1 m,p 2 op 1 4 0 4 E

 112144 Ωm,p ≤ E V -1 m,p -V -Ωm,p ≤ s V m,p -V m,p 2 op .

2

 2

7. 10 .

 10 Proof of Proposition 7. Let g j := Γ ϕ j be such that Lg j = ϕ j -ϕ j ρ, as dened in Proposition 3, which is such that by (13)[START_REF] Molginer | A non-local model for a swarm[END_REF] g j (x) = 2 σ 2 ρ(x) x -∞ (ϕ j (y) -ϕ j ρ)ρ(y)dy = -2 σ 2 ρ(x) +∞ x (ϕ j (y) -ϕ j ρ)ρ(y)dy, X(s)) -ϕ j ρ ds = σ T 0g j (X(s))dW (s) + g j (X(0)) -g j (X(T )).

First

  by Pythogoras Theorem, the equality (50)E ρ m -ρ 2 = ρ -ρ m 2 + E ρ m -ρ mholds. The result will follow from a bound on the second rhs term,E ρ m -ρ m 2 = m-1 j=0 E[( a ja j ) 2 ]which is the variance term. By formula[START_REF] Nickl | Nonparametric statistical inference for drift vector elds of multi-dimensional diusions[END_REF], and Inequality[START_REF] Comte | Penalized nonparametric mean square estimation of the coecients of diusion processes[END_REF] in Proposition 3, we getE( ρ m -ρ m2k (y)ρ(y)dy -( ϕ j ρ) y)ρ(y)dy -( ϕ j ρ)

  p . For this term, we have

	Lemma 4. Under [H5], it holds
	(34)

(k, ). This brings no new diculty and we can prove that this sum is bounded by a constant × Lϕ(m) T 2 . 2

with θ ,u (z) = θ (u -z), and Γ θ ,u (y) = 2 σ 2 ρ(y)

Therefore

with

Then, (

We work on B 1 (u) and split it into in order to use the two forms of Γ ϕ ,u (y). We have:

Now, we use [START_REF] Kessler | Statistical methods for stochastic dierential equations[END_REF] for the rst stochastic integral:

By the stochastic Fubini theorem (see e.g. [START_REF] Hutton | Interchanging the order of dierentiation and stochastic integration[END_REF], Lemma 2.1), we interchange the ordinary and the stochastic integral. This yields: ρ(X (2) (s)) dW 2 (s) dz.


Anne Raingeard
email: a.raingeard@brgm.frpa.reninger@brgm.fr

Pierre-Alexandre Reninger

Aurélie Peyrefitte
email: a.peyrefitte@brgm.fr

Guillaume Martelet
email: g.martelet@brgm.fr

Bertrand Aunay
email: b.aunay@brgm.frarnauld.malard@isska.ch

Arnauld Malard

AEM imagery down to one kilometer depth: New constraints for geological and hydrogeological modeling in volcanic contexts

Keywords: Airborne electromagnetics, Time Domain ElectroMagnetics, aeromagnetic, volcanic island, geological modelling, hydrogeological modelling

INTRODUCTION

In order to bring fresh water to the northern part of La Réunion volcanic island (Indian Ocean), the County Council of La Réunion (CD974) requested an assessment of the geological and hydrogeological risks related to the drilling of a new water drainage gallery at great depths (300 to 1000 m) below the Plaine des Fougères. This new "GANOR" gallery will be connected to the existing deep GSAM gallery, which is covered by up to 1000 meters of volcanic materials. To assess the risks, BRGM has established a 3D geological model. This model then supported a 3D hydrogeological model in order to identify potential hydrogeological risks of drilling in the surrounding of highly permeable aquifers under high pressure. [START_REF] Peyrefitte | 71628-28: Aide à l'implantation d'une gallerie souterraine (GANOR): acquisition et valorisation de données géophysiques héliportées pour la caractérisation profonde géologique et hydrogéologique de la Plaine des Fougères[END_REF] present this work in a technical report. The local geological map exhibits very monotonous surface information, mainly with lava flows and tuffs PN3 and PN4 (Figure 1). The few existing soundings in the area never go deeper than a few dozen of meters.

DATA ACQUISITION

In 2014, an airborne electromagnetics (AEM) survey was conducted using the SkyTEM 304 system over the entire La Réunion island. This survey gave an image of the first 300 meters (line path in Figure 1). A first study creating a 3D geological model combining this AEM dataset with magnetic modelling was performed. To allow a preliminary geological/ hydrogeological interpretation, the lack of resistivity information at great depths, but also the lack of rock magnetization data, resulted in a preliminary understanding of the deep hydrogeology. This result was promising, but still highly hill-constrained.

Two surveys were then flown in 2021 over the Plaine des Fougères area, using SkyTEM 306HP and SkyTEM 312HP systems in order to extend the imagery, down to 1 000 m depth. These two new AEM datasets and the previous one are presented in Table 1. Figure 4 is an example of South-North resistivity profiles, illustrating the increase of investigation depth: SkyTEM 304 images resistivity up to 200-400 m depth, whereas SkyTEM 306 HP and 312 HP system reach 800-1000 m depth in this area, down to a maximum of 1200 m. Combining the three AEM datasets, a high-resolution / deep

SUMMARY

We present the integration of airborne magnetic data and five different airborne electromagnetics data sets spaning from 3 000 NIA up to 1 000 000 NIA magnetic moments (three different AEM systems were used) in La Réunion volcanic island. Subsequently, a 3D geological model of the first kilometer beneath the Plaine des Fougères was built, in order to constrain 3D hydrogeological modeling. This approach allowed for the correlation of different datasets, providing a comprehensive image of the subsurface and enabling a greater hydrogeological understanding. It was used to position the route of a deep water drainage gallery and has great potential for applications in other areas. penetration 3D resistivity block of the Plaine des Fougères was generated from the surface down to 1000 m depth. The careful preprocessing of each dataset, their inversion and their combination into a single resistivity model is explained below.

DATA PROCESSING

BRGM processing methodology of AEM data is presented in [START_REF] Reninger | Processing methodology for regional AEM surveys and local implications[END_REF], it includes the use of the singular value decomposition [START_REF] Reninger | Singular value decomposition as a denoising tool for airborne time domain electromagnetic data[END_REF]. Also, to increase the signal/noise ratio, this processing includes an adaptative stack that changes its span according to the noise level, optimizing resolution and depth of investigation. A manual editing of the residual noise was also achieved to complete the processing of each EM dataset. Such a processing strategy has proven to be essential for imaging resistivity contrasts up to one kilometer.

A joint smooth SCI inversion [START_REF] Viezzoli | Spatially constrained inversion for quasi 3D modelling of aem data[END_REF] was run considering all AEM datasets. This inversion resulted at the location of each AEM data in a 1D-resistivity model, made of 40 layers, having a fixed thickness and variable resistivity with depth.

A 3D resistivity model was then generated by the interpolation of the 304 and 306HP 1D resistivity models for the near surface (the first 200 m) and the interpolation of 306HP and 312HP 1D resistivity models for greater depths (from 200 to 1000 meters); the estimated depth of investigation was taken into account for interpolating data. The resulting 3D resistivity model takes full advantage of the resolution of the 304 and 306 HP systems for the 200 upper meters and of the depth of investigation of the 306HP and 312HP systems for mapping deeper geological bodies, allowing for greater modeling of the geological structure and understanding of the hydrogeology of the Plaine des Fougères.

During the AEM surveys, airborne magnetic data were also acquired. They were all combined to obtain a high-resolution magnetic map. Given that the surveys were conducted with the same magnetic instrumentation and parameters (acquisition height, diurnal removal, the assembly was quite simple: data points measured along flight lines of the surveys were checked for consistency and subsequently gathered in a common database. A high resolution 100 m magnetic anomaly grid was derived using a simple interpolation of this new combined dataset (Figure 2).

Oriented rock sampling was made by geologists over the area from May to June 2021. A sample magnetization library was compiled, combining laboratory measurements (performed on oriented samples with a JR5-A magnetometer) and field magnetic susceptibility (with a kappameter). This library provides rock magnetization constraints for deposits of each activity period of the Piton des Neiges (the volcano on which is located the study area). A total of 75 rock samples were used as reference for further magnetic modelling.

Taking into account the volcanic history, 2D magnetic modeling was conducted, using the geological map, the AEM resistivity model as background preliminary geometry and rock sample magnetic measurements to constrain the geological magnetic responses. The obtained modelling generated geologically realistic cross-sections (example in Figure 5). The resistivity model provided the overall geometry of the volcanic layers and their degree of alteration. The magnetic data refined the geological model in depth, thereby highlighting distinct geological units within relatively homogenous resistivity layers.

The integration of the different datasets had the added benefit of providing a more detailed geological understanding of the area and allowing for more accurate interpretations of the hydrogeology.

RESULTS

Combining all the available geological information with series of 2D geophysically constrained cross-sections, a 3D geological model was derived all over the Plaine des Fougères on a depth range of -1000 m below the topography (1500 meters above sea level). We used GeoModeller® software which previously proved to allow the successful integration of geophysics in complex 3-D geological models (e.g. [START_REF] Martelet | Integrated 3D geophysical and geological modelling of the Hercynian Suture Zone in the Champtoceaux area (South Brittany, France)[END_REF][START_REF] Calcagno | Geological modelling from field data and geological knowledge: Part I. Modelling method coupling 3D potentialfield interpolation and geological rules[END_REF]. Figure 3 displays the resulting 3D geological model, built up upon geological units classified by their increasing age, ranging from PN0 to PN4. Based on the interpretation of resistivity signatures, some of the modelled units are differentiated based on their alteration or water saturation. The 3D geological model was developed to assess the hydrogeological risk of different scenarios of gallery route. The 3D model was therefore implemented in a hydrogeological model performed in Visual KARSYS, a free software designed for assessing hydrogeological risk [START_REF] Malard | Visual KARSYS, a web-platform for the documentation of karst aquifers including online geological modelling[END_REF]. The 3D geological geometries were used as input (as well as springs, geological sounding…) to simulate the hydrogeological conditions within the 3D model. The analysis of this hydrogeological model was focused along several possible gallery routes. Taking into account the expected lithological and hydrological characteristics of the modelled geological layers the hydrogeological risk was assessed, along the different gallery routes, depending on the proximity of overlying aquifers and hydraulic head.

CONCLUSIONS

Combining high-resolution near surface and deep AEM surveys with aeromagnetic data, geological knowledge and rock sample measurements, it was possible to generate a comprehensive 3D geological model down to a depth of 1000 meters below the Plaine des Fougères. The deep AEM investigation was possible with the combination of a high power emission and a resistive underground (mainly above 300 Ohm.m in the first 400 meters). The generated 3D geological model allowed for the characterisation of the subsurface materials, as well as a better understanding of the geological structure and hydrogeological system of the area. This made it possible to map geometry of hydrogeological bodies within the first kilometer below the surface, and thus reduce the risk associated with the drilling of a deep gallery. Since, previous deep tunnelling in the area faced enormous difficulties (and loss of money), anticipating the geological and hydrogeological environment of tunnelling is of utmost importance to mitigate the geological risks prior committing heavy construction operations.

Figure 1 :

 1 Figure 1: Geological map over the Plaine des Fougères. GSAM gallery is represented with a white line, the possible area of the future gallery in a white polygon and in blue, the aquifers crossed during the GSAM drilling. Airborne acquisitions are represented in black (304 system), light orange (306 HP system) and green (312 HP system).

Figure 2 :

 2 Figure 2: Magnetic anomaly map made with the fusion of 2014 and 2021 surveys (green lines show 2021 flight path).

Figure 3 :

 3 Figure 3: 3D geological model of Plaine des Fougères.

Figure 4 :

 4 Figure 4: South-North resistivity profiles acquired with 304 (a), 306HP (b) or 312HP (c) systems in the Plaine des Fougères area. (a) is oriented S/N, (b) and (c) are overlapping, SW/NE. Dotted blue line shows a point where lines are intersecting.

Figure 5 :

 5 Figure 5: Example of an interpolated resistivity profile along the existing underground Salazie amont gallery, with geological interpretation. PN0 to PN4 represent old to recent volcanic formations.

th International Airborne Electromagnetics Workshop, 3-7 th September 2023, Fitzroy Island

th International Airborne Electromagnetics Workshop, 3-7 th September 2023, Fitzroy Island

ACKNOWLEDGMENTS

This project is funded by the County Council of Réunion (CD974) and the Water Office of La Réunion. We thank the University of La Réunion for cooperation in sampling volcanic rocks.