Infill asymptotics for logistic regression estimators for spatio-temporal
point processes

M.N.M. van Lieshout1,2 and C. Lu2
1Centrum Wiskunde & Informatica
P.O. Box 94079, NL-1090 GB, Amsterdam, The Netherlands
2Department of Applied Mathematics, University of Twente
P.O. Box 217, NL-7500 AE, Enschede, The Netherlands

> Abstract
> This paper discusses infill asymptotics for logistic regression estimators
> for spatio-temporal point processes whose intensity functions are of log-
> linear form. We establish strong consistency and asymptotic normality for
> the parameters of a Poisson point process model and demonstrate how these
> results can be extended to general point process models. Additionally, under
> proper conditions, we also extend our central limit theorem to other
> unbiased estimating equations that are based on the Campbell–Mecke theorem.
>
>
> Keywords & Phrases: Campbell–Mecke theorem, infill asymptotics, logistic
> regression estimator, spatio-temporal point process, unbiased estimating
> equation.
>  2010 Mathematics Subject Classification: 60G55, 62M30.

##  1 Introduction

Spatial and spatio-temporal point processes have been widely used to model,
for example, earthquakes (Bray and Schoenberg, 2013), fires (Møller and Díaz-
Avalos, 2010; Lu et al., 2021) and tropical rain forests (Waagepetersen,
2008).

Usually, the modelling procedure of a spatial or spatio-temporal point pattern
starts with the intensity function which characterizes the probability of a
point occurring in an infinitesimal ball centred at a given location and time.
In many applications, the intensity function is defined as a log-linear
parametric function of certain covariates (e.g., Møller and Waagepetersen,
2007; Waagepetersen, 2008; Guan and Shen, 2010; Coeurjolly and Møller, 2014).
To estimate the parameters, classical options include maximum likelihood
estimation (e.g., Ogata and Tanemura, 1981; Møller and Waagepetersen, 2004),
maximum pseudo-likelihood estimation (e.g., Besag, 1977; Baddeley and Turner,
2000), logistic regression estimation (e.g., Baddeley et al., 2010; Baddeley
et al., 2014) and minimum contrast estimation (e.g., Guyon, 1995).

Among them, maximum likelihood estimation is the most computationally
intensive method, since for many point process models the likelihood function
involves an intractable normalizing constant that must be approximated by
Markov chain Monte Carlo simulations (e.g., Møller and Waagepetersen, 2004).
Instead, maximum pseudo-likelihood estimation and logistic regression
estimation are based on the well-known Campbell–Mecke and Nguyen–Zessin
theorems (see, e.g., Daley and Vere-Jones, 2009) and are easy to implement
using standard software for generalized linear models (e.g., Baddeley et al.,
2015).

From a theoretical perspective, it is important to analyze the asymptotic
properties of the parameter estimators mentioned above. For spatial and
spatio-temporal point processes, two asymptotic regimes can be formulated:
increasing-domain and infill asymptotics (Ripley, 1988). In the former regime,
the observation window grows; whereas in the latter, the window remains fixed
but contains more and more points.

In the literature, Waagepetersen (2008) studied infill asymptotics for maximum
pseudo-likelihood estimators for Poisson and cluster point process models with
log-linear parametric intensity functions using various weighting schemes to
approximate the integral involved. Again, assuming log-linear parametric
intensity functions, Thurman et al. (2015) and Choiruddin et al. (2018)
developed increasing-domain asymptotics for regularized versions of pseudo-
likelihood and logistic regression estimators. In addition, Baddeley et al.
(2014) proved increasing-domain asymptotics for logistic regression estimators
for stationary Gibbs point processes. In this paper, we study consistency and
asymptotic normality for logistic regression estimators for spatio-temporal
point processes in the infill asymptotic regime and extend the central limit
theorem to other unbiased estimating equations that are based on the
Campbell–Mecke theorem.

The remainder of the paper is organized as follows. Section 2 introduces some
background and notation regarding spatio-temporal point processes: parametric
intensity function, logistic regression estimation, computations of the first
two moments of the first-order $U$-statistics and infill asymptotic regime.
Section 3 provides proofs for strong consistency and asymptotic normality for
logistic regression estimators in the case of Poisson point process models.
Section 4 discusses the extension of the asymptotic results to general point
process models and general unbiased estimating equations. Finally, the paper
finishes with some conclusions.

##  2 Background and notation

Let $X$ be a spatio-temporal point process on a bounded non-empty open set
$W\times T\subset\mathbb{R}^{2}\times\mathbb{R}$ equipped with the Borel
$\sigma$-algebra (see, e.g., Daley and Vere-Jones, 2009). Suppose that the
first-order moment measure of $X$, which is defined as

$\Lambda\left(B\right)=\mathbb{E}\left\\{\sum_{x\in X}1\left(x\in
B\right)\right\\}$

for any Borel subset $B\subset W\times T$, exists as a $\sigma$-finite measure
that is absolutely continuous with respect to Lebesgue measure with
Radon–Nikodym derivative $\lambda$. Here, $1(\cdot)$ denotes the indicator
function. Then $\lambda:W\times T\to[0,\infty)$ is called the intensity
function of $X$.

Higher order moment measures are defined analogously. For instance, the
second-order factorial moment measure of $X$ is

$\Lambda^{(2)}\left(B_{1}\times
B_{2}\right)=\mathbb{E}\left\\{\mathop{\sum\nolimits^{\neq}}_{x,y\in
X}1\left(x\in B_{1},y\in B_{2}\right)\right\\}$

for any Borel subsets $B_{1},B_{2}\subset W\times T$. Suppose this measure
also exists as a $\sigma$-finite measure that is absolutely continuous with
respect to Lebesgue measure with Radon–Nikodym derivative $\lambda^{(2)}$.
Then $\lambda^{(2)}:(W\times T)^{2}\to[0,\infty)$ is called the second-order
product density function of $X$. Upon scaling, one obtains the pair
correlation function
$g(x,y)=\lambda^{(2)}\left(x,y\right)/(\lambda\left(x\right)\lambda\left(y\right))$
of $X$, provided that $\lambda\left(x\right)\lambda\left(y\right)>0$.

###  2.1 Parametric intensity function

In this paper, we assume that $X$ admits an intensity function $\lambda$ that
is of log-linear form offset by a measurable function $b$ and parameterized by
a vector $\boldsymbol{\theta}$ in some parameter space
$\boldsymbol{\Theta}\subset\mathbb{R}^{m}$. Specifically,

$\lambda\left(u;\boldsymbol{\theta}\right)=b\left(u\right)\exp\left(\boldsymbol{\theta}^{\top}\bm{z}\left(u\right)\right),$
(1)

where $u=(w,t)\in W\times T$ denotes a location and time combination,
$b:W\times T\to[0,\infty)$ is a measurable function that serves as the
baseline or reference intensity,
$\bm{z}=[z^{(1)},\dots,z^{(m)}]^{\top}:W\times T\to\mathbb{R}^{m}$ is an
$m$-dimensional measurable vector of spatio-temporal covariates and
$\boldsymbol{\theta}=[\theta^{(1)},\dots,\theta^{(m)}]$ is the parameter
vector. The gradient vector of $\lambda\left(u;\boldsymbol{\theta}\right)$
with respect to $\boldsymbol{\theta}$ then takes the form

$\nabla\lambda\left(u;\boldsymbol{\theta}\right)=\lambda\left(u;\boldsymbol{\theta}\right)\bm{z}\left(u\right).$

Conditions must be imposed on $b$ and $\bm{z}$ to ensure that $\lambda$ is
absolutely integrable. In the sequel, it will sometimes also be necessary to
assume that $b$, and thus $\lambda$, is strictly positive. When this is the
case, we will state it explicitly.

###  2.2 Logistic regression estimation

Estimating equations for the parameters of a spatio-temporal point process
model in general and the logistic regression estimator in particular are based
on the Campbell–Mecke theorem (see, e.g., Daley and Vere-Jones, 2009).

Consider a spatio-temporal point process $X$ on $W\times T$ with intensity
function $\lambda$. For any real-valued measurable function $f$ defined on
$W\times T$ such that $f\lambda$ is absolutely integrable, the Campbell–Mecke
theorem reads

$\mathbb{E}\left\\{\sum_{x\in X}f\left(x\right)\right\\}=\int_{W\times
T}f\left(u\right)\lambda\left(u\right)du,$ (2)

where $x$ runs through the points of $X$. When $\lambda$ is parameterized by a
vector $\boldsymbol{\theta}$ as $\lambda\left(u;\boldsymbol{\theta}\right)$,
(2) provides a basis for estimating $\boldsymbol{\theta}$.

The logistic regression estimator is based on the components of the vector
function

$\boldsymbol{f}\left(u;\boldsymbol{\theta}\right)=\nabla\log\left[\frac{\lambda\left(u;\boldsymbol{\theta}\right)}{\lambda\left(u;\boldsymbol{\theta}\right)+\rho\left(u\right)}\right]=\frac{\rho\left(u\right)/\lambda\left(u;\boldsymbol{\theta}\right)}{\lambda\left(u;\boldsymbol{\theta}\right)+\rho\left(u\right)}\nabla\lambda\left(u;\boldsymbol{\theta}\right).$
(3)

Here, one assumes that $\lambda(u;\boldsymbol{\theta})$ is a positive-valued
differentiable function such that its gradient vector
$\nabla\lambda(u;\boldsymbol{\theta})$ is absolutely integrable, and that
$\rho(u)$ is a positive-valued measurable function also defined on $W\times
T$. The idea is then to estimate both sides of (2) and solve the equations for
$\boldsymbol{\theta}$. In order to approximate the right-hand side, one may
use a ‘dummy’ point process $D$ on $W\times T$ which is independent of $X$ and
has absolutely integrable intensity function $\rho$. Applying the
Campbell–Mecke theorem to $D$, one finds that

$\sum_{x\in
D}\frac{1}{\lambda(x;\boldsymbol{\theta})+\rho(x)}\nabla\lambda(x;\boldsymbol{\theta})$

is an unbiased estimator for the right-hand side of (2) with $\boldsymbol{f}$
as in (3). Hence,

${\boldsymbol{s}}(X,D;\boldsymbol{\theta})=\sum_{x\in
X}\frac{\rho(x)/\lambda(x;\boldsymbol{\theta})}{\lambda(x;\boldsymbol{\theta})+\rho(x)}\nabla\lambda(x;\boldsymbol{\theta})-\sum_{x\in
D}\frac{1}{\lambda(x;\boldsymbol{\theta})+\rho(x)}\nabla\lambda(x;\boldsymbol{\theta})=\boldsymbol{0}$
(4)

is an unbiased estimating equation. It is interesting to observe that the
middle part of (4) is exactly the gradient of a logistic log-likelihood
function

$l(X,D;\boldsymbol{\theta})=\sum_{x\in
X}\log\left[\frac{\lambda(x;\boldsymbol{\theta})}{\lambda(x;\boldsymbol{\theta})+\rho(x)}\right]+\sum_{x\in
D}\log\left[\frac{\rho(x)}{\lambda(x;\boldsymbol{\theta})+\rho(x)}\right].$
(5)

Thus, (4) can also be interpreted as a score function and solved using
standard software such as the R-package stats (Venables and Ripley, 2002). The
existence and uniqueness of a maximizer of the log-likelihood function (5) are
ensured under proper conditions (Silvapulle, 1981).

For simplicity of writing, in the remainder of this paper, we write
$\boldsymbol{s}(\boldsymbol{\theta})$ and $l(\boldsymbol{\theta})$ for
$\boldsymbol{s}(X,D;\boldsymbol{\theta})$ and $l(X,D;\boldsymbol{\theta})$ to
suppress the dependence on the point patterns $X$ and $D$. Moreover, we use
$\boldsymbol{\theta_{0}}$ and $\widehat{\boldsymbol{\theta}}$ to respectively
denote the true value and the estimator of $\boldsymbol{\theta}$.

###  2.3 First two moments of first-order $U$-statistics

Following Reitzner and Schulte (2013), we call random vectors of the form

$\bm{H}=\sum_{x\in X}\left[h^{(1)}(x),\dots,h^{(m)}(x)\right]^{\top}$

first-order $U$-statistics of $X$ if $h^{(l)}\lambda$ is absolutely integrable
for every $1\leq l\leq m$.

Suppose that $X$ has intensity function $\lambda^{(1)}=\lambda$ and second-
order product density function $\lambda^{(2)}$. By the Campbell–Mecke theorem
and in analogy to (2),

$\mathbb{E}\left\\{\mathop{\sum\nolimits^{\neq}}_{x_{1},\dots,x_{r}\in
X}f\left(x_{1},\dots,x_{r}\right)\right\\}=\int_{(W\times
T)^{r}}f\left(u_{1},\dots,u_{r}\right)\lambda^{(r)}\left(u_{1},\dots,u_{r}\right)du_{1}\dots
du_{r}$

for $r=1,2$ and any real-valued measurable function $f$ defined on $W\times T$
such that $f\lambda$ and $f\lambda^{(2)}$ are absolutely integrable. Then, it
can be used to derive the first two moments of $\bm{H}$. Indeed,

$\mathbb{E}\left\\{\bm{H}\right\\}=\left[\int_{W\times
T}h^{(l)}(u)\lambda\left(u\right)du\right]_{l=1}^{m}$

and

$\mathbb{E}\left\\{\bm{H}\bm{H}^{\top}\right\\}=\left[\int_{W\times
T}h^{(k)}\left(u\right)h^{(l)}\left(u\right)\lambda\left(u\right)du+\int_{W\times
T}\int_{W\times
T}h^{(k)}(u)h^{(l)}(v)\lambda^{(2)}\left(u,v\right)dudv\right]_{k,l=1}^{m}.$

If $\mathbb{E}\left\\{\bm{H}\bm{H}^{\top}\right\\}$ is finite and $\lambda>0$,
the entries in the covariance matrix of $\bm{H}$ are finite and can be
expressed in terms of the pair correlation function $g(u,v)$ as

$\begin{split}\operatorname{\bf{Cov}}\left\\{\bm{H}\right\\}&=\left[\int_{W\times
T}h^{(k)}\left(u\right)h^{(l)}\left(u\right)\lambda\left(u\right)du\right]_{k,l=1}^{m}\\\
&+\left[\int_{W\times T}\int_{W\times
T}h^{(k)}(u)h^{(l)}(v)(g(u,v)-1)\lambda\left(u\right)\lambda\left(v\right)dudv\right]_{k,l=1}^{m}.\end{split}$

###  2.4 Infill asymptotic regime

The infill asymptotic regime considered in this paper is as defined in
Waagepetersen (2008) and Lieshout (2020).

Let $\\{Y_{i}\\}$ and $\\{E_{i}\\}$, with $i\in\mathbb{N}^{+}$, be two
sequences of independent and identically distributed spatio-temporal point
processes with intensity functions $\lambda$ and $\rho$, respectively. Set

$X_{n}=\bigcup_{i=1}^{n}Y_{i},\quad D_{n}=\bigcup_{i=1}^{n}E_{i}.$

Write $\lambda_{n}$ for the intensity function of $X_{n}$ and $\rho_{n}$ for
that of the ‘dummy’ point process $D_{n}$. Thus, $\lambda_{n}=n\lambda$ and
$\rho_{n}=n\rho$, and we assume that the intensity functions of $X_{n}$ and
$D_{n}$ increase at the same rate. For any $n\in\mathbb{N}^{+}$, the score
function (4) based on $X_{n}$ and $D_{n}$ then becomes

$\begin{split}\boldsymbol{s}_{n}(\boldsymbol{\theta})&=\sum_{x\in
X_{n}}\frac{\rho_{n}(x)/\lambda_{n}(x;\boldsymbol{\theta})}{\lambda_{n}(x;\boldsymbol{\theta})+\rho_{n}(x)}\nabla\lambda_{n}(x;\boldsymbol{\theta})-\sum_{x\in
D_{n}}\frac{1}{\lambda_{n}(x;\boldsymbol{\theta})+\rho_{n}(x)}\nabla\lambda_{n}(x;\boldsymbol{\theta})\\\
&=\sum_{x\in
X_{n}}\frac{\rho(x)/\lambda(x;\boldsymbol{\theta})}{\lambda(x;\boldsymbol{\theta})+\rho(x)}\nabla\lambda(x;\boldsymbol{\theta})-\sum_{x\in
D_{n}}\frac{1}{\lambda(x;\boldsymbol{\theta})+\rho(x)}\nabla\lambda(x;\boldsymbol{\theta})=\boldsymbol{0}.\end{split}$
(6)

Note that the terms in the two sums above do not depend on $n$ after
simplification while the subscripts $x$ run through the points of $X_{n}$ and
$D_{n}$. Taking the limit as $n\to\infty$, one obtains an asymptotic regime
that Ripley (1988) calls ‘infill asymptotics’.

Intuitively, under this regime, the estimate for the parameter vector
$\boldsymbol{\theta}$ will become more precise when the points observed in the
fixed domain $W\times T$ become more dense.

Our aim in the remainder of this paper is to analyze the asymptotic behaviour
of the estimator $\widehat{\boldsymbol{\theta}_{n}}$ defined by (6) as
$n\rightarrow\infty$.

##  3 Infill asymptotics for Poisson point processes

For the sake of completeness, we first recall the definition of Poisson point
processes (see, e.g., Lieshout, 2000; Daley and Vere-Jones, 2009).

######  Definition 3.1.

A spatio-temporal point process $X$ defined as in the beginning of Section 2
is a Poisson point process if it satisfies the following properties:

  * (i)

for any bounded Borel set $B\subset W\times T$, the number of points that fall
in $B$ is Poisson distributed with mean $\Lambda(B)$;

  * (ii)

in disjoint bounded Borel sets $B_{1},B_{2}\subset W\times T$, the numbers of
points that fall in $B_{1},B_{2}$ are independent.

The probability distribution of a Poisson point process $X$ is completely
specified by its first-order moment measure $\Lambda$. Thus, the assumption of
a log-linear intensity function parameterized as (1) is quite natural and
constitutes an exponential family with the components of $\sum_{x\in
X}\bm{z}(x)$ as sufficient statistics. Moreover, the pair correlation function
of $X$ is always equal to one. The infill asymptotic regime (cf., Section 2.4)
is then particularly appropriate for a Poisson point process model because of
the conditional independence of the points given a fixed number.

In this section, we establish strong consistency and asymptotic normality for
logistic regression estimators in the case of spatio-temporal Poisson point
processes. For ease of referencing, we list here the conditions
($\mathcal{C}.1$)–($\mathcal{C}.8$) that are required to derive the asymptotic
results:

$(\mathcal{C}.1)$

$\\{Y_{i}\\}$ and $\\{E_{i}\\}$ with $i\in\mathbb{N}^{+}$ are two independent
sequences of independent and identically distributed spatio-temporal point
processes on some bounded open set $W\times
T\subset\mathbb{R}^{2}\times\mathbb{R}$ and defined on the same underlying
probability space $\Omega$. Set $X_{n}=\cup_{i=1}^{n}Y_{i}$ and
$D_{n}=\cup_{i=1}^{n}E_{i}$.

$(\mathcal{C}.2)$

$Y_{i}$ is a Poisson point process with intensity function
$\lambda(u;\boldsymbol{\theta})$ as (1), where $b>0$ is an absolutely
integrable function, $\bm{z}$ is a measurable vector of covariates and the
parameter vector $\boldsymbol{\theta}$ lies in an open set
$\boldsymbol{\Theta}\subset\mathbb{R}^{m}$; $E_{i}$ has absolutely integrable
intensity function $\rho(u)>0$.

$(\mathcal{C}.3)$

$E_{i}$ has bounded pair correlation function $g$, that is,
$\sup_{(u,v)\in(W\times T)^{2}}g(u,v)$ $<\infty$.

$(\mathcal{C}.4)$

For every $\boldsymbol{\theta}\in\boldsymbol{\Theta}$, there exist
$\epsilon_{1}(\boldsymbol{\theta}),\epsilon_{2}(\boldsymbol{\theta})>0$ such
that $\epsilon_{1}(\boldsymbol{\theta})<\inf_{u\in W\times
T}\rho(u)/\lambda(u;\boldsymbol{\theta})$ and $\sup_{u\in W\times
T}\rho(u)/\lambda(u;\boldsymbol{\theta})<\epsilon_{2}(\boldsymbol{\theta})$.

$(\mathcal{C}.5)$

The elements of the measurable covariate vector $\bm{z}$ are bounded, that is,
$\sup_{u\in W\times T}\|\bm{z}(u)\|$ $<\infty$.

$(\mathcal{C}.6)$

The parameter space $\boldsymbol{\Theta}$ is convex.

$(\mathcal{C}.7)$

The parametric model for $\lambda$ is identifiable, that is,
$\lambda(u;\boldsymbol{\theta})=\lambda(u;\boldsymbol{\tilde{\theta}})$ almost
everywhere on $W\times T$ implies
$\boldsymbol{\theta}=\boldsymbol{\tilde{\theta}}$.

$(\mathcal{C}.8)$

The matrix $\bm{U}$, whose ($k,l$)-th entry reads $\int_{W\times
T}\lambda(u;\boldsymbol{\theta_{0}})\rho(u)z^{(k)}(u)z^{(l)}(u)/(\lambda(u;\boldsymbol{\theta_{0}})+\rho(u))du$
with $\boldsymbol{\theta_{0}}\in\boldsymbol{\Theta}$ and $1\leq k,l\leq m$, is
positive definite.

A few remarks on some of the conditions are appropriate. In condition
($\mathcal{C}.2$), $Y_{i}$ is assumed to be a Poisson point process, however,
it may be relaxed (cf., Section 4.1). In condition ($\mathcal{C}.3$), $E_{i}$
is not required to be a Poisson point process, as its realisations are only
used to approximate an integral. Condition ($\mathcal{C}.4$) is reasonable,
recalling the rule of thumb recommended by Baddeley et al. (2014) for
selecting $\rho(u)$ that $\rho=4\lambda$. Condition ($\mathcal{C}.7$) is
necessary for strong consistency (cf., Theorem 3.3) and condition
($\mathcal{C}.8$) helps ensure the attainment of the logistic regression
estimator $\widehat{\boldsymbol{\theta}_{n}}$ as $n\to\infty$ (cf., Theorem
3.11).

In the remainder of this paper, we will use $P_{\boldsymbol{\theta_{0}}}$ to
denote the distribution of $(X_{n},D_{n})$ under the true parameter value
$\boldsymbol{\theta_{0}}$.

###  3.1 Strong consistency

To establish strong consistency, we start our investigations with the
asymptotic behaviour of the scaled log-likelihood function.

######  Lemma 3.2.

Assume that the conditions ($\mathcal{C}.1$)–($\mathcal{C}.2$) and
($\mathcal{C}.4$)–($\mathcal{C}.5$) hold. Define
$l_{n}(\boldsymbol{\theta})=l(X_{n},D_{n};\boldsymbol{\theta})$ by (5) with
$\boldsymbol{\theta}\in\boldsymbol{\Theta}$. Then, as $n\to\infty$,
$l_{n}(\boldsymbol{\theta})/n$ converges $P_{\boldsymbol{\theta_{0}}}$-almost
surely to

$\int_{W\times
T}\left\\{\lambda(u;\boldsymbol{\theta_{0}})\log\left[\frac{\lambda(u;\boldsymbol{\theta})}{\lambda(u;\boldsymbol{\theta})+\rho(u)}\right]+\rho(u)\log\left[\frac{\rho(u)}{\lambda(u;\boldsymbol{\theta})+\rho(u)}\right]\right\\}du.$

###### Proof.

Under conditions ($\mathcal{C}.1$)–($\mathcal{C}.2$) and recalling the
logistic log-likelihood function (5),

$\frac{l_{n}(\boldsymbol{\theta})}{n}=\frac{1}{n}\sum_{x\in
X_{n}}\log\left[\frac{\lambda(x;\boldsymbol{\theta})}{\lambda(x;\boldsymbol{\theta})+\rho(x)}\right]+\frac{1}{n}\sum_{x\in
D_{n}}\log\left[\frac{\rho(x)}{\lambda(x;\boldsymbol{\theta})+\rho(x)}\right],$

which consists of two first-order $U$-statistics defined on $X_{n}$ and
$D_{n}$. We shall derive their strong convergence separately.

Write the first-order $U$-statistic defined on $X_{n}$ as

$\frac{1}{n}\sum_{x\in
X_{n}}\log\left[\frac{\lambda(x;\boldsymbol{\theta})}{\lambda(x;\boldsymbol{\theta})+\rho(x)}\right]=\frac{1}{n}\sum_{i=1}^{n}\left\\{\sum_{x\in
Y_{i}}\log\left[\frac{\lambda(x;\boldsymbol{\theta})}{\lambda(x;\boldsymbol{\theta})+\rho(x)}\right]\right\\}.$

The sum is the average of independent and identically distributed real-valued
random variables. By the Campbell–Mecke theorem,

$\mathbb{E}_{\boldsymbol{\theta_{0}}}\left\\{\sum_{x\in
Y_{i}}\log\left[\frac{\lambda(x;\boldsymbol{\theta})}{\lambda(x;\boldsymbol{\theta})+\rho(x)}\right]\right\\}=\int_{W\times
T}\log\left[\frac{\lambda(u;\boldsymbol{\theta})}{\lambda(u;\boldsymbol{\theta})+\rho(u)}\right]\lambda(u;\boldsymbol{\theta_{0}})\,du.$
(7)

Conditions ($\mathcal{C}.2$) and ($\mathcal{C}.5$) imply that the intensity
function $\lambda(u;\boldsymbol{\theta}_{0})$ is absolutely integrable on
$W\times T$. Furthermore, by conditions ($\mathcal{C}.2$) and
($\mathcal{C}.4$),

$\log\left[\frac{1}{1+\epsilon_{2}(\boldsymbol{\theta})}\right]<\log\left[\frac{1}{1+\rho(u)/\lambda(u;\boldsymbol{\theta})}\right]=\log\left[\frac{\lambda(u;\boldsymbol{\theta})}{\lambda(u;\boldsymbol{\theta})+\rho(u)}\right]<\log\left[\frac{1}{1+\epsilon_{1}(\boldsymbol{\theta})}\right].$

Thus, the $P_{\boldsymbol{\theta_{0}}}$-mean in (7) is finite for all
$\boldsymbol{\theta}$. Then, Kolmogorov’s strong law of large numbers implies
that the first-order $U$-statistic defined on $X_{n}$ converges
$P_{\boldsymbol{\theta_{0}}}$-almost surely to the integral in the right-hand
side of equation (7).

Similarly, for the first-order $U$-statistic defined on the dummy point
process $D_{n}$, as $n\to\infty$,

$\begin{split}\frac{1}{n}\sum_{x\in
D_{n}}\log\left[\frac{\rho(x)}{\lambda(x;\boldsymbol{\theta})+\rho(x)}\right]&\stackrel{{\scriptstyle
P_{\boldsymbol{\theta_{0}}}-a.s.}}{{\to}}\mathbb{E}_{\boldsymbol{\theta_{0}}}\left\\{\sum_{x\in
E_{i}}\log\left[\frac{\rho(x)}{\lambda(x;\boldsymbol{\theta})+\rho(x)}\right]\right\\}\\\
&=\int_{W\times
T}\log\left[\frac{\rho(u)}{\lambda(u;\boldsymbol{\theta})+\rho(u)}\right]\rho(u)du\end{split}$

because $\rho(u)$ is absolutely integrable by condition ($\mathcal{C}.2$) and

$\log\left[\frac{\epsilon_{1}(\boldsymbol{\theta})}{1+\epsilon_{1}(\boldsymbol{\theta})}\right]<\log\left[\frac{\rho(u)/\lambda(u;\boldsymbol{\theta})}{1+\rho(u)/\lambda(u;\boldsymbol{\theta})}\right]=\log\left[\frac{\rho(u)}{\lambda(u;\boldsymbol{\theta})+\rho(u)}\right]<\log\left[\frac{\epsilon_{2}(\boldsymbol{\theta})}{1+\epsilon_{2}(\boldsymbol{\theta})}\right]$

by conditions ($\mathcal{C}.2$) and ($\mathcal{C}.4$).

The proof is then completed by combining these two strong convergence results.
∎

Recall that the logistic regression estimator minimizes
$U_{n}(\boldsymbol{\theta})=-l_{n}(\boldsymbol{\theta})/{n}$ and thus can be
considered as a minimum contrast estimator. The next theorem is then concerned
with strong consistency.

######  Theorem 3.3.

Assume that the conditions ($\mathcal{C}.1$)–($\mathcal{C}.2$) and
($\mathcal{C}.4$)–($\mathcal{C}.7$) hold. Define
$l_{n}(\boldsymbol{\theta})=l(X_{n},D_{n};\boldsymbol{\theta})$ by (5) with
$\boldsymbol{\theta}\in\boldsymbol{\Theta}$ and set
$\widehat{\boldsymbol{\theta}_{n}}=\operatorname*{arg\,max}_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}l_{n}(\boldsymbol{\theta})$.
If $\widehat{\boldsymbol{\theta}_{n}}$ is attained, then as
$n\rightarrow\infty$, $\widehat{\boldsymbol{\theta}_{n}}$ converges
$P_{\boldsymbol{\theta_{0}}}$-almost surely to $\boldsymbol{\theta_{0}}$.

###### Proof.

Conditions ($\mathcal{C}.2$) and ($\mathcal{C}.6$) ensure that the parameter
space $\boldsymbol{\Theta}$ is open and convex.

Firstly, we prove that, for every $\omega$ in the underlying probability space
$\Omega$, the realisations of the function $\boldsymbol{\theta}\mapsto
U_{n}(\boldsymbol{\theta})$ are convex. By condition ($\mathcal{C}.2$), the
intensity function $\lambda(x;\boldsymbol{\theta})$ has the log-linear form
(1) and is thus twice differentiable with respect to the parameter vector
$\boldsymbol{\theta}\in\boldsymbol{\Theta}$. Then, the Hessian matrix of
$-l_{n}(\boldsymbol{\theta})$ reads

$\left[\sum_{x\in X_{n}\cup
D_{n}}\frac{\lambda(x;\boldsymbol{\theta})\rho(x)}{\left(\lambda(x;\boldsymbol{\theta})+\rho(x)\right)^{2}}z^{(k)}(x)z^{(l)}(x)\right]_{k,l=1}^{m}.$

By decomposition, it can be written into the product of a matrix $\bm{M}$ and
its transpose as $\bm{M}^{\top}\bm{M}$, where $\bm{M}$ is a [$|X_{n}\cup
D_{n}|\times m$]-matrix given by

$\bm{M}=\left[\frac{(\lambda(x_{i};\boldsymbol{\theta})\rho(x_{i}))^{1/2}}{\lambda(x_{i};\boldsymbol{\theta})+\rho(x_{i})}z^{(k)}(x_{i})\right]_{i=1,k=1}^{\left|X_{n}\cup
D_{n}\right|,m}.$

Here, $|X_{n}\cup D_{n}|$ denotes the number of points in $X_{n}\cup D_{n}$,
the subscript $i$ runs through all points in $X_{n}\cup D_{n}$ and $k$ runs
through the $m$ covariates. One can readily obtain that the Hessian matrix of
$-l_{n}(\boldsymbol{\theta})$ is positive semi-definite for every
$\boldsymbol{\theta}\in\boldsymbol{\Theta}$, which implies that the function
$\boldsymbol{\theta}\mapsto-l_{n}(\boldsymbol{\theta})$, and thus
$-l_{n}(\boldsymbol{\theta})/n$, is convex.

Secondly, we prove that, as $n\rightarrow\infty$,
$U_{n}(\boldsymbol{\theta})-U_{n}(\boldsymbol{\theta_{0}})$ converges
$P_{\boldsymbol{\theta_{0}}}$-almost surely to a function
$K(\boldsymbol{\theta},\boldsymbol{\theta_{0}})$ that is non-negative and
vanishes only at $\boldsymbol{\theta}=\boldsymbol{\theta_{0}}$. By Lemma 3.2,
under conditions ($\mathcal{C}.1$)–($\mathcal{C}.2$) and
($\mathcal{C}.4$)–($\mathcal{C}.5$),
$U_{n}(\boldsymbol{\theta})-U_{n}(\boldsymbol{\theta_{0}})$ converges
$P_{\boldsymbol{\theta_{0}}}$-almost surely to

$\int_{W\times
T}\left\\{\lambda(u;\boldsymbol{\theta_{0}})\log\left[\frac{\lambda(u;\boldsymbol{\theta_{0}})}{\lambda(u;\boldsymbol{\theta})}\right]-\left(\lambda(u;\boldsymbol{\theta_{0}})+\rho(u)\right)\log\left[\frac{\lambda(u;\boldsymbol{\theta_{0}})+\rho(u)}{\lambda(u;\boldsymbol{\theta})+\rho(u)}\right]\right\\}du.$

Denote this limit by $K(\boldsymbol{\theta},\boldsymbol{\theta_{0}})$ and the
integrand by $k(u;\boldsymbol{\theta},\boldsymbol{\theta_{0}})$. Then, clearly
$K(\boldsymbol{\theta_{0}},\boldsymbol{\theta_{0}})=0$. Furthermore,

$k(u;\boldsymbol{\theta},\boldsymbol{\theta_{0}})=\lambda(u;\boldsymbol{\theta})\left\\{\frac{\lambda(u;\boldsymbol{\theta_{0}})}{\lambda(u;\boldsymbol{\theta})}\log\left[\frac{\lambda(u;\boldsymbol{\theta_{0}})}{\lambda(u;\boldsymbol{\theta})}\right]-\left(\frac{\lambda(u;\boldsymbol{\theta_{0}})}{\lambda(u;\boldsymbol{\theta})}+\frac{\rho(u)}{\lambda(u;\boldsymbol{\theta})}\right)\log\left[\frac{\frac{\lambda(u;\boldsymbol{\theta_{0}})}{\lambda(u;\boldsymbol{\theta})}+\frac{\rho(u)}{\lambda(u;\boldsymbol{\theta})}}{1+\frac{\rho(u)}{\lambda(u;\boldsymbol{\theta})}}\right]\right\\}.$

By condition ($\mathcal{C}.2$), $\lambda(u;\boldsymbol{\theta})>0$ as
$b(u)>0$. Consider the function $a\mapsto a\log a-(a+b)\log[(a+b)/(1+b)]$ with
$a,b>0$. Its derivative with respect to $a$ is $\log[a(1+b)/(a+b)]$. Hence,
the function is strictly decreasing when $a\in(0,1)$ and strictly increasing
when $a\in(1,+\infty)$. For $a=1$, $a\log a-(a+b)\log[(a+b)/(1+b)]=0$. Thus,
$k(u;\boldsymbol{\theta},\boldsymbol{\theta_{0}})$ is non-negative, and is
strictly positive when
$\lambda(u;\boldsymbol{\theta})\neq\lambda(u;\boldsymbol{\theta_{0}})$.

Under condition ($\mathcal{C}.7$), strong consistency follows from an appeal
to the Proposition below Guyon (1995, Theorem 3.4.4). ∎

######  Remark 3.4.

It is important to note that the proofs of Lemma 3.2 and Theorem 3.3 do not
depend on the assumption that $Y_{i}$, and thus $X_{n}$, is a Poisson point
process. However, strong consistency relies on the concavity of the log-
likelihood function of logistic regression estimation.

###  3.2 Asymptotic normality

To establish asymptotic normality, we start our investigations with the Taylor
series of the score function (6).

For every component of the score function, denoted by
$\boldsymbol{s}_{n}^{(i)}(\boldsymbol{\theta})$ with $1\leq i\leq m$, the
second-order Taylor expansion of
$\boldsymbol{s}_{n}^{(i)}(\widehat{\boldsymbol{\theta}_{n}})$ with respect to
$\boldsymbol{\theta}$ at $\boldsymbol{\theta_{0}}$ reads

$0=\boldsymbol{s}_{n}^{(i)}(\widehat{\boldsymbol{\theta}_{n}})=\boldsymbol{s}_{n}^{(i)}(\boldsymbol{\theta_{0}})+\nabla\boldsymbol{s}_{n}^{(i)}(\boldsymbol{\theta_{0}})(\widehat{\boldsymbol{\theta}_{n}}-\boldsymbol{\theta_{0}})+\frac{1}{2}(\widehat{\boldsymbol{\theta}_{n}}-\boldsymbol{\theta_{0}})^{\top}\nabla^{2}\boldsymbol{s}_{n}^{(i)}(\boldsymbol{\theta^{\prime,i}})(\widehat{\boldsymbol{\theta}_{n}}-\boldsymbol{\theta_{0}}),$
(8)

where $\nabla\boldsymbol{s}_{n}^{(i)}(\boldsymbol{\theta})$ is the $[1\times
m]$-vector containing the first-order partial derivatives of
$\boldsymbol{s}_{n}^{(i)}(\boldsymbol{\theta})$ with respect to
$\boldsymbol{\theta}$ and
$\nabla^{2}\boldsymbol{s}_{n}^{(i)}(\boldsymbol{\theta})$ is the $[m\times
m]$-matrix containing the second-order partial derivatives of
$\boldsymbol{s}_{n}^{(i)}(\boldsymbol{\theta})$. Moreover,
$\boldsymbol{\theta^{\prime,i}}$ is a convex combination of
$\widehat{\boldsymbol{\theta}_{n}}$ and $\boldsymbol{\theta_{0}}$ which, by
condition ($\mathcal{C}.6$), lies in $\boldsymbol{\Theta}$ as well.

Write $\nabla\boldsymbol{s}_{n}(\boldsymbol{\theta_{0}})$ for the matrix whose
$i$-th row is $\boldsymbol{s}_{n}^{(i)}(\boldsymbol{\theta_{0}})$ with $1\leq
i\leq m$ and assume that its inverse is well-defined. Heuristically, the idea
is then to ignore the error term which is the quadratic form in (8) and
rearrange the remaining terms to obtain

$n^{1/2}(\widehat{\boldsymbol{\theta}_{n}}-\boldsymbol{\theta_{0}})\approx\left[-\frac{\nabla\boldsymbol{s}_{n}(\boldsymbol{\theta_{0}})}{n}\right]^{-1}\frac{\boldsymbol{s}_{n}(\boldsymbol{\theta_{0}})}{n^{1/2}}.$
(9)

In the next two lemmas, we first study the asymptotic behaviours of the terms
$-\nabla\boldsymbol{s}_{n}(\boldsymbol{\theta_{0}})/n$ and
$\boldsymbol{s}_{n}(\boldsymbol{\theta_{0}})/n^{1/2}$, respectively.

######  Lemma 3.5.

Assume that the conditions ($\mathcal{C}.1$)–($\mathcal{C}.2$) and
($\mathcal{C}.5$) hold. Define $\boldsymbol{s}_{n}(\boldsymbol{\theta})$ $=$
$\boldsymbol{s}(X_{n},D_{n};$ $\boldsymbol{\theta})$ by (6) with
$\boldsymbol{\theta}\in\boldsymbol{\Theta}$. Then, as $n\rightarrow\infty$,
$-\nabla\boldsymbol{s}_{n}(\boldsymbol{\theta_{0}})/n$ converges
$P_{\boldsymbol{\theta_{0}}}$-almost surely to the matrix $\bm{U}$ given by

$\bm{U}=\left[\int_{W\times
T}\frac{\lambda(u;\boldsymbol{\theta_{0}})\rho(u)}{\lambda(u;\boldsymbol{\theta_{0}})+\rho(u)}z^{(k)}(u)z^{(l)}(u)du\right]_{k,l=1}^{m}.$

###### Proof.

Under conditions ($\mathcal{C}.1$)–($\mathcal{C}.2$) and recalling the score
function (6),

$-\frac{\nabla\boldsymbol{s}_{n}(\boldsymbol{\theta_{0}})}{n}=\left[\sum_{x\in
X_{n}\cup
D_{n}}\frac{\lambda(x;\boldsymbol{\theta_{0}})\rho(x)}{n\left(\lambda(x;\boldsymbol{\theta_{0}})+\rho(x)\right)^{2}}z^{(k)}(x)z^{(l)}(x)\right]_{k,l=1}^{m}.$

We shall prove component-wise strong convergence.

Consider the $(k,l)$-th entry of the matrix above which, recalling condition
($\mathcal{C}.1$), is given by

$\frac{1}{n}\sum_{i=1}^{n}\left\\{\sum_{x\in Y_{i}\cup
E_{i}}\frac{\lambda(x;\boldsymbol{\theta_{0}})\rho(x)}{\left(\lambda(x;\boldsymbol{\theta_{0}})+\rho(x)\right)^{2}}z^{(k)}(x)z^{(l)}(x)\right\\}.$

The sum is the average of independent and identically distributed real-valued
random variables. By the Campbell–Mecke theorem,

$\mathbb{E}_{\boldsymbol{\theta_{0}}}\left\\{\sum_{x\in Y_{i}\cup
E_{i}}\frac{\lambda(x;\boldsymbol{\theta_{0}})\rho(x)}{\left(\lambda(x;\boldsymbol{\theta_{0}})+\rho(x)\right)^{2}}z^{(k)}(x)z^{(l)}(x)\right\\}=\int_{W\times
T}\frac{\lambda(u;\boldsymbol{\theta_{0}})\rho(u)}{\lambda(u;\boldsymbol{\theta_{0}})+\rho(u)}z^{(k)}(u)z^{(l)}(u)du.$

By condition ($\mathcal{C}.5$), the covariate terms in the integrand above are
all bounded. Furthermore, by condition ($\mathcal{C}.2$),

$0<\frac{\lambda(u;\boldsymbol{\theta_{0}})\rho(u)}{\lambda(u;\boldsymbol{\theta_{0}})+\rho(u)}<\rho(u)$

and $\rho(u)$ is absolutely integrable. Kolmogorov’s strong law of large
numbers implies the claimed $P_{\boldsymbol{\theta_{0}}}$-almost sure
convergence. ∎

######  Lemma 3.6.

Assume that the conditions ($\mathcal{C}.1$)–($\mathcal{C}.3$) and
($\mathcal{C}.5$) hold. Define
$\boldsymbol{s}_{n}(\boldsymbol{\theta})=\boldsymbol{s}(X_{n},D_{n};$
$\boldsymbol{\theta})$ by (6) with
$\boldsymbol{\theta}\in\boldsymbol{\Theta}$. Then, as $n\rightarrow\infty$,
$\boldsymbol{s}_{n}(\boldsymbol{\theta_{0}})/n^{1/2}$ converges under
$P_{\boldsymbol{\theta_{0}}}$ in distribution to an $m$-dimensional normally
distributed random vector with mean zero and covariance matrix

$\begin{split}\bm{V}&=\left[\int_{W\times
T}\frac{\lambda(u;\boldsymbol{\theta_{0}})\rho(u)}{\lambda(u;\boldsymbol{\theta_{0}})+\rho(u)}z^{(k)}(u)z^{(l)}(u)du\right]_{k,l=1}^{m}\\\
&+\left[\int_{W\times T}\int_{W\times
T}\frac{\lambda(u;\boldsymbol{\theta_{0}})\lambda(v;\boldsymbol{\theta_{0}})\rho(u)\rho(v)z^{(k)}(u)z^{(l)}(v)}{(\lambda(u;\boldsymbol{\theta_{0}})+\rho(u))(\lambda(v;\boldsymbol{\theta_{0}})+\rho(v))}\left(g(u,v)-1\right)dudv\right]_{k,l=1}^{m}.\end{split}$

###### Proof.

Under conditions ($\mathcal{C}.1$)–($\mathcal{C}.2$) and recalling the score
function (6),

$\frac{\boldsymbol{s}_{n}(\boldsymbol{\theta_{0}})}{n^{1/2}}=\sum_{x\in
X_{n}}\frac{\rho(x)}{n^{1/2}\left(\lambda(x;\boldsymbol{\theta_{0}})+\rho(x)\right)}\boldsymbol{z}(x)-\sum_{x\in
D_{n}}\frac{\lambda(x;\boldsymbol{\theta_{0}})}{n^{1/2}\left(\lambda(x;\boldsymbol{\theta_{0}})+\rho(x)\right)}\boldsymbol{z}(x).$

It consists of two first-order $U$-statistics defined on $X_{n}$ and $D_{n}$.
We write it as

$\begin{split}\frac{\boldsymbol{s}_{n}(\boldsymbol{\theta_{0}})}{n^{1/2}}&=\frac{\boldsymbol{s}_{n}(X_{n};\boldsymbol{\theta_{0}})}{n^{1/2}}-\frac{\boldsymbol{s}_{n}(D_{n};\boldsymbol{\theta_{0}})}{n^{1/2}}\\\
&=n^{1/2}\left\\{\frac{1}{n}\sum_{i=1}^{n}\sum_{x\in
Y_{i}}\frac{\rho(x)}{\lambda(x;\boldsymbol{\theta_{0}})+\rho(x)}\boldsymbol{z}(x)-\frac{1}{n}\sum_{i=1}^{n}\sum_{x\in
E_{i}}\frac{\lambda(x;\boldsymbol{\theta_{0}})}{\lambda(x;\boldsymbol{\theta_{0}})+\rho(x)}\boldsymbol{z}(x)\right\\}\end{split}$

and discuss the two first-order $U$-statistics separately.

Consider $n^{-1/2}\boldsymbol{s}_{n}(X_{n};\boldsymbol{\theta_{0}})$. Note
that the first term between the curly brackets above is again the average of
independent and identically distributed real-valued random vectors. By the
Campbell–Mecke theorem and recalling Section 2.3,

$\text{\bf{M}}_{\boldsymbol{\theta_{0}}}=\mathbb{E}_{\boldsymbol{\theta_{0}}}\left\\{\sum_{x\in
Y_{i}}\frac{\rho(x)\boldsymbol{z}(x)}{\lambda(x;\boldsymbol{\theta_{0}})+\rho(x)}\right\\}=\left[\int_{W\times
T}\frac{\rho(u)\lambda(u;\boldsymbol{\theta_{0}})z^{(l)}(u)}{\lambda(u;\boldsymbol{\theta_{0}})+\rho(u)}du\right]_{l=1}^{m}$

and, because $Y_{i}$ is a Poisson point process under condition
($\mathcal{C}.2$),

$\operatorname{\bf{Cov}}_{\boldsymbol{\theta_{0}}}\left\\{\sum_{x\in
Y_{i}}\frac{\rho(x)\boldsymbol{z}(x)}{\lambda(x;\boldsymbol{\theta_{0}})+\rho(x)}\right\\}=\left[\int_{W\times
T}\frac{\lambda(u;\boldsymbol{\theta_{0}})\rho^{2}(u)z^{(k)}(u)z^{(l)}(u)}{(\lambda(u;\boldsymbol{\theta_{0}})+\rho(u))^{2}}du\right]_{k,l=1}^{m}.$
(10)

By condition ($\mathcal{C}.5$), the covariate terms in the integrand above are
all bounded. Furthermore, by condition ($\mathcal{C}.2$),

$0<\frac{\lambda(u;\boldsymbol{\theta_{0}})\rho^{2}(u)}{(\lambda(u;\boldsymbol{\theta_{0}})+\rho(u))^{2}}<\lambda(u;\boldsymbol{\theta_{0}}),\quad
0<\frac{\lambda(u;\boldsymbol{\theta_{0}})\lambda(v;\boldsymbol{\theta_{0}})\rho(u)\rho(v)}{(\lambda(u;\boldsymbol{\theta_{0}})+\rho(u))(\lambda(v;\boldsymbol{\theta_{0}})+\rho(v))}<\lambda(u;\boldsymbol{\theta_{0}})\lambda(v;\boldsymbol{\theta_{0}})$

and $\lambda(u;\boldsymbol{\theta_{0}}),\lambda(v;\boldsymbol{\theta_{0}})$
are absolutely integrable. Then, the multi-variate Lindeberg–Lévy central
limit theorem implies that
$n^{-1/2}(\boldsymbol{s}_{n}(X_{n};\boldsymbol{\theta_{0}})-n\text{\bf{M}}_{\boldsymbol{\theta_{0}}})$
converges under $P_{\boldsymbol{\theta_{0}}}$ in distribution to an
$m$-dimensional normally distributed random vector with mean zero and
covariance matrix given by (10).

Analogously, by condition ($\mathcal{C}.3$),
$n^{-1/2}(\boldsymbol{s}_{n}(D_{n},\boldsymbol{\theta_{0}})-n\text{\bf{M}}_{\boldsymbol{\theta_{0}}})$
converges under $P_{\boldsymbol{\theta_{0}}}$ in distribution to an
$m$-dimensional normally distributed random vector with mean zero and
covariance matrix given by

$\begin{split}&\operatorname{\bf{Cov}}_{\boldsymbol{\theta_{0}}}\left\\{\sum_{x\in
E_{i}}\frac{\lambda(x;\boldsymbol{\theta_{0}})\boldsymbol{z}(x)}{\lambda(x;\boldsymbol{\theta_{0}})+\rho(x)}\right\\}=\left[\int_{W\times
T}\frac{\lambda^{2}(u;\boldsymbol{\theta_{0}})\rho(u)z^{(k)}(u)z^{(l)}(u)}{(\lambda(u;\boldsymbol{\theta_{0}})+\rho(u))^{2}}du\right]_{k,l=1}^{m}\\\
&+\left[\int_{W\times T}\int_{W\times
T}\frac{\lambda(u;\boldsymbol{\theta_{0}})\lambda(v;\boldsymbol{\theta_{0}})\rho(u)\rho(v)z^{(k)}(u)z^{(l)}(v)}{(\lambda(u;\boldsymbol{\theta_{0}})+\rho(u))(\lambda(v;\boldsymbol{\theta_{0}})+\rho(v))}(g(u,v)-1)dudv\right]_{k,l=1}^{m}.\end{split}$

Applying Lévy’s continuity theorem and using the independence of $X_{n}$ and
$D_{n}$ under condition ($\mathcal{C}.1$), the weak limits for the two first-
order $U$-statistics can be combined, which completes the proof. ∎

######  Remark 3.7.

The pair correlation function $g(u,v)$ of $E_{i}$ can be tuned to control the
covariance matrix $\bm{V}$. When $E_{i}$, and thus $D_{n}$, is a Poisson point
process, $g(u,v)\equiv 1$ and the second term in $\bm{V}$ vanishes. Moreover,
regular point processes of $E_{i}$, whose $g(u,v)<1$, may be useful to reduce
the variance (Gautier et al., 2019).

Recalling (9), we are now in a position to conjecture a central limit theorem
for the logistic regression estimator. For a formal proof, though, we need to
analyze the error term in the Taylor series (8). The next theorem is then
concerned with asymptotic normality.

######  Theorem 3.8.

Assume that the conditions ($\mathcal{C}.1$)–($\mathcal{C}.8$) hold. Define
$\boldsymbol{s}_{n}(\boldsymbol{\theta})=\boldsymbol{s}(X_{n},D_{n};\boldsymbol{\theta})$
by (6) with $\boldsymbol{\theta}\in\boldsymbol{\Theta}$ and let
$\widehat{\boldsymbol{\theta}_{n}}$ be the logistic regression estimator for
which $\boldsymbol{s}_{n}(\widehat{\boldsymbol{\theta}_{n}})=\bm{0}$. If
$\widehat{\boldsymbol{\theta}_{n}}$ is attained, then as $n\to\infty$,
$n^{1/2}(\widehat{\boldsymbol{\theta}_{n}}-\boldsymbol{\theta_{0}})$ converges
under $P_{\boldsymbol{\theta_{0}}}$ in distribution to an $m$-dimensional
normally distributed random vector with mean zero and covariance matrix
$\bm{U}^{-1}\bm{V}(\bm{U}^{-1})^{\top}$, where $\bm{U},\bm{V}$ are as defined
in Lemma 3.5 and Lemma 3.6.

###### Proof.

Consider the error term in the Taylor series (8) which, under conditions
($\mathcal{C}.1$) and ($\mathcal{C}.2$) and recalling the score function (6),
reads

$\begin{split}&\frac{1}{2}(\widehat{\boldsymbol{\theta}_{n}}-\boldsymbol{\theta_{0}})^{\top}\nabla^{2}\boldsymbol{s}_{n}^{(i)}(\boldsymbol{\theta^{\prime,i}})(\widehat{\boldsymbol{\theta}_{n}}-\boldsymbol{\theta_{0}})=\frac{1}{2}\sum_{k=1}^{m}\sum_{l=1}^{m}(\widehat{\boldsymbol{\theta}_{n}}-\boldsymbol{\theta_{0}})^{(k)}(\widehat{\boldsymbol{\theta}_{n}}-\boldsymbol{\theta_{0}})^{(l)}\times\\\
&\times\sum_{x\in X_{n}\cup
D_{n}}\frac{\left(\lambda(x;\boldsymbol{\theta^{\prime,i}})-\rho(x)\right)\lambda(x;\boldsymbol{\theta^{\prime,i}})\rho(x)z^{(i)}(x)z^{(k)}(x)z^{(l)}(x)}{\left(\lambda(x;\boldsymbol{\theta^{\prime,i}})+\rho(x)\right)^{3}}\end{split}$
(11)

for all choices of $\boldsymbol{\theta^{\prime,i}}$ among the convex
combinations of $\widehat{\boldsymbol{\theta}_{n}}$ and
$\boldsymbol{\theta_{0}}$.

Firstly, note that

$\sum_{x\in X_{n}\cup
D_{n}}\frac{\left(\lambda(x;\boldsymbol{\theta^{\prime,i}})-\rho(x)\right)\lambda(x;\boldsymbol{\theta^{\prime,i}})\rho(x)z^{(i)}(x)z^{(k)}(x)z^{(l)}(x)}{n\left(\lambda(x;\boldsymbol{\theta^{\prime,i}})+\rho(x)\right)^{3}}$
(12)

is bounded in absolute value by

$G^{(n)}_{i,k,l}=\sum_{x\in X_{n}\cup
D_{n}}\frac{\left|z^{(i)}(x)z^{(k)}(x)z^{(l)}(x)\right|}{n}.$

This bound does not depend on $\boldsymbol{\theta^{\prime,i}}$ and thus does
not depend on $\widehat{\boldsymbol{\theta}_{n}}$. By the Campbell–Mecke
theorem,

$\mathbb{E}_{\boldsymbol{\theta_{0}}}\left\\{G^{(n)}_{i,k,l}\right\\}=\int_{W\times
T}(\lambda(u;\boldsymbol{\theta_{0}})+\rho(u))\left|z^{(i)}(u)z^{(k)}(u)z^{(l)}(u)\right|du.$

Under conditions ($\mathcal{C.}2$) and ($\mathcal{C.}5$),
$\mathbb{E}_{\boldsymbol{\theta_{0}}}\\{G^{(n)}_{i,k,l}\\}$ is a non-negative
constant. If it is zero, the product of the covariate terms inside has to be
zero almost everywhere on $W\times T$, which contradicts condition
($\mathcal{C.}8$). Thus,
$\mathbb{E}_{\boldsymbol{\theta_{0}}}\\{G^{(n)}_{i,k,l}\\}$ is strictly
positive. Then, by Markov’s inequality, for any $\delta>0$, there always
exists a finite
$H_{i,k,l}(\delta)=\mathbb{E}_{\boldsymbol{\theta_{0}}}\\{G^{(n)}_{i,k,l}\\}/\delta$
such that

$\mathbb{P_{\boldsymbol{\theta_{0}}}}\left\\{\left|G^{(n)}_{i,k,l}\right|\geq
H_{i,k,l}(\delta)\right\\}\leq\delta.$

Since $H_{i,k,l}(\delta)$ depends only on $\delta$ but not on $n$,
$G^{(n)}_{i,k,l}$, and thus (12), converges under
$P_{\boldsymbol{\theta_{0}}}$ in probability to $O(1)$. In addition, note that
this result applies to every component of (12) with $1\leq i,k,l\leq m$.

Now, move back to the Taylor expansion (8). By Lemma 3.5,
$-\nabla\boldsymbol{s}_{n}(\boldsymbol{\theta_{0}})/n$ converges almost
surely, and thus in probability, under $P_{\boldsymbol{\theta_{0}}}$ to
$\boldsymbol{U}$. Then, collecting all $i$ with $1\leq i\leq m$ and recalling
the results obtained above, (8) can be rewritten as

$\left[\boldsymbol{U}+o_{P}(1)-\frac{1}{2}(\widehat{\boldsymbol{\theta}_{n}}-\boldsymbol{\theta_{0}})^{\top}O_{P}(1)\right](\widehat{\boldsymbol{\theta}_{n}}-\boldsymbol{\theta_{0}})=\frac{\boldsymbol{s}_{n}(\boldsymbol{\theta_{0}})}{n}.$

Next, we appeal to strong consistency (cf., Theorem 3.3) to conclude that
$(\widehat{\boldsymbol{\theta}_{n}}-\boldsymbol{\theta_{0}})^{\top}O_{P}(1)$
converges under $P_{\boldsymbol{\theta_{0}}}$ in probability to zero.
Furthermore, by condition $\mathcal{C}.8$, the matrix $\boldsymbol{U}$ is
invertible. Thus, $\boldsymbol{U}+o_{P}(1)$ is also invertible with a
probability tending to one as $n\to\infty$. Multiplication by this inverse
then yields that

$n^{1/2}(\widehat{\boldsymbol{\theta}_{n}}-\boldsymbol{\theta_{0}})=\left[\boldsymbol{U}+o_{P}(1)\right]^{-1}\frac{\boldsymbol{s}_{n}(\boldsymbol{\theta_{0}})}{n^{1/2}}.$

The remainder of the proof is then a straightforward application of Slutsky’s
theorem. Obviously, $[\boldsymbol{U}+o_{P}(1)]^{-1}$ converges under
$P_{\boldsymbol{\theta_{0}}}$ in probability to $\bm{U}^{-1}$. By Lemma 3.6,
$\boldsymbol{s}_{n}(\boldsymbol{\theta_{0}})/n^{1/2}$ converges under
$P_{\boldsymbol{\theta_{0}}}$ in distribution to an $m$-dimensional normally
distributed random vector with mean zero and covariance matrix $\bm{V}$. Thus,
as $n\rightarrow\infty$,
$n^{1/2}(\widehat{\boldsymbol{\theta}_{n}}-\boldsymbol{\theta_{0}})$ converges
under $P_{\boldsymbol{\theta_{0}}}$ in distribution to an $m$-dimensional
normally distributed random vector with mean zero and covariance matrix
$\bm{U}^{-1}\bm{V}(\bm{U}^{-1})^{\top}$. ∎

######  Remark 3.9.

In Theorem 3.3, we have proved that the logistic regression estimator
$\widehat{\boldsymbol{\theta}_{n}}$ is a strongly consistent estimator. Here,
observing that $\bm{U}$ and $\bm{V}$ are constant matrices, we thus further
provide the convergence rate, which is
$(\widehat{\boldsymbol{\theta}_{n}}-\boldsymbol{\theta_{0}})^{(l)}=O_{P}(n^{-1/2})$
for every $1\leq l\leq m$.

######  Remark 3.10.

The approximate variance $n^{-1}\bm{U}^{-1}\bm{V}(\bm{U}^{-1})^{\top}$ of
${\widehat{\boldsymbol{\theta}_{n}}}$ decreases, as $n$ increases. The
covariance matrix depends on the unknown intensity function
$\lambda(u;\boldsymbol{\theta_{0}})$. In practice, a plug-in estimator of
$\lambda(u;\boldsymbol{\theta_{0}})$ would be used and the dummy point process
$D_{n}$ can be employed to approximate the integrals inside based on the
Campbell–Mecke theorem. For instance,

$\begin{split}n\bm{V}&\approx\left[\sum_{x\in
D_{n}}\frac{\lambda(x;\widehat{\boldsymbol{\theta}_{n}})z^{(k)}(x)z^{(l)}(x)}{\lambda(x;\widehat{\boldsymbol{\theta}_{n}})+\rho(x)}\right]_{k,l=1}^{m}+\left[\mathop{\sum\nolimits^{\neq}}_{x,y\in
D_{n}}\frac{\lambda(x;\widehat{\boldsymbol{\theta}_{n}})z^{(k)}(x)}{\lambda(x;\widehat{\boldsymbol{\theta}_{n}})+\rho(x)}\frac{\lambda(y;\widehat{\boldsymbol{\theta}_{n}})z^{(l)}(y)}{\lambda(y;\widehat{\boldsymbol{\theta}_{n}})+\rho(y)}\right]_{k,l=1}^{m}\\\
&-\left[\sum_{x\in
D_{n}}\frac{\lambda(x;\widehat{\boldsymbol{\theta}_{n}})z^{(k)}(x)}{\lambda(x;\widehat{\boldsymbol{\theta}_{n}})+\rho(x)}\sum_{y\in
D_{n}}\frac{\lambda(y;\widehat{\boldsymbol{\theta}_{n}})z^{(l)}(y)}{\lambda(y;\widehat{\boldsymbol{\theta}_{n}})+\rho(y)}\right]_{k,l=1}^{m}.\end{split}$

Similarly, $n\bm{U}$ can be approximated by the first of the three terms in
the right hand side of the equation above. When $E_{i}$, and thus $D_{n}$, is
a Poisson point process, $\bm{U}=\bm{V}$ and
$\bm{U}^{-1}\bm{V}(\bm{U}^{-1})^{\top}$ reduces to $\bm{U}^{-1}$. The
construction of approximate confidence intervals for the components of
$\bm{\theta_{0}}$ is now straightforward.

Note that both strong consistency and asymptotic normality of the logistic
regression estimator rely on the attainment of
$\widehat{\boldsymbol{\theta}_{n}}$. Hence, for the sake of completeness, we
also formulate the asymptotic existence of the logistic regression estimator
in the next theorem.

######  Theorem 3.11.

Assume that the conditions ($\mathcal{C}.1$)–($\mathcal{C}.8$) hold. Define
$\boldsymbol{s}_{n}(\boldsymbol{\theta})=\boldsymbol{s}(X_{n},D_{n};\boldsymbol{\theta})$
by (6) with $\boldsymbol{\theta}\in\boldsymbol{\Theta}$. Then, for every $n$,
an estimator $\widehat{\boldsymbol{\theta}_{n}}$ exists that solves
$\boldsymbol{s}_{n}(\widehat{\boldsymbol{\theta}_{n}})=\bm{0}$ with a
probability tending to one as $n\rightarrow\infty$.

###### Proof.

The proof follows from an appeal to Sørensen (1999, Corollary 2.6). By
condition ($\mathcal{C}.2$), $\boldsymbol{s}_{n}(\boldsymbol{\theta})$ is
continuously differentiable with respect to $\boldsymbol{\theta}$ for all
$\boldsymbol{\theta}\in\boldsymbol{\Theta}$. Thus, we only need to verify
Sørensen (1999, Condition 2.5). For elegance of writing, we omit the proof
here but provide a proof in more general cases in Lemma 4.5. ∎

##  4 Generalized infill asymptotic results

In this section, we discuss the extensions of our infill asymptotic results to
general point processes and to other unbiased estimating equations that are
based on the Campbell–Mecke theorem.

###  4.1 Extensions to general point processes

By releasing the Poisson assumption of $X_{n}$ and modifying corresponding
conditions, the infill asymptotic results for logistic regression estimators
obtained in Section 3 can be extended to other point processes, although a
log-linear form (1) for the intensity function may sometimes not be that
natural.

In the remainder of this section, we first demonstrate the extended infill
asymptotics for general spatio-temporal point processes. Afterwards, we
discuss the applicability for some specific families of point processes.

To extend the asymptotic results to general point processes, we shall need the
following modified conditions:

$(\mathcal{C}.9)$

$Y_{i}$ is a point process with intensity function
$\lambda(u;\boldsymbol{\theta})$ given by (1), where $b>0$ is an absolutely
integrable function, $\bm{z}$ is a measurable vector of covariates and the
parameter vector $\boldsymbol{\theta}$ lies in an open set
$\boldsymbol{\Theta}\subset\mathbb{R}^{m}$; $E_{i}$ has absolutely integrable
intensity function $\rho(u)>0$.

$(\mathcal{C}.10)$

Both $E_{i}$ and $Y_{i}$ have bounded pair correlation functions, $g$ and $h$,
that is, $\sup_{(u,v)\in(W\times T)^{2}}$ $g(u,v)$ $<\infty$ and
$\sup_{(u,v)\in(W\times T)^{2}}h(u,v)<\infty$ .

The theorems concerned with strong consistency, asymptotic normality and
existence of the logistic regression estimator are then as follows.

######  Theorem 4.1.

Assume that the conditions ($\mathcal{C}.1$),
($\mathcal{C}.4$)–($\mathcal{C}.7$) and ($\mathcal{C}.9$) hold. Define
$l_{n}(\boldsymbol{\theta})=l(X_{n},D_{n};\boldsymbol{\theta})$ by (5) with
$\boldsymbol{\theta}\in\boldsymbol{\Theta}$ and set
$\widehat{\boldsymbol{\theta}_{n}}=\operatorname*{arg\,max}_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}l_{n}(\boldsymbol{\theta})$.
If $\widehat{\boldsymbol{\theta}_{n}}$ is attained, then as
$n\rightarrow\infty$, $\widehat{\boldsymbol{\theta}_{n}}$ converges
$P_{\boldsymbol{\theta_{0}}}$-almost surely to $\boldsymbol{\theta_{0}}$.

###### Proof.

The same proofs as for Lemma 3.2 and Theorem 3.3 can be applied
straightforwardly. ∎

######  Theorem 4.2.

Assume that the conditions ($\mathcal{C}.1$),
($\mathcal{C}.4$)–($\mathcal{C}.10$) hold. Define
$\boldsymbol{s}_{n}(\boldsymbol{\theta})$ $=\boldsymbol{s}(X_{n},D_{n};$
$\boldsymbol{\theta})$ by (6) with $\boldsymbol{\theta}\in\boldsymbol{\Theta}$
and let $\widehat{\boldsymbol{\theta}_{n}}$ be the logistic regression
estimator for which
$\boldsymbol{s}_{n}(\widehat{\boldsymbol{\theta}_{n}})=\bm{0}$. If
$\widehat{\boldsymbol{\theta}_{n}}$ is attained, then as $n\rightarrow\infty$,
$n^{1/2}(\widehat{\boldsymbol{\theta}_{n}}-\boldsymbol{\theta_{0}})$ converges
under $P_{\boldsymbol{\theta_{0}}}$ in distribution to an $m$-dimensional
normally distributed random vector with mean zero and covariance matrix
$\bm{U}^{-1}\bm{V}(\bm{U}^{-1})^{\top}$, where $\bm{U}$ is as defined in Lemma
3.5 and $\bm{V}$ is now given by

$\begin{split}\bm{V}&=\left[\int_{W\times
T}\frac{\lambda(u;\boldsymbol{\theta_{0}})\rho(u)}{\lambda(u;\boldsymbol{\theta_{0}})+\rho(u)}z^{(k)}(u)z^{(l)}(u)du\right]_{k,l=1}^{m}\\\
&+\left[\int_{W\times T}\int_{W\times
T}\frac{\lambda(u;\boldsymbol{\theta_{0}})\lambda(v;\boldsymbol{\theta_{0}})\rho(u)\rho(v)z^{(k)}(u)z^{(l)}(v)}{(\lambda(u;\boldsymbol{\theta_{0}})+\rho(u))(\lambda(v;\boldsymbol{\theta_{0}})+\rho(v))}\left(g(u,v)+h(u,v)-2\right)dudv\right]_{k,l=1}^{m}.\end{split}$

###### Proof.

Similar proofs as for Lemma 3.5, Lemma 3.6 and Theorem 3.8 can be used. Note
that, by condition ($\mathcal{C}.10$), the pair correlation function of
$Y_{i}$, denoted by $h(u,v)$, is bounded, which ensures the entries of
$\bm{V}$ to be finite. ∎

######  Remark 4.3.

The covariance matrix now depends not only on the unknown intensity function
$\lambda(u;\boldsymbol{\theta_{0}})$ but also on the pair correlation function
$h(u,v)$. Similar to Remark 3.10, in practice, a plug-in estimator of
$\lambda(u;\boldsymbol{\theta_{0}})$ would be used and the point processes
$X_{n}$ and $D_{n}$ can be employed to approximate the integrals. Compared to
Remark 3.10, the approximation of $n\bm{V}$ needs the two extra terms below

$\left[\mathop{\sum\nolimits^{\neq}}_{x,y\in
X_{n}}\frac{\rho(x)z^{(k)}(x)}{\lambda(x;\widehat{\boldsymbol{\theta}_{n}})+\rho(x)}\frac{\rho(y)z^{(l)}(y)}{\lambda(y;\widehat{\boldsymbol{\theta}_{n}})+\rho(y)}\right]_{k,l=1}^{m}-\left[\sum_{x\in
D_{n}}\frac{\lambda(x;\widehat{\boldsymbol{\theta}_{n}})z^{(k)}(x)}{\lambda(x;\widehat{\boldsymbol{\theta}_{n}})+\rho(x)}\sum_{y\in
D_{n}}\frac{\lambda(y;\widehat{\boldsymbol{\theta}_{n}})z^{(l)}(y)}{\lambda(y;\widehat{\boldsymbol{\theta}_{n}})+\rho(y)}\right]_{k,l=1}^{m},$

while the approximation of $n\bm{U}$ stays the same. The construction of
approximate confidence intervals for the components of $\bm{\theta_{0}}$ is
then straightforward again.

######  Theorem 4.4.

Assume that the conditions ($\mathcal{C}.1$),
($\mathcal{C}.4$)–($\mathcal{C}.10$) hold. Define
$\boldsymbol{s}_{n}(\boldsymbol{\theta})=\boldsymbol{s}(X_{n},D_{n};\boldsymbol{\theta})$
by (6) with $\boldsymbol{\theta}\in\boldsymbol{\Theta}$. Then, for every $n$,
an estimator $\widehat{\boldsymbol{\theta}_{n}}$ exists that solves
$\boldsymbol{s}_{n}(\widehat{\boldsymbol{\theta}_{n}})=\bm{0}$ with a
probability tending to one as $n\rightarrow\infty$.

###### Proof.

A similar proof as for Theorem 3.11 can be applied as well. ∎

In application, for general point processes, the strength of interactions can
sometimes enter into the intensity functions.

For log-Gaussian Cox processes (Coles and Jones, 1991; Møller et al., 1998),
their intensity functions are the exponentials of Gaussian random variables
driven by associated Gaussian random fields with mean function
$\left(\boldsymbol{\theta}^{\top}\bm{z}\left(u\right)\right)\log
b\left(u\right)$ and covariance function $\sigma^{2}c(u,v)$. Here, $c(u,v)$ is
some correlation function and $\sigma^{2}>0$. Then, the intensity functions
read $\lambda(u;\boldsymbol{\theta})\exp{(\sigma^{2}/2)}$, which conforms with
the log-linear form (1) upon adding a component with the entry one to
$\boldsymbol{z}$. Should $c(u,v)$ depend on further parameters, e.g. the decay
rate of interactions, additional estimating equations are required.

For Markov point processes (Lieshout, 2000), in general, their intensity
functions are not known explicitly. However, in the vector function (3), one
can replace the intensity function $\lambda(u;\boldsymbol{\theta})$ by the
Papangelou conditional intensity function
$\lambda(u;\boldsymbol{\theta}|\bm{x})$ (Baddeley et al., 2014). The latter
characterizes the probability of a new point occurring in an infinitesimal
ball centred at a location and time combination $u\in W\times T$ given the
current point realisations $\bm{x}$ of $X$. Similar infill asymptotics in such
cases may be obtained upon adding appropriate conditions on the Papangelou
conditional intensity functions.

###  4.2 Extensions to other unbiased estimating equations

The Campbell–Mecke theorem provides the theoretical foundation for a wide
range of unbiased estimating equations. Although strong consistency of these
estimators may not hold due to the loss of the likelihood interpretation,
asymptotic normality and existence hold under appropriate conditions.

In the remainder of this section, we extend the central limit theorem (cf.,
Theorems 4.2 and 4.4) to more general unbiased estimating equations by
introducing the necessary background and providing full proofs.

Recalling (2), a more general unbiased estimating function based on the
Campbell–Mecke theorem can take the form

$\boldsymbol{s}(\boldsymbol{\theta})=\sum_{x\in
X}\boldsymbol{f}(x;\boldsymbol{\theta})-\sum_{x\in
D}\boldsymbol{f}(x;\boldsymbol{\theta})\frac{\lambda(x;\boldsymbol{\theta})}{\rho(x)}=\boldsymbol{0},$
(13)

for some test vector function $\boldsymbol{f}:W\times T\to\mathbb{R}^{m}$ such
that every component of $\boldsymbol{f}$, denoted by $f^{(i)}$ with $1\leq
i\leq m$, has the property that $f^{(i)}\lambda$ is absolutely integrable.
Note that, the test function $\boldsymbol{f}$ can contain some non-negative
quadrature weight functions (e.g., Waagepetersen, 2008; Guan and Shen, 2010)
depending on the first- and second-order characteristics of $X$ to reduce the
bias caused by deterministic numerical approximation.

To extend the central limit theorem to such unbiased estimating equations, we
shall need the following additional conditions:

$(\mathcal{C}.11)$

The component functions $f^{(i)}(u,\boldsymbol{\theta})$ are twice
continuously differentiable with respect to $\boldsymbol{\theta}$. Its
components and the first- and second-order partial derivatives, denoted by
$\frac{\partial}{\partial\theta^{(l)}}f^{(i)}(u;\boldsymbol{\theta})$ and
$\frac{\partial^{2}}{\partial\theta^{(k)}\partial\theta^{(l)}}f^{(i)}(u;\boldsymbol{\theta})$,
in a neighbourhood of $\boldsymbol{\theta_{0}}$ are bounded in absolute value
by some functions $d^{(i)}(u),d^{(i)}_{l}(u)$ and $d^{(i)}_{k,l}(u)$,
respectively, that are absolutely integrable with respect to
$\lambda(u;\boldsymbol{\theta_{0}})$.

$(\mathcal{C}.12)$

The second-order partial derivatives satisfy an adapted Hőlder condition in a
neighbourhood of $\boldsymbol{\theta_{0}}$: there exist some $\alpha\in(0,1]$
and some functions $e^{(i)}_{k,l}(u)$ which are absolutely integrable with
respect to $\lambda(u;\boldsymbol{\theta_{0}})$ such that

$\left|\frac{\partial^{2}}{\partial\theta^{(k)}\partial\theta^{(l)}}f^{(i)}(u;\boldsymbol{\theta})-\frac{\partial^{2}}{\partial\theta^{(k)}\partial\theta^{(l)}}f^{(i)}(u;\boldsymbol{\theta_{0}})\right|\leq
e^{(i)}_{k,l}(u)\|\boldsymbol{\theta}-\boldsymbol{\theta_{0}}\|^{\alpha}.$

$(\mathcal{C}.13)$

The matrix $\widetilde{\bm{U}}$, whose ($k,l$)-th entry reads $\int_{W\times
T}f^{(k)}(u;\boldsymbol{\theta_{0}})\lambda(u;\boldsymbol{\theta_{0}})z^{(l)}(u)du$
with $\boldsymbol{\theta_{0}}\in\boldsymbol{\Theta}$ and $1\leq k,l\leq m$, is
invertible.

$(\mathcal{C}.14)$

The product of any two component functions, denoted by
$f^{(k)}(u;\boldsymbol{\theta_{0}})f^{(l)}(u;\boldsymbol{\theta_{0}})$ with
$1\leq k,l\leq m$, is absolutely integrable with respect to
$\lambda(u;\boldsymbol{\theta_{0}})$.

A few remarks on the conditions above are appropriate. Condition
($\mathcal{C}.11$) and ($\mathcal{C}.12$) restrict the smoothness of $\bm{f}$
and its first-order and second-order derivatives. Condition ($\mathcal{C}.13$)
provides a statement corresponding to condition ($\mathcal{C}.8$) but in the
case of general unbiased estimating equations. Condition ($\mathcal{C}.14$)
complements necessary technical constraints to ensure that the entries in the
variances of these estimators are finite.

To establish the central limit theorem, we start our investigations with the
following lemma.

######  Lemma 4.5.

Assume that the conditions ($\mathcal{C}.1$), ($\mathcal{C}.5$),
($\mathcal{C}.9$) and ($\mathcal{C}.11$) hold. Define
$\boldsymbol{s}_{n}(\boldsymbol{\theta})=\boldsymbol{s}(X_{n},D_{n};\boldsymbol{\theta})$
by (13) with $\boldsymbol{\theta}\in\boldsymbol{\Theta}$. Then, as
$n\to\infty$, for all $\beta>0$,

$\sup_{\boldsymbol{\theta^{\prime,k}}\in
C_{\beta}^{(n)}}\left\|-\frac{1}{n}\nabla\boldsymbol{s}_{n}(\boldsymbol{\theta^{\prime,1}},\dots,\boldsymbol{\theta^{\prime,m}})-\widetilde{\bm{U}}\right\|$

converges under $P_{\boldsymbol{\theta_{0}}}$ in probability to zero, where
$C_{\beta}^{(n)}=\\{\boldsymbol{\theta}\in\boldsymbol{\Theta}:\|\boldsymbol{\theta}-\boldsymbol{\theta_{0}}\|\leq\beta/n^{1/2}\\}$
and

$\widetilde{\bm{U}}=\left[\int_{W\times
T}f^{(k)}(u;\boldsymbol{\theta_{0}})\lambda(u;\boldsymbol{\theta_{0}})z^{(l)}(u)du\right]_{k,l=1}^{m}.$

###### Proof.

Under conditions ($\mathcal{C}.1$), ($\mathcal{C}.9$) and ($\mathcal{C}.11$)
and recalling the unbiased estimating equation (13),
$-\nabla\boldsymbol{s}_{n}(\boldsymbol{\theta^{\prime,1}},\dots,\boldsymbol{\theta^{\prime,m}})/n$
reads

$\begin{split}&\frac{1}{n}\left[\sum_{x\in
X_{n}}-\frac{\partial}{\partial\theta^{(l)}}f^{(k)}(x;\boldsymbol{\theta^{\prime,k}})+\sum_{x\in
D_{n}}\frac{\lambda(x;\boldsymbol{\theta^{\prime,k}})}{\rho(x)}\left(\frac{\partial}{\partial\theta^{(l)}}f^{(k)}(x;\boldsymbol{\theta^{\prime,k}})+f^{(k)}(x;\boldsymbol{\theta^{\prime,k}})z^{(l)}(x)\right)\right]_{k,l=1}^{m}.\end{split}$

Fix $\beta>0$. We then prove the convergence of the supremum component-wisely.

Firstly, by condition ($\mathcal{C}.11$), one can obtain that, for large
enough $n$ and every $1\leq i,k,l\leq m$,

$\left|f^{(i)}(u;\boldsymbol{\theta})\right|\leq
d^{(i)}(u),\quad\left|\frac{\partial}{\partial\theta^{(l)}}f^{(i)}(u;\boldsymbol{\theta})\right|\leq
d_{l}^{(i)}(u),\quad\left|\frac{\partial^{2}}{\partial\theta^{(k)}\partial\theta^{(l)}}f^{(i)}(u;\boldsymbol{\theta})\right|\leq
d^{(i)}_{k,l}(u),$ (14)

for all $u\in W\times T$ and $\boldsymbol{\theta}\in C_{\beta}^{(n)}$. Note
that, taking $n$ large enough, $C_{\beta}^{(n)}$ lies entirely within
$\boldsymbol{\Theta}$, as the parameter space $\boldsymbol{\Theta}$ is open by
condition ($\mathcal{C}.9$).

Furthermore, under condition ($\mathcal{C}.9$), the first-order Taylor
expansion of $\lambda(u;\boldsymbol{\theta})$ with respect to
$\boldsymbol{\theta}$ at $\boldsymbol{\theta_{0}}$ reads

$\lambda(u;\boldsymbol{\theta})-\lambda(u;\boldsymbol{\theta_{0}})=b\left(u\right)\exp\left(\boldsymbol{\underline{\theta}}^{\top}\bm{z}\right)\sum_{i=1}^{m}(\boldsymbol{\theta}-\boldsymbol{\theta_{0}})^{(i)}z^{(i)}(u),$

where $\boldsymbol{\underline{\theta}}$ is a convex combination of
$\boldsymbol{\theta}$ and $\boldsymbol{\theta_{0}}$. By condition
($\mathcal{C}.5$), the covariate terms in the inner product
$\boldsymbol{\underline{\theta}}^{\top}\bm{z}$ and the sum above are all
bounded. Thus, recalling the definition of $C_{\beta}^{(n)}$, when
$n\to\infty$, one obtains that
$\lambda(u;\boldsymbol{\theta})=\lambda(u;\boldsymbol{\theta_{0}})+b(u)o(1)$
for all $u\in W\times T$ and $\boldsymbol{\theta}\in C_{\beta}^{(n)}$.

Now, consider
$-\nabla\boldsymbol{s}_{n}(\boldsymbol{\theta^{\prime,1}},\dots,\boldsymbol{\theta^{\prime,m}})/n$.
Denote its $(k,l)$-th entry by
$[-\nabla\boldsymbol{s}_{n}^{(k)}(\boldsymbol{\theta^{\prime,k}})/n]_{l}$ in
the remainder of the proof. The difference between
$[-\nabla\boldsymbol{s}_{n}^{(k)}(\boldsymbol{\theta^{\prime,k}})/n]_{l}$ and
$[-\nabla\boldsymbol{s}_{n}(\boldsymbol{\theta_{0}})/n]_{k,l}$ reads

$\begin{split}&-\frac{1}{n}\sum_{x\in
X_{n}}\frac{\partial}{\partial\theta^{(l)}}f^{(k)}(x;\boldsymbol{\theta^{\prime,k}})+\frac{1}{n}\sum_{x\in
X_{n}}\frac{\partial}{\partial\theta^{(l)}}f^{(k)}(x;\boldsymbol{\theta_{0}})\\\
&+\frac{1}{n}\sum_{x\in
D_{n}}\frac{\lambda(x;\boldsymbol{\theta^{\prime,k}})}{\rho(x)}\frac{\partial}{\partial\theta^{(l)}}f^{(k)}(x;\boldsymbol{\theta^{\prime,k}})-\frac{1}{n}\sum_{x\in
D_{n}}\frac{\lambda(x;\boldsymbol{\theta_{0}})}{\rho(x)}\frac{\partial}{\partial\theta^{(l)}}f^{(k)}(x;\boldsymbol{\theta_{0}})\\\
&+\frac{1}{n}\sum_{x\in
D_{n}}\frac{\lambda(x;\boldsymbol{\theta^{\prime,k}})}{\rho(x)}f^{(k)}(x;\boldsymbol{\theta^{\prime,k}})z^{(l)}(x)-\frac{1}{n}\sum_{x\in
D_{n}}\frac{\lambda(x;\boldsymbol{\theta_{0}})}{\rho(x)}f^{(k)}(x;\boldsymbol{\theta_{0}})z^{(l)}(x).\end{split}$
(15)

Below, we specifically analyze the asymptotic behavior of the second line in
(15). The analysis of the other two lines proceeds along similar lines.

Under condition ($\mathcal{C}.11$), the first-order Taylor expansion of
$\frac{\partial}{\partial\theta^{(l)}}f^{(k)}(x;\boldsymbol{\theta^{\prime,k}})$
with respect to $\boldsymbol{\theta}$ at $\boldsymbol{\theta_{0}}$ reads

$\frac{\partial}{\partial\theta^{(l)}}f^{(k)}(x;\boldsymbol{\theta^{\prime,k}})-\frac{\partial}{\partial\theta^{(l)}}f^{(k)}(x;\boldsymbol{\theta_{0}})=\sum_{i=1}^{m}(\boldsymbol{\theta^{\prime,k}}-\boldsymbol{\theta_{0}})^{(i)}\frac{\partial^{2}}{\partial\theta^{(i)}\partial\theta^{(l)}}f^{(k)}(x;\boldsymbol{\theta^{\prime\prime,k}}),$

where $\boldsymbol{\theta^{\prime\prime,k}}$ is a convex combination of
$\boldsymbol{\theta^{\prime,k}}$ and $\boldsymbol{\theta_{0}}$, and thus
$\boldsymbol{\theta^{\prime\prime,k}}\in C_{\beta}^{(n)}$. Then, recalling
that
$\lambda(x;\boldsymbol{\theta^{\prime,k}})=\lambda(x;\boldsymbol{\theta_{0}})+b(x)o(1)$,
the second line in (15) is bounded in absolute value by

$\frac{1}{n}\sum_{x\in
D_{n}}\left\\{\sum_{i=1}^{m}\frac{\lambda(x;\boldsymbol{\theta_{0}})}{\rho(x)}\left|(\boldsymbol{\theta^{\prime,k}}-\boldsymbol{\theta_{0}})^{(i)}\right|\left|\frac{\partial^{2}}{\partial\theta^{(i)}\partial\theta^{(l)}}f^{(k)}(x;\boldsymbol{\theta^{\prime\prime,k}})\right|+\left|\frac{b(x)o(1)}{\rho(x)}\right|\left|\frac{\partial}{\partial\theta^{(l)}}f^{(k)}(x;\boldsymbol{\theta^{\prime,k}})\right|\right\\},$

which, recalling the definition of $C_{\beta}^{(n)}$ and the bound functions
in (14), is further bounded by

$\sum_{i=1}^{m}\sum_{x\in
D_{n}}\frac{\beta}{n^{3/2}}\frac{\lambda(x;\boldsymbol{\theta_{0}})}{\rho(x)}d_{i,l}^{(k)}(x)+\sum_{x\in
D_{n}}\frac{1}{n}\frac{b(x)}{\rho(x)}|o(1)|d_{l}^{(k)}(x)$ (16)

for large enough $n$. Note that this bound depends only on $n$ and not on
$\boldsymbol{\theta^{\prime,k}}$. Furthermore, by condition
($\mathcal{C}.11$), one can apply the Campbell-Mecke theorem and Markov’s
inequality as in the proof of Theorem 3.8 to obtain that the inner sum in the
first term of (16) converges under $P_{\boldsymbol{\theta_{0}}}$ in
probability to zero. Since this result applies to every component of the inner
sum with $1\leq i\leq m$ and to the second term, (16) converges under
$P_{\boldsymbol{\theta_{0}}}$ in probability to zero. Hence,

$\sup_{\boldsymbol{\theta^{\prime,k}}\in
C_{\beta}^{(n)}}\left|\frac{1}{n}\sum_{x\in
D_{n}}\frac{\lambda(x;\boldsymbol{\theta^{\prime,k}})}{\rho(x)}\frac{\partial}{\partial\theta^{(l)}}f^{(k)}(x;\boldsymbol{\theta^{\prime,k}})-\frac{1}{n}\sum_{x\in
D_{n}}\frac{\lambda(x;\boldsymbol{\theta_{0}})}{\rho(x)}\frac{\partial}{\partial\theta^{(l)}}f^{(k)}(x;\boldsymbol{\theta_{0}})\right|=o_{P}(1).$

Applying similar proofs to the other two lines in (15) and combining the
obtained results as above, by the triangle inequality, one obtains that

$\sup_{\boldsymbol{\theta^{\prime,k}}\in
C_{\beta}^{(n)}}\left|\left[-\frac{1}{n}\nabla\boldsymbol{s}_{n}^{(k)}(\boldsymbol{\theta^{\prime,k}})\right]_{l}-\left[-\frac{1}{n}\nabla\boldsymbol{s}_{n}(\boldsymbol{\theta_{0}})\right]_{k,l}\right|=o_{P}(1).$
(17)

Finally, consider the analogue of Lemma 3.5 for the unbiased estimating
equation (13). Under conditions ($\mathcal{C}.1$), ($\mathcal{C}.9$) and
($\mathcal{C}.11$), $-\nabla\boldsymbol{s}_{n}(\boldsymbol{\theta_{0}})/n$ now
reads

$\begin{split}&\frac{1}{n}\sum_{i=1}^{n}\left[\sum_{x\in
Y_{i}}-\frac{\partial}{\partial\theta^{(l)}}f^{(k)}(x;\boldsymbol{\theta_{0}})+\sum_{x\in
E_{i}}\frac{\lambda(x;\boldsymbol{\theta_{0}})}{\rho(x)}\left(\frac{\partial}{\partial\theta^{(l)}}f^{(k)}(x;\boldsymbol{\theta_{0}})+f^{(k)}(x;\boldsymbol{\theta_{0}})z^{(l)}(x)\right)\right]_{k,l=1}^{m}.\end{split}$

Again, by conditions ($\mathcal{C}.5$) and ($\mathcal{C}.11$), one can apply
the Campbell–Mecke theorem and Kolmogorov’s strong law of large numbers to
obtain that $-\nabla\boldsymbol{s}_{n}(\boldsymbol{\theta_{0}})/n$ converges
almost surely, and thus in probability, under $P_{\boldsymbol{\theta_{0}}}$ to
$\widetilde{\bm{U}}$. Accordingly,

$\left|\left[-\frac{1}{n}\nabla\boldsymbol{s}_{n}(\boldsymbol{\theta_{0}})\right]_{k,l}-\widetilde{\bm{U}}_{k,l}\right|=o_{P}(1).$
(18)

The proof is completed by applying the triangle inequality to combine (17) and
(18). ∎

The next theorem is then concerned with the central limit theorem.

######  Theorem 4.6.

Assume that the conditions ($\mathcal{C}.1$),
($\mathcal{C}.4$)–($\mathcal{C}.5$) and ($\mathcal{C}.9$)–($\mathcal{C}.14$)
hold. Define $\boldsymbol{s}_{n}(\boldsymbol{\theta})$ by (13) with
$\boldsymbol{\theta}\in\boldsymbol{\Theta}$. Then, as $n\rightarrow\infty$,
for every $n$, an estimator $\widehat{\boldsymbol{\theta}_{n}}$ exists that
solves $\boldsymbol{s}_{n}(\boldsymbol{\theta})=0$ with a probability tending
to one. Moreover, $\widehat{\boldsymbol{\theta}_{n}}$ converges under
$P_{\boldsymbol{\theta_{0}}}$ in probability to $\boldsymbol{\theta_{0}}$ and
$n^{1/2}(\widehat{\boldsymbol{\theta}_{n}}-\boldsymbol{\theta_{0}})$ converges
under $P_{\boldsymbol{\theta_{0}}}$ in distribution to an $m$-dimensional
normally distributed random vector with mean $\boldsymbol{0}$ and covariance
matrix
$\widetilde{\bm{U}}^{-1}\widetilde{\bm{V}}(\widetilde{\bm{U}}^{-1})^{\top}$,
where $\widetilde{\bm{U}}$ is as defined in Lemma 4.5 and $\widetilde{\bm{V}}$
is given by

$\begin{split}\widetilde{\bm{V}}&=\left[\int_{W\times
T}f^{(k)}(u;\boldsymbol{\theta_{0}})f^{(l)}(u;\boldsymbol{\theta_{0}})\lambda(u;\boldsymbol{\theta_{0}})\left(1+\frac{\lambda(u;\boldsymbol{\theta_{0}})}{\rho(u)}\right)du\right]_{k,l=1}^{m}\\\
&+\left[\int_{W\times T}\int_{W\times
T}f^{(k)}(u;\boldsymbol{\theta_{0}})\lambda(u;\boldsymbol{\theta_{0}})f^{(l)}(v;\boldsymbol{\theta_{0}})\lambda(v;\boldsymbol{\theta_{0}})(g(u,v)+h(u,v)-2)dudv\right]_{k,l=1}^{m}.\end{split}$

###### Proof.

First of all, by condition ($\mathcal{C}.11$),
$\boldsymbol{s}_{n}(\boldsymbol{\theta})$ is twice continuously differentiable
with respect to $\boldsymbol{\theta}$. Then, we shall verify the analogues of
Lemmas 3.5 and 3.6 for the unbiased estimating equation (13).

The analogue of Lemma 3.5 was established in Lemma 4.5. For Lemma 3.6, under
conditions ($\mathcal{C}.1$) and ($\mathcal{C}.9$),
$\boldsymbol{s}_{n}(\boldsymbol{\theta})/n^{1/2}$ now reads

$\frac{\boldsymbol{s}_{n}(\boldsymbol{\theta_{0}})}{n^{1/2}}=n^{1/2}\left\\{\frac{1}{n}\sum_{i=1}^{n}\sum_{x\in
Y_{i}}\boldsymbol{f}(x;\boldsymbol{\theta_{0}})-\frac{1}{n}\sum_{i=1}^{n}\sum_{x\in
E_{i}}\boldsymbol{f}(x;\boldsymbol{\theta_{0}})\frac{\lambda(x;\boldsymbol{\theta_{0}})}{\rho(x)}\right\\}.$

The two terms between the curly brackets above are both averages of
independent and identically distributed real-valued random vectors. By
conditions ($\mathcal{C}.9$)–($\mathcal{C}.11$), their means are finite and
identical. Additionally by conditions ($\mathcal{C}.4$) and
($\mathcal{C}.14$), their variances are also finite and given by the entries
of $\widetilde{\bm{V}}$. Then, similar to the proof of Lemma 3.6, one can
apply the Campbell–Mecke theorem and multi-variate Lindeberg–Lévy’s central
limit theorem to obtain that
$\boldsymbol{s}_{n}(\boldsymbol{\theta_{0}})/n^{1/2}$ converges under
$P_{\boldsymbol{\theta_{0}}}$ in distribution to an $m$-dimensional normally
distributed random vector with mean zero and covariance matrix
$\widetilde{\bm{V}}$.

From Sørensen (1999, Condition 2.7, Corollary 2.8 and Theorem 2.9) and by
condition ($\mathcal{C}.13$), we thus only need to prove that

$\sup_{\boldsymbol{\theta}\in
C_{\beta}^{(n)}}\left\|-\frac{1}{n}\boldsymbol{s}_{n}(\boldsymbol{\theta})\right\|,\quad\sup_{\boldsymbol{\theta^{\prime,k}}\in
C_{\beta}^{(n)}}\left\|-\frac{1}{n}\nabla\boldsymbol{s}_{n}(\boldsymbol{\theta^{\prime,1}},\dots,\boldsymbol{\theta^{\prime,m}})-\widetilde{\bm{U}}\right\|$

and

$\sup_{\boldsymbol{\theta^{\prime,k}}\in
C_{\beta}^{(n)}}\left\|-\frac{1}{n}\nabla^{2}\boldsymbol{s}^{(i)}_{n}(\boldsymbol{\theta^{\prime,1}},\dots,\boldsymbol{\theta^{\prime,m}})-\widetilde{\bm{Q}}^{(i)}\right\|$

all converge under $P_{\boldsymbol{\theta_{0}}}$ in probability to zero, where
$\widetilde{\bm{Q}}^{(i)}$ is given by

$\begin{split}\widetilde{\bm{Q}}^{(i)}&=\left[\int_{W\times
T}f^{(i)}(u,\boldsymbol{\theta_{0}})\lambda(u,\boldsymbol{\theta_{0}})z^{(k)}(u)z^{(l)}(u)du\right]_{k,l=1}^{m}\\\
&+\left[\int_{W\times
T}\left(\frac{\partial}{\partial\theta^{(k)}}f^{(i)}(u,\boldsymbol{\theta_{0}})z^{(l)}(u)+\frac{\partial}{\partial\theta^{(l)}}f^{(i)}(u,\boldsymbol{\theta_{0}})z^{(k)}(u)\right)\lambda(u,\boldsymbol{\theta_{0}})du\right]_{k,l=1}^{m}.\end{split}$

In Lemma 4.5, we already proved that

$\sup_{\boldsymbol{\theta^{\prime,k}}\in
C_{\beta}^{(n)}}\left\|-\frac{1}{n}\nabla\boldsymbol{s}_{n}(\boldsymbol{\theta^{\prime,1}},\dots,\boldsymbol{\theta^{\prime,m}})-\widetilde{\bm{U}}\right\|\stackrel{{\scriptstyle
P_{\boldsymbol{\theta_{0}}}}}{{\to}}0.$

The remainder of the proof mostly proceeds along similar lines for the other
two cases.

Note that condition ($\mathcal{C}.12$) is needed in proving the convergence of
some terms in the $(k,l)$-th entry of
$-\nabla^{2}\boldsymbol{s}^{(i)}_{n}(\boldsymbol{\theta^{\prime,1}},\dots,\boldsymbol{\theta^{\prime,m}})/n$.
For instance, we need to prove that

$\sup_{\boldsymbol{\theta^{\prime,k}}\in
C_{\beta}^{(n)}}\left|\frac{1}{n}\sum_{x\in
D_{n}}\frac{\lambda(x;\boldsymbol{\theta^{\prime,k}})}{\rho(x)}\frac{\partial^{2}}{\partial\theta^{(j)}\theta^{(l)}}f^{(k)}(x;\boldsymbol{\theta^{\prime,k}})-\frac{1}{n}\sum_{x\in
D_{n}}\frac{\lambda(x;\boldsymbol{\theta_{0}})}{\rho(x)}\frac{\partial^{2}}{\partial\theta^{(j)}\theta^{(l)}}f^{(k)}(x;\boldsymbol{\theta_{0}})\right|=o_{P}(1).$

Since $f^{(k)}(u,\boldsymbol{\theta})$ is only twice continuously
differentiable with respect to $\boldsymbol{\theta}$, the proof in Lemma 4.5
does not apply here. However, recalling from the proof of Lemma 4.5 that
$\lambda(x;\boldsymbol{\theta^{\prime,k}})=\lambda(x;\boldsymbol{\theta_{0}})+b(x)o(1)$,
the absolute value term above is bounded by

$\frac{1}{n}\sum_{x\in
D_{n}}\left\\{\frac{\lambda(x;\boldsymbol{\theta_{0}})}{\rho(x)}\left|\frac{\partial^{2}f^{(k)}(x;\boldsymbol{\theta^{\prime,k}})}{\partial\theta^{(j)}\theta^{(l)}}-\frac{\partial^{2}f^{(k)}(x;\boldsymbol{\theta_{0}})}{\partial\theta^{(j)}\theta^{(l)}}\right|+\left|\frac{b(x)o(1)}{\rho(x)}\right|\left|\frac{\partial^{2}f^{(k)}(x;\boldsymbol{\theta^{\prime,k}})}{\partial\theta^{(j)}\theta^{(l)}}\right|\right\\},$

which, under the adapted Hőlder condition in condition ($\mathcal{C}.12$), is
further bounded by

$\sum_{x\in
D_{n}}\left\\{\frac{\beta}{n^{\alpha/2+1}}\frac{\lambda(x;\boldsymbol{\theta_{0}})}{\rho(x)}e_{i,l}^{(k)}(x)+\frac{1}{n}\frac{b(x)}{\rho(x)}\left|o(1)\right|d_{i,l}^{(k)}(x)\right\\}$

for large enough $n$. This bound depends only on $n$ and not on
$\boldsymbol{\theta}^{\prime,k}$. By conditions ($\mathcal{C}.11$) and
($\mathcal{C}.12$), one can again apply the Campbell-Mecke theorem and
Markov’s inequality as in the proof of Theorem 3.8 to obtain that the bound
converges under $P_{\boldsymbol{\theta_{0}}}$ in probability to zero. ∎

######  Remark 4.7.

Looking back on the logistic regression estimator given by (3), its component
functions are smooth with respect to $\boldsymbol{\theta}$ and thus the
second-order derivatives satisfy the adapted Hőlder condition. Moreover, its
component functions and the first- and second-order derivatives are bounded
under condition ($\mathcal{C}.5$).

##  5 Conclusion

In this paper, we established strong consistency and asymptotic normality for
logistic regression estimators for spatio-temporal point processes under the
regime of infill asymptotics. We also extended the asymptotic results to other
unbiased estimating functions that are based on the Campbell–Mecke theorem. We
demonstrated how to construct approximate confidence intervals in practice.
Compared to increasing-domain asymptotics, infill asymptotics provide the
theoretical foundation on the asymptotic behaviour of parameter estimates
where identical point patterns are only valid within a fixed region.

In addition, we assumed that spatio-temporal point processes of interest have
log-linear parametric intensity functions, which simplifies the computation.
However, we believe that the asymptotic results can be proved for a wide class
of intensity functions. Moreover, motivated by Baddeley et al. (2014), it
would also be interesting to extend infill asymptotics to the Gibbs domain by
studying associated background and conditions.

## Acknowledgements

This research was funded by the Dutch Research Council (NWO) for the project
‘Data driven risk management for fire services’ (18004). We thank Maurits de
Graaf and Clément Lezane for valuable discussions and comments.

## References

  * Baddeley and Turner (2000) A. Baddeley and R. Turner.  Practical maximum pseudolikelihood for spatial point patterns.  _Australian & New Zealand Journal of Statistics_, 42(3):283–322, 2000.
  * Baddeley et al. (2010) A. Baddeley, M. Berman, N.I. Fisher, A. Hardegen, R.K. Milne, D. Schuhmacher, R. Shah, and R. Turner.  Spatial logistic regression and change-of-support in poisson point processes.  _Electronic Journal of Statistics_ , 4:1151–1201, 2010\.
  * Baddeley et al. (2014) A. Baddeley, J. Coeurjolly, E. Rubak, and R. Waagepetersen.  Logistic regression for spatial gibbs point processes.  _Biometrika_ , 101:377–392, 2014.
  * Baddeley et al. (2015) A. Baddeley, E. Rubak, and R. Turner.  _Spatial Point Patterns: Methodology and Applications with R_.  Chapman and Hall/CRC, 2015.
  * Besag (1977) J. Besag.  Some methods of statistical analysis for spatial data.  _Bulletin of the International Statistical Institute_ , 47:77–92, 1977.
  * Bray and Schoenberg (2013) A. Bray and F. Schoenberg.  Assessment of point process models for earthquake forecasting.  _Statistical Science_ , 28(4):520–520, 2013.
  * Choiruddin et al. (2018) A. Choiruddin, J. Coeurjolly, and F. Letué.  Convex and non-convex regularization methods for spatial point processes intensity estimation.  _Electronic Journal of Statistics_ , 12(1):1210–1255, 2018.
  * Coeurjolly and Møller (2014) J. Coeurjolly and Jesper Møller.  Variational approach for spatial point process intensity estimation.  _Bernoulli_ , 20(3):1097–1125, 2014.
  * Coles and Jones (1991) P. Coles and B. Jones.  A lognormal model for the cosmological mass distribution.  _Monthly Notices of the Royal Astronomical Society_ , 248:1–13, 1991.
  * Daley and Vere-Jones (2009) D.J. Daley and D. Vere-Jones.  _An Introduction to the Theory of Point Processes_.  Springer, New York, 2009.
  * Gautier et al. (2019) Guillaume Gautier, Rémi Bardenet, and Michal Valko.  On two ways to use determinantal point processes for monte carlo integration.  In _NeurIPS_ , 2019.
  * Guan and Shen (2010) Y. Guan and Y. Shen.  A weighted estimating equation approach for inhomogeneous spatial point processes.  _Biometrika_ , 97(4):867–880, 2010.
  * Guyon (1995) X. Guyon.  _Random Fields on a Network: Modeling, Statistics and Applications_.  New York: Springer, 1995.
  * Lieshout (2000) M.N.M. van Lieshout.  _Markov point processes and their applications_.  Imperial College Press, 2000.
  * Lieshout (2020) M.N.M. van Lieshout.  Infill asymptotics and bandwidth selection for kernel estimators of spatial intensity functions.  _Methodology and Computing in Applied Probability_ , 22(3):995–1008, 2020.
  * Lu et al. (2021) C. Lu, M.N.M. van Lieshout, M. de Graaf, and P. Visscher.  Data-driven chimney fire risk prediction using machine learning and point process tools, 2021.
  * Møller and Díaz-Avalos (2010) J. Møller and C. Díaz-Avalos.  Structured spatio-temporal shot-noise Cox point process models, with a view to modelling forest fires.  _Scandinavian Journal of Statistics_ , 37(1):2–25, 2010.
  * Møller and Waagepetersen (2007) J. Møller and R. P. Waagepetersen.  Modern statistics for spatial point processes.  _Scandinavian Journal of Statistics_ , 34:643–684, 2007\.
  * Møller and Waagepetersen (2004) J. Møller and R.P. Waagepetersen.  _Statistical Inference and Simulation for Spatial Point Processes_.  Chapman and Hall/CRC, 2004.
  * Møller et al. (1998) J. Møller, A.R. Syversveen, and R.P. Waagepetersen.  Log gaussian cox processes.  _Scandinavian Journal of Statistics_ , 25(3):451–482, 1998.
  * Ogata and Tanemura (1981) Y. Ogata and M. Tanemura.  Estimation of interaction potentials of spatial point patterns through the maximum likelihood procedure.  _Annals of the Institute of Statistical Mathematics_ , 33:315–338, 1981.
  * Reitzner and Schulte (2013) M. Reitzner and M. Schulte.  Central limit theorems for $U$-statistics of poisson point processes.  _Annals of Probability_ , 41(6):3879–3909, 2013\.
  * Ripley (1988) B.D. Ripley.  _Statistical Inference for Spatial Processes._ Cambridge University Press, 1988.
  * Silvapulle (1981) M.J. Silvapulle.  On the existence of maximum likelihood estimators for the binomial response models.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 43(3):310–313, 1981.
  * Sørensen (1999) M. Sørensen.  On asymptotics of estimating functions.  _Brazilian Journal of Probability and Statistics_ , 13(2):111–136, 1999.
  * Thurman et al. (2015) A.L. Thurman, R. Fu, Y. Guan, and J. Zhu.  Regularized estimating equations for model selection of clustered spatial point processes.  _Statistica Sinica_ , 25(1):173–188, 2015.
  * Venables and Ripley (2002) W.N. Venables and B.D. Ripley.  _Modern Applied Statistics With S_.  New York: Springer, 2002.
  * Waagepetersen (2008) R.P. Waagepetersen.  Estimating functions for inhomogeneous spatial point processes with incomplete covariate data.  _Biometrika_ , 95:351–363, 2008.