Jean-François Coeurjolly
email: jean-francois.coeurjolly@upmf-grenoble.fr.

Frédéric Lavancier
email: frederic.lavancier@univ-nantes.fr.

Parametric estimation of pairwise Gibbs point processes with infinite range interaction

Keywords: Lennard-Jones potential, pseudolikelihood, central limit theorem

Parametric estimation of pairwise Gibbs point processes with infinite range interaction

Introduction

Spatial Gibbs point processes are an important class of models used in spatial point pattern analysis [START_REF] Van Lieshout | Markov point processes and their applications[END_REF][START_REF] Møller | Statistical inference and simulation for spatial point processes[END_REF][START_REF] Illian | Statistical analysis and modelling of spatial point patterns[END_REF]. Gibbs point processes can be viewed as modifications of the Poisson point process in order to introduce dependencies, such as attraction or repulsion, between points. These models initially arise from statistical physics to approximate the interaction between pairs of particles [START_REF] Ruelle | Statistical Mechanics[END_REF][START_REF] Preston | Random fields[END_REF][START_REF] Georgii | Gibbs measures and phase transitions[END_REF]. The most well-known example is the Lennard-Jones model [START_REF] Lennard-Jones | On the determination of molecular fields[END_REF] which yields repulsion at short scales and attraction at long scales.

Assuming that the Gibbs model has a parametric form, an important question concerns the estimation of the parameters from a realization of the point process observed on a finite subset of R d . Popular solutions include likelihood (e.g. [START_REF] Ogata | Estimation of interaction potentials of spatial point patterns through the maximum likelihood procedure[END_REF][START_REF] Huang | Improvements of the maximum pseudo-likelihood estimators in various spatial statistical models[END_REF], pseudolikelihood (e.g. [START_REF] Besag | Some methods of statistical analysis for spatial data[END_REF][START_REF] Jensen | Pseudolikelihood for exponential family models of spatial point processes[END_REF]Baddeley and Turner, 2000b) and logistic regression [START_REF] Baddeley | Logistic regression for spatial Gibbs point processes[END_REF]. The two latter methods are more interesting from a practical point of view as they avoid the computation of the normalizing constant in the likelihood, which is in most cases inaccessible for Gibbs point processes and must be approximated by simulation-based methods. We focus in this paper on the pseudolikelihood and logistic regression methods.

When the Gibbs model is assumed to have a finite range interaction, consistency and asymptotic normality of the pseudolikelihood and logistic regression estimators are established in [START_REF] Jensen | Pseudolikelihood for exponential family models of spatial point processes[END_REF]; [START_REF] Jensen | On asymptotic normality of pseudolikelihood estimates of pairwise interaction processes[END_REF]; [START_REF] Billiot | Maximum pseudolikelihood estimator for exponential family models of marked Gibbs point processes[END_REF]; [START_REF] Dereudre | Campbell equilibrium equation and pseudo-likelihood estimation for non-hereditary Gibbs point processes[END_REF]; [START_REF] Coeurjolly | Asymptotic properties of the maximum pseudolikelihood estimator for stationary Gibbs point processes including the Lennard-Jones model[END_REF]; [START_REF] Baddeley | Logistic regression for spatial Gibbs point processes[END_REF], for large families of Gibbs models. The finite range assumption means that there exists R ą 0 such that the particles do not interact at all if they are at a distance greater than R ą 0 apart. For the two aforementioned inference methods, this assumption turns out to be crucial from both a practical point of view and a theoretical point of view, as explained below. However this assumption may imply an artificial discontinuity of the interactions between particles, where two particles at a distance R ´ apart interact while they do not at a distance R ` , for any small ą 0. This is for instance the case for the widely used Strauss model, see e.g. [START_REF] Møller | Statistical inference and simulation for spatial point processes[END_REF]. In fact, this assumption rules out many interesting Gibbs models from statistical physics like the Lennard-Jones model. The purpose of this work is to extend the pseudolikelihood and logistic regression methods to infinite range interaction Gibbs models.

From a practical point of view, an important issue is edge effects. Assume we observe a Gibbs point process with finite range interaction R ą 0 on a window W Ă R d . Then the pseudolikelihood computed on W actually depends on the point process on W ' R, where W ' R denotes the dilation of W by a ball with radius R. Some border correction is often used to make the pseudolikelihood score unbiased. An obvious solution is to compute the pseudolikelihood on the eroded set W a R, and in view of pW aRq'R Ď W (see [START_REF] Chiu | Stochastic geometry and its applications[END_REF]) the observation of the point process on W is sufficient for the computation. From a theoretical point of view, standard technical tools for unbiased estimating equations are available to derive the asymptotic properties of the associated estimator. If the Gibbs point process has infinite range interaction, then the pseudolikelihood computed on W depends on the point process over the whole space R d . It is in general impossible to apply a border correction that preserves unbiasedness of the pseudolikelihood score function. We propose in Section 2 a family of contrast functions that involve an eroded set, following the previous border correction, and a truncated range of interaction. The details are exposed in Section 2. However these contrast functions still lead to biased score functions and the standard ingredients to derive consistency and asymptotic normality of the estimators do not apply.

The strong consistency of the maximum pseudolikelihood estimator was studied by [START_REF] Mase | Consistency of the maximum pseudo-likelihood estimator of continuous state space Gibbs processes[END_REF] for pairwise interaction Gibbs point processes, including the infinite range interaction case, but under the assumption that the configuration of points outside W is known. Under the more realistic setting where the point process is observed only on W , we prove the strong consistency of our pseudolikelihood estimator in Proposition 3.1. Our result is valid for a large family of pairwise Gibbs models, namely the class of Ruelle superstable and lower regular models. The asymptotic normality is more challenging to establish. When the pseudolikelihood score function is unbiased, the main ingredient is a central limit theorem for conditionally centered random fields proved and generalized in [START_REF] Guyon | Asymptotic comparison of estimators in the Ising model[END_REF]; [START_REF] Jensen | On asymptotic normality of pseudolikelihood estimates of pairwise interaction processes[END_REF]; [START_REF] Comets | A central limit theorem for conditionally centred random fields with an application to Markov fields[END_REF]; [START_REF] Dedecker | A central limit theorem for stationary random fields[END_REF]; [START_REF] Coeurjolly | Residuals for stationary marked Gibbs point processes[END_REF]. It allows in particular to avoid mixing assumptions for Gibbs point process that are only known in restrictive frameworks (see for instance [START_REF] Heinrich | Mixing properties of Gibbsian point processes and asymptotic normality of Takacs-Fiksel estimates[END_REF]; [START_REF] Jensen | Asymptotic normality of estimates in spatial point processes[END_REF]). In our infinite range setting where the score function is biased, a new ingredient is needed. We establish in Appendix A a new central limit theorem for triangular arrays of almost conditionally centered random fields. This allows us to derive in Theorem 3.3 the asymptotic normality of our pseudolikelihood estimator. Assume the potential decreases with a hyperbolic rate with exponent γ 2 . While γ 2 ą d is (in general) required to ensure the existence of a Gibbs measure and the consistency of the pseudolikelihood estimator, we require the condition γ 2 ą 2d to ensure a central limit theorem. The asymptotic normality when d ă γ 2 ă 2d remains an open question. Nonetheless our result includes the important example of the Lennard-Jones model in dimension d " 2, for which γ 2 " 6. Proposition 3.4 discusses similar asymptotic results for the logistic regression estimator.

The remainder of this paper is organized as follows. In Section 2 we recall some basic facts about Gibbs point processes and we explain how to generalize inference methods for Gibbs models with infinite range interaction. We derive in Section 3 the asymptotic properties of our estimators. Section 4 contains a simulation study concerning the estimation of the Lennard-Jones potential, where some recommendations are derived for the practical choice of tuning parameters in the pseudolikelihood contrast function. Appendix A contains our main technical tool, namely a central limit theorem for almost conditionally centered random fields, and Appendix B gathers auxiliary lemmas.

Background and statistical methodology

Notation

We write Λ Ť R d for a bounded set Λ in R d . A configuration of points x is a locally finite subset of R d , which means that the set x Λ :" x X Λ is finite for any set Λ Ť R d . We use the notation x Λ c " xzx Λ and denote by Ω 0 the space of all locally finite configurations of points in R d . For a pp, qq matrix M with real entries, we denote by }M } " trpM J M q 1{2 its Frobenius norm where tr is the trace operator and M J is the transpose matrix of M . For a vector z P R p , }z} reduces to its Euclidean norm. For a bounded set E Ă Z d , |E| denotes the number of elements of E, while for z P R p or i P Z p , |z| and |i| stand for the uniform norm.

At many places in the document, we use the notation c to denote a generic positive constant which may vary from line to line.

Pairwise interaction Gibbs point processes

We briefly recall the needed background material on point processes and we refer to [START_REF] Daley | An Introduction to the Theory of Point Processes[END_REF] for more details. A point process is a probability measure on Ω 0 . The reference distribution on Ω 0 is the homogeneous Poisson point process with intensity β ą 0, denoted by π β . For Λ Ť R d , we write π β Λ for the restriction of π β to Λ. For any ∆ Ť R d and x P Ω 0 , N ∆ pxq denotes the number of elements of x X ∆. Let ∆ i be the unit cube centered at i P Z d . We consider the following space of tempered configurations.

Ω T " tx P Ω 0 ; Dt ą 0, @n ě 1, ÿ iPZ d ,|i|ďn N 2 ∆ i pxq ď tp2n `1q d u.
From the ergodic theorem (see [START_REF] Guyon | Random fields on a network[END_REF]), any second order stationary measure on Ω 0 is supported on Ω T . We denote by Φ : R d Ñ R Y t`8u a pair potential function, to which we associate the pairwise energy function H Λ : Ω T Ñ R Y t`8u, indexed by Borel sets Λ Ť R d and defined by

H Λ pxq " 1 2 ÿ u,v Px, u‰v, tu,vuXx Λ ‰H Φpu ´vq (2.1)
and we let

Ω " tx P Ω T , @Λ Ť R d H Λ pxq ă 8u.
Following the Dobrushin-Lanford-Ruelle formalism, see [START_REF] Preston | Random fields[END_REF], we say that P is a Gibbs measure with activity parameter β ą 0 for the pair potential function Φ if P pΩq " 1 and for P -almost every configuration x and any Λ Ť R d , the conditional law of P given x Λ c is absolutely continuous with respect to π β Λ with the density expt´H Λ pxqu{Z Λ px Λ c q, where Z

Λ px Λ c q " ş Ω expt´H Λ px Λ Y x Λ c qu π β Λ p dx Λ c q is the normalizing constant.
We use at many places in this paper the GNZ equation, after [START_REF] Georgii | Canonical and grand canonical Gibbs states for continuum systems[END_REF] and Nguyen and Zessin (1979b) (2.2)

This quantity does not depend on Λ, provided u P Λ, and can be viewed as the conditional probability to have a point in a vicinity of u, given that the configuration elsewhere is x.

Theorem 2.1 (GNZ formula). A probability measure P on Ω is a Gibbs measure with activity parameter β ą 0 for the pair potential function Φ if for any measurable function f : Ω ˆRd Ñ R such that the following expectations are finite,

E # ÿ uPX f pu, Xzuq + " E "ż f pu, Xqλpu, Xq du * (2.3)
where E denotes the expectation with respect to P .

This result can be refined by a conditional version stated in the following lemma. Its proof is actually part of the initial proof of (2.3), see also Billiot et al. (2008, Proof of Theorem 2) for a particular case. We reproduce the demonstration below.

Lemma 2.2 (Conditional GNZ formula). Let P be a Gibbs measure with activity parameter β ą 0, with pair potential Φ and Papangelou conditional intensity λ. Then for any measurable function f : Ω ˆRd Ñ R and for any Λ Ť R d such that the following expectations are finite

E # ÿ uPX Λ f pu, Xzuq | X Λ c + " E "ż Λ f pu, Xqλpu, Xq | X Λ c * (2.4)
where E denotes the expectation with respect to P .

Proof. By definition of the conditional law of P given x 1

Λ c E # ÿ uPX Λ f pu, Xzuq | X Λ c " x 1 Λ c + " ż Ω ÿ uPx Λ f pu, x Λ Y x 1 Λ c zuq e ´HΛ px Λ Yx 1 Λ c q Z Λ px 1 Λ c q π β Λ p dx Λ q " E π β Λ # ÿ uPX Λ f pu, X Λ Y x 1 Λ c zuq e ´HΛ pX Λ Yx 1 Λ c q Z Λ px 1 Λ c q + where E π β Λ
denotes the expectation with respect to π β Λ . From the Slivnyak-Mecke formula [START_REF] Slivnyak | Some properties of stationary flows of homogeneous random events[END_REF]; [START_REF] Mecke | Eine charakteristische eigenschaft der doppelt stochastischen poissonschen prozesse[END_REF]), we know that for any admissible measurable function h

E π β Λ # ÿ uPX hpu, Xzuq + " β E π β Λ "ż R d hpu, Xq du * .
By definition of the Papangelou conditional intensity (2.2), we also have for any u P Λ, β e ´HΛ pxYuq " e ´HΛ pxq λpu, xq. Using these two facts, we conclude by

E " ÿ uPX Λ f pu, Xzuq | X Λ c " x 1 Λ c * " β E π β Λ "ż Λ f pu, X Λ Y x 1 Λ c q e ´HΛ pX Λ Yx 1 Λ c Yuq Z Λ px 1 Λ c q du * " E π β Λ "ż Λ f pu, X Λ Y x 1 Λ c qλpu, X Λ Y x 1 Λ c q e ´HΛ pX Λ Yx 1 Λ c q Z Λ px 1 Λ c q du * " ż Ω ż Λ f pu, x Λ Y x 1 Λ c qλpu, x Λ Y x 1 Λ c q e ´HΛ px Λ Yx 1 Λ c q Z Λ px 1 Λ c q du π β Λ p dx Λ q " E "ż Λ f pu, Xqλpu, Xq | X Λ c " x 1 Λ c * .
The existence of a Gibbs measure P satisfying the above definition and characterization is a difficult question. Sufficient conditions on the pair potential Φ can be found in [START_REF] Ruelle | Statistical Mechanics[END_REF] and are also discussed in [START_REF] Preston | Random fields[END_REF]. The special case of finite range potentials, i.e. compactly supported functions Φ, is treated in [START_REF] Bertin | Existence of "nearest-neighbour" spatial Gibbs models[END_REF]. As we are mainly interested in this paper by infinite range potentials, we introduce the following assumption, that leads to the existence of at least one stationary Gibbs measure, as proved in [START_REF] Ruelle | Statistical Mechanics[END_REF].

[Φ] The potential Φ is bounded from below and there exist 0 ă r 1 ă r 2 ă 8, c ą 0 and γ 1 , γ 2 ą d such that Φpuq ě c}u} ´γ1 for }u} ď r 1 and |Φpuq| ď c}u} ´γ2 for }u} ě r 2 .

Examples of potentials satisfying [Φ] are Φpuq " }u} ´γ with γ ą d and Φpuq " e ´}u} }u} ´γ with γ ą d, in which cases the assumption is satisfied with γ 1 " γ 2 " γ. Another important example is the general Lennard-Jones pair potential defined for some d ă γ 2 ă γ 1 and some A, B ą 0 by Φpuq " A}u} ´γ1 ´B}u} ´γ2 . The standard Lennard-Jones model corresponds to d " 2, γ 1 " 12 and γ 2 " 6. The main interest of this model is that it can model repulsion at small scales and attraction at large scales.

Inference for infinite range Gibbs point processes

In this section, we extend the usual statistical methodologies available for finite range Gibbs point processes to the infinite range case. We assume that the Gibbs measure is parametric, in that the explicit expression of the associated Papangelou conditional intensity (2.2) is entirely determined by the knowledge of some parameter θ P Θ, including the activity parameter β ą 0, where Θ is an open bounded set of R p . We stress this assumption by writing λ θ instead of λ and Φ θ instead of Φ. For brevity, assumption [Φ] now means that Φ θ fulfills this assumption for any θ P Θ.

Assume that we observe the point process X in W n where pW n q ně1 is a sequence of bounded domains which converges to R d as n Ñ 8. As outlined in the introduction, the pseudolikelihood and the logistic regression methods are popular alternatives to the maximum likelihood as they do not involve the normalizing constant. The associated estimators are respectively defined as the maximum of where ρ is some fixed positive real number.

A problem however occurs. The integrals in (2.5) and (2.6) are not computable in practice because for values of u close to the boundary of W n , λ θ pu, Xq depends on X W c n which is not observed. When X has a finite range 0 ă R ă 8, meaning that Φ θ is compactly supported on the euclidean ball Bp0, Rq or equivalently that for any u P R d and any x P Ω, λ θ pu, xq " λ θ pu, x Bpu,Rq q, we can simply substitute W n by W n a R in (2.5) and (2.6), where for Λ Ť R d and some κ ě 0 the notation Λ a κ stands for the domain Λ eroded by the ball Bp0, κq. Using this border correction λ θ pu, Xq can be indeed computed for any u P W n aR. As a remaining practical issue, the integrals have to be approximated by some numerical scheme or by Monte-Carlo, see [START_REF] Baddeley | Logistic regression for spatial Gibbs point processes[END_REF] for an efficient solution.

The asymptotic properties of the pseudolikelihood and the logistic regression estimators are well understood in this finite range setting, see the references in introduction. Maximizing the log-pseudolikelihood (or the logistic regression likelihood) on W n a R is equivalent to cancel the score, i.e. the gradient of LPL WnaR pX; θq (or LRL WnaR pX; θq) with respect to θ. The key-ingredient is that both scores constitute unbiased estimating functions, since by application of the GNZ formula (2.3) their expectation vanishes when θ corresponds to the true parameter of the underlying Gibbs measure. Standard theoretical tools for unbiased estimating equations (see e.g. [START_REF] Guyon | Random fields on a network[END_REF]) can therefore be used to study the consistency and asymptotic normality of the associated estimators.

In the infinite range setting, the situation becomes more delicate since for any u, λ θ pu, Xq depends on X Λ for any Λ Ă R d . In this case, we introduce the following modifications of (2.5) and (2.6) that depend on two sequences of positive real numbers pα n q and pR n q Ą LPL Wnaαn,Rn pX; θq "

ÿ uPX Wnaαn log λ θ pu, X u,Rn zuq ´żWnaαn λ θ pu, X u,Rn q du (2.7) Ą LRL Wnaαn,Rn pX; θq " ÿ uPX Wnaαn log λ θ pu, X u,Rn zuq λ θ pu, X u,Rn zuq `ρ ´żWnaαn ρ log λ θ pu, X u,Rn q `ρ ρ du (2.8)
where we denote X u,Rn " X Bpu,RnqXWn . These expressions are computable from the single observation of X on W n , provided the integrals are approximated as usual by numerical scheme or by Monte-Carlo. Since they depend on two sequences α n and R n , (2.7) and (2.8) actually form a general family of contrast functions, important particular cases being the choices

• α n " R n , which agrees with the classical border correction for finite range interaction models with range R taking R n " R;

• R n " 8, accounting for the maximal possible range of interaction;

• R n " 8 and α n " 0, which is a particular case of the previous choice where in addition no erosion is considered.

We study in the next section the asymptotic properties of estimators derived from (2.7) and (2.8) for a wide class of sequences α n and R n , and based on a simulation study in Section 4, we give some recommendations for the choice of these sequences in practice. From a theoretical point view, these contrast functions introduce new challenges since the gradients of Ą LPL Wnaαn,Rn pX; θq and Ą LRL Wnaαn,Rn pX; θq are no longer unbiased estimating equations in the infinite range case. To overcome this difficulty we prove a new central limit theorem in Appendix A for almost conditionally centered random fields.

Asymptotic properties

We present asymptotic properties of the maximum pseudolikelihood estimate, derived from (2.7), for infinite range Gibbs point process. Similar results for the max-imum logistic regression derived from (2.8) are presented at the end of this section without proof. We assume the window of observation expands to R d as follows.

[W n ] pW n q is an increasing sequence of convex compact sets, such that W n Ñ R d as n Ñ 8.

We focus on exponential family models of pairwise interaction Gibbs point processes and rewrite the model (2.2) for any u P R d and x P Ω as

λ θ pu, xq " βe ´řvPx Φ θ pv´uq " e ´θJ tpu,xq (3.1)
with θ 1 " ´log β and t " pt 1 , . . . , t p q J where t 1 pu, xq " 1 and

t m pu, xq " ÿ vPx g m pv ´uq, m " 2, . . . , p. (3.2)
In that connection, our framework amounts to assume that Φ " ř p m"2 θ m g m . For convenience we let g 1 " 0 and we denote by g the p-dimensional vector g " p0, g 2 , . . . , g p q J . We make the following assumption on g.

[g] For all m ě 2, g m is bounded from below and there exist γ 1 , γ 2 ą d and c g , r 0 ą 0 such that (i) @}x} ă r 0 and @θ P Θ, θ 2 g 2 pxq ě c g }x} ´γ1 (ii) @m ě 3, g m pxq " op}x} ´γ1 q as }x} Ñ 0 (iii) @m ě 2 and @}x} ě r 0 , |g m pxq| ď c}x} ´γ2 .

Since Θ is bounded, [g] implies [Φ] which yields that for any θ P Θ there exists a Gibbs measure P θ . Assumption [g] allows us to specify which function g m is responsible for the behavior at the origin of Φ θ , namely g 2 . Note that the Lennard-Jones model defined in Section 2.2 (and the other examples presented in this section) fits this setting with θ 2 " A, θ 3 " ´B, g 2 puq " }u} ´γ1 and g 3 puq " }u} ´γ2 . In the sequel, θ ‹ stands for the true parameter vector to estimate. In other words, we assume observing a realization of a spatial point process X with Gibbs measure P θ ‹ on W n .

For exponential family models (3.1) the score function of the log-pseudolikelihood defined by (2.7) writes s Wnaαn,Rn pX; θq where for any

∆ Ă R d s ∆,Rn pX; θq " ż ∆ tpu, X u,Rn qλ θ pu, X u,Rn q du ´ÿ uPX ∆ tpu, X u,Rn zuq. (3.3)
Our first result establishes the strong consistency of the maximum pseudolikelihood based on (2.7) for infinite range Gibbs point processes and for a wide class of sequences pα n , R n q. In close relation, [START_REF] Mase | Consistency of the maximum pseudo-likelihood estimator of continuous state space Gibbs processes[END_REF] proved the strong consistency of estimators derived from (2.5). As pointed out in Section 2.3, the form (2.5) of log-pseudolikelihood is however unusable as it can only be computed if X is observed on R d . We obtain the same result but for estimators derived from the computable pseudolikelihood given by (2.7). [START_REF] Preston | Random fields[END_REF], it is sufficient to prove consistency for ergodic measures. So, we assume here that P θ ‹ is ergodic. Since Θ is an open bounded set, and by convexity of θ Ñ ´Ą LPL Wnaαn,Rn px; θq, then from Guyon (1995, Theorem 3.4.4) we only need to

prove that K n pθ, θ ‹ q " |W n a α n | ´1 ! Ą LPL Wnaαn,Rn pX; θ ‹ q ´Ą LPL Wnaαn,Rn pX; θq ) Ñ Kpθ, θ ‹ q almost surely as n Ñ 8
, where θ Ñ Kpθ, θ ‹ q is a nonnegative function which vanishes at θ " θ ‹ only. We decompose K n pθ, θ ‹ q as the sum of the three terms T 1 `T2 pθ ‹ q ´T2 pθq where for any θ P Θ

T 1 " |W n a α n | ´1 tLPL Wnaαn pX; θ ‹ q ´LPL Wnaαn pX; θqu T 2 pθq " |W n a α n | ´1 ! Ą LPL Wnaαn,Rn pX; θq ´LPL Wnaαn pX; θq ) .
Lemma B.2 shows in particular that λ θ p0, Xq and |θ J tp0, Xq|λ θ ‹ p0, Xq have finite expectation under P θ ‹ . Hence, using the ergodic theorem for spatial processes of Nguyen and Zessin (1979a), we can follow the proof of [START_REF] Mase | Consistency of the maximum pseudo-likelihood estimator of continuous state space Gibbs processes[END_REF] or the proof of Billiot et al. (2008, Theorem 1) to prove that T 1 Ñ Kpθ, θ ‹ q almost surely as n Ñ 8 where

Kpθ, θ ‹ q " E ´λθ ‹ p0, Xq " e pθ ‹ ´θq J tp0,Xq ´t1 `pθ ‹ ´θq J tp0, Xqu
ıw hich is a nonnegative function that vanishes at θ " θ ‹ only, under the identifiability condition (3.5). So the rest of the proof consists in proving that T 2 pθq Ñ 0 almost surely for any θ P Θ. We have Concerning the remaining term T 1 2 , we have

T 2 pθq " T 1 1 `T 1 2 where T 1 1 " |W n a α n | ´1 ÿ uPX Wnaαn θ J ttpu, X u,Rn zuq ´tpu, Xzuqu T 1 2 " |W n a α n | ´1 ż Wnaαn tλ θ pu, Xq ´λθ pu, X u,
T 1 2 " |W n a α n | ´1 ż Wnaαn λ θ pu, Xq ! 1 ´e´θ J ř vPX pu,Rnq c gpv´uq ) du.
We can use exactly the same decomposition as in (3.6) by introducing δ n , then use Lemmas B.1 and B.2 to apply the ergodic theorem, leading to T 1 2 Ñ 0 almost surely. These details are omitted. Hence T 2 pθq Ñ 0 almost surely for any θ P Θ and the proof of Proposition 3.1 is completed.

As a preliminary result towards the asymptotic normality of p θ Ą LPL , we state in the next lemma general conditions on the sequences α n and R n leading to the equivalence in probability of the score functions, up to |W n | 1{2 . Lemma 3.2. Assume that [W n ] and [g] hold. Let D be the set of all sequences pα n , R n q such that α n |W n | ´1{d Ñ 0 and such that there exists

0 ă γ 1 ă γ 2 ´d such that α ´γ1 n |W n | 1{2 Ñ 0 and R ´γ1 n |W n | 1{2 Ñ 0.
Then, for any pα n , R n q P D and pα 1 n , R 1 n q P D we have

s Wnaαn,Rn pX; θ ˚q ´sWnaα 1 n ,R 1 n pX; θ ˚q " o P p|W n | 1{2 q. (3.8)
Proof. For ∆ Ť R d , we denote for short s ∆,Rn :" s ∆,Rn pX; θ ˚q and set

s 1 ∆,Rn " ż ∆ tpu, X u,Rn qλ θ ‹ pu, Xq du ´ÿ uPX ∆ tpu, X u,Rn zuq. (3.9)
We prove below that for any pα n , R n q P D, s where ∆ j is the unit cube centered at j P Z d , ∆ n,j " ∆ j X pW n zpW n a α n qq and J n Ă Z d is the set such that W n zpW n a α n q " Y jPJn ∆ n,j . We have

P p|W n | 1{2 q. First, we have A " A 1 `A2 with A 1 " ż Wnaαn ttpu, X u,Rn qλ θ ‹ pu, X u,Rn q ´tpu, X Wn qλ θ ‹ pu, X Wn qu du, A 2 " ÿ uPX Wnaαn ttpu, X Wn zuq ´tpu, X u,
VarpCq " ÿ j,kPJn Covps 1 ∆ n,j ,8 , s 1 ∆ n,k ,8 q.
It is not difficult to check that all results of Lemma B.4 (i)-(iii) remain true if the ball B u,n " Bpu, α n q therein is replaced by W n , or in other words s 1 ∆ is replaced by s 1 ∆,8 . Therefore, from (iii) of Lemma B.4

} VarpCq} " Op|J n |q " Op|W n zpW n a α n q|q " Op|W n | 1´1{d α n q showing that |W n | ´1} VarpCq} Ñ 0. Hereby |W n | ´1{2 C " o P p1q and the proof is completed.
The next result establishes the asymptotic normality of the score function associated to the modified pseudolikelihood Ą LPL Wnaαn,Rn pX; θq at the true value of the parameter θ " θ ‹ , whenever pα n , R n q belongs to the set D introduced in Lemma 3.2. The proof relies on a new central limit theorem stated in Appendix A. As a consequence we deduce the asymptotic normality of p θ Ą LPL . These results require the following notation: let Σ 8 and U 8 the pp, pq matrices Theorem 3.3. Under the assumptions of Lemma 3.2 with γ 2 ą 2d and the assumption that Σ 8 is a positive definite matrix, then we have the two following convergences in distribution for any pα n , R n q P D, as n Ñ 8, (i)

Σ 8 " E tp0, Xqtp0, Xq J λ θ ‹ p0, Xq ( `żR d E tp0, Xqtpv, Xq J λ θ ‹ p0, Xqλ θ ‹ pv, Xq ( t1 ´e´Φ θ ‹ pvq u dv `żR d E λ θ ‹ p0, Xqλ θ ‹ pv, Xq ( gpvqgpvq J e ´Φθ ‹ pvq dv (3.11) U 8 " E tp0, Xqtp0, Xq J λ θ ‹ p0, Xq ( . ( 3
|W n | ´1{2 s Wnaαn,Rn pX; θ ‹ q d Ñ N p0, Σ 8 q, (ii) |W n | 1{2 ´p θ Ą LPL ´θ‹ ¯d Ñ N `0, U ´1 8 Σ 8 U ´1 8 ˘.
Some remarks on this theorem are in order. The condition γ 2 ą 2d is clearly the most restrictive one. Nonetheless it includes the standard Lennard-Jones model in dimension d " 2 for which γ 2 " 6. Under [g], existence of the model is ensured if γ 2 ą d but it remains an open problem to prove the asymptotic normality of the pseudolikelihood estimator when d ă γ 2 ă 2d. Concerning the set D of possible sequences pα n , R n q, it includes the natural choices pα n , R n q " pα n , α n q and pα n , R n q " pα n , 8q discussed in Section 2.3, provided α n tends to infinity at a good rate. However, D does not include the particular case pα n , R n q " p0, 8q, whereas this choice leads to a consistent estimator as proved in Proposition 3.1. In fact, when the erosion parameter α n does not tend to infinity, some edge effects occur due to the infinite range of the process. These edge effects are negligible with respect to |W n | but not with respect to |W n | 1{2 . Finally, following [START_REF] Coeurjolly | Fast covariance estimation for innovations computed from a spatial Gibbs point process[END_REF], it is possible to construct a fast estimator of the asymptotic covariance matrices Σ 8 and U 8 , but its asymptotic properties are out of the scope of the present paper.

Proof. (i) Since γ 2 ą 2d, there exists ε ą 0 and a ą 0 such that

d 2pγ 2 ´d ´εq ă a ă 1 2 `ε. (3.13)
For such ε, a, we let γ 1 " γ 2 ´d ´ε and consider the particular case R n " α n "

|W n | a{d . Then, α n |W n | ´1{d Ñ 0 and α ´γ1 n |W n | 1{2 " R ´γ1 n |W n | 1{2 " |W n | d{2´apγ 2 ´d´εq d Ñ 0.
This particular choice of pα n , R n q thus belongs to D. From Lemma 3.2, if we prove (i) for this choice, then the same convergence holds true for all sequences pα n , R n q in D, completing the proof of (i). Henceforth we let R n " α n " |W n | a{d where a is such that (3.13) holds.

Denoting by ∆ j the unit cube centered at j P Z d , we let ∆ n,j " ∆ j X pW n a α n q and I n Ă Z d the set such that W n a α n " Y jPIn ∆ n,j . At several places in the proof the sequence ρ n " |I n | 1{2 {α γ 1 n is involved. Then,

ρ n " |I n | 1{2 α γ 1 n " O ! |W n | d{2´apγ 2 ´d´εq d ) (3.14)
tends to 0. We write for short s ∆ n,j " s ∆ n,j ,Rn pX; θ ‹ q " s ∆ n,j ,αn pX; θ ‹ q and B u,n " Bpu, α n q. Note that for any u P W n a α n , tpu, X u,Rn q " tpu, X u,αn q " tpu, X Bu,n q and similarly λ θ ‹ pu, X u,Rn q " λ θ ‹ pu, X u,αn q " λ θ ‹ pu, X Bu,n q. Therefore for any j P I n ,

s ∆ n,j " ż ∆ n,j
tpu, X Bu,n qλ θ pu, X Bu,n q du ´ÿ uPX ∆ n,j tpu, X Bu,n zuq.

Letting Z n,j " s ∆ n,j ´Eps ∆ n,j q, we have

s Wnaαn " S n `Eps Wnaαn q
where S n " ř jPIn Z n,j . Define

p Σ n " ÿ jPIn ÿ kPIn |k´j|ďαn Z n,j Z J n,k and Σ n " E p Σ n .
The proof of (i) is completed if we show that " op|I n |q which, following the remark after Theorem A.1, satisfies the assumption of Theorem A.1 if Assumption (a) of this theorem is satisfied for any q ě 1. And the latter holds by definition of Z n,j and Lemma B.3.

Σ ´1{2 n S n d Ñ N p0, I p q, Σ ´1{2 n Eps Wnaαn q Ñ 0 and |W n | ´1Σ n Ñ Σ 8 . Let
Concerning assumption (b), we use for short the notation s 1 ∆ " s 1 ∆,αn for any ∆ Ť R d , see (3.9), namely

s 1 ∆ " ż ∆ tpu, X Bu,n qλ θ ‹ pu, Xq du ´ÿ uPX ∆ tpu, X Bu,n zuq.
Note that from the GNZ formula E s 1 ∆ " 0. We have from Lemma B.4, for any sequence Lemma B.3 and from (3.14).

J n Ă I n such that |J n | Ñ 8, ÿ j,kPJn › › EpZ n,j Z J n,k q › › " ÿ j,kPJn › › Covps ∆ n,j , s ∆ n,k q › › ď ÿ j,kPJn › › Covps ∆ n,j , s ∆ n,k q ´Covps 1 ∆ n,j , s 1 ∆ n,k q › › `ÿ j,kPJn › › Covps 1 ∆ n,j , s 1 ∆ n,k q › › ď ÿ j,kPJn ˜c α γ 1 n p1 `|k ´j| γ 2 q `c α 2γ 1 n ¸`ÿ j,kPJn |j´k|ď2r 0 › › Covps 1 ∆ n,j , s 1 ∆ n,k q › › `ÿ j,kPJn |j´k|ą2r 0 c|k ´j| ´γ2 ď c α ´γ1 n |J n | `c|J n |ρ 2 n `c|J n |} Varps 1 ∆ n,j q} `c|J n | which is Op|J n |q by
Since Σ 8 is assumed to be a positive definite matrix, assumption (c) holds if we prove that |I n | ´1Σ n Ñ Σ 8 as n Ñ 8. For this, let

Σ 1 n " Varps 1 Wnaαn q " ÿ j,kPIn E ! s 1 ∆ n,j ps 1 ∆ n,k q J ) and p Σ 1 n " ÿ j,kPIn |k´j|ďαn s 1 ∆ n,j ps 1 ∆ n,k q J .
We have

} |I n | ´1Σ n ´Σ8 } ď T 1 `T2 `T3
where

T 1 " |I n | ´1} Ep p Σ n q ´Ep p Σ 1 n q}, T 2 " |I n | ´1} Ep p Σ 1 n q ´Σ1 n }, T 3 " }|I n | ´1Σ 1 n ´Σ8 }.
First applying Lemma B.4

T 1 ď |I n | ´1 ÿ jPIn ÿ kPIn |k´j|ďαn } EpZ n,j Z J n,k q ´Ets 1 ∆ n,j ps 1 ∆ n,k q J u} " |I n | ´1 ÿ jPIn ÿ kPIn |k´j|ďαn › › Covps ∆ n,j , s ∆ n,k q ´Covps 1 ∆ n,j , s 1 ∆ n,k q › › " Opα ´γ1 n q `Opα d´2γ 1 n q
and T 1 Ñ 0 since γ 2 ą 2d implies 2γ 1 ´d ą 0 as soon as ε ă d{2, which can be assumed without loss of generality. Second, from (iii) in Lemma B.4

T 2 ď |I n | ´1 ÿ jPIn ÿ kPIn |k´j|ąαn } Ets 1 ∆ n,j ps 1 ∆ n,k q J u} ď c |I n | ´1 ÿ jPIn ÿ kPIn |k´j|ąαn 1 |k ´j| γ 2 ď c ÿ |i|ąαn 1 |i| γ 2 " op1q.
Finally T 3 Ñ 0 from (ii) in Lemma B.4, which concludes the proof of condition (c) of Theorem A.1.

To prove assumption (d), we apply the conditional GNZ formula (2.4) to write, for any j P I n ,

E `Zn,j | X ∆ n,k , k ‰ j " E « ż ∆ n,j tpu, X Bu,n q λ θ ‹ pu, Xq ´λθ ‹ pu, X Bu,n q ( du | X ∆ n,k , k ‰ j ff ´E ż ∆ n,j
tpu, X Bu,n q λ θ ‹ pu, Xq ´λθ ‹ pu, X Bu,n q ( du.

From Lemma

E › › E `Zn,j | X ∆ n,k , k ‰ j ˘› › ď c α γ 1 n EtY p0, Xqu " Opα ´γ1 n q. Hence |I n | ´1{2 ÿ jPIn E } E `Zn,j | X ∆ n,k , k ‰ j ˘} " Opρ n q
tends to 0 from (3.14). All conditions of Theorem A.1 are therefore satisfied, which yields that Σ

´1{2 n S n d Ñ N p0, I p q. The convergence |W n | ´1Σ n Ñ Σ 8 is an immediate consequence of assumption (c) checked above. It remains to prove that Σ ´1{2 n Eps Wnaαn q Ñ 0. This is a consequence of the GNZ formula, (A.3), Lemma B.1 and the condition (3.13) since

}Σ ´1{2 n Eps Wnaαn q} ď }Σ ´1{2 n } } Eps Wnaαn q} ď c |I n | ´1{2 α ´γ1 n E ż Wnaαn Y pu, Xq du " O pρ n q " op1q.
(ii) It is worth repeating that θ Ñ ´Ą LPL Wnaαn,Rn px; θq is a convex function with Hessian matrix given by (3.4). Following Lemmas B.1-B.2 and arguments developed in the proof of Proposition 3.1, we leave the reader to check that almost surely d dθ dθ J LPL Wnaαn pX; θq Ñ E tp0, Xqtp0, Xq J λ θ p0, Xq ( as n Ñ 8, which equals to U 8 when θ " θ ‹ . We also note that (3.5) implies that U 8 is a positive definite matrix. These facts and (i) allow us to apply Guyon (1995, Theorem 3.4.5) to deduce the result.

|W n a α n | ´1 " d dθ dθ J Ą LPL Wnaαn,
The following proposition focuses on the maximum logistic regression and states its strong consistency and asymptotic normality. The result is given without proof, but we claim that it follows by the same arguments as those involved in the proofs of Proposition 3.1 and Theorem 3. (4.1) We chose β " 100 and σ " 0.1 and considered three cases where ε takes the values 0.1, 0.5 and 1 respectively, which, following [START_REF] Baddeley | Variational estimators for the parameters of Gibbs point process models[END_REF] we call low, moderate and high rigidity models. The realizations are generated using the Metropolis-Hastings algorithm, implemented in the R package spatstat [START_REF] Baddeley | Spatstat: an R package for analyzing spatial point patterns[END_REF][START_REF] Baddeley | Spatial Point Patterns: Methodology and Applications with R[END_REF], on W n " r´n, ns 2 and for n " 1{2, 1, 2. To take into account the infinite range characteristic of the Lennard-Jones model, the processes are simulated on r´n ´2, n `2s 2 and then clipped to W n . Figure 1 depicts some typical realizations on r´1, 1s 2 .

|W n | 1{2 ´p θ Ą LRL ´θ‹ ¯d Ñ N `0, V ´1 8 Γ 8 V ´1 8 where denoting hpu, xq " ρ tpu, xq{tλ θ ‹ pu, xq `ρu for any u P R d , x P Ω, Γ 8 " E hp0, Xqhp0, Xq J λ θ ‹ p0, Xq ( `żR d E hp0, Xqhpv, Xq J λ θ ‹ p0, Xqλ θ ‹ pv, Xq ( t1 ´e´Φ θ ‹ pvq u dv `żR d E λ θ ‹ p0, Xqλ θ ‹ pv, Xq∆ v hp0, Xq∆ 0 hpv, Xq J ( e ´Φθ ‹ pvq dv V 8 " 1 ρ E hp0, Xqhp0, Xq J λ θ ‹ p0,
For each model and each observation window, we considered three versions of maximum pseudolikelihood estimators given by (2.7) of the parameter vector tlogpβq, σ, εu J : (i) α n " R n P r0.05, 0.3s, (ii) α n P r0.05, 0.3s and R n " 8, (iii) α n " 0, R n " 8. We remind that the values α n " 0 and R n " 8 respectively mean that no border erosion is considered (i.e. W n a α n " W n ) and the maximal possible range of interaction in W n is taken into account (i.e. λpu, x u,Rn q " λpu, x Wn q). Writing α n and/or R n P r0.05, 0.3s means that we evaluated the estimates for 30 values regularly sampled in r0.05, 0.3s.

We computed the pseudolikelihood estimator by using a 100ˆ100 grid of quadrature points to discretize the integral involved in (2.7). We did not use the Berman-Turner approximation, implemented in spatstat for a large class of models excluding (4.1) (see Baddeley and Turner (2000a)), because the latter may artificially lead to biased estimates for very repulsive patterns. As suggested by [START_REF] Baddeley | Variational estimators for the parameters of Gibbs point process models[END_REF], to minimise numerical problems (overflow, instability, slow convergence) we rescaled the interpoint distances to a unit equal to the true value of σ.

We define the weighted mean squared error WMSE by

WMSE " p E " ´z log β ´log β ¯2* plog βq 2 `p E pp σ ´σq 2 ( σ 2 `p E pp ε ´εq 2 ( ε 2
and we consider in the following its root RWMSE " ? WMSE. Similarly, we define the root-weighted squared bias and the root-weighted variance respectively denoted by RWSQ and RWV.

Tables 1 and2 summarize the simulation study based on 100 replications, where we report the values of RWMSE, RWSB and RWV. When α n and/or R n vary, we report in Table 1 the smallest value of RWMSE and the associated value α opt of α n between brackets. To be consistent, we report in Table 2 the values of RWSB and RWV associated to α opt . We observe that the three versions of the estimates have a RWMSE decreasing with n for the three Lennard-Jones models. In case (i) where α n " R n , the optimal value seems to be around α n " 0.15. A closer look at the estimates showed us that their average behavior (sample mean and standard deviation) fluctuate quite a lot with α n . In case (ii) where R n " 8, we observed that the biases of the estimates do not fluctuate that much with α n . Since the estimates had smaller standard deviation when the amount of information is maximal, i.e. when α n is low, this explains why the smallest value of α n led in almost all cases to the smallest RWMSE. Surprisingly, the third situation corresponding to α n " 0 and R n " 8 produced very interesting results which are optimal or close to the optimal ones in all cases considered. This estimator may be very time consuming to evaluate for very large datasets since all the points are involved in the evaluations of the Papangelou conditional intensity. Nonetheless, for the setting considered in this simulation study the computational time differences were negligible. The situation α n " 0 and R n " 8 is supported by Proposition 3.1 (consistency) but not by Theorem 3.3 (asymptotic normality). However, the normal QQ-plots in Figures 234seem to show a convergence to a Gaussian behavior for all our choices of α n and R n , with approximatively the same rate of convergence, i.e. |W n | ´1{2 , if we refer to the decreasing rate of the slopes in each QQ-plot. Note that the Gaussian behavior is less clear in the low rigidity Lennard-Jones model than in the moderate and high rigidity cases, but this seems specific to the model rather than to the estimators. In conclusion, to estimate the parameters of a Lennard-Jones model using the pseudolikelihood method, we recommend to use no erosion and no finite range correction. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Figure 1: Typical realizations on r´1, 1s 2 of a Lennard-Jones model with parameters β " 100, σ " 0.1 and ε " 0.1 (left), 0.5 (middle) and 1 (right).

A A new central limit theorem

When the Gibbs point process has a finite range, the asymptotic normality of the pseudolikelihood or the logistic regression estimators are essentially derived from a central limit theorem for conditionally centered random fields, see the references in introduction. This connection comes from the fact that in the finite range case, the score function of the pseudolikelihood (or the logistic regression) is conditionally centered, by application of the conditional GNZ formula (2.4). In the infinite range case, the score functions of the log-pseudolikelihood and the logistic regression are neither centered, nor conditionally centered. In the following theorem, the RWMSE r´1{2, 1{2s 2 r´1, 1s 2 r´2, 2s 2 Low (ε " 0.1) α n " R n P r0.05, 0.3s 3.26 (0.13) 1.25 (0.13) 0.62 (0.12) α n P r0.05, 0.3s, R n " 8 3.72 (0.05) 1.79 (0.05) 0.63 (0.06) α n " 0, R n " 8 3.5

1.66 0.69 Moderate (ε " 0.5) α n " R n P r0.05, 0.3s 0.65 (0.12) 0.34 (0.14) 0.2 (0.15) α n P r0.05, 0.3s, R n " 8 0.68 (0.05) 0.38 (0.05) 0.19 (0.05) α n " 0, R n " 8 0.59 0.33 0.18 High (ε " 1) α n " R n P r0.05, 0.3s 1.04 (0.08) 0.42 (0.16) 0.13 (0.16) α n P r0.05, 0.3s, R n " 8 1.34 (0.05) 0.36 (0.05) 0.16 (0.05) α n " 0, R n " 8 1.23 0.27 0.17

Table 1: Root-weighted mean squared errors (RWMSE) of parameters estimates for different Lennard-Jones models. The results are based on 100 replications. The realizations are generated on r´n ´2, n `2s 2 for n " 1{2, 1, 2 and the window of observation corresponds to r´n, ns 2 . When it makes sense, we indicate between brackets the value α opt of α n leading to the minimal value of RWMSE.

RWSB and RWV r´1{2, 1{2s 2 r´1, 1s 2 r´2, 2s 2 Low (ε " 0.1) α n " R n P r0.05, 0.3s 1.82 2.70 0.57 1.11 0.07 0.62 α n P r0.05, 0.3s, R n " 8 2.49 2.76 0.82 1.59 0.03 0.63 α n " 0, R n " 8 2.36 2.59 0.78 1.46 0.20 0.66 Moderate (ε " 0.5) α n " R n P r0.05, 0.3s 0.23 0.60 0.16 0.30 0.07 0.19 α n P r0.05, 0.3s, R n " 8 0.04 0.66 0.10 0.37 0.02 0.19 α n " 0, R n " 8 0.07 0.58 0.02 0.33 0.02 0.18 High (ε " 1) α n " R n P r0.05, 0.3s 0.43 0.71 0.16 0.39 0.07 0.12 α n P r0.05, 0.3s, R n " 8 0.11 1.27 0.13 0.34 0.12 0.10 α n " 0, R n " 8 0.06 1.23 0.05 0.26 0.12 0.11

Table 2: Root-weighted squared biases (RWSB) and variances (RWS) of parameters estimates for different Lennard-Jones models. The setting is the same as in Table 1. When α n and R n vary, we report the values leading to the minimal RWMSE, i.e. the values associated to α n " α opt as given in Table 1. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -2 -1 0 1 2 0.09 0.10 0.11 0.12 0.13

q [-1/2,1/2] 2 [-1,1] 2 [-2
,2] 2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -2 -1 0 1 2 0.0

0.2 0.4 0.6 0.8 1.0 q [-1/2,1/2] 2 [-1,1] 2 [-2
,2] 2 q q q q q q q q qq q q q q q q q q q q q q qq q q q q q q q q -1 0 1 2 0.10 0.12 0.14 0.16 0.18

0.20 q [-1/2,1/2] 2 [-1,1] 2 [-2
,2] 2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -2 -1 0 1 2 0.0

0.2 0.4 0.6 0.8 q [-1/2,1/2] 2 [-1,1] 2 [-2,2] 2
q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -2 -1 0 1 2 0.09 0.10 0.11 0.12 0.13

q [-1/2,1/2] 2 [-1,1] 2 [-2
,2] 2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -2 -1 0 1 2 0.0 The first row (resp. second and third) corresponds to estimates obtained with α n " R n " α opt (resp. with α n " α opt , R n " 8 and with α n " 0, R n " 8). The optimal values α opt are given in Table 1. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -2 -1 0 1 2 0.100 0.105 0.110 0.115

0.2 0.4 0.6 0.8 q [-1/2,1/2] 2 [-1,1] 2 [-2,2] 2
0.120 q [-1/2,1/2] 2 [-1,1] 2 [-2,2] 2
q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -2 -1 0

1 2 0.2 0.4 0.6 0.8 1.0 1.2 q [-1/2,1/2] 2 [-1,1] 2 [-2
,2] 2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -2 -1 0 1 2 0.095 0.100 0.105 0.110 0.115 0.120 0.125 0.130

q [-1/2,1/2] 2 [-1,1] 2 [-2
,2] 2 q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -2 -1 0 1 2 0.0 0.5 1.0

1.5 q [-1/2,1/2] 2 [-1,1] 2 [-2,2] 2
q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -2 -1 0 1 2 0.095 0.100 0.105 0.110 0.115

q [-1/2,1/2] 2 [-1,1] 2 [-2,2] 2
q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -2 -1 0 1 2 0.0 0.5 1.0

1.5 q [-1/2,1/2] 2 [-1,1] 2 [-2,2] 2
Figure 3: Normal QQ-plots for estimates of the parameter ε (left) and σ (right) for the moderate interaction Lennard-Jones model, i.e. plogpβq, ε, σq " plogp100q, 0.5, 0.1q. The first row (resp. the second and third) corresponds to estimates obtained with α n " R n " α opt (resp. with α n " α opt , R n " 8 and with α n " 0, R n " 8). The optimal values α opt are given in Table 1.

q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -2 -1 0 1 2 0.10 0.11 0.12 0.13 0.14 0.15

q [-1/2,1/2] 2 [-1,1] 2 [-2
,2] 2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -2 -1 0 1 2 0.0 0.5 1.0

1.5 q [-1/2,1/2] 2 [-1,1] 2 [-2
,2] 2 q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q -2 -1 0 1 2 0.095 0.100 0.105 0.110 0.115

0.120 q [-1/2,1/2] 2 [-1,1] 2 [-2
,2] 2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -2 -1 0 1 2

0 1 2 3 4 5 6 q [-1/2,1/2] 2 [-1,1] 2 [-2,2] 2
q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -2 -1 0 1 2 0.100 0.105 0.110

0.115 q [-1/2,1/2] 2 [-1,1] 2 [-2
,2] 2 q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -2 -1 0 1 2 0.0 0.5 The first row (resp. the second and third) corresponds to estimates obtained with α n " R n " α opt (resp. with α n " α opt , R n " 8 and with α n " 0, R n " 8). The optimal values α opt are given in Table 1.

1.0 1.5 2.0 2.5 q [-1/2,1/2] 2 [-1,1] 2 [-2,2] 2
conditional centering condition is replaced by condition (d), which turns out to be sufficient for our application to s Wnaαn,Rn pX; θ ‹ q in Theorem 3.3. The other conditions are mainly due to the non-stationary setting induced by the presence of α n and R n . They allow in particular to control the asymptotic behavior of the empirical covariance matrix in (A.1). For two square matrices A, B we write A ě B when A ´B is a positive semi-definite matrix.

Theorem A.1. For n P N and j P Z d , let X n,j be a triangular array field in a measurable space S. For n P N, let I n Ă Z d and α n P R `such that |I n | Ñ 8 and α n Ñ 8 as n Ñ 8. Define S n " ř jPIn Z n,j where Z n,j " f n,j pX n,k , k P K n,j q with K n,j " tk P Z d , |k ´j| ď α n u and where f n,j : S K n,j Ñ R p is a measurable function. We define p Σ n and Σ n by

p Σ n " ÿ jPIn ÿ kPIn |k´j|ďαn Z n,j Z J n,k and Σ n " E p Σ n .
We assume that (a) E Z n,j " 0 and there exists q ě 1 such that sup ně1 sup jPIn E }Z n,j } 4q ă 8, (b) for any sequence

J n Ă I n such that |J n | Ñ 8 as n Ñ 8, |J n | ´1 ÿ j,kPJn › › EpZ n,j Z J n,k q › › " Op1q.
E }E pZ n,j | X n,k , k ‰ jq} Ñ 0, then Σ ´1{2 n S n d Ñ N p0, I p q. (A.2)
Before detailing the proof, let us remark that if assumption (a) is valid for any q ě 1 then the result remains true if α p2`εqd n " op|I n |q for any ε ą 0.

Proof. For m, m 1 " 1, . . . , p, let ∆ mm 1 " p p Σ n ´Σn q mm 1 . Let q ě 1 be as in assumption (a), the assertion (A.1) will be proved if we prove that E `|I n | ´2q ∆ 2q mm 1 ˘" op1q. We have ∆ mm 1 " ř jPIn U n,j where

U n,j " ÿ kPI n,j
Z n,j Z J n,k ´EpZ n,j Z J n,k q ( mm 1 and I n,j " tk P I n : |k ´j| ď α n u. Let j 1 , . . . , j 2q P I n such that |j k ´j1 | ą 4α n for k " 2, . . . , 2q. It is clear that for any j P I n , U n,j depends only on X n,k for |k ´j| ď 2α n . So,

EpU n,j 1 . . . U n,j 2q q " E E `Un,j 1 . . . U n,j 2q | X n,k , |k ´j1 | ą 2α n ˘( " E U n,j 2 . . . U n,j 2q E pU n,j 1 | X n,k , |k ´j1 | ą 2α n q ( " E U n,j 2 . . . U n,j 2q E pU n,j 1 q ( " 0
whereby we deduce that Ep∆ 2q mm 1 q " ÿ j 1 ,...,j 2q PIn |j k ´j1 |ď4αn,k"2,...,2q

EpU n,j 1 . . . U n,j 2q q. Now, by condition (a) and Hölder's inequality, we have for any

j P I n E U 2q n,j " ÿ k 1 ,...,k 2q PI n,j E " Z n,j Z J n,k 1 ´EpZ n,j Z J n,k 1 q ( mm 1 . . . ! Z n,j Z J n,k 2q ´EpZ n,j Z J n,k 2q q ) mm 1 ı ď c ÿ k 1 ,...,k 2q PI n,j µ 4q ď c ˆsup jPIn # k P I n , |k ´j| ď α n ( ˙2q " Opα 2qd n q.
From Hölder's inequality, we continue with Ep∆ 2q mm 1 q ď ÿ j 1 ,...,j 2q PIn |j k ´j1 |ď4αn,k"2,...,2q

EpU 2q

n,j 1 q 1{2q EpU 2q n,j 2q q 1{2q " O α p4q´1qd

n |I n | ( leading to Ep|I n | ´2q ∆ 2q mm 1 q " O ˜αp4q´1qd n |I n | 2q´1 ¸" op1q
by assumption on α n , which completes the proof of (A.1). We now focus on (A.2) and we let

S n " Σ ´1{2 n S n , S n,j " ÿ kPI n,j Z n,k and S n,j " Σ ´1{2 n S n,j
where we recall the notation I n,j " K n,j X I n . According to Stein's method (see [START_REF] Bolthausen | On the central limit theorem for stationary mixing random fields[END_REF], in order to show (A.2) it suffices to prove that for all u P R p such that }u} " 1 and for all ω P R E piω ´uJ S n qe iωu J Sn ( Ñ 0 as n Ñ 8 where i " ? ´1. Letting v " ωu, this is equivalent to show that for all v P R p , } EpA n q} Ñ 0 where A n " piv ´Sn qe iv J Sn . We decompose the term A n in (i) For any q ě 0, Etλ θ p0, Xq q u ă 8.

(ii) Let f : R d Ñ R be a measurable function such that |f puq| ď cp1 `}u}q ´γ with γ ą d, then for any q ě 0 E ! e q| ř uPX f p}u}q| ) ă 8.

(iii) For any q ě 0, q 1 ą 0 and θ P Θ, Et|t m |p0, Xq q λ θ p0, Xq q 1 u ă 8.

(iv) Let f 1 and f 2 be two functions as in (ii), then for any q 1 , q 2 , q 3 ě 0 and q 1 ą 0,

E # |t m |p0, Xq q 1 ˇˇˇˇÿ uPX f 1 p}u}q ˇˇˇˇq 2 e q 3| ř uPX f 2 p}u}q| λ θ p0, Xq q 1 + ă 8.
Proof. The first statement is a consequence of Proposition 5.2 (a) in [START_REF] Ruelle | Superstable interactions in classical statistical mechanics[END_REF]. It relies on the following property, see also Mase (1995, Lemma 2). If ψ : R `Ñ R ìs a decreasing function with ş 8 0 ψptqt d´1 dt ă 8, then for any q ě 0, Epe q ř uPX ψp}u}q q ă 8.

The proof of (ii) is an easy consequence of this property. We deduce in particular that all moments of ř uPX f p}u}q exist and are finite. Assuming (iii) is true, then (iv) is a straightforward consequence of the previous properties and Hölder's inequality. Let us prove (iii). For any ε ą 0, using the fact that for any q ě 0, κ ą 0, x Þ Ñ x q e ´κx is bounded on r0, 8q, we have

|t m |p0, xq q λ θ p0, xq q 1 " |t m |p0, xq q e ´q1 ř p k"2 θ k t k p0,xq " |t m |p0, xq q e ´q1 ε|θmtmp0,xq| e q 1 ε|θmtmp0,xq|´q 1 ř p k"2 θ k t k p0,xq ď c e ´q1 ř uPx Φθ puq ,
where Φθ puq " ř p k"2 θ k g k puq ´ε|θ m g m puq|. The proof of (iii) is completed in view of (i) if we show that Φθ satisfies [Φ] for any θ. Write Φθ puq " Φ1 puq `Φ 2 puq with

Φ1 puq " θ 2 2 g 2 puq `p ÿ k"3 θ k g k puq, Φ2 puq " θ 2 2 g 2 puq ´ε|θ m g m puq|.
From [g], we deduce that there exists r ą 0 such that }u} ă r implies Φ1 puq ą c}u} ´γ1 . Moreover if m " 2, Φ2 puq ą 0 for all }u} ă r 0 , provided ε ă 1{2. If m ě 3, there exists r 1 such that }u} ă r 1 implies |θ m g m puq| ă c g }u} ´γ1 {p4εq where c g is the constant in [g], yielding Φ2 puq ą pc g {4q}u} ´γ1 . In all cases, we obtain that for some r 1 ą 0, }u} ă r 1 implies Φpuq ą c}u} ´γ1 . On the other hand, it is clear that if }u} ą r 0 then | Φpuq| ď c}u} ´γ2 and that Φθ is bounded from below, proving that it satisfies [Φ].

Lemma B.3. Let j P I n and s m " ps ∆ n,j q m , respectively s 1 m " ps 1 ∆ n,j q m , be the m-th coordinate of s ∆ n,j given by (B.1), respectively s 1 ∆ n,j given by (B.2). Under [g], if α n ě r 0 then, for any q P N, Ep|s m | q q ă 8 and Ep|s 1 m | q q ă 8.

Proof. The proof being similar for s m and s 1 m , we only give the details concerning s m . From (B.1) and the binomial formula, the statement is a consequence of

E $ & % ˇˇˇˇˇÿ uPX ∆ n,j t m pu, X Bu,n zuq ˇˇˇˇˇp 1 ˇˇˇˇż ∆ n,j t m pu, X Bu,n qλ θ ‹ pu, X Bu,n q du ˇˇˇˇp 2 , .
ă 8

for any p 1 , p 2 P N. Applying the Cauchy-Schwarz's inequality, we consider each term above separately. First, for any p P N, by Hölder's inequality and using Lemma B.1

we get

E "ˇˇˇˇż ∆ n,j t m pu, X Bu,n qλ θ ‹ pu, X Bu,n q du ˇˇˇp * ď c ż ∆ n,j E ˇˇt m pu, X Bu,n q ˇˇp λ θ ‹ pu, X Bu,n q p ( du ď c ż ∆ n,j E " |t m pu, Xq| p λ θ ‹ pu, Xq p `|λ θ ‹ pu, X Bu,n q ´λθ ‹ pu, Xq| p (‰ du ď c ż ∆ n,j E " |t m pu, Xq| p λ θ ‹ pu, Xq p 1 `Gp pu, xqe c p Gpu,xq (‰ du
which is finite by Lemma B.2 and the stationarity of X.

Second, we can prove by induction and successive application of the GNZ formula, see Corollary 3.1 in [START_REF] Decreusefond | Moment formulae for general point processes[END_REF], that

E » - $ & % ÿ uPX ∆ n,j t m pu, X Bu,n zuq , . - p fi fl " p ÿ k"1 ÿ pP 1 ,...,P k qPT k p E ż ∆ k n,j λ θ ‹ ptu 1 , . . . , u k u, Xq k ź "1 t |P | m pu , X Bu,n Y tuzu uq du
where T k p is the set of all partitions of t1, . . . , pu into k subsets, |P| is the cardinality of P, u " pu 1 , . . . , u k q and uzu " pu 1 , . . . , u ´1, u `1, . . . , u k q. Since

λ θ ‹ ptu 1 , . . . , u k u, Xq " k ź "1 λ θ ‹ pu , Xq k ź i"1,i‰
e ´Φθ ‹ pu i ´u q , we obtain by application of Hölder's inequality,

E $ & % ˇˇˇˇˇÿ uPX ∆ n,j t m pu, X Bu,n zuq ˇˇˇˇˇp , . - ď p ÿ k"1 ÿ PPT k p k ź "1 E 1{k ż ∆ k n,j λ θ ‹ pu , Xq k |t m pu , X Bu,n Y tuzu uq| k|P| k ź i"1,i‰ e ´kΦ θ ‹ pu i ´u q du.
The proof is completed if we show that all expectations above are finite. To that end, note that

t m pu , X Bu,n Y tuzu uq " t m pu , X Bu,n q `k ÿ h"1,h‰ g m pu h ´u q whereby, denoting q " k|P| E ż ∆ k n,j λ θ ‹ pu , Xq k |t m pu , X Bu,n Y tuzu uq| q k ź i"1,i‰
e ´kΦ θ ‹ pu i ´u q du ď q ÿ r"0 ˆq r ˙ż∆ k n,j ˇˇˇˇk ÿ

h"1,h‰ g m pu h ´u q ˇˇˇˇr k ź i"1,i‰ e ´kΦ θ ‹ pu i ´u q E |t m pu , X Bu,n q| q´r λ θ ‹ pu , Xq k ( du

ď c q ÿ r"0 k ÿ h"1,h‰ ż ∆ k n,j
|g m pu h ´u q| r k ź i"1,i‰ e ´kΦ θ ‹ pu i ´u q E |t m pu , X Bu,n q| q´r λ θ ‹ pu , Xq k ( du.

The last expectation is finite in view of Lemma B.2, so the above expression is lower than (i) For Λ 1 , Λ 2 two bounded Borel sets of R d

Covps 1 Λ 1 , s 1 Λ 2 q " E ż Λ 1 XΛ 2
tpu, X Bu,n qtpu, X Bu,n q J λ θ ‹ pu, Xq du

`E ż Λ 1 ż Λ 2
tpu, X Bu,n qtpv, X Bv,n q J tλ θ ‹ pu, Xqλ θ ‹ pv, Xq ´λθ ‹ ptu, vu, Xqu du dv

`E ż Λ 1 ż Λ 2
∆ v tpu, X Bu,n qt∆ u tpv, X Bv,n qu J λ θ ‹ ptu, vu, Xq du dv where for any u, v P R d , x P Ω and any measurable function f : R d ˆΩ Ñ R p , the difference operator ∆ v is defined by ∆ v f pu, xq " f pu, x Y vq ´f pu, xq.

(ii) Let p∆ n q be a sequence of increasing domains such that ∆ n Ñ R d as n Ñ 8, then |∆ n | ´1 Varps 1 ∆n q Ñ Σ 8 where Σ 8 is defined by (3.11).

(iii) Let j, k P I n . Then if |k ´j| ą 2r 0 , › › ›Covps 1 ∆ n,j , s 1 ∆ n,k q › › › ď c|k ´j| ´γ2 .

(iv) For any j, k P I n , if α n ě r 0 , then › › ›Covps ∆ n,j , s ∆ n,k q ´Covps 1 ∆ n,j , s 1 ∆ n,k q › › › ď c α γ 1 n p1 `|k ´j| γ 2 q `c α 2γ 1 n as n Ñ 8, where we recall that γ 1 " γ 2 ´d ´ε with 0 ă ε ă γ 2 ´d.

Proof. (i) is a slight extension of Coeurjolly and Rubak (2013, Lemma 3.1) where the case Λ 1 " Λ 2 was considered. The proof is omitted. For (ii), we note that for any u, v P R d , m ě 1 and x P Ω }tpu, X Bu,n qtpv, X Bv,n q J } ď } |t|pu, Xq|t|pv, Xq J }.

λ θ ‹ pu
The result is derived using the dominated convergence theorem, the stationarity of X and since from Lemma B.2 the random variables } |t|p0, Xq|t|p0, Xq J }λ θ ‹ p0, Xq and } |t|p0, Xq|t|pv, Xq J }λ θ ‹ p0, Xqλ θ ‹ pv, Xq have expectation uniformly bounded in v while by Plugging these inequalities in (B.5) shows (iii), as the remaining terms have finite expectations from Lemma B.2. We now focus on (iv). Let us write s ∆ n,j " s 1 ∆ n,j `Ij where I j " ş ∆ n,j π n pu, xq du and π n pu, xq " tpu, x Bu,n qtλ θ ‹ pu, x Bu,n q ´λθ ‹ pu, xqu.

We have

Covps ∆ n,j ,s ∆ n,k q ´Covps 1 ∆ n,j , s 1 ∆ n,k q " Eps 1 ∆ n,j I J k q `Eps 1 ∆ n,k I J j q `EpI j I J k q ´Eps ∆ n,j q E ps ∆ n,k q J . (B.7)

Let us control each term in (B.7). From the GNZ formula tpu, Xqλ θ ‹ pu, Xqe ´Φθ ‹ pv´uq gpv ´uq J tλ θ ‹ pv, X Bv,n q ´λθ ‹ pv, Xqu du dv.

Eps 1 ∆ n,j I J k q " E ż ∆ n,j ż ∆ n,
(B.8)

[g] implies [Φ] which in turn yields |1 ´e´Φ θ ‹ pv´uq | ď 1 `e´Φ θ ‹ pv´uq ď c since Φ θ ‹ is bounded from below. On the other hand, for any m ě 2, denoting Φθ ‹ " Φ θ ‹ ´ε|θ ‹ m g m | for some ε ą 0, we have |g m |e ´Φθ ‹ " |g m |e ´ε|θ ‹ m gm| e ´Φ θ ‹ ď c since x Þ Ñ xe ´κx is bounded on r0, 8q for any κ ą 0 and Φθ ‹ satisfies [Φ] as seen in the proof of Lemma B.2. This proves that for any u, v, }e ´Φθ ‹ pv´uq gpv ´uq J } is bounded. Moreover, from (B.6), we know that if |k ´j| ą 2r 0 , then |1´e ´Φθ ‹ pv´uq | ď c|k ´j| ´γ2 and similarly }e ´Φθ ‹ pv´uq gpv´uq J } ď c|k´j| ´γ2 . We deduce that for any u P ∆ n,j , any v P ∆ n,k and any j, k, |1 ´e´Φ θ ‹ pv´uq | ď cp1 `|k ´j|q ´γ2 and }e ´Φθ ‹ pv´uq gpv ´uq J } ď cp1 `|k ´j|q ´γ2 . Plugging these inequalities in (B.8) and applying Lemmas B.1-B.2 to the remaining terms shows that for any j, k } Eps 1 ∆ n,j I J k q} ď c α γ 1 n p1 `|k ´j| γ 2 q .

(B.9)

The same inequality obviously holds for } Eps 1 ∆ n,k I J j q}. For the two last terms in the right hand side of (B.7), namely

EpI j I J k q " E ż ∆ n,j ż ∆ n,k
π n pu, Xqπ n pv, Xq J du dv and, after application of the GNZ formula, Eps ∆ n,j q Eps ∆ n,k q J " ż ∆ n,j ż ∆ n,k E π n pu, Xq E π n pv, Xq J du dv, we deduce from Lemmas B.1-B.2 that their norm is bounded by α ´2γ 1 n for any j, k, up to a positive constant. The latter and (B.9) prove (iv).

  , recalled below, which is a characterization of a Gibbs measure. It is given in terms of the Papangelou conditional intensity λ : R d ˆΩ Ñ R defined for any Λ Q u by λpu, xq " β e ´HΛ pxYuq e ´HΛ pxq " β e ´řvPx Φpv´uq .

  .12) These matrices are indeed correctly defined, as [g] implies on the one hand that all the expectations involved are uniformly bounded in v by Lemmas B.1-B.2, and on the other hand that ż R d |1 ´e´Φ θ ‹ pvq | dv ă 8 and ż R d }gpvqgpvq J }e ´Φθ ‹ pvq dv ă 8. We denote by d Ñ the convergence in distribution.

3 .

 3 Proposition 3.4. Under the assumptions of Proposition 3.1, the maximum logistic regression estimator defined by p θ Ą LRL " argmax θPΘ Ą LRL Wnaαn pX; θq converges almost surely to θ ‹ as n Ñ 8 and under the assumptions of Theorem 3.3 it satisfies the following convergence in distribution

Figure 2 :

 2 Figure2: Normal QQ-plots for estimates of the parameter ε (left) and σ (right) for the low interaction Lennard-Jones model, i.e. plogpβq, ε, σq " plogp100q, 0.1, 0.1q. The first row (resp. second and third) corresponds to estimates obtained with α n " R n " α opt (resp. with α n " α opt , R n " 8 and with α n " 0, R n " 8). The optimal values α opt are given in Table1.

Figure 4 :

 4 Figure4: Normal QQ-plots for estimates of the parameter ε (left) and σ (right) for the high interaction Lennard-Jones model, i.e. plogpβq, ε, σq " plogp100q, 1, 0.1q. The first row (resp. the second and third) corresponds to estimates obtained with α n " R n " α opt (resp. with α n " α opt , R n " 8 and with α n " 0, R n " 8). The optimal values α opt are given in Table1.

  (i) |t m |pu, x Bu,n q| ď |t m |pu, xq (ii) |t m pu, xq ´tm pu, x Bu,n q| ď c mintGpu, xq, α ´γ1 n Hpu, xqu (iii) @θ P Θ, |λ θ pu, xq ´λθ pu, x Bu,n q| ď c e c Gpu,xq λ θ pu, xq mintGpu, xq, α ´γ1 n Hpu, xqu. Proof. The first statement is straightforward from the definition. For the second one, from [g] and since α n ě r 0 , |t m pu, xq ´tm pu, x Bu,n q| " | ÿ vPx g m pv ´uq1p}v ´u} ě α n q| ď c ÿ vPx }v ´u} ´γ2 1p}v ´u} ě α n q, which is clearly lower than c Gpu, xq. Pushing one step further, we get |t m pu, xq ´tm pu, x Bu,n q| ď c α ´γ1 n ÿ vPx }v ´u} ´d´ε 1p}v ´u} ě α n q ď c α ´γ1 n Hpu, xq, which proves (ii). For the third statement, since for all x, |1 ´ex | ă |x|e |x| , we have |λ θ pu, xq ´λθ pu, x Bu,n q| " λ θ pu, xq ˇˇ1 ´eř ´uq|1p}v ´u} ě α n q where c " pp ´1q sup θPΘ sup m |θ m | ă 8, since Θ is bounded.Lemma B.2. Under the assumption [g], then for any θ P Θ we have the following statements where E denotes the expectation with respect to P θ ‹ .

  m pvq| r e ´kΦ θ ‹ pvq dv, which is finite from [g]. Lemma B.4. The following properties hold under the assumption [g].

  θ ‹ pvq | dv ă 8 and ż R d }gpvqgpvq J }e ´Φθ ‹ pvq dv ă 8.To prove (iii), we apply (i) to the disjoint sets ∆ n,j , ∆ n,k and relations(B.3)-(B.4Xq|t|pv, Xq J }λ θ ‹ pu, Xqλ θ ‹ pv, Xq|1 ´e´Φ θ ‹ pv´uq | du dv J } λ θ ‹ pu, Xqλ θ ‹ pv, Xqe ´Φθ ‹ pv´uq du dv. (B.5)Since |k ´j| ą 2r 0 , we deduce from [g] that for any pu, vq P ∆ n,j ˆ∆n,k and any m ě 2, |g m pv ´uq| ď c|k ´j| ´γ2 . This leads to }gpv ´uqgpv ´uq J } ď c|k ´j| ´γ2 .Similarly since Φ θ ‹ " ř p m"2 θ ‹ m g m ,for any pu, vq P ∆ n,j ˆ∆n,k , e ´Φθ ‹ pv´uq ď e |Φ θ ‹ pv´uq| ď c and |1 ´e´Φ θ ‹ pv´uq | ď |Φ θ ‹ pv ´uq|e |Φ θ ‹ pv´uq| ď c|k ´j| ´γ2 . (B.6)

LPL Wn pX; θq " ÿ uPX Wn

			log λ θ pu, Xzuq	´żWn	λ θ pu, Xq du		(2.5)
	LRL Wn pX; θq "	ÿ uPX Wn	log	λ θ pu, Xzuq λ θ pu, Xzuq `ρ	´żWn	ρ log	λ θ pu, Xq ρ	`ρ	du	(2.6)

  Proposition 3.1. Assume that [W n ] and [g] hold. Then for any x P Ω the function θ Ñ ´Ą LPL Wnaαn,Rn px; θq is a convex function with Hessian matrix given by ´d dθ dθ J Moreover if α n |W n | ´1{d Ñ 0 and R n Ñ 8 as n Ñ 8, and if for any y P R p zt0u

	P y J tp0, Xq ‰ 0 (	ą 0,	(3.5)
	then the maximum pseudolikelihood estimator
	p θ

Ą LPL Wnaαn,Rn px; θq " ´d dθ J s Wnaαn,Rn px; θq " ż Wnaαn tpu, x u,Rn qtpu, x u,Rn q J λ θ pu, x u,Rn q du. (3.4) Ą LPL " argmax θPΘ Ą LPL Wnaαn,Rn pX; θq converges almost surely to θ ‹ as n Ñ 8. Proof. By [W n ] and the basic assumption on α n , we are ensured that W n a α n is a sequence of regular bounded domains of R d and that |W n aα n | Ñ 8 as n Ñ 8. Since any stationary Gibbs measure can be represented as a mixture of ergodic measures

  Let us first look at T 1 1 . By boundedness of Θ, it is sufficient to prove that |W n a α n | ´1 ř uPX Wnaαn tt m pu, X u,Rn zuq ´tm pu, Xzuqu tends to 0 almost surely. Let δ n be a sequence of real numbers such that δ n Ñ 8 and δ n |W n | ´1{d Ñ 0 as n Ñ 8. For brevity, let X pu,Rnq c :" X R d ztWnXBpu,Rnqu .In the first sum above, }u ´v} ě R n ^pα n `δn q and using the same arguments and the same notation as in (ii) of Lemma B.1 we get that the absolute value of this sum is lower than c pR n ^pα n `δn qq ´γ1 ř uPX Wnapαn`δnq Hpu, Xzuq for some γ 1 ą 0. Hence|W n a α n | ´1 ˇˇˇˇˇÿ

		ÿ
		tt m pu, X u,Rn zuq ´tm pu, Xzuqu
	uPX Wnaαn
		ÿ		ÿ
	"				g m pv ´uq
		uPX Wnaαn	vPX pu,Rnq c
		ÿ			ÿ			ÿ
	"	uPX Wnapαn`δnq	vPX pu,Rnq c	g m pv ´uq	`ÿ uPX WnaαnzWnapαn`δnq	vPX pu,Rnq c	g m pv ´uq (3.6)
								ˇˇˇˇˇÑ
								ÿ
								g m pv ´uq	0	(3.7)
						uPX Wnapαn`δnq	vPX pu,Rnq c
	almost surely. For the second sum in (3.6), using the notation |t m |pu, xq "	ř vPx |g m pvú
	q| we have
		|W n a α n |	´1 ˇˇˇˇˇÿ	ÿ	g m pv ´uq	ˇˇˇˇď
					uPX WnaαnzWnapαn`δnq	vPX pu,Rnq c
		|W n a α n | ´1	ÿ	|t m |pu, Xq
					uPX WnaαnzWnapαn`δnq
		" |W n a α n |

Rn qu du. uPX Wnapαn`δnq ÿ vPX pu,Rnq c g m pv ´uq ˇˇˇˇď c pR n ^pα n `δn qq ´γ1 |W n a pα n `δn q| |W n a α n | |W n a pα n `δn q| ´1 ÿ uPX Wnapαn`δnq Hpu, Xzuq. By Lemma B.2, the random variable |Hp0, Xq|λ θ ‹ p0, Xq has finite expectation under P θ ‹ . Moreover our assumptions ensure that W n a pα n `δn q is a sequence of regular bounded domains of R d with |W n a pα n `δn q| Ñ 8 as n Ñ 8. So by the ergodic theorem |W n a pα n `δn q| ´1 ř uPX Wnapαn`δnq Hpu, Xzuq Ñ E tHp0, Xqλ θ ‹ p0, Xqu almost surely whereby |W n a α n | ´1 ˇˇˇˇˇÿ ´1 ÿ uPX Wnaαn |t m |pu, Xq ´|W n a α n | ´1 ÿ uPX Wnapαn`δnq |t m |pu, Xq which tends to 0 almost surely by application of the ergodic theorem, Lemma B.2 and since |W n a pα `δn q|{|W n a α n | Ñ 1 as n Ñ 8. Combining this result with (3.7) in (3.6) shows that T 1 1 Ñ 0 almost surely.

  Rn zuqu.Both terms above are handled similarly and we give the details for A 2 only. Denoting A 2,m the m-th coordinate of A 2 , we obtain using the same arguments and the same notation as in Lemma B.1 (ii) that for any m ´1{2 A 2 " o P p1q. The same result holds for A 1 by the arguments developed in Lemma B.1 (ii)-(iii). Wn qtλ θ ‹ pu, X Wn q ´λθ ‹ pu, Xqu du. that |λ θ pu, xq ´λθ pu, x Wn q| " λ θ pu, xq ˇˇ1 ´eř last sum }v ´u} ě α n since u P W n a α n . The same arguments as in the proof ofLemma B.1 (iii) thus apply and we get for all u P W n a α n , |λ θ pu, xq ´λθ pu, x Wn q| ď c α ´γ1 n e c Gpu,xq λ θ pu, xqHpu, xq. From Lemma B.1 (i) and Lemma B.2, we obtain E |B| " Opα ´γ1 n |W n a α n |q " Opα ´γ1 n |W n |q and thus |W n | ´1{2 B " o P p1q. Wn qλ θ ‹ pu, Xq du `ÿ uPX WnzpWnaαnq tpu, X Wn zuq

	Second,
		ż
	B " tpu, X Note vPx W c Wnaαn n	Φ θ pv´uq	ˇw
	here in the Third,
	C " tpu, X " ´żWnzpWnaαnq ´ÿ jPJn s 1 ∆ n,j ,8
	|A 2,m | ď	ÿ	ÿ	|g m pv ´uq| ď R n ´γ1	ÿ	Hpu, Xq.
		uPX Wnaαn	vPX WnzBpu,Rnq	uPX Wn

Applying the GNZ formula and Lemma B.2, we deduce that E |A 2,m | " OpR ´γ1 n |W n |q showing that |W n |

  us prove the first convergence by application of Theorem A.1. By [W n ] and the definition of I n , we have |I n | " Op|W n |q, see e.g. Coeurjolly and Møller (2014, Lemma A.1). From (3.13), α

	p2`εqd n

  Bu,n q} |λ θ ‹ pu, xq ´λθ ‹ pu, x Bu,n q| ď c xq " }t|t m |pu, xqu mě1 }Hpu, xqe c Gpu,xq λ θ ‹ pu, xq using the notation of the lemma. Since |∆ n,j | ď 1, we deduce from the stationarity of X and Lemma B.2 that

	α γ 1

B.1, we have for any u P R d and x P Ω }tpu, x n Y pu, xq where Y pu,

  Rn pX; θq ´d dθ dθ J LPL Wnaαn pX; θq

		*
		Ñ 0
	and
	´|W n a α n |	´1

  , xqλ θ ‹ pv, xq ´λθ ‹ ptu, vu, xq " λ θ ‹ pu, xqλ θ ‹ pv, xqt1 ´e´Φ θ ‹ pv´uq u (B.3) and ∆ v t m pu, xq " t m pu, x Y vq ´tm pu, xq " " 0 if m " 1 g m pvq if m ě 2 (B.4) which leads to ∆ u tpv, xq " gpvq. Letting |t|pu, xq " t|t m |pu, xqu mě1 for any u P R d and x P Ω, we have for any u, v P R d

  k tpu, Xqλ θ ‹ pu, Xqtπ n pv, Xq ´πn pv, X Y uqu J du dv.By definition of λ θ ‹ and t (see (2.2) and (3.2)), we have for any u, v P R d and x P Ω π n pv, x Y uq " e ´Φθ ‹ pv´uq rπ n pv, xq `gpv ´uqtλ θ ‹ pv, x Bv,n q ´λθ ‹ pv, xqus, Xqλ θ ‹ pu, Xq 1 ´e´Φ θ ‹ pv´uq ( π n pv, Xq J du dv

	whereby
		ż	ż
	Eps 1 ∆ n,j I J k q " E tpu, ´E ż ∆ n,j ∆ n,k ż
	∆ n,j	∆ n,k

Acknowledgements

The authors would like to thank sincerely the referee for providing many interesting remarks and relevant suggestions. Her/his work allowed us to have a new look on a previous version of this paper and allowed us to improve our results in several directions. The research of J.-F. Coeurjolly is partially funded by Persyval-lab EA Oculo-Nimbus.

the same spirit as [START_REF] Bolthausen | On the central limit theorem for stationary mixing random fields[END_REF] : A n " A n,1 ´An,2 ´An,3 where

A n,2 " e iv J Sn Σ ´1{2 n ÿ jPIn Z n,j p1 ´iv J S n,j ´e´iv J S n,j q

A n,3 " Σ ´1{2 n ÿ jPIn Z n,j e iv J pSn´S n,j q and prove in the following that } E A n,r } Ñ 0 for r " 1, 2, 3 as n Ñ 8. First, assumption (c) implies that |I n | ´1Σ n is a positive definite matrix for n sufficiently large, which is now assumed in the following. By we denote the constant p{λ min pQq where λ min pM q stands for the smallest eigenvalue of a positive definite squared matrix M . For n sufficiently large, λ min p|I n | ´1Σ n q ě λ min pQq ą 0 whereby we deduce

Using this result, Jensen's inequality and the sub-multiplicative property of the Frobenius norm, we get for q ě 1 satisfying (a) and the assumption on

whereby we deduce that } E A n,1 } Ñ 0 from (A.1). Second, since |1 ´e´iz ´iz| ď z 2 {2 for any z P R, we have

Let us decompose B n,j " B p1q n,j `Bp2q n,j where B p1q n,j " B n,j 1p}Z n,j } ď |I n | τ q and B p2q n,j " B n,j 1p}Z n,j } ą |I n | τ q with 1{τ " 2p4p ´1q. By assumption (b), we have

By assumption (a), using Hölder and Bienaymé-Chebyshev inequalities, we continue with

Combining (A.4)-(A.5), we deduce that as n Ñ 8

by definition of τ and α n . Third, for any j P I n , S n ´Sn,j does not depend on X n,j . This yields

whereby we deduce, in view of (A.3), that

which tends to 0 by assumption (d).

B Auxiliary results

We gather in this section several auxiliary results. They are established under the setting, assumptions and notation of Section 3. In particular, we recall that ∆ j is the cube centered at j P Z d with volume 1, ∆ n,j " ∆ j X pW n a α n q, I n Ă Z d is the set such that W n a α n " Y jPIn ∆ n,j , B u,n " Bpu, α n q and for any Then, if α n ě r 0