# Inverse regression for spatially distributed functional data

F.Suneel Babu<EMAIL_ADDRESS>[ S.Ruiqi
<EMAIL_ADDRESS>[ Department of Mathematical Sciences, University
of Texas at El Paso, TX 79968, USApresep=, ]e1  Department of Mathematics and
Statistics, Texas Tech University, TX 79409, USApresep=, ]e2

# Supplement to “Inverse regression for spatially distributed functional data”

F.Suneel Babu<EMAIL_ADDRESS>[ S.Ruiqi
<EMAIL_ADDRESS>[ Department of Mathematical Sciences, University
of Texas at El Paso, TX 79968, USApresep=, ]e1  Department of Mathematics and
Statistics, Texas Tech University, TX 79409, USApresep=, ]e2

###### Abstract

Spatially distributed functional data are prevalent in many statistical
applications such as meteorology, energy forecasting, census data, disease
mapping, and neurological studies. Given their complex and high-dimensional
nature, functional data often require dimension reduction methods to extract
meaningful information. Inverse regression is one such approach that has
become very popular in the past two decades. We study the inverse regression
in the framework of functional data observed at irregularly positioned spatial
sites. The functional predictor is the sum of a spatially dependent functional
effect and a spatially independent functional nugget effect, while the
relation between the scalar response and the functional predictor is modeled
using the inverse regression framework. For estimation, we consider local
linear smoothing with a general weighting scheme, which includes as special
cases the schemes under which equal weights are assigned to each observation
or to each subject. This framework enables us to present the asymptotic
results for different types of sampling plans over time such as non-dense,
dense, and ultra-dense. We discuss the domain-expanding infill (DEI) framework
for spatial asymptotics, which is a mix of the traditional expanding domain
and infill frameworks. The DEI framework overcomes the limitations of
traditional spatial asymptotics in the existing literature. Under this unified
framework, we develop asymptotic theory and identify conditions that are
necessary for the estimated eigen-directions to achieve optimal rates of
convergence. Our asymptotic results include pointwise and $L_{2}$ convergence
rates. Simulation studies using synthetic data and an application to a real-
world dataset confirm the effectiveness of our methods.

Covariance operator,

domain-expanding infill asymptotics,

irregularly positioned,

local linear smoothing,

nugget effect,

unified framework,

###### keywords:

##  1 Introduction

Sufficient dimension reduction methods have become popular owing to their
usefulness in extracting useful information from high-dimensional data and
have found a wide range of applications. For a brief overview of sufficient
dimension reduction methods, we refer to Li (2018) and Cook (2018). For recent
advancements on post-reduction inference, we refer to Kim et al. (2020). In
this proposal, we present a method for estimating the sufficient dimension
reduction space where the scalar response and the functional predictors are
spatially distributed and related through linear indices as assumed in Müller
and Stadtmüller (2005), Ferré and Yao (2003), Ferré and Yao (2005), Hsing and
Ren (2009) and Jiang, Yu and Wang (2014). Spatially distributed functional
data are prevalent in many statistical applications such as meteorology,
energy forecasting, census data, disease mapping, and neurological studies.
One well-known application is the spatial prediction of the Canadian
temperature and precipitation data (Ramsay and Silverman, 2005; Giraldo,
Delicado and Mateu, 2010). The dataset includes the daily average temperatures
at 35 weather stations over 30 years. The locations of these 35 stations serve
as spatial coordinates for the spatial prediction of weather data.

Although there are model-based and model-free sufficient dimension reduction
methods, the latter is more robust and flexible as it relies on weaker
assumptions. The inverse regression or sliced inverse regression (SIR)
proposed by Li (1991) and Duan and Li (1991) is one important approach in this
paradigm. The SIR approach assumes that the scalar response variable in a
regression model depends on the multivariate predictor through an unknown
number of linear projections. Formally,

$\displaystyle Y\perp\\!\\!\\!\perp\bm{X}~{}|~{}\bm{\beta}^{T}\bm{X}$

where $\bm{X}$ is a $p$-dimensional random vector, $Y$ is a random variable,
and the coefficient matrix $\bm{\beta}$ is a $p\times K$ matrix such that $K$
is much smaller than $p$. The matrix $\bm{\beta}$ is not identifiable.
However, the space spanned by the columns of $\bm{\beta}$, the effective
dimension reduction (e.d.r.) space, is of interest. Li (1991) showed that the
e.d.r. directions $\\{b_{j}\\}_{j=1}^{K}$ are estimated from the generalized
eigen-analysis

$\displaystyle\text{cov}[E(\bm{X}|Y)]b_{j}$
$\displaystyle=\lambda_{j}\text{cov}(\bm{X})b_{j},$ (1.1)

where the $\lambda_{j}$’s are eigenvalues, under the following design
condition:

$\displaystyle E(\bm{c}^{T}\bm{X}|\bm{\beta}^{T}\bm{X})\text{ is linear in
}\bm{\beta}^{T}\bm{X},\text{ for any direction }\bm{c}\in\mathbb{R}^{p}.$

For more variants of the above framework, we refer to, for example, Cook and
Weisberg (1991), Cook and Li (2002), Cook (2007), and Cook, Forzani and Yao
(2010).

Proposing the SIR approach to functional data $X(\cdot)$, which are stochastic
processes observed over a time interval $\mathcal{I}\subset\mathbb{R}$, is not
simple due to the complication of inverting a covariance operator on
$L_{2}(\mathcal{I})$. The first papers to take functional SIR into
consideration were Ferré and Yao (2003, 2005). They showed that it is feasible
to estimate the e.d.r. space under regularity conditions as long as the true
dimension of the space is known (Forzani and Cook, 2007; Ferré and Yao, 2007).
For some further refinements and alternatives, we refer to Hsing and Ren
(2009), Cook, Forzani and Yao (2010), Chen, Hall and Müller (2011), Li and
Song (2017); Song (2019), Li and Song (2022) and Lee and Li (2022). An
assumption commonly made in functional data analysis is that the complete
trajectories for a sample of $n$ random functions are fully observed. This
assumption, however, may not hold in practice because the measurements are
obtained at discrete and scattered time points. As mentioned in Jiang, Yu and
Wang (2014), one way to overcome this limitation is to borrow information
across all sample functions by applying smoothing to obtain estimators.
Unfortunately, the spatial points are often irregularly positioned in space
which complicate the smoothing of functional data.

One important issue with spatial data is that they require a $d$-dimensional
index space $\\{\bm{s}\in\mathbb{R}^{d},d\geq 2\\}$ and one needs to take the
dependence across all directions into consideration. Instead of defining a
random field through specific model equations, the dependence conditions are
often imposed by the decay of covariances or the use of mixing conditions.
Different asymptotics holds in the presence of clusters and for sparsely
distributed points. Traditionally, there exist two different approaches for
spatial data: expanding domain and infill asymptotics (Cressie, 2015). The
expanding domain asymptotics requires that the domain of spatial measurement
locations tends to infinity. It is similar to time series and more suited for
a situation where the measurements are on a grid. In the infill asymptotics,
the total domain is kept fixed, but the density of measurements in that domain
is allowed to increase. By proposing a general sampling framework, Hörmann and
Kokoszka (2013) established conditions for which sample mean and sample
covariance of the spatially distributed functional data are consistent. Their
framework, however, assumes that the complete trajectories of the functions
are known.

This study proposes a framework of inverse regression for functional data
observed at irregularly positioned spatial locations as well as sparsely
observed in time. The data consist of a scalar response $Y(\bm{s}_{j})$ and a
functional covariate $X(\bm{s}_{j})=\\{X(\bm{s}_{j};t),t\in\mathcal{I}\\}$,
observed at spatial points $\bm{s}_{j}\in\mathbb{R}^{2}$, $j=1,\ldots,n$. We
assume that each observed functional covariate has a location-specific random
process which is interpreted as the “nugget” effect. The nugget effects are
uncorrelated across functions at different spatial locations. Our contribution
is twofold, which is summarized as follows.

  1. [label=()]

  2. 1.

Methodologically, we consider local linear smoothing (Fan and Gijbels, 2018)
to estimate the inverse regression function $E(X(\bm{s};t)|Y(\bm{s}))$ and the
covariance function of $X(\bm{s};t)$ by borrowing the information across all
the subjects. The method is applicable to the unified weighting scheme
proposed in Zhang and Wang (2016), including non-dense (cannot attain the
$\sqrt{n}$ rate), dense (attain $\sqrt{n}$ rate with non-negligible bias) and
ultra-dense (attain $\sqrt{n}$ rate without asymptotic bias) categories of
functions. The proposed approach allows both the density of observations and
the domain to increase, and therefore retains the advantages of both spatial
data frameworks.

  3. 2.

Theoretically, we identify the sufficient conditions under which the spatial-
temporal covariance, the inverse regression function and the covariance
function are consistent in Theorems 4.1, 4.5 and 4.9, respectively. As they
are derived for irregular spatial locations, each of these results is new to
the literature and hence is of independent interest. Our asymptotic results
include both pointwise rates and $L_{2}$ rates which are new to the
literature. It is also worth mentioning that the theoretical results of
spatially correlated data are not trivial extensions from the independent
case. For example, Lemma S.1.5 in the supplement Chatla and Liu (2024) proves
a uniform convergence for kernel estimators with spatially correlated
observations, which makes use Bradley’s coupling lemma (Lemma 1.2 in Bosq,
2012). Applying the uniform convergence results, we give the convergence rate
for the e.d.r. directions in Theorem 4.10. The rate for e.d.r. direction
depends on the convergence rates of spatio-temporal covariance, inverse
regression function, and the number of basis functions used for the truncation
of covariance operator.

The remainder of the paper is organized as follows. We outline the inverse
regression procedure for spatially distributed functional data in Section 2.
The formulation of the domain-expanding infill framework is included in
Section 3, which also covers the estimation of covariance operators and e.d.r.
directions under local linear smoothing framework. In Section 4, we outline
the conditions that are required for our asymptotic study and derive the
convergence rates for covariance operators and e.d.r. directions. Due to space
constraints, the proofs of the main results, a validation of the theoretical
findings through simulations, and a real data analysis are deferred to the
supplement Chatla and Liu (2024).

Notation: For functions $f(x)$ and $g(x,y)$ defined on a compact set, let
$\|f\|_{L_{2}}=\left[\int f^{2}(x)dx\right]^{1/2}$ and
$\|g\|_{HS}=\left[\int\int g^{2}(x,y)dxdy\right]^{1/2}$. For any vector
$\bm{x}$, denote its $L_{2}$ norm as $\|\bm{x}\|$. For any positive sequences
$\\{a_{n}\\}$ and $\\{b_{n}\\}$, we write $a_{n}\lesssim b_{n}$ if
$a_{n}/b_{n}$ is bounded by a constant, and write $a_{n}\asymp b_{n}$ if
$C_{1}\leq a_{n}/b_{n}\leq C_{2}$ for some $C_{1},C_{2}>0$.

##  2 Inverse regression for spatial function

In this section, we formulate the inverse regression framework for spatially
distributed functional data. We assume
$\\{(Y(\bm{s}),X(\bm{s})),\bm{s}\in\mathcal{D}_{n}\subseteq\mathbb{R}^{2}\\}$
is a stationary random field with $Y(\bm{s})$ being a scalar response variable
and a functional covariate $X(\bm{s})$ taking values in $L^{2}(\mathcal{I})$
where $\mathcal{I}\subset\mathbb{R}$. The space $L^{2}(\mathcal{I})$ is a
collection of square-integrable functions on $\mathcal{I}$. The value of the
function $X(\bm{s})$ at $t\in\mathcal{I}$ is denoted by $X(\bm{s};t)$. We
consider the following model

$\displaystyle Y(\bm{s})$
$\displaystyle=f(\langle\beta_{1}(\bm{s}),X(\bm{s})\rangle,\ldots,\langle\beta_{K}(\bm{s}),X(\bm{s})\rangle,\epsilon(\bm{s})),$
(2.1)

where the $K$ unknown functions $\beta_{1}(\cdot),\ldots,\beta_{K}(\cdot)$ are
in $L^{2}(\mathcal{I})$, $f$ is an arbitrary unknown function on
$\mathbb{R}^{K+1}$, $\epsilon(\bm{s})$ is a random error that is independent
of $X(\bm{s})$, and the inner product is defined as

$\displaystyle\langle\beta_{K}(\bm{s}),X(\bm{s})\rangle$
$\displaystyle=\int_{\mathcal{I}}\beta_{K}(\bm{s};t)X(\bm{s};t)dt.$

Assuming that the spatial dependence is second-order stationary and isotropic,
we write the covariance function of $X(\bm{s};t)$ as

$\displaystyle R(\|\bm{s}_{1}-\bm{s}_{2}\|;t_{1},t_{2})$
$\displaystyle=\text{cov}\\{X(\bm{s}_{1};t_{1}),X(\bm{s}_{2};t_{2})\\}$ (2.2)

for any
$(\bm{s}_{1},t_{1}),~{}(\bm{s}_{2},t_{2})\in\mathcal{D}_{n}\times\mathcal{I}$.
The covariance operator $R(\cdot,\cdot,\cdot)$ exists and bounded given that

$\displaystyle E\|X(\bm{s})\|_{L_{2}}^{4}<\infty\qquad\text{ for all
}\bm{s}\in\mathcal{D}_{n},$ (2.3)

where the $L_{2}$ norm is defined as

$\displaystyle\|X(\bm{s})\|_{L_{2}}$
$\displaystyle:=\left\\{\int_{\mathcal{I}}X^{2}(\bm{s};t)dt\right\\}^{1/2}.$

Similarly, we define the conditional covariance operator

$\displaystyle
R_{e}(\|\bm{s}_{1}-\bm{s}_{2}\|;t_{1},t_{2})=\text{cov}\\{E(X(\bm{s}_{1};t_{1})|Y(\bm{s}_{1})),E(X(\bm{s}_{2};t_{2})|Y(\bm{s}_{2}))\\},$
(2.4)

for any
$(\bm{s}_{1},t_{1}),~{}(\bm{s}_{2},t_{2})\in\mathcal{D}_{n}\times\mathcal{I}$.
We note that, $R(0;t_{1},t_{2})$ and $R_{e}(0;t_{1},t_{2})$ are the variance
operators of $X(\bm{s})$ and $E(X(\bm{s})|Y(\bm{s}))$, respectively.

Following Li (1991) and Ferré and Yao (2003), we obtain that

$\displaystyle E(\langle
b(\bm{s}),X(\bm{s})\rangle|\langle\beta_{1}(\bm{s}),X(\bm{s})\rangle,\ldots,\langle\beta_{K}(\bm{s}),X(\bm{s})\rangle)\text{
is linear in }$
$\displaystyle\qquad\langle\beta_{1}(\bm{s}),X(\bm{s})\rangle,\ldots,\langle\beta_{K}(\bm{s}),X(\bm{s})\rangle\text{
for any direction $b(\bm{s})$ in }L^{2}(\mathcal{I}).$ (2.5)

In other words, under model (2.1), the functions
$\langle\beta_{j}(\bm{s}),X(\bm{s})\rangle$, for $j=1,\ldots,K$, summarize all
the information contained in the functional variable $X(\bm{s})$ for
predicting $Y(\bm{s})$. The individual coefficient functions
$\beta_{j}(\bm{s})$, for $j=1,\ldots,K$, are not identifiable since the link
function $f(\cdot)$ in (2.1) is unknown. However, our interest lies in
estimating the effective dimension reduction (e.d.r.) space which is spanned
by the individual coefficient functions and is identifiable. Moreover, the
operator $R_{e}(0)$ is degenerate in any direction $R(0)-$ orthogonal to the
e.d.r. space. This implies we can recover the basis of e.d.r. space through
the $R(0)-$ orthogonal eigenvectors of $R_{e}(0)$ associated with the $K$
largest eigenvalues

$\displaystyle R_{e}(0)\beta_{j}(\bm{s})$
$\displaystyle=\lambda_{j}(\bm{s})R(0)\beta_{j}(\bm{s}),$ (2.6)

where $\beta_{i}^{{}^{\prime}}(\bm{s})R_{e}(0)\beta_{j}(\bm{s})=1$, if $i=j$,
and $0$ otherwise.

Based on (2.6), we compute the coefficient functions by performing a spectral
decomposition of the operator $R^{-1}(0)R_{e}(0)$ or the operator
$R^{-1/2}(0)R_{e}(0)R^{-1/2}(0)$, provided the orthogonal eigenvectors have
norm one. However, the compact covariance operator $R(0)$ is not invertible
for functional data. To overcome this limitation, we restrict the operator to
a smaller domain as in Jiang, Yu and Wang (2014). Note that under the
assumption that $E(\|X(\bm{s})\|_{L_{2}}^{4})<\infty$, $R(0)$ is a self-
adjoint, positive semidefinite, and Hilbert-Schmidt operator. Therefore, there
exists an orthonormal basis $\\{\pi_{i}(\cdot)\\}_{i=1}^{\infty}$ in
$L^{2}(\mathcal{I})$ such that $X(\bm{s};t)$ can be expressed as

$\displaystyle X(\bm{s};t)$
$\displaystyle=\mu(t)+\sum_{j=1}^{\infty}A_{j}(\bm{s})\pi_{j}(t),$

where $\mu(t)=E(X(\bm{s};t))$ is the mean function, and $\\{A_{j}(\bm{s})\\}$
are assumed to be second-order stationary and uncorrelated random fields with
$E(A_{j}(\bm{s}))=0$ and $E(A_{j}^{2}(\bm{s}))=\xi_{j}$. Here,
$\xi_{j}^{\prime}s$ are the eigenvalues of

$\displaystyle R(0;t_{1},t_{2})$
$\displaystyle=\sum_{i=1}^{\infty}\xi_{i}\pi_{i}(t_{1})\pi_{i}(t_{2}),$

satisfying $\xi_{1}\geq\xi_{2}\geq\ldots\geq\xi_{j}\geq\cdots\geq 0$ with
corresponding eigenfunctions $\pi_{j}(\cdot)$. By employing the arguments
analogous to Jiang, Yu and Wang (2014), we can show that when the following
condition is satisfied,

$\displaystyle\sum_{i,j=1}^{\infty}\frac{E^{2}\\{E(A_{i}(\bm{s})|Y(\bm{s}))E(A_{j}(\bm{s})|Y(\bm{s}))\\}}{\xi_{i}^{2}\xi_{j}}<\infty,$
(2.7)

then $R^{-1/2}(0)R_{e}(0)R^{-1/2}(0)$ is well-defined on the range space of
$R^{1/2}(0)$. This means for any $\mathfrak{f}\in\mathfrak{R}_{R^{1/2}(0)}$
with

$\displaystyle\mathfrak{R}_{R^{1/2}(0)}:=\left\\{\mathfrak{f}\in
L_{2}(\mathcal{I}):\sum_{i}\xi^{-1}|\langle\mathfrak{f},\pi_{i}\rangle|^{2}<\infty,\mathfrak{f}\perp\text{ker}(R(0))\right\\},$

it follows that

$\displaystyle
R^{-1/2}(0)\mathfrak{f}=\sum_{i=1}^{\infty}\xi^{-1/2}\langle\mathfrak{f},\pi_{i}\rangle\pi_{i}.$

While this solves the theoretical difficulty, often the estimation of
$R^{-1/2}(0)$ requires some regularization (Ferré and Yao, 2003; Jiang, Yu and
Wang, 2014). Let
$\widetilde{R}^{-1/2}(0)=R^{-1/2}(0)|_{\mathfrak{R}_{R^{1/2}(0)}}$ be the
inverse operator defined on the range space of $R^{1/2}(0)$. The directions we
obtain from $\widetilde{R}^{-1/2}(0)R_{e}(0)\widetilde{R}^{-1/2}(0)$ are still
in the e.d.r. space because $\mathfrak{R}_{R^{1/2}(0)}$ is a subspace of
$L_{2}(\mathcal{I})$.

Let
$R^{-1/2}(0)R_{e}(0)R^{-1/2}(0)(t_{1},t_{2})=\sum_{i=1}^{\infty}\lambda_{i}\eta_{i}(t_{1})\eta_{i}(t_{2})$.
The arguments analogous to Jiang, Yu and Wang (2014) can be used to show that
this is a well-defined operator given (2.7). Furthermore, the e.d.r.
directions $\beta_{j}=R^{-1/2}\eta_{j}$, for $j=1,\ldots,K$, are well-defined.

##  3 Sampling plan and methodology

This section discusses the sampling plan for spatial sites and for functional
data before outlining the procedure for estimating the coefficient functions
in model (2.1).

###  3.1 Domain-expanding infill framework

In general, there are two different paradigms for asymptotics (or the data
collection) in spatial statistics: expanding domain asymptotics, where the
domain can be expanded while maintaining at least roughly constant the
distance between two neighboring observations (see, e.g., Quenouille (1949);
Matérn (2013); Dalenius, Hájek and Zubrzycki (1961)), and infill or fixed
domain asymptotics, where the domain is fixed but the distance between
neighboring observations may tend to zero (see, e.g., Traub and Werschulz
(1998); Novak (2006)). The estimators may achieve consistency under expanding
domain asymptotics, but they may not do so under fixed domain asymptotics,
which is a noteworthy distinction between these two paradigms. Despite these
differences, each of these frameworks for spatial statistics has a large body
of research; for references, see Lu and Tjøstheim (2014). Similar sampling
frameworks are also taken into account for functional data (see, for example,
Delicado et al. (2010); Chouaf and Laksaci (2012); Hörmann and Kokoszka
(2013); Laksaci, Rachdi and Rahmani (2013)).

In this study, the domain-expanding infill (DEI) framework by Lu and Tjøstheim
(2014) is taken into consideration. The traditional infill and domain-
expanding asymptotics have disadvantages. While the former cannot guarantee
the estimators to be consistent as showed by Lahiri (1996) and Zhang (2004)
the latter has a less attractive assumption, at least in some applications,
that the distance between neighboring observations does not tend to zero. The
DEI asymptotics framework overcomes these drawbacks and is defined as

$\displaystyle\delta_{n}=\underset{1\leq j\leq n}{\max}\delta_{j,n}\rightarrow
0,\text{ with }\delta_{j,n}=\min\\{\|\bm{s}_{i}-\bm{s}_{j}\|:1\leq i\leq
n,i\neq j\\},$ (3.1)

which means that the distance between neighboring observations tends to $0$,
as $n\rightarrow\infty$, and

$\displaystyle\Delta_{n}=\underset{1\leq j\leq
n}{\min}\Delta_{j,n}\rightarrow\infty,\text{ with
}\Delta_{j,n}=\max\\{\|\bm{s}_{i}-\bm{s}_{j}\|:1\leq i\leq n,i\neq j\\},$
(3.2)

which implies that the domain at each location is expanding to $\infty$, as
$n\rightarrow\infty$, where $\|\cdot\|$ is the usual Euclidean norm in
$\mathbb{R}^{2}$. The existing frameworks may be described using the above
terminology. The traditional infill asymptotics require (3.1) but $\max_{1\leq
j\leq n}\Delta_{j,n}\leq C_{0}<\infty$, while the traditional domain-expanding
asymptotics require (3.2) but $\min_{1\leq j\leq n}\delta_{j,n}\geq c_{0}>0$
for all $n$. By (3.1), it is clear that the sampling locations $\bm{s}_{i}$
are non-stochastic. There are also some related works of stochastic sampling
designs (Zhang and Li (2022) and Kurisu (2022a) and the references therein)
combining both infill and increasing domain.

In many situations, the DEI framework may be natural as a result of the data
structure. For example, socio-economic data are often collected over
individual cities and suburbs: on one hand, more cities or suburbs are taken
into account which means expanding the domain; on the other hand, more
observations are collected within each city or suburb, necessitating the use
of an infill asymptotic framework.

###  3.2 Methodology

Assume we observe data at locations
$\\{\bm{s}_{i}\in\mathcal{D}_{n},i=1,\ldots,n\\}$ that are allowed to be
irregularly positioned. Suppose we observe

$\displaystyle Z(\bm{s}_{i};T_{ij})$
$\displaystyle=X(\bm{s}_{i};T_{ij})+U_{i}(T_{ij})+e_{ij},\qquad
i=1,\ldots,n,j=1,\ldots,N_{i},$ (3.3)

where

$\displaystyle X(\bm{s};t)$ $\displaystyle=\mu(t)+V(\bm{s},t),$

with $V(\bm{s},t)$ being the stochastic part of $X(\bm{s};t)$ with
$EV(\bm{s},t)=0$ and $e_{ij}$’s are i.i.d. copies of a random error $e$ with
mean zero and finite variance $\sigma^{2}$, and $U_{i}(t)$ is a mean zero
temporal process called as functional nugget effect. Assume that
$V(\cdot,\cdot)$, $U_{i}(\cdot)$, and $e$ are mutually independent. As
mentioned in Zhang and Li (2022), the functional nugget effects $U_{i}(\cdot)$
specify local variations that are not correlated with neighbor functions. We
denote the covariance between nugget effects as
$\Lambda(t_{1},t_{2})=\text{cov}\\{U_{i}(t_{1}),U_{i}(t_{2})\\}$. The
observation time points $\\{T_{ij}\\}$ are i.i.d. copies of a random variable
$T$ which is defined on the compact interval $\mathcal{I}$. We consider the
local linear smoothing (Fan and Gijbels, 2018) approach for the estimation of
e.d.r. directions. For convenience in notation, we denote
$Z_{ij}:=Z(\bm{s}_{i};T_{ij})$, $X_{ij}:=X(\bm{s}_{i};T_{ij})$,
$Y_{i}=Y(\bm{s}_{i})$, $V_{ij}=V(\bm{s}_{i};T_{ij})$, and
$U_{ij}:=U_{i}(T_{ij})$.

In the following, we discuss the estimation of the required covariance
matrices.

Estimation of $R(0;t_{1},t_{2})$:  Denote $K_{h}(u)=h^{-1}K_{1}(u)$ where
$K_{1}(\cdot)$ is a univariate kernel function defined in Condition (III)
below. Note that the covariance
$R(0;t_{1},t_{2})=\text{cov}\\{X(\bm{s};t_{1}),~{}X(\bm{s};t_{2})\\}$. We
first estimate the mean function $\mu(t):=E(X(\bm{s};t))$ using a one-
dimensional local linear smoother applied to the pooled data
$(Z_{ij},T_{ij})$, $1\leq i\leq n$, $1\leq j\leq N_{i}$. Therefore, the mean
estimator is $\widehat{\mu}(t)=\widehat{\alpha}_{0}$ with

$\displaystyle(\widehat{\alpha}_{0},\widehat{\alpha}_{1})$
$\displaystyle=\underset{(\alpha_{0},\alpha_{1})\in\mathbb{R}^{2}}{\text{argmin}}\sum_{i=1}^{n}w_{i}\sum_{j=1}^{N_{i}}\\{Z_{ij}-\alpha_{0}-\alpha_{1}(T_{ij}-t)\\}^{2}K_{h_{\mu}}\left(T_{ij}-t\right).$
(3.4)

Here $w_{1},\ldots,w_{n}$ are weights satisfying $\sum_{i=1}^{n}w_{i}N_{i}=1$.
The following two weighting schemes are commonly used in the literature (Yao,
Müller and Wang, 2005; Li and Hsing, 2010; Zhang and Wang, 2016).

  1. [label=()]

  2. 1.

Equal-weight-per-observation (OBS): $w_{i}=1/N_{S,0}$ with
$N_{S,0}=\sum_{i=1}^{n}N_{i}$;

  3. 2.

Equal-weight-per-subject (SUBJ): $w_{i}=1/nN_{i}$.

As suggested by Zhang and Wang (2016), using a weighing scheme that is neither
OBS nor SUBJ could achieve a better rate. Related works include a mixed scheme
of OBS and SUBJ in Zhang and Wang (2016) and a bandwidth-dependent scheme in
Zhang and Wang (2018). For simplicity, we focus on the OBS and SUBJ schemes.

We assume $N_{i}\geq 2$ to estimate $R(0;t_{1},t_{2})$. Let
$L_{b}(u)=b^{-1}L(u/b)$ where $L(\cdot)$ is a kernel function defined similar
to $K_{1}(\cdot)$ with bandwidth $b$. We apply a three-dimensional local
linear smoother to the cross-products
$C_{ij,i^{\prime}j^{\prime}}:=[Z_{ij}-\widehat{\mu}(T_{ij})][Z_{i^{\prime}j^{\prime}}-\widehat{\mu}(T_{i^{\prime}j^{\prime}})]$,
$1\leq i,i^{\prime}\leq n$, $1\leq j\leq N_{i}$, $1\leq j^{\prime}\leq
N_{i^{\prime}}$, and attach weight $v_{i,i^{\prime}}$ to each
$C_{ij,i^{\prime}j^{\prime}}$ for the $(i,i^{\prime})$th subject pair such
that $\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}N_{i}N_{i^{\prime}}v_{i,i^{\prime}}=1$. Specifically, we consider the
estimator $\widehat{R}(\|\bm{s}_{0}\|;t_{1},t_{2})=\widehat{\alpha}_{0}$ where

$\displaystyle(\widehat{\alpha}_{0},\widehat{\alpha}_{1},\widehat{\alpha}_{2},\widehat{\alpha}_{3})$
(3.5) $\displaystyle=$
$\displaystyle\underset{(\alpha_{0},\alpha_{1},\alpha_{2},\alpha_{3})\in\mathbb{R}^{4}}{\text{arg
min}}\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v_{i,i^{\prime}}\sum_{j=1}^{N_{i}}\sum_{j^{\prime}=1}^{N_{i^{\prime}}}$
$\displaystyle\qquad\bigg{\\{}C_{ij,i^{\prime}j^{\prime}}-\alpha_{0}-\alpha_{1}(T_{ij}-t_{1})-\alpha_{2}(T_{i^{\prime}j^{\prime}}-t_{2})-\alpha_{3}(\|\bm{s}_{i}-\bm{s}_{i^{\prime}}\|-\|\bm{s}_{0}\|)\bigg{\\}}^{2}$
$\displaystyle\qquad\qquad\times
K_{h_{c}}\left(T_{ij}-t_{1}\right)K_{h_{c}}\left(T_{i^{\prime}j^{\prime}}-t_{2}\right)L_{b}(\|\bm{s}_{i}-\bm{s}_{i^{\prime}}\|-\|\bm{s}_{0}\|).$

Similar to Zhang and Wang (2016), we consider weights
$v_{i,i^{\prime}}=1/\sum_{i=1}^{n}\sum_{i^{\prime}\neq i}N_{i}N_{i^{\prime}}$
for the OBS scheme and weights
$v_{i,i^{\prime}}=1/[n(n-1)N_{i}N_{i^{\prime}}]$ for the SUBJ scheme. Our
asymptotic theory regards $N_{i}$, $w_{i}$ and $v_{i,i^{\prime}}$ as fixed
quantities that are allowed to vary over $n$. For random $N_{i}$, the theory
proceeds as if conditional on the values of $N_{i}$. The asymptotic properties
of $\widehat{R}(\|\bm{s}_{0}\|;t_{1},t_{2})$ are discussed in Theorem 4.1 and
Lemma 4.2 below.

Estimation of $\Gamma(t_{1},t_{2})$ and $\Lambda(t_{1},t_{2})$: Because of the
independence between $X(\bm{s};t)$ and the functional nugget effect $U(t)$, we
obtain

$\displaystyle\Gamma(t_{1},t_{2}):=\text{cov}\\{X(\bm{s};t_{1}),X(\bm{s};t_{2})\\}=R(0;t_{1},t_{2})+\Lambda(t_{1},t_{2}),$

for $t_{1}\neq t_{2}$. This motivates us to apply a two-dimensional local
linear smoother to the cross-products
$C_{ijk}:=[Z_{ij}-\widehat{\mu}(T_{ij})][Z_{ik}-\widehat{\mu}(T_{ik})]$,
$1\leq i\leq n$, $1\leq j\neq k\leq N_{i}$ to estimate $\Gamma(t_{1},t_{2})$
and attach weight $v_{i}$ to each $C_{ijk}$ for the $i$th subject such that
$\sum_{i=1}^{n}N_{i}(N_{i}-1)v_{i}=1$. Specifically, the estimator
$\widehat{\Gamma}(t_{1},t_{2})=\widehat{\alpha}_{0}$ is considered where

$\displaystyle(\widehat{\alpha}_{0},\widehat{\alpha}_{1},\widehat{\alpha}_{2})=\underset{(\alpha_{0},\alpha_{1},\alpha_{2})\in\mathbb{R}^{3}}{\text{arg
min}}\sum_{i=1}^{n}v_{i}\sum_{1\leq k\neq j\leq
N_{i}}\left\\{C_{ijk}-\alpha_{0}-\alpha_{1}(T_{ij}-t_{1})-\alpha_{2}(T_{ik}-t_{2})\right\\}^{2}$
$\displaystyle\times
K_{h_{c}}\left(T_{ij}-t_{1}\right)K_{h_{c}}\left(T_{ik}-t_{2}\right).$ (3.6)

As in Zhang and Wang (2016), weights $v_{i}=1/\sum_{i=1}^{n}N_{i}(N_{i}-1)$
yield the OBS scheme and weights $v_{i}=1/[nN_{i}(N_{i}-1)]$ yield the SUBJ
scheme. Here, a natural estimator for the functional nugget effect is

$\displaystyle\widehat{\Lambda}(t_{1},t_{2})$
$\displaystyle=\widehat{\Gamma}(t_{1},t_{2})-\widehat{R}(0;t_{1},t_{2}).$

The asymptotic properties of $\widehat{\Gamma}(t_{1},t_{2})$ are discussed in
Theorem 4.9 below. Given $\widehat{\Gamma}(t_{1},t_{2})$, we can estimate
$\sigma^{2}$ by proceeding similar to Yao, Müller and Wang (2005). We do not
discuss the properties of $\widehat{\sigma}^{2}$ here as it is not the main
focus of our study.

Estimation of $R_{e}(0,t_{1},t_{2})$:  Let us introduce
$K_{H}(\bm{u})=|\bm{H}|^{-1/2}K_{2}(\bm{H}^{-1/2}\bm{u})$ where $K_{2}(\cdot)$
is a two-dimensional kernel defined in Condition (III) and
$\bm{H}=\text{diag}(h_{t}^{2},h_{y}^{2})$ with $h_{y}$ and $h_{t}$ are the
bandwidths for $Y_{i}$ and $T_{ij}$, respectively. Let
$m(t,y)=E(X(\bm{s};T_{ij}=t)|Y(\bm{s})=y)$. The function $m(t,y)$ is assumed
to be smooth since we employ local linear smoothing for estimation. For
estimation of $m(t,y)$, we apply a two-dimensional smoothing method to the
pooled sample $\\{Z_{ij}\\}$ over $\\{T_{ij},Y_{i}\\}$ for $i=1,\ldots,n$,
$j=1,\ldots,N_{i}$. Therefore, we consider the estimator
$\widehat{m}(t,y)=\widehat{\alpha}_{0}$ where

$\displaystyle(\widehat{\alpha}_{0},\widehat{\alpha}_{1},\widehat{\alpha}_{2})$
$\displaystyle=\underset{(\alpha_{0},\alpha_{1},\alpha_{2})\in\mathbb{R}^{3}}{\text{argmin}}\sum_{i=1}^{n}w_{i}\sum_{j=1}^{N_{i}}\left\\{Z_{ij}-\alpha_{0}-\alpha_{1}(T_{ij}-t)-\alpha_{2}(Y_{i}-y)\right\\}^{2}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\times
K_{H}\left((T_{ij}-t,Y_{i}-y)^{T}\right),$ (3.7)

where the weight, $w_{i}$ is attached to each observation for the $i$-th
subject for $i=1,\ldots,n$, with $\sum_{i=1}^{n}N_{i}w_{i}=1$. Recall that
$R_{e}(0)=\text{cov}(E(X(t)|Y))$. Let $\mathfrak{B}$ be a compact subset of
$\mathbb{R}$. We estimate $R_{e}(0)$ empirically as follows:

$\displaystyle\widehat{R}_{e}(0,t_{1},t_{2})$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\widehat{m}(t_{1},Y_{i})\widehat{m}(t_{2},Y_{i})\mathbbm{1}_{\\{Y_{i}\in\mathfrak{B}\\}}-\frac{1}{n}\sum_{i=1}^{n}\widehat{m}(t_{1},Y_{i})\mathbbm{1}_{\\{Y_{i}\in\mathfrak{B}\\}}\frac{1}{n}\sum_{j=1}^{n}\widehat{m}(t_{2},Y_{j})\mathbbm{1}_{\\{Y_{j}\in\mathfrak{B}\\}},$

where $\mathbbm{1}_{\\{Y_{j}\in\mathfrak{B}\\}}$ is an indicator function that
takes value one if $Y_{j}\in\mathfrak{B}$ and zero otherwise. In practice, we
may choose $\mathfrak{B}=[y_{\alpha},y_{1-\alpha}]$ where $y_{\alpha}$ is the
$\alpha$-th percentile of $Y_{i}$ values. In general, we can take $\alpha$
equal to 5% or 10%. The asymptotic properties of the estimator
$\widehat{m}(t,y)$ and $\widehat{R}_{e}(0)$ are discussed in Theorem 4.5 and
Lemma 4.6 below.

Estimation of the e.d.r. directions $\beta_{j}$ for $j=1,\ldots,K$ and the
link function $f(\cdot)$:  We calculate the e.d.r. directions using the eigen-
analysis of $R^{-1/2}(0)R_{e}(0)R^{-1/2}(0)$. Let $\widehat{\eta}_{j}$ be the
estimator of the standardized e.d.r. direction $\eta$ and
$\widehat{\beta}_{j}$ be the estimator of $\beta_{j}$. According to Jiang, Yu
and Wang (2014), the following procedures are used to perform an eigen-
analysis on a discrete grid in practice:

  1. [label=()]

  2. 1.

on an equally-spaced grid $\\{t_{1},\ldots,t_{p}\\}$, compute the $p\times p$
matrices $\widehat{\bm{R}}_{e}$ and $\widehat{\bm{R}}$,

  3. 2.

obtain $\widehat{\bm{R}}^{-1/2}$ from the singular value decomposition of
$\widehat{\bm{R}}$,

  4. 3.

obtain the first $K$ eigenfunctions
$\widehat{\bm{\eta}}_{1},\ldots,\widehat{\bm{\eta}}_{K}$ corresponding to the
largest $k$ eigenvalues $\widehat{\lambda}_{1},\ldots,\widehat{\lambda}_{K}$
from the eigen-analysis of
$\widehat{\bm{R}}^{-1/2}\widehat{\bm{R}}_{e}\widehat{\bm{R}}^{-1/2}$,

  5. 4.

compute the e.d.r. directions
$\widehat{\bm{\beta}}_{j,p}=\widehat{\bm{R}}^{-1/2}\widehat{\bm{\eta}}_{j}$,
for $j=1,\ldots,K$.

The link function $f(\cdot)$ is estimated, assuming $K$ is significantly
smaller than $p$, by using standard nonparametric smoothing techniques to the
indices $\langle\beta_{j}(t),Z_{i}(t)\rangle$ for $j=1,\ldots,K$, where
$Z_{i}(\cdot)$ is the observed function for the $i$th location,
$i=1,\ldots,n$. The traditional nonparametric methods, however, suffer from
the curse of dimensionality when $K$ is not small. In this situation, we may
need to use methods like Multi-layer Neural Network (Ferré and Yao, 2005),
which are not affected by the curse of dimensionality.

##  4 Asymptotic results

In this section, we derive the convergence rates for the covariance operators
and the e.d.r. directions. Standard smoothing assumptions are made on the
second derivatives of $R$ and $R_{e}$ since their estimators $\widehat{R}$ and
$\widehat{R}_{e}$ are obtained by employing local linear smoothing.

We denote $\bm{T}_{i}=(T_{i1},\ldots,T_{iN_{i}})^{T}$ and
$\bm{Z}_{i}=(Z_{i1},\ldots,Z_{iN_{i}})^{T}$ for $i=1,\ldots,n$. Let
$Y_{i}=Y(\bm{s}_{i})$ for $i=1,\ldots,n$. For convenience, we summarize the
main assumptions here. Conditions (I) and (II) are related to the random field
itself.

Condition (I) (spatial process):

  1. [label=()]

  2. 1.

The random field
$\\{(Y_{i},\bm{T}_{i},\bm{Z}_{i}):\bm{s}_{i}\in\mathcal{D}_{n},i=1,\ldots,n\\}$
is strictly stationary.

  3. 2.

The time points $\\{T_{ij}:i=1,\ldots,n,j=1,\ldots,N_{i}\\}$, are i.i.d.
copies of a random variable $T$ defined on $\mathcal{I}$, say $[0,1]$. The
density $f_{T}(\cdot)$ of $T$ is bounded away from zero and infinity. The
second derivative of $f_{T}(\cdot)$, denoted by $f^{(2)}_{T}(\cdot)$, is
bounded.

    1. [label=()]

    2. (a)

$V$, $U_{i}$ and $e$ are mutually independent, and are independent of $T$.

  4. 3.

Assume that $\mu^{(2)}(t)$ exists and bounded on $\mathcal{I}$.

  5. 4.

    1. [label=()]

    2. (a)

For all $1\leq i\neq i^{\prime}\leq n$, the vectors $(T_{ij},Y_{i})$ and
$(T_{i^{\prime}j^{\prime}},Y_{i^{\prime}})$, $1\leq j\leq N_{i}$, $1\leq
j^{\prime}\leq N_{i^{\prime}}$, admit a joint density $g_{i,i^{\prime}}$ and a
marginal density $g$; further,
$|g_{i,i^{\prime}}(t,y,t^{\prime},y^{\prime})-g(t,y)g(t^{\prime},y^{\prime})|\leq
C$ for $t,t^{\prime},y,y^{\prime}\in\mathbb{R}$, where $C>0$ is some constant.
Let $g_{3}(t,t^{\prime},y)$ be the joint density of
$(T_{ij}=t,T_{ij^{\prime}}=t^{\prime},Y_{i}=y)$. We assume $g$, $g_{3}$ and
$m(t,y)$ have continuous and bounded second derivatives.

    3. (b)

The data points $Y_{i}$’s are independent of $U_{ij}$ and $e_{ij}$. Let
$\Psi(t,y)$ be the variance of $X_{ij}$ given $T_{ij}=t$ and $Y_{i}=y$ and
$\Psi(t,y,t^{\prime})$ be the covariance of $X_{ij}$ and $X_{ij^{\prime}}$
given $T_{ij}=t$, $T_{ij^{\prime}}=t^{\prime}$ and $Y_{i}=y$. We assume that
$\Psi$ has continuous and bounded second derivatives. Let $\Lambda(t)$ be the
variance of $U_{ij}$ given $T_{ij}$ and $\Lambda(t_{1},t_{2})$ be the
covariance of $U_{ij}$ and $U_{ik}$ given $T_{ij}$ and $T_{ik}$. We assume
$\Lambda$ has continuous and bounded second derivatives.

    4. (c)

$\sup_{\bm{s}}E[\sup_{t}|V(\bm{s},t)|^{2+\varrho}]<\infty$,
$E[\sup_{t}|U(t)|^{2+\varrho}]<\infty$, and $E|e|^{2+\varrho}<\infty$ for some
$\varrho>0$.

  6. 5.

    1. [label=()]

    2. (a)

Let $\Xi(t_{1},t_{2},t_{1}^{\prime},t_{2}^{\prime})$ be the covariance of
$[X_{ij}-\mu(T_{ij})]$ $[X_{ik}-\mu(T_{ik})]$ and
$[X_{ij^{\prime}}-\mu(T_{ij^{\prime}})]$
$[X_{ik^{\prime}}-\mu(T_{ik^{\prime}})]$ and
$\Lambda_{1}(t_{1},t_{2},t_{1}^{\prime},t_{2}^{\prime})$ be the covariance of
$U_{ij}U_{ik}$ and $U_{ij^{\prime}}U_{ik^{\prime}}$ given $T_{ij}=t_{1}$,
$T_{ik}=t_{2}$, $T_{ij^{\prime}}=t_{1}^{\prime}$ and
$T_{ik^{\prime}}=t_{2}^{\prime}$. We assume $\Xi$, $\Lambda_{1}$ and
$\Gamma(t_{1},t_{2})$ have continuous and bounded second derivatives.

    3. (b)

$\sup_{\bm{s}}E[\sup_{t}|V(\bm{s},t)|^{4+\varrho_{c}}]<\infty$,
$E[\sup_{t}|U(t)|^{4+\varrho_{c}}]<\infty$, and $E|e|^{4+\varrho_{c}}<\infty$
with $\varrho_{c}/2=\varrho>0$.

  7. 6.

    1. (a)

Let
$\Xi_{R}(t_{1},t_{2},t_{1}^{\prime},t_{2}^{\prime};\|\bm{s}_{i}-\bm{s}_{i^{\prime}}\|)$
be the covariance of $[Z_{ij}-\mu(T_{ij})]$
$[Z_{i^{\prime}j^{\prime}}-\mu(T_{i^{\prime}j^{\prime}})]$ and
$[Z_{ik}-\mu(T_{ik})]$
$[Z_{i^{\prime}k^{\prime}}-\mu(T_{i^{\prime}k^{\prime}})]$ given
$T_{ij}=t_{1}$, $T_{i^{\prime}j^{\prime}}=t_{2}$, $T_{ik}=t_{1}^{\prime}$ and
$T_{i^{\prime}k^{\prime}}=t_{2}^{\prime}$ for $1\leq i\neq i^{\prime}\neq n$
and $1\leq j,k\leq N_{i}$, $1\leq j^{\prime},k^{\prime}\leq N_{i^{\prime}}$.
We assume that $\Xi_{R}$ and $R(s;t_{1},t_{2})$ have continuous and bounded
second derivatives.

Conditions (I)[(i), (ii)(a), (iii)(a)] are standard for the spatial data
generating process in the context of the problem under study. For example,
similar conditions are used by Masry (1986) in time series, by Tran (1990) and
Hallin, Lu and Tran (2004) in the gridded spatial data context, and by Lu and
Tjøstheim (2014) for irregular spatial measurements. These conditions are
required to ensure uniform consistency for joint probability densities across
different spatial sites. Conditions (I)[(ii)(b), (iii)(b)] are mild
assumptions on the covariance of the function. They are similar to Assumption
(A.6) in Jiang, Yu and Wang (2014) and needed for local linear smoothing.
Conditions (I)[(ii)(c), (iii)(c)] are typical in the spatial regression setup.
For example, they are similar to (A2) of Hallin, Lu and Tran (2004),
Assumption (3.2) in Gao, Lu and Tjøstheim (2006) and Conditions (C5)-(C7) in
Li and Hsing (2010). These conditions are required for Lemma S.1.1 in the
supplement Chatla and Liu (2024). Condition (iv) is similar to Assumption
(I)(iii) in Lu and Tjøstheim (2014) which is required for the spatio-temporal
covariance.

Condition (II) provides the assumptions on spatial mixing as well as spatial
sites, which include irregularly positioned observations. Let
$\mathcal{B}(\mathcal{S})$ be the Borel $\sigma$-field generated by
$\\{(Y_{i},\bm{T}_{i},\bm{Z}_{i}):\bm{s}_{i}\in\mathcal{S},i=1,\ldots,n\\}$
for any collection of sites $\mathcal{S}\subset\mathbb{R}^{2}$. For each
couple $\mathcal{S}^{\prime}$, $\mathcal{S}^{\prime\prime}$, let
$d(\mathcal{S}^{\prime},\mathcal{S}^{\prime\prime}):=\text{min}\\{\|\bm{s}^{\prime}-\bm{s}^{\prime\prime}\|:\bm{s}^{\prime}\in\mathcal{S}^{\prime},\bm{s}^{\prime\prime}\in\mathcal{S}^{\prime\prime}\\}$
be the distance between $\mathcal{S}^{\prime}$ and
$\mathcal{S}^{\prime\prime}$. Finally, let Card$(\mathcal{S})$ be the
cardinality of $\mathcal{S}$.

Condition (II) (mixing and sampling sites):

  1. [label=()]

  2. 1.

There exist a function $\phi$ such that $\phi(t)\rightarrow 0$ as
$t\rightarrow\infty$, and a symmetric function
$\psi:\mathbb{N}^{2}\rightarrow\mathbb{R}^{+}$ decreasing in each of its two
arguments such that

$\displaystyle\alpha\left(\mathcal{B}(\mathcal{S^{\prime}}),\mathcal{B}(\mathcal{S^{\prime\prime}})\right)$
$\displaystyle:=\sup\left\\{|P(AB)-P(A)P(B)|,A\in\mathcal{B}(\mathcal{S^{\prime}}),B\in\mathcal{B}(\mathcal{S^{\prime\prime}})\right\\}$
$\displaystyle\leq\psi(\text{Card}(\mathcal{S^{\prime}}),\text{Card}(\mathcal{S^{\prime\prime}}))\phi(d(\mathcal{S^{\prime}},\mathcal{S^{\prime\prime}})),$
(4.1)

for any
$\mathcal{S^{\prime}},\mathcal{S^{\prime\prime}}\subset\mathbb{R}^{2}$.
Moreover, the function $\phi$ is such that

$\displaystyle\underset{m\rightarrow\infty}{\lim}m^{\kappa}\sum_{j=m}^{\infty}j^{2}\\{\phi(j)\\}^{\varrho/(2+\varrho)}$
$\displaystyle=0$ (4.2)

for some constant $\kappa>$ and for the constant $\varrho$ defined in
Condition (I)[(iv), (v)]. In this study, we consider $\phi(j)\lesssim
j^{-\rho}$ for some constant $\rho>3(2+\varrho)/\varrho$ and
$\psi(i,j)\lesssim\min(i,j)$.

The observations are positioned as $\\{\bm{s}_{i},i=1,\ldots,n\\}$ for which
(3.1) and (3.2) hold true with $\underset{1\leq j\leq
n}{\min}\delta_{j,n}/\delta_{n}\geq c_{1}>0$ and $\underset{1\leq j\leq
n}{\max}\Delta_{j,n}/\Delta_{n}\leq C_{1}<\infty$ for all $n$. Suppose
$g_{\bm{S}}$ is the spatial density function (or sampling intensity function)
defined on $\mathbb{R}^{2}$ such that:

  1. [label=()]

  2. 2.

$\sum_{i=1}^{n}\mathbb{W}_{i}\bm{I}(\bm{s}_{i}\in
A)\rightarrow\int_{A}g_{\bm{S}}(\bm{s})d\bm{s}$ as $n\rightarrow\infty$, for
any measurable set $A\subset\mathbb{R}^{2}$ and for some weights
$\mathbb{W}_{i}>0$ and $\sum_{i=1}^{n}\mathbb{W}_{i}=1$.

  3. 3.

The spatial density $g_{\bm{S}}(\bm{s})$ is bounded and has second derivatives
that are continuous on $\mathbb{R}^{2}$.

  4. 4.

The following two integrals exist:

$\displaystyle A_{0}^{*}(\|\bm{s}_{0}\|)$
$\displaystyle=\|\bm{s}_{0}\|\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\bigg{[}\int_{0}^{2\pi}g_{\bm{S}}(u+\|\bm{s}_{0}\|\cos(\vartheta),v+\|\bm{s}_{0}\|\sin(\vartheta))d\vartheta\bigg{]}$
$\displaystyle\qquad\qquad\times g_{\bm{S}}(u,v)dudv>0,$

and

$\displaystyle A_{1}^{*}(\|\bm{s}_{0}\|)$
$\displaystyle=\|\bm{s}_{0}\|\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\bigg{[}\int_{0}^{2\pi}(\cos(\vartheta),\sin(\vartheta))$
$\displaystyle\qquad\times\bigg{[}\frac{\partial
g_{\bm{S}}(u+\|\bm{s}_{0}\|\cos(\vartheta),v+\|\bm{s}_{0}\|\sin(\vartheta))}{\partial\bm{s}}d\vartheta\bigg{]}g_{\bm{S}}(u,v)dudv>0.$

The first part of Condition (II)(i), the $\alpha$-mixing condition, is
standard for spatial data. For example, it is similar to (A4) in Hallin, Lu
and Tran (2004) and Assumption I(i) in Lu and Tjøstheim (2014). This condition
characterizes spatial dependence. If the condition (4.1) holds with $\psi=1$,
then the random field $\\{Y_{i},\bm{T}_{i},\bm{Z}_{i}\\}$ is called strongly
mixing. Many time series and stochastic processes are shown to be strongly
mixing, at least, in the serial case $d=1$. The last part of Condition (II)(i)
is similar to Assumption 2 in Hansen (2008) and Assumption E2 in Kurisu
(2022b) and needed for uniform convergence. Conditions (II)[(ii), (iii)] are
from Lu and Tjøstheim (2014). These are mild conditions on the spatial sites
and are needed for the derivation of asymptotic bias and variances of the
concerned estimators. Condition (II)(iv) is also from Lu and Tjøstheim (2014)
which assumes that the joint probability density is isotropic.

Condition (III) deals with the kernel functions that are used for estimating
both covariance operators. These are similar to Assumption (A.3) of Jiang, Yu
and Wang (2014).

Condition (III) (kernel functions):

  1. [label=()]

  2. 1.

Let $K_{2}(\cdot,\cdot)$ be a bivariate kernel function that is compactly
supported, symmetric and B-Lipschitz. Further, it is a kernel of order
$(|b_{1}|,|b_{2}|)$, that is,

$\displaystyle\sum_{\ell_{1}+\ell_{2}=\ell}\int\int
u^{\ell_{1}}v^{\ell_{2}}K_{2}(u,v)dudv=\Bigg{\\{}\begin{array}[]{ll}0,&0\leq\ell<|b_{2}|,\ell\neq|b_{1}|,\\\
(-1)^{|b_{1}|}|b_{1}|!,&\ell=|b_{1}|,\\\ \neq 0,&\ell=|b_{2}|.\end{array}$

Similarly, the univariate kernel function $K_{1}$ of order $(c_{1},c_{2})$ is
compactly supported, symmetric and B-Lipschitz, and defined as

$\displaystyle\int
u^{\ell}K_{1}(u)du=\Bigg{\\{}\begin{array}[]{ll}0,&0\leq\ell<c_{2},\ell\neq
c_{1},\\\ (-1)^{c_{1}}c_{1}!,&\ell=c_{1},\\\ \neq 0,&\ell=c_{2}.\end{array}$

We assume that $\int|K_{1}(u)|^{2+\varrho}du<\infty$,
$\int|uK_{1}(u)|^{2+\varrho}du<\infty$,
$\int|K_{2}(u,v)|^{2+\varrho}dudv<\infty$,
$\int|u^{2}K_{2}(u,v)|^{2+\varrho}dudv<\infty$, and
$\int|v^{2}K_{2}(u,v)|^{2+\varrho}dudv<\infty$. In this study, we set
$|b_{1}|=0$, $|b_{2}|=2$ for $K_{2}$, and $c_{1}=0,c_{2}=2$ for $K_{1}$.

  3. 2.

The kernel function $L(\cdot)$ satisfies that $\int L(s)ds=1$, $\int
sL(s)ds=0$, and $\mu_{L,2}=\int s^{2}L(s)ds<\infty$ and $\nu_{L,2}=\int
L^{2}(s)ds<\infty$. Moreover, it is Lipschitz continuous.

Before we state the conditions on bandwidths, define

$\displaystyle\gamma_{n1}$
$\displaystyle=nh^{-2}\text{max}^{2}\\{w_{i}N_{i}\\}^{2},$ (4.5)
$\displaystyle\gamma_{n2}$
$\displaystyle=nh^{-2}\text{max}^{2}\\{v_{i}N_{i}(N_{i}-1)\\},\text{ and }$
(4.6) $\displaystyle\gamma_{n3}$
$\displaystyle=n(n-1)^{2}h^{-2}\text{max}^{2}\\{v_{i,i^{\prime}}N_{i}N_{i^{\prime}}\\}.$
(4.7)

Let $\overline{N}_{1}=n^{-1}\sum_{i=1}^{n}N_{i}(N_{i}-1)$ and
$\overline{N}_{2}=n^{-1}(n-1)^{-1}\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}N_{i}N_{i^{\prime}}$. For the SUBJ and OBS weighting schemes, we obtain
$\gamma_{n1}=1/nh^{2}$ and
$\gamma_{n1}=(nh^{2})^{-1}\text{max}^{2}\\{N_{i}\\}/\overline{N}^{2}$.
Similarly, we obtain $\gamma_{n2}=1/nh^{2}$ for SUBJ scheme and
$\gamma_{n2}=(nh^{2})^{-1}\text{max}^{2}\\{N_{i}(N_{i}-1)\\}/\overline{N}_{1}^{2}$
for OBS scheme. Moreover, for SUBJ scheme $\gamma_{n3}=1/nh^{2}$ and for OBS
scheme
$\gamma_{n3}=(nh^{2})^{-1}\text{max}^{2}\\{N_{i}N_{i^{\prime}}\\}/\overline{N}_{2}^{2}$.
We now provide conditions related to bandwidths.

Condition (IV) (bandwidth):

  1. [label=()]

  2. 1.

Assume that the bandwidths $h_{t}=O(h)$ and $h_{y}=O(h)$. The bandwidth $h$
satisfies $h\rightarrow 0$, $nh^{2}\delta_{n}^{2}\rightarrow\infty$,
$\sum_{i=1}^{n}N_{i}w_{i}^{2}/h^{2}\rightarrow 0$,
$\sum_{i=1}^{n}N_{i}(N_{i}-1)w_{i}^{2}/h\rightarrow 0$ and

$\displaystyle\begin{split}\underset{n\rightarrow\infty}{\lim}\gamma_{n1}^{1+2/\kappa}\delta_{n}^{-2(1+2/\kappa)}h^{2\left(1-\frac{2\varrho}{\kappa(2+\varrho)}\right)}&=0,\text{
for }\kappa>\max\\{1,2\varrho/(2+\varrho)\\},\text{ and }\\\
\underset{n\rightarrow\infty}{\lim}\gamma_{n1}^{1+2/\kappa}\delta_{n}^{-2(1+2/\kappa)}h^{\frac{2\varrho}{4+\varrho}\left(1-\frac{2(4+\varrho)}{\kappa(2+\varrho)}\right)}&=0,\text{
for }\kappa>2(4+\varrho)/(2+\varrho),\end{split}$ (4.8)

where $\kappa$ and $\varrho>0$ defined in (4.2), $\gamma_{n1}$ is defined in
(4.5), and $\delta_{n}$ is defined in (3.1).

  3. 2.

The bandwidth $h_{\mu}$ satisfies $h_{\mu}\rightarrow 0$,
$nh_{\mu}\delta_{n}^{2}\rightarrow\infty$,
$\sum_{i=1}^{n}N_{i}w_{i}^{2}/h_{\mu}\rightarrow 0$,
$\sum_{i=1}^{n}N_{i}(N_{i}-1)w_{i}^{2}\rightarrow 0$, and (4.8) with $h^{2}$
replaced by $h_{\mu}$.

  4. 3.

The bandwidth $h_{c}$ satisfies $h_{c}\rightarrow 0$,
$nh_{c}^{2}\delta_{n}^{2}\rightarrow\infty$,
$\sum_{i=1}^{n}N_{i}(N_{i}-1)v_{i}^{2}/h_{c}^{2}\rightarrow 0$,
$\sum_{i=1}^{n}N_{i}(N_{i}-1)(N_{i}-2)v_{i}^{2}/h_{c}\rightarrow 0$,
$\sum_{i=1}^{n}N_{i}(N_{i}-1)(N_{i}-2)(N_{i}-3)v_{i}^{2}\rightarrow 0$, and
(4.8) with $h$ replaced by $h_{c}$ and $\gamma_{n1}$ replaced by
$\gamma_{n2}$.

  5. 4.

The bandwidths $h_{c}$ and $b$ satisfy $h_{c},b\rightarrow 0$,
$nh_{c}^{2}\delta_{n}^{2}\rightarrow\infty$,
$\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v^{2}_{i,i^{\prime}}N_{i}N_{i^{\prime}}/h_{c}^{2}b\rightarrow 0$,
$\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v_{i,i^{\prime}}^{2}N_{i}N_{i^{\prime}}(N_{i}-1)(N_{i^{\prime}}-1)b^{-1}\rightarrow
0$, $\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v_{i,i^{\prime}}^{2}N_{i}N_{i^{\prime}}(N_{i}-1)h_{c}^{-1}b^{-1}\rightarrow
0$, $\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v_{i,i^{\prime}}^{2}N_{i}N_{i^{\prime}}(N_{i^{\prime}}-1)h_{c}^{-1}b^{-1}\rightarrow
0$, and (4.8) with $h$ replaced by $h_{c}$ and $\gamma_{n1}$ replaced by
$\gamma_{n3}$.

The first part of each of the Condition (IV) is standard in the nonparametric
smoothing literature and is needed to obtain the consistency of estimators.
For example, they are similar to (C1a) and (D1a) of Zhang and Wang (2016). The
conditions involving $\delta_{n}$ and $\kappa$ in (4.8) are mild and similar
to the condition in Assumption (IV) (a) in Lu and Tjøstheim (2014) and (A4) in
Hallin, Lu and Tran (2004).

Condition (V) (uniform convergence)

  1. [label=()]

  2. 1.

Assume $\sup_{i}nw_{i}N_{i}<\infty$,

$\displaystyle
nQ_{n}^{\frac{2+\varrho}{1+\varrho}}h_{\mu}^{\frac{2+\varrho}{1+\varrho}}\delta_{n}^{2}$
$\displaystyle\to\infty,\quad\frac{n^{\rho+1}Q_{n}^{\frac{(2+\varrho)\rho+4+3\varrho}{1+\varrho}}h_{\mu}^{\frac{(2+\varrho)\rho+6+5\varrho}{1+\varrho}}\delta_{n}^{2}\log(n)}{\Delta_{n}^{6}}\to\infty,$
(4.9)

with

$\displaystyle Q_{n}$
$\displaystyle=\left\\{\sum_{i=1}^{n}w_{i}^{2}N_{i}h_{\mu}^{-1}+\sum_{i=1}^{n}w_{i}^{2}N_{i}(N_{i}-1)\right\\}^{1/2}.$

  3. 2.

Assume $\sup_{i}nw_{i}N_{i}=O(1)$, $nh^{2}\delta_{n}^{2}\to\infty$,
$h/\delta_{n}\rightarrow 0$, and

$\displaystyle\frac{n^{1+\rho}h^{14+2\rho}\delta_{n}^{5+2\rho}}{\Delta_{n}^{6}\log^{4}(n)}$
$\displaystyle\to$ $\displaystyle\infty.$ (4.10)

  4. 3.

Assume $v_{i,i^{\prime}}=v_{i^{\prime},i}$,
$\sup_{i,i^{\prime}}n^{2}v_{i,i^{\prime}}N_{i}N_{i^{\prime}}\leq C$ for some
constant $C$, and $\frac{nh_{c}^{2}b}{\log(n)}\to\infty$.

The above conditions are mild and sufficient for uniform convergence.
Condition (V) illustrates the roles of $\rho,\delta_{n},\Delta_{n}$, namely,
the spatial correlation and DEI parameters. A larger $\rho$ means weaker
spatial dependence, and the choice of the bandwidths $h_{\mu},h$ becomes more
flexible. In the case without spatial correlation, namely $\rho=\infty$, the
conditions in (4.9) and (4.10) are almost equivalent (up to logarithm factors)
to

$\displaystyle nQ_{n}h_{\mu}\delta_{n}^{2}\to\infty,\quad
nQ_{n}^{\frac{2+\varrho}{1+\varrho}}h_{\mu}^{\frac{2+\varrho}{1+\varrho}}\to\infty,\quad
nh^{2}\delta_{n}^{2}\to\infty.$

The above conditions are similar to Conditions (C2c) and (C3c) from Zhang and
Wang (2016), except the DEI parameter $\delta_{n}$ is involved.

We now provide the asymptotic results for the covariance estimators proposed
in the study. These results require the uniform convergence of
$\widehat{\mu}(t)$ which is provided in Lemma S.1.3 in the supplement Chatla
and Liu (2024).

###  4.1 Results for $\widehat{R}(0;t_{1},t_{2})$

In the following theorem, we provide the asymptotic bias and variance of the
spatial covariance estimator $\widehat{R}(\|\bm{s}_{0}\|;t_{1},t_{2})$.

######  Theorem 4.1.

Under conditions (I)[(i), (ii) (iii), (v)(b), (vi)], (II), (III), (IV)[(ii),
(iv)] and (V)(i), for a fixed interior point $(\|\bm{s}_{0}\|,t_{1},t_{2})$,
we have that

$\displaystyle E\\{\widehat{R}$
$\displaystyle(\|\bm{s}_{0}\|;t_{1},t_{2})-R(\|\bm{s}_{0}\|;t_{1},t_{2})~{}|~{}\mathcal{X}_{T}\\}$
$\displaystyle=\frac{1}{2}h_{c}^{2}\left\\{\frac{\partial^{2}R(\|\bm{s}_{0}\|;t_{1},t_{2})}{\partial
t_{1}^{2}}+\frac{\partial^{2}R(\|\bm{s}_{0}\|;t_{1},t_{2})}{\partial
t_{2}^{2}}\right\\}+\frac{1}{2}b^{2}\mu_{L,2}\frac{\partial^{2}R(\|\bm{s}_{0}\|;t_{1},t_{2})}{\partial\|\bm{s}_{0}\|^{2}}+o_{p}(h_{c}^{2}+b^{2})$

and

var
$\displaystyle\\{\widehat{R}(\|\bm{s}_{0}\|;t_{1},t_{2})~{}|\mathcal{X}_{T}\\}$
$\displaystyle\qquad=\frac{\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v^{2}_{i,i^{\prime}}N_{i}N_{i^{\prime}}}{h_{c}^{2}bA_{0}^{*}f_{T}(t_{1})f_{T}(t_{2})}\nu_{L,2}\nu_{K_{1},2}^{2}~{}\Xi_{R}(t_{1},t_{2};\|\bm{s}_{0}\|)$
$\displaystyle\qquad\qquad+\frac{\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v^{2}_{i,i^{\prime}}N_{i}N_{i^{\prime}}(N_{i^{\prime}}-1)}{h_{c}bA_{0}^{*}f_{T}(t_{1})}\nu_{L,2}\nu_{K_{1},2}~{}\Xi_{R}(t_{1},t_{2},t_{2};\|\bm{s}_{0}\|)$
$\displaystyle\qquad\qquad+\frac{\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v^{2}_{i,i^{\prime}}N_{i}N_{i^{\prime}}(N_{i}-1)}{h_{c}bA_{0}^{*}f_{T}(t_{2})}\nu_{L,2}\nu_{K_{1},2}~{}\Xi_{R}(t_{1},t_{2},t_{1};\|\bm{s}_{0}\|)$
$\displaystyle\qquad\qquad+\frac{\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v^{2}_{i,i^{\prime}}N_{i}(N_{i}-1)N_{i^{\prime}}(N_{i^{\prime}}-1)}{bA_{0}^{*}}\nu_{L,2}\Xi_{R}(t_{1},t_{2},t_{1},t_{2};\|\bm{s}_{0}\|)+o_{p}(1),$

where $\mathcal{X}_{T}=\\{T_{ij},i=1,\ldots,n,j=1,\ldots,N_{i}\\}$,
$\nu_{L,2}=\int L^{2}(u)du$, $\nu_{K_{1},2}=\int K_{1}^{2}(u)du$ and $\Xi_{R}$
is defined in Condition (I)(vi), $A_{0}^{*}$ is defined in Condition (II)(iv),
and $\mu_{L,2}$ is defined in Condition (III)(ii).

The proof of Theorem 4.1 is available in the supplement Chatla and Liu (2024).

######  Remark 1.

The asymptotic bias of the estimator $\widehat{R}(\|\bm{s}_{0}\|;t_{1},t_{2})$
consists of the biases from two kernel smoothing functions: over the time,
$T_{ij}$, and over the values of spatial locations, $\bm{s}_{i}$. It is
different from the rate of standard two-dimensional smoothing, which includes
only the bias from kernel smoothing over time. As discussed in Lu and
Tjøstheim (2014), under the condition $b/h_{c}\rightarrow 0$ as
$n\rightarrow\infty$, the asymptotic bias will simplify to the standard two-
dimensional rate $O(h_{c}^{2})$. Similarly, the asymptotic variance is
different from the standard rate due to the extra bandwidth $b$. In Remark 4,
we give sufficient conditions under which the asymptotic variance will
simplify to the standard rate.

######  Remark 2.

It is a common practice in the literature to provide asymptotic normality of
estimator instead of pointwise bias and variance as in Theorem 4.1 when the
observations are dependent (Fan and Yao, 2008). However, the asymptotic
normality of $\widehat{R}(\|\bm{s}_{0}\|;t_{1},t_{2})$ is very challenging
especially under the broad framework such as the one adopted in this
manuscript and requires stronger conditions. We leave this for future
research.

We define the Hilbert-Schmidt norm by
$\|\Phi\|_{HS}=[\int\int\Phi(s,t)^{2}dsdt]^{1/2}$. The following lemma
provides the convergence rate of $\widehat{R}(0)$ in Hilbert-Schmidt norm.

######  Lemma 4.2.

Suppose conditions (I)[(i), (ii) (iii), (v)(b), (vi)], (II), (III), (IV)[(ii),
(iv)], and (V)[(i), (iii)] hold. Then, we have that

$\displaystyle\|\widehat{R}(0)$
$\displaystyle-R(0)\|_{HS}=O_{p}\bigg{(}h_{c}^{2}+b^{2}+\bigg{[}\frac{\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v^{2}_{i,i^{\prime}}N_{i}N_{i^{\prime}}}{h_{c}^{2}b}+\frac{\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v^{2}_{i,i^{\prime}}N_{i}N_{i^{\prime}}(N_{i^{\prime}}-1)}{h_{c}b}$
$\displaystyle+\frac{\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v^{2}_{i,i^{\prime}}N_{i}N_{i^{\prime}}(N_{i}-1)}{h_{c}b}+\frac{\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v^{2}_{i,i^{\prime}}N_{i}(N_{i}-1)N_{i^{\prime}}(N_{i^{\prime}}-1)}{b}\bigg{]}^{1/2}\bigg{)}.$
(4.11)

The proof of the above lemma is available in the supplement Chatla and Liu
(2024). In the following corollary, we provide the corresponding rates for the
OBS and SUBJ weighting schemes. Prior to that, we introduce some notation. Let
$N_{S}=\sum_{i=1}^{n}N_{i}$, $\overline{N}=N_{S}/n$,
$\overline{N}_{H}=[n^{-1}\sum_{i=1}^{n}N_{i}^{-1}]^{-1}$, and

$\displaystyle\overline{N}_{2,12}$
$\displaystyle=n^{-1}(n-1)^{-1}\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}N_{i}^{2}N_{i^{\prime}}^{2},\qquad\overline{N}_{H,2}=\big{(}n^{-1}(n-1)^{-1}\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}N_{i}^{-1}N_{i^{\prime}}^{-1}\big{)}^{-1},$
$\displaystyle\overline{N}_{2,2}$
$\displaystyle=n^{-1}(n-1)^{-1}\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}N_{i}N_{i^{\prime}}^{2},\qquad\overline{N}_{2,1}=n^{-1}(n-1)^{-1}\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}N_{i}^{2}N_{i^{\prime}}.$

Let $\widehat{R}_{obs}$ and $\widehat{R}_{subj}$ be the specific estimators of
$R(\cdot)$ for OBS and SUBJ weighting schemes, respectively. The following
corollary provides the convergence rate of both these estimators.

######  Corollary 4.3.

Suppose that conditions (I)[(i), (ii) (iii), (v)(b), (vi)], (II), (III),
(IV)[(ii), (iv)], and (V)[(i), (iii)] hold. Then, we have the following:

  1. (1)

OBS:

$\displaystyle\|\widehat{R}_{obs}(0)$
$\displaystyle-R(0)\|_{HS}=O_{p}\bigg{(}h_{c}^{2}+b^{2}+\sqrt{\left(\frac{1}{\overline{N}_{2}h_{c}^{2}b}+\frac{\overline{N}_{2,1}}{\overline{N}_{2}^{2}h_{c}b}+\frac{\overline{N}_{2,12}}{\overline{N}_{2}^{2}b}\right)\frac{1}{n(n-1)}}\bigg{)}.$

  2. (2)

SUBJ:

$\displaystyle\|\widehat{R}_{subj}(0)$
$\displaystyle-R(0)\|_{HS}=O_{p}\bigg{(}h_{c}^{2}+b^{2}+\sqrt{\left(\frac{1}{\overline{N}_{H,2}h_{c}^{2}b}+\frac{1}{\overline{N}_{H}h_{c}b}+\frac{1}{b}\right)\frac{1}{n(n-1)}}\bigg{)}.$

######  Remark 3.

We note that $\overline{N}\geq\overline{N}_{H}$ and
$\overline{N}_{2}\geq\overline{N}_{H2}$. Define
$\overline{N}_{S2}=n^{-1}\sum_{i=1}^{n}N_{i}^{2}$. Analogously to the
arguments used in Zhang and Wang (2016), the rate for $\widehat{R}_{obs}$ is
as good as $\widehat{R}_{subj}$ if
$\lim\sup_{n}\overline{N}_{S2}/(\overline{N})^{2}<\infty$,
$\lim\sup_{n}(\overline{N}\cdot\overline{N}_{2,1})/(\overline{N}_{2})^{2}<\infty$,
and $\lim\sup_{n}(\overline{N}_{2,12})/(\overline{N}_{2}^{2})<\infty$.

######  Remark 4.

When $\overline{N}_{2}<\infty$ and $\overline{N}_{H,2}<\infty$, the rate of
convergence of $\widehat{R}(0)$ without spatial smoothing is
$O_{p}(1/\sqrt{N(\Delta)}h_{c}^{2})$ (Liu, Ray and Hooker, 2017) where
$N(\Delta)$ is the collection of location pairs that share the same spatial
correlation structure. Comparing this with the asymptotic order in Corollary
4.3, it appears reasonable to take $b=O(1/\sqrt{n-1})$. Then under the
conditions in Remark 3 the convergence rates in Corollary 4.3 will be of
similar order as that of two-dimensional smoothing as in Corollary 4.3 of
Zhang and Wang (2016). Moreover, $b=O(1/\sqrt{n-1})$ satisfy
$b/h_{c}\rightarrow 0$ under the condition that $nh_{c}^{2}\rightarrow\infty$.

######  Corollary 4.4.

Suppose that conditions (I)[(i), (ii) (iii), (v)(b), (vi)], (II), (III),
(IV)[(ii), (iv)], and (V)[(i), (iii)] hold. Then, we have the following:

  1. (1)

OBS: Suppose the bandwidth satisfies $b=O(1/\sqrt{n-1})$ and
$nh_{c}^{2}\rightarrow\infty$. Moreover, assume that
$\lim\sup_{n}\overline{N}_{S2}/(\overline{N})^{2}<\infty$,
$\lim\sup_{n}(\overline{N}\cdot\overline{N}_{2,1})/(\overline{N}_{2})^{2}<\infty$,
and $\lim\sup_{n}(\overline{N}_{2,12})/(\overline{N}_{2}^{2})<\infty$.

    1. (a)

If $\overline{N}_{2}/n^{1/2}\rightarrow 0$ and
$h_{c}\asymp(n\overline{N}_{2})^{-1/6}$, then

$\displaystyle\|\widehat{R}_{obs}(0)-R(0)\|_{HS}$
$\displaystyle=O_{p}(h_{c}^{2})+o_{p}\left(\frac{1}{\sqrt{n\overline{N}_{2}h_{c}^{2}}}\right).$

    2. (b)

If $\overline{N}_{2}/n^{1/2}\rightarrow C$ for some $C<\infty$ and
$h_{c}\asymp n^{-1/4}$, then

$\displaystyle\|\widehat{R}_{obs}(0)-R(0)\|_{HS}$
$\displaystyle=O_{p}\left(\frac{1}{\sqrt{n}}\right)+o_{p}\left(\frac{1}{\sqrt{n}}\right).$

    3. (c)

If $\overline{N}_{2}/n^{1/2}\rightarrow\infty$, $h_{c}=o(n^{-1/4})$ and
$h_{c}\overline{N}_{2}\rightarrow\infty$, then

$\displaystyle\|\widehat{R}_{obs}(0)-R(0)\|_{HS}$
$\displaystyle=O_{p}\left(\frac{1}{\sqrt{n}}\right)+o_{p}\left(\frac{1}{\sqrt{n}}\right).$

  2. (2)

SUBJ: Replacing $\overline{N}$ and $\overline{N}_{2}$ in (1) with
$\overline{N}_{H}$ and $\overline{N}_{H,2}$, respectively, leads to the
results for $\widehat{R}_{subj}$.

######  Remark 5.

As in Zhang and Wang (2016), under either the OBS or SUBJ scheme, Corollary
4.4 yields a systematic partition of functional data into three categories:
non-dense, dense and ultra-dense. Cases (b) and (c) correspond to dense data
as root-$n$ convergence can be achieved, and case (a) corresponds to non-dense
data where root-$n$ convergence can never be achieved.

######  Remark 6.

As suggested in Zhang and Wang (2016), one may use a mixture of the OBS and
SUBJ schemes which attain at least the better rate between
$\widehat{R}_{obs}(0)$ and $\widehat{R}_{subj}(0)$. One could use
$v_{i,i^{\prime}}=\theta/[\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}N_{i}N_{i^{\prime}}]+(1-\theta)/[n(n-1)N_{i}N_{i^{\prime}}]$ for some
$0\leq\theta\leq 1$. Then the rate of the corresponding estimator, denoted by
$\widehat{R}_{\theta}(0)$, is

$\displaystyle\|\widehat{R}_{\theta}(0)-R(0)\|_{HS}$
$\displaystyle=O_{p}\left(h_{c}^{2}+\sqrt{\theta^{2}d_{n1}+(1-\theta)^{2}d_{n2}}\right),$

where

$\displaystyle
d_{n1}=\left(\frac{1}{\overline{N}_{2}h_{c}^{2}}+\frac{\overline{N}_{2,1}}{\overline{N}_{2}^{2}h_{c}}+\frac{\overline{N}_{2,12}}{\overline{N}_{2}^{2}}\right)\frac{1}{n}\quad\text{
and }\quad
d_{n2}=\left(\frac{1}{\overline{N}_{H,2}h_{c}^{2}}+\frac{1}{\overline{N}_{H}h_{c}}+1\right)\frac{1}{n}.$

It is easy to note that the rate is minimized for
$\theta^{*}=d_{n2}/(d_{n1}+d_{n2})$. Recently, Zhang and Wang (2018) proposed
an optimal weighting scheme for longitudinal and functional data. This scheme
is designed to minimize the convergence rate given a bandwidth. One may adopt
a similar approach in our setting.

###  4.2 Results for $\widehat{R}_{e}(0,t_{1},t_{2})$

The following theorem provides the asymptotic bias and variance for the
inverse regression function $m(t,y)$ in (3.7). It extends Lemma 2.1 of Jiang,
Yu and Wang (2014) to spatially correlated functional data.

######  Theorem 4.5.

Suppose conditions (I)[(i), (ii), (iv)], (II)(i), (III)(i) and (IV)(i) hold.
For a fixed interior point $(t,y)$, if
$h_{t}\sum_{i=1}^{n}N_{i}(N_{i}-1)w_{i}^{2}/\sum_{i=1}^{n}N_{i}w_{i}^{2}\rightarrow
C_{0}$, for some constant $C_{0}$, then

$\displaystyle E\\{\widehat{m}(t,y)-m(t,y)~{}|~{}\mathcal{X}_{y}\\}$
$\displaystyle=$
$\displaystyle\frac{c(K_{2})}{2}tr\left(\bm{H}.\bm{\mathcal{H}}_{m}(t,y)\right)+o_{p}(tr(\bm{H}))$
(4.12)

and

$\displaystyle\text{var}\\{\widehat{m}(t,y)~{}|~{}\mathcal{X}_{y}\\}$
$\displaystyle=$
$\displaystyle\frac{\sum_{i=1}^{n}N_{i}w_{i}^{2}}{|\bm{H}|^{1/2}g^{2}(t,y)}(\Psi(t,y)+\Lambda(t)+\sigma^{2})\nu_{K_{2},2}g(t,y)$
(4.13)
$\displaystyle+\frac{\sum_{i=1}^{n}N_{i}(N_{i}-1)w_{i}^{2}}{h_{y}g^{2}(t,y)}\times\\{(\Psi(t,y,t)+\Lambda(t,t))g_{3}(t,y,t)\\}+o_{p}(1),$

where $\mathcal{X}_{y}=\\{(T_{ij},Y_{i}),i=1,\ldots,n,j=1,\ldots,N_{i}\\}$ is
the set of design points, $\bm{H}=\text{diag}\\{h_{t}^{2},h_{y}^{2}\\}$,
$\int\bm{u}\bm{u}^{T}K_{2}(\bm{u})d\bm{u}=c(K_{2})\bm{I}$ with an identity
matrix $\bm{I}$, $\bm{\mathcal{H}}_{m}(t,y)$ is the second derivative of
$m(t,y)$, and $\nu_{K_{2},2}=\int K_{2}^{2}(\bm{u})d\bm{u}$.

The proof of Theorem 4.5 is available in the supplement Chatla and Liu (2024).
By letting $w_{i}=1/N_{S,0}$ and $w_{i}=1/nN_{i}$, respectively, in Theorem
4.5, we obtain the corresponding rates for the OBS and SUBJ estimators. Let
$\mathfrak{B}\subset\mathbb{R}$ be a compact set. We now estimate the
conditional covariance $R_{e}(0,t_{1},t_{2})$ empirically as,

$\displaystyle\widehat{R}_{e}(0,t_{1},t_{2})$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\widehat{m}(t_{1},Y_{i})\widehat{m}(t_{2},Y_{i})\mathbbm{1}_{\\{Y_{i}\in\mathfrak{B}\\}}-\frac{1}{n}\sum_{i=1}^{n}\widehat{m}(t_{1},Y_{i})\mathbbm{1}_{\\{Y_{i}\in\mathfrak{B}\\}}\frac{1}{n}\sum_{j=1}^{n}\widehat{m}(t_{2},Y_{j})\mathbbm{1}_{\\{Y_{j}\in\mathfrak{B}\\}},$

where $\mathbbm{1}_{\\{Y_{j}\in\mathfrak{B}\\}}$ is an indicator function that
takes value one if $Y_{j}\in\mathfrak{B}$ and zero otherwise. With the help of
Theorem 4.5, in the following lemma, we obtain the convergence rate of
covariance operator $\widehat{R}_{e}(0,t_{1},t_{2})$ in Hilbert-Schmidt norm.

######  Lemma 4.6.

Assume the conditions (I)[(i),(ii), (iv)], (II)(i), (III)(i), (IV)(i) and
(V)(ii) hold. If $h_{t}\asymp h_{y}\asymp h\asymp n^{-1/d}$ where
$d\in(2,14)$, then

$\displaystyle\|\widehat{R}_{e}(0)-R_{e}(0)\|_{HS}$
$\displaystyle=O_{p}\left(h^{2}+\sqrt{\frac{\sum_{i=1}^{n}N_{i}w_{i}^{2}}{h^{2}}+\frac{\sum_{i=1}^{n}N_{i}(N_{i}-1)w_{i}^{2}}{h}}\right).$
(4.14)

######  Remark 7.

The proof of the above lemma is available in the supplement Chatla and Liu
(2024). It requires uniform convergence of the cross-product terms $u_{l,k}$
in (S.36) in the supplement Chatla and Liu (2024). For this reason, we
restrict both uniform convergence and Hilbert-Schmidt norm computations to the
compact set $(t,y)\in\mathcal{I}\times\mathfrak{B}$. In practice, we can trim
both ends of the support of $y$ and consider the middle 90% or 95% of $y$
values to compute $\widehat{R}_{e}(0)$. This in turn improves the quality of
the estimator by alleviating boundary effects associated with smoothing.

Let $\widehat{R}_{e,obs}(0)$ and $\widehat{R}_{e,subj}(0)$ be the specific
estimators of $R_{e}(\cdot)$ for the OBS and SUBJ weighting schemes,
respectively. The following corollary provides the convergence rates of these
estimators.

######  Corollary 4.7.

Suppose conditions (I)[(i), (ii), (iv)], (II)(i),(III)(i), (IV)(i) and (V)(ii)
hold. If the bandwidths satisfy $h_{t}\asymp h_{y}\asymp h\asymp n^{-1/d}$
where $d\in(2,14)$, then

$\displaystyle\textrm{(1)
OBS:}\quad\|\widehat{R}_{e,obs}(0)-R_{e}(0)\|_{HS}=O_{p}\left(h^{2}+\sqrt{\left(\frac{1}{\overline{N}h^{2}}+\frac{\overline{N}_{S2}}{(\overline{N})^{2}h}\right)\frac{1}{n}}\right);$
$\displaystyle\textrm{(2)
SUBJ:}\quad\|\widehat{R}_{e,subj}(0)-R_{e}(0)\|_{HS}=O_{p}\left(h^{2}+\sqrt{\left(\frac{1}{\overline{N}_{H}h^{2}}+\frac{1}{h}\right)\frac{1}{n}}\right).$

The above rates are consistent with the rates of two-dimensional smoothing.
The following corollary provides the convergence rates for sparse, dense and
ultra-dense functional data.

######  Corollary 4.8.

Suppose conditions (I)[(i), (ii), (iv)], (II)(i), (III)(i), (IV)(i) and
(V)(ii) hold. Assume $h_{t}\asymp h_{y}\asymp h\asymp n^{-1/d}$ where
$d\in(2,14)$. Then, we have the following:

  1. (1)

OBS: Assume $\lim\sup_{n}\overline{N}_{S2}/(\overline{N})^{2}<\infty$.

    1. (a)

If $\overline{N}/n^{1/2}\rightarrow 0$ and $h\asymp(n\overline{N})^{-1/6}$,
then
$\|\widehat{R}_{e,obs}(0)-R_{e}(0)\|_{HS}=O_{P}(h^{2}+1/\sqrt{n\overline{N}h^{2}})$.

    2. (b)

If $\overline{N}/n^{1/2}\rightarrow C$ for some $C<\infty$ and $h\asymp
n^{-1/4}$, then $\|\widehat{R}_{e,obs}(0)-R_{e}(0)\|_{HS}=O_{P}(1/\sqrt{n})$.

    3. (c)

If $\overline{N}/n^{1/2}\rightarrow\infty$, $h=o(n^{-1/4})$, and
$h\overline{N}\rightarrow\infty$, then
$\|\widehat{R}_{e,obs}(0)-R_{e}(0)\|_{HS}=O_{P}(1/\sqrt{n})$.

  2. (2)

SUBJ: Replacing $\overline{N}$ in (1) with $\overline{N}_{H}$ leads to the
results for $\widehat{R}_{e,subj}$.

###  4.3 Results for $\widehat{\Gamma}(t_{1},t_{2})$

We now derive the asymptotic bias and variance expression for the covariance
operator $\widehat{\Gamma}(t_{1},t_{2})$ in (3.5) using general weight $v_{i}$
for the $i$th observation. This result requires the uniform convergence of
$\widehat{\mu}(t)$ which is provided in Lemma S.1.3 in the supplement Chatla
and Liu (2024). Let $\Xi_{1,2}=\Xi(t_{1},t_{2})$,
$\Xi_{1,2,2}=\Xi(t_{1},t_{2},t_{2})$, $\Xi_{1,2,1}=\Xi(t_{1},t_{2},t_{1})$ and
$\Xi_{1,2,1,2}=\Xi(t_{1},t_{2},t_{1},t_{2})$ for $\Xi$ defined in Condition
(I)(v). Similarly, let $\Lambda_{1,2}=\Lambda(t_{1},t_{2})$,
$\Lambda_{1,2,2}=\Lambda(t_{1},t_{2},t_{2})$,
$\Lambda_{1,2,1}=\Lambda(t_{1},t_{2},t_{1})$ and
$\Lambda_{1,2,1,2}=\Lambda(t_{1},t_{2},t_{1},t_{2})$ for $\Lambda$ defined in
condition (I)(v). The following theorem extends Theorem 3.2 of Zhang and Wang
(2016) to spatially correlated functional data except the asymptotic normality
result, which is not pursued here given the complicated setting.

######  Theorem 4.9.

Under the conditions (I)[(i), (ii), (iii), (v)], (II)(i), (III)(i), (IV)(iii)
and (V)(i), for a fixed interior point $(t_{1},t_{2})$, we have that

$\displaystyle
E\\{\widehat{\Gamma}(t_{1},t_{2})-\Gamma(t_{1},t_{2})~{}|~{}\mathcal{X}_{T}\\}$
$\displaystyle=\frac{c(K_{1})}{2}tr\left(\bm{H}_{c}.\bm{\mathcal{H}}_{c}(t_{1},t_{2})\right)+o_{p}(tr(\bm{H}_{c}))$

and

$\displaystyle\text{var}\\{\widehat{\Gamma}(t_{1},t_{2})~{}|~{}\mathcal{X}_{T}\\}$
$\displaystyle=\frac{\sum_{i=1}^{n}N_{i}(N_{i}-1)v_{i}^{2}}{h_{c}^{2}}\bigg{(}\frac{\Xi_{1,2}+\Lambda_{1,2}+\sigma^{4}}{f_{T}(t_{1})f_{T}(t_{2})}\left[\nu^{2}_{K_{1},2}+K_{1}^{*^{2}}\left(\frac{t_{1}-t_{2}}{h_{c}}\right)\right]\bigg{)}$
$\displaystyle+\frac{\sum_{i=1}^{n}N_{i}(N_{i}-1)(N_{i}-2)v_{i}^{2}}{h_{c}}$
$\displaystyle\times\bigg{(}\frac{[\Xi_{1,2,2}+\Lambda_{1,2,2}]f_{T}(t_{1})+[\Xi_{1,2,1}+\Lambda_{1,2,1}]f_{T}(t_{2})}{f_{T}(t_{1})f_{T}(t_{2})}\left[\nu_{K_{1},2}+K_{1}^{*}\left(\frac{t_{1}-t_{2}}{h_{c}}\right)\right]\bigg{)}$
$\displaystyle\qquad\qquad\qquad+o_{p}(1),$

where $\mathcal{X}_{T}$ is defined in Theorem 4.1, $\nu_{K_{1},2}=\int
K_{1}^{2}(u)du$, $c(K_{1})=\int u^{2}K_{1}(u)du$,
$\bm{H}_{c}=\text{diag}\\{h_{c}^{2},h_{c}^{2}\\}$,
$\bm{\mathcal{H}}_{c}(t_{1},t_{2})$ is the second derivatives of
$\Gamma(t_{1},t_{2})$, and $K_{1}^{*}(\cdot)$ is a convolution of kernel
$K_{1}(\cdot)$ with itself.

######  Remark 8.

We omit the proof, as it proceeds similarly to Theorem 4.1. The asymptotic
bias and variance expressions are similar to those from Theorem 3.2 of Zhang
and Wang (2016), which assumes that the functions are independent. As
mentioned in Zhang and Wang (2016), the covariance estimator exhibits the
“discontinuity of the asymptotic variance” because the asymptotic variance in
Theorem 4.9 is different for $t_{1}=t_{2}$ and $t_{1}\neq t_{2}$. We refer to
Corollaries 3.4 and 3.5 in Zhang and Wang (2016) for a comprehensive
discussion.

###  4.4 Results for e.d.r. directions

To provide convergence rates for e.d.r. directions, we consider the following
condition on the decay rate of $\\{\xi_{i}\\}$ and an additional assumption on
the coefficients $a_{k,i}=\langle\pi_{i},\eta_{k}\rangle$. Recall that
$R(0,t,t^{\prime})=\sum_{i=1}^{\infty}\xi_{i}\pi_{i}(t)\pi_{i}(t^{\prime})$,
where $\\{\xi_{i}\\}_{1}^{\infty}$ are the eigenvalues and
$\\{\pi_{i}\\}_{1}^{\infty}$ are the corresponding eigenfunctions.

Condition (VI) (eigenvalues):

  1. [label=()]

  2. 1.

Assume that $\xi_{i}-\xi_{i+1}>\overline{M}_{0}i^{-\overline{\alpha}_{1}-1}$
for some constant $\overline{M}_{0}>0$, $i\geq 1$ and
$\overline{\alpha}_{1}>1$. This implies
$\xi_{i}\geq\overline{M}_{1}i^{-\overline{\alpha}_{1}}$ for some constant
$\overline{M}_{1}>0$.

  3. 2.

The coefficient $|a_{k,j}|<\overline{M}_{2}j^{-\overline{\alpha}_{2}}$ for
some positive constants $\overline{M}_{2}$ and $\overline{\alpha}_{2}$.

The above condition (a) ensures that the gaps between the eigenvalues will not
be too small. It also provides a lower bound on the rate at which $\xi_{i}$
decreases. The condition (b) assumes that the coefficients $a_{k,j}$ do not
decrease too quickly. Similar conditions exist in the literature for
functional data, for example, we refer to Hall and Horowitz (2007) and Jiang
and Wang (2015). Because $\eta_{k}=\sum_{i=1}^{\infty}a_{k,i}\pi_{i}(t)$ and
$\eta_{k}\in\mathfrak{R}_{R^{1/2}(0)}$, we obtain that
$\sum_{i=1}^{\infty}a_{k,i}^{2}\xi_{i}^{-1}<\infty$. Under Condition (VI),
this implies that $2\overline{\alpha}_{2}-\overline{\alpha}_{1}>1$. A
sufficient condition for $2\overline{\alpha}_{2}-\overline{\alpha}_{1}>1$ is
$\overline{\alpha}_{2}\geq\overline{\alpha}_{1}$ and it means that the space
spanned by $\\{\eta_{k}\\}_{k=1}^{K}$ is at least as smooth as operator
$R(0)$.

Let $R_{L}^{-1/2}(0)=\sum_{i=1}^{L}\xi_{i}^{-1/2}\pi_{i}(s)\pi_{i}(t)$ be the
truncated version of $R^{-1/2}(0)$. Let

$\displaystyle r_{n1}$
$\displaystyle=h^{4}+\frac{\sum_{i=1}^{n}N_{i}w_{i}^{2}}{h^{2}}+\frac{\sum_{i=1}^{n}N_{i}(N_{i}-1)w_{i}^{2}}{h}~{}\text{
and}$ $\displaystyle r_{n2}$
$\displaystyle=h_{c}^{4}+b^{4}+\frac{\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v^{2}_{i,i^{\prime}}N_{i}N_{i^{\prime}}}{h_{c}^{2}b}+\frac{\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v^{2}_{i,i^{\prime}}N_{i}N_{i^{\prime}}(N_{i^{\prime}}-1)}{h_{c}b}$
$\displaystyle+\frac{\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v^{2}_{i,i^{\prime}}N_{i}N_{i^{\prime}}(N_{i}-1)}{h_{c}b}+\frac{\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v^{2}_{i,i^{\prime}}N_{i}(N_{i}-1)N_{i^{\prime}}(N_{i^{\prime}}-1)}{b}.$

In the following theorem, we provide the convergence rate for the e.d.r.
direction $\widehat{\beta}_{j}$.

######  Theorem 4.10.

Under the conditions of Lemma 4.2 and Lemma 4.6 and condition (VI), we have
that

$\displaystyle\|\widehat{\beta}_{j}-\beta_{j}\|_{L_{2}}^{2}$
$\displaystyle=O_{p}\left(L^{(-2\overline{\alpha}_{2}+\overline{\alpha}_{1}+1)}+L^{(3\overline{\alpha}_{1}+2)}r_{n1}+L^{(3\overline{\alpha}_{1}-2\overline{\alpha}_{2}+4)}r_{n2}\right)\text{
for }j=1,\ldots,K.$

######  Remark 9.

The proof uses the results from Hall and Horowitz (2007) and proceeds
similarly to Theorem 2.1 in Jiang and Wang (2015). Therefore, it is omitted
here for brevity. The convergence rate of e.d.r. direction
$\widehat{\beta}_{j}$ depends on the rates of involved covariance operators
$\widehat{R}_{e}(0)$ and $\widehat{R}(0)$. The second and third terms in the
above rate are equal for the choice of
$L=(r_{n2}/r_{n1})^{1/2(\overline{\alpha}_{2}-1)}$. Corollary 4.11 provides
conditions under which $r_{n2}\asymp r_{n1}$.

######  Remark 10.

The rates obtained in Theorem 4.10 may not be optimal for the e.d.r.
directions $\widehat{\beta}_{j}$ which, as discussed in Jiang, Yu and Wang
(2014), can be estimated at the optimal rate for a one-dimensional smoother.
In the following Corollary 4.11, under some regularity conditions, we show
that $\widehat{\beta}_{j}$ achieve the optimal convergence rate.

The following corollary provides the convergence rate of $\widehat{\beta}_{j}$
for sparse, dense and ultra-dense functional data similar to Corollaries 4.4
and 4.8. Let $\widehat{\beta}_{obs}$ and $\widehat{\beta}_{subj}$ be the
specific estimators of $\beta$ for the OBS and SUBJ weighting schemes,
respectively.

######  Corollary 4.11.

Suppose all the conditions in corollaries 4.4 and 4.8 including the ones for
OBS and SUBJ schemes hold. Additionally, assume $h_{c}\asymp h$ and
$\overline{\alpha}_{2}>>\overline{\alpha}_{1}>1$. Then, we have the following:

  1. (1)

OBS:

    1. (a)

If $\overline{N}/n^{1/2}\rightarrow 0$ and $h\asymp(n\overline{N})^{-1/6}$,
then
$\|\widehat{\beta}_{obs}-\beta\|_{L_{2}}=O_{P}(h^{2}+1/\sqrt{n\overline{N}h^{2}})$.

    2. (b)

If $\overline{N}/n^{1/2}\rightarrow C$ for some $C<\infty$ and $h\asymp
n^{-1/4}$, then $\|\widehat{\beta}_{obs}-\beta\|_{L_{2}}=O_{p}(1/\sqrt{n}).$

    3. (c)

If $\overline{N}/n^{1/2}\rightarrow\infty$, $h=o(n^{-1/4})$, and
$h\overline{N}\rightarrow\infty$, then
$\|\widehat{\beta}_{obs}-\beta\|_{L_{2}}=O_{p}(1/\sqrt{n}).$

  2. (2)

SUBJ: Replacing $\overline{N}$ and $\overline{N}_{2}$ in (1) with
$\overline{N}_{H}$ and $\overline{N}_{H,2}$, respectively, leads to the
results for $\widehat{\beta}_{subj}$.

##  5 Conclusion

In this study, we propose a new inverse regression framework for spatially
correlated functional data. Our framework adopts the general weighting scheme
for functional data proposed in Zhang and Wang (2016). Similarly, it extends
the mixed asymptotic framework for the spatial sampling proposed by Lu and
Tjøstheim (2014) to the functional data setting. Due to both of these
frameworks, our methodology provides a unified framework for spatially
correlated functional data which includes the existing results as special
cases. In terms of theoretical results, we provide the point-wise convergence
results, $L_{2}$ convergence results for the covariance operators and uniform
convergence result for the mean. These results are new to the literature and
hence of independent interest. Finally, we show that the estimated e.d.r.
directions achieve the $\sqrt{n}$-convergence rate for the dense and ultra-
dense functional data and achieve the standard smoothing rate for the sparse
functional data. The simulations in the supplement Chatla and Liu (2024)
validate our theoretical findings. The proposed method shows a good
performance. We find that it improves upon the existing methods when the
functional data include a nugget effect. At the same time, it also does not
lose much efficiency when the data do not include a nugget effect.

The present study considers the estimation of the e.d.r. directions and does
not discuss the issues related to inference. An interesting direction for
future research is to provide an inference framework for the estimated e.d.r.
directions.

## Acknowledgements

The authors would like to thank an anonymous referee, an Associate Editor and
the Editor for their constructive comments which led to significant
improvements in the paper. The authors would also like to thank Dr. Ci-Ren
Jiang for sharing the code for Sliced Inverse Regression. Supplement to
“Inverse regression for spatially distributed functional data” The
supplementary material contains the technical proofs of the main results, a
simulation study, and a real data analysis.

## References

  * Bosq (2012) [author] Bosq, DenisD. (2012). Nonparametric Statistics for Stochastic Processes: Estimation and Prediction. Second Edition, Springer-Verlag, New York.
  * Bowman and Azzalini (2021) [author] Bowman, A. W.A. W. and Azzalini, A.A. (2021). R package sm: nonparametric smoothing methods (version 2.2-5.7),
  * Chamberlain and Hocking (2023) [author] Chamberlain, ScottS. and Hocking, DanielD. (2023). rnoaa: ’NOAA’ Weather Data from R. R package version 1.4.0.
  * Chatla and Liu (2024) [author] Chatla, Suneel BabuS. B. and Liu, RuiqiR. (2024). Supplement to “Inverse regression for spatially distributed functional data”.
  * Chen, Hall and Müller (2011) [author] Chen, DongD., Hall, PeterP. and Müller, Hans-GeorgH.-G. (2011). Single and multiple index functional regression models with nonparametric link. Ann. Statist. 39 1720–1747.
  * Chouaf and Laksaci (2012) [author] Chouaf, AbdelhakA. and Laksaci, AliA. (2012). On the functional local linear estimate for spatial regression. Stat. Risk Model. 29 189–214.
  * Cook (2007) [author] Cook, R DennisR. D. (2007). Fisher lecture: Dimension reduction in regression. Statist. Sci. 22 1–26.
  * Cook (2018) [author] Cook, R DennisR. D. (2018). An Introduction to Envelopes: Dimension Reduction for Efficient Estimation in Multivariate Statistics. John Wiley & Sons.
  * Cook, Forzani and Yao (2010) [author] Cook, RDR., Forzani, LL. and Yao, AFA. (2010). Necessary and sufficient conditions for consistency of a method for smoothed functional inverse regression. Statist. Sinica 20 235–238.
  * Cook and Li (2002) [author] Cook, R DennisR. D. and Li, BingB. (2002). Dimension reduction for conditional mean in regression. Ann. Statist. 30 455–474.
  * Cook and Weisberg (1991) [author] Cook, R DennisR. D. and Weisberg, SanfordS. (1991). Sliced inverse regression for dimension reduction: Comment. J. Amer. Statist. Assoc. 86 328–332.
  * Cressie (2015) [author] Cressie, NoelN. (2015). Statistics for Spatial Data. John Wiley & Sons.
  * Dalenius, Hájek and Zubrzycki (1961) [author] Dalenius, ToreT., Hájek, JaroslavJ. and Zubrzycki, StefanS. (1961). On plane sampling and related geometrical problems. In Proceedings of the 4th Berkeley symposium on probability and mathematical statistics 1 125–150.
  * Delicado et al. (2010) [author] Delicado, PedroP., Giraldo, RamónR., Comas, CarlosC. and Mateu, JorgeJ. (2010). Statistics for spatial functional data: some recent contributions. Environmetrics 21 224–239.
  * Deo (1973) [author] Deo, Chandrakant MC. M. (1973). A note on empirical processes of strong-mixing sequences. Ann. Statist. 1 870–875.
  * Duan and Li (1991) [author] Duan, NaihuaN. and Li, Ker-ChauK.-C. (1991). Slicing regression: a link-free regression method. Ann. Statist. 19 505–530.
  * Fan and Gijbels (2018) [author] Fan, JianqingJ. and Gijbels, IreneI. (2018). Local Polynomial Modelling and its Applications: Monographs on Statistics and Applied Probability 66. Routledge.
  * Fan and Yao (2008) [author] Fan, JianqingJ. and Yao, QiweiQ. (2008). Nonlinear Time Series: Nonparametric and Parametric Methods. Springer Ser. Statist. Springer New York, NY.
  * Ferré and Yao (2003) [author] Ferré, LouisL. and Yao, Anne-FrançoiseA.-F. (2003). Functional sliced inverse regression analysis. Statistics 37 475–488.
  * Ferré and Yao (2005) [author] Ferré, LouisL. and Yao, Anne-FrançoiseA.-F. (2005). Smoothed functional inverse regression. Statist. Sinica 15 665–683.
  * Ferré and Yao (2007) [author] Ferré, LL. and Yao, AFA. (2007). Reply to the paper by Liliana Forzani and R. Dennis Cook:”a note on smoothed functional inverse regression”. Statist. Sinica 17 1683–1687.
  * Forzani and Cook (2007) [author] Forzani, LilianaL. and Cook, R DennisR. D. (2007). A note on smoothed functional inverse regression. Statist. Sinica 17 1677–1681.
  * Gao, Lu and Tjøstheim (2006) [author] Gao, JitiJ., Lu, ZudiZ. and Tjøstheim, DagD. (2006). Estimation in semiparametric spatial regression. Ann. Statist. 34 1395–1435.
  * Giraldo, Delicado and Mateu (2010) [author] Giraldo, RamónR., Delicado, PedroP. and Mateu, JorgeJ. (2010). Continuous time-varying kriging for spatial prediction of functional data: An environmental application. J. Agric. Biol. Environ. Stat. 15 66–82.
  * Hall and Horowitz (2007) [author] Hall, PeterP. and Horowitz, Joel LJ. L. (2007). Methodology and convergence rates for functional linear regression. Ann. Statist. 35 70–91.
  * Hallin, Lu and Tran (2004) [author] Hallin, MarcM., Lu, ZudiZ. and Tran, Lanh TL. T. (2004). Local linear spatial regression. Ann. Statist. 32 2469–2500.
  * Hansen (2008) [author] Hansen, Bruce EB. E. (2008). Uniform convergence rates for kernel estimation with dependent data. Econometric Theory 24 726–748.
  * Hörmann and Kokoszka (2013) [author] Hörmann, SiegfriedS. and Kokoszka, PiotrP. (2013). Consistency of the mean and the principal components of spatially distributed functional data. Bernoulli 19 1535–1558.
  * Hsing and Ren (2009) [author] Hsing, TailenT. and Ren, HaoboH. (2009). An RKHS formulation of the inverse regression dimension-reduction problem. Ann. Statist. 37 726–755.
  * Ibragimov (1975) [author] Ibragimov, II. (1975). Independent and Stationary Sequences of Random Variables. Wolters, Noordhoff Pub.
  * Jiang and Wang (2015) [author] Jiang, Ci-RenC.-R. and Wang, Jane-LingJ.-L. (2015). Correction: Inverse regression for longitudinal data: Ann. Statist. 42 (2014) 563-591. Ann. Statist. 43 2326–2329.
  * Jiang, Yu and Wang (2014) [author] Jiang, Ci-RenC.-R., Yu, WeiW. and Wang, Jane-LingJ.-L. (2014). Inverse regression for longitudinal data. Ann. Statist. 42 563–591.
  * Kim et al. (2020) [author] Kim, KyongwonK., Li, BingB., Yu, ZhouZ. and Li, LexinL. (2020). On post dimension reduction statistical inference. Ann. Statist. 48 1567–1592.
  * Kurisu (2022a) [author] Kurisu, DaisukeD. (2022a). Nonparametric regression for locally stationary random fields under stochastic sampling design. Bernoulli 28 1250–1275.
  * Kurisu (2022b) [author] Kurisu, DaisukeD. (2022b). Nonparametric regression for locally stationary functional time series. Electron. J. Stat. 16 3973–3995.
  * Lahiri (1996) [author] Lahiri, Soumendra NathS. N. (1996). On inconsistency of estimators based on spatial data under infill asymptotics. Sankhya A 58 403–417.
  * Laksaci, Rachdi and Rahmani (2013) [author] Laksaci, AliA., Rachdi, MustaphaM. and Rahmani, SaâdiaS. (2013). Spatial modelization: local linear estimation of the conditional distribution for functional data. Spat. Stat. 6 1–23.
  * Lee and Li (2022) [author] Lee, Kuang-YaoK.-Y. and Li, LexinL. (2022). Functional sufficient dimension reduction through average Fréchet derivatives. Ann. Statist. 50 904–929.
  * Li (1991) [author] Li, Ker-ChauK.-C. (1991). Sliced inverse regression for dimension reduction. J. Amer. Statist. Assoc. 86 316–327.
  * Li (2018) [author] Li, BingB. (2018). Sufficient Dimension Reduction: Methods and Applications with R. Chapman & Hall/CRC.
  * Li and Hsing (2010) [author] Li, YehuaY. and Hsing, TailenT. (2010). Uniform convergence rates for nonparametric regression and principal component analysis in functional/longitudinal data. Ann. Statist. 38 3321–3351.
  * Li and Song (2017) [author] Li, BingB. and Song, JunJ. (2017). Nonlinear sufficient dimension reduction for functional data. Ann. Statist. 45 1059–1095.
  * Li and Song (2022) [author] Li, BingB. and Song, JunJ. (2022). Dimension reduction for functional data based on weak conditional moments. Ann. Statist. 50 107–128.
  * Liu, Ray and Hooker (2017) [author] Liu, ChongC., Ray, SurajitS. and Hooker, GilesG. (2017). Functional principal component analysis of spatially correlated data. Stat. Comput. 27 1639–1654.
  * Lu and Tjøstheim (2014) [author] Lu, ZudiZ. and Tjøstheim, DagD. (2014). Nonparametric estimation of probability density functions for irregularly observed spatial data. J. Amer. Statist. Assoc. 109 1546–1564.
  * Masry (1986) [author] Masry, EliasE. (1986). Recursive probability density estimation for weakly dependent stationary processes. IEEE Trans. Inform. Theory 32 254–267.
  * Matérn (2013) [author] Matérn, BertilB. (2013). Spatial Variation 36\. Springer New York, NY.
  * Müller and Stadtmüller (2005) [author] Müller, Hans-GeorgH.-G. and Stadtmüller, UlrichU. (2005). Generalized functional linear models. Ann. Statist. 33 774–805.
  * Novak (2006) [author] Novak, ErichE. (2006). Deterministic and Stochastic Error Bounds in Numerical Analysis 1349\. Springer Berlin, Heidelberg.
  * Quenouille (1949) [author] Quenouille, Maurice HM. H. (1949). Approximate tests of correlation in time-series. J. R. Stat. Soc. Ser. B. Stat. Methodol. 11 68–84.
  * Ramsay and Silverman (2005) [author] Ramsay, Jim OJ. O. and Silverman, Bernard WB. W. (2005). Functional Data Analysis. Springer Ser. Statist. Springer-Verlag New York.
  * Schlather et al. (2015) [author] Schlather, MartinM., Malinowski, AlexanderA., Menck, Peter J.P. J., Oesting, MarcoM. and Strokorb, KirstinK. (2015). Analysis, simulation and prediction of multivariate random fields with package RandomFields. J. Stat. Softw. 63 1–25.
  * Song (2019) [author] Song, JunJ. (2019). On sufficient dimension reduction for functional data: Inverse moment-based methods. Wiley Interdiscip. Rev. Comput. Stat. 11 e1459.
  * R Core Team (2022) [author] R Core Team (2022). R: a Language and environment for statistical computing.
  * Tran (1990) [author] Tran, Lanh TatL. T. (1990). Kernel density estimation on random fields. J. Multivariate Anal. 34 37–53.
  * Traub and Werschulz (1998) [author] Traub, Joseph FrederickJ. F. and Werschulz, Arthur GA. G. (1998). Complexity and Information 26862\. Cambridge University Press.
  * Wainwright (2019) [author] Wainwright, Martin JM. J. (2019). High-Dimensional Statistics: A Non-Asymptotic Viewpoint 48\. Cambridge University Press.
  * Yao, Müller and Wang (2005) [author] Yao, FangF., Müller, Hans-GeorgH.-G. and Wang, Jane-LingJ.-L. (2005). Functional data analysis for sparse longitudinal data. J. Amer. Statist. Assoc. 100 577–590.
  * Zhang (2004) [author] Zhang, HaoH. (2004). Inconsistent estimation and asymptotically equal interpolations in model-based geostatistics. J. Amer. Statist. Assoc. 99 250–261.
  * Zhang and Li (2022) [author] Zhang, HaozheH. and Li, YehuaY. (2022). Unified principal component analysis for sparse and dense functional data under spatial dependency. J. Bus. Econom. Statist. 40 1523–1537.
  * Zhang and Wang (2016) [author] Zhang, XiaokeX. and Wang, Jane-LingJ.-L. (2016). From sparse to dense functional data and beyond. Ann. Statist. 44 2281–2321.
  * Zhang and Wang (2018) [author] Zhang, XiaokeX. and Wang, Jane-LingJ.-L. (2018). Optimal weighting schemes for longitudinal and functional data. Statist. Probab. Lett. 138 165–170.
  * Zhou et al. (2022) [author] Zhou, YidongY., Bhattacharjee, SatarupaS., Carroll, CodyC., Chen, YaqingY., Dai, XiongtaoX., Fan, JianingJ., Gajardo, AlvaroA., Hadjipantelis, Pantelis Z.P. Z., Han, KyungheeK., Ji, HaoH., Zhu, ChangboC., Müller, Hans-GeorgH.-G. and Wang, Jane-LingJ.-L. (2022). Fdapace: functional data analysis and empirical dynamics. R package version 0.5.9.

##  S.1 Proofs

In this supplement, we present the proofs of our main results in Section 4 of
the manuscript and provide the results from a simulation study, and the
weather data analysis. The asymptotic results for the covariance operators
$\widehat{R}$ and $\widehat{R}_{e}$ require uniform convergence result for the
mean estimator. Therefore, in Section S.1.1 we provide the uniform convergence
result for the mean estimator $\widehat{\mu}$. Section S.1.2 presents the
pointwise and $L_{2}$ convergence results for the spatial covariance matrix
$\widehat{R}$. Similarly, Section S.1.3 provides the pointwise and $L_{2}$
convergence results for the conditional mean function $E(X(t)|Y)$ and the
conditional covariance matrix $\widehat{R}_{e}$, respectively. Finally, we
present the results from the simulation study and the data analysis in
Sections S.2 and S.3, respectively.

The following two lemmas are required for our theoretical arguments. For
completeness, we present them here.

######  Lemma S.1.1.

Let $\mathcal{L}_{r}(\mathcal{F})$ denote the class of $\mathcal{F}$\-
measurable random variables $W$ satisfying
$\|W\|_{r}:=(E|W|^{r})^{1/{r}}<\infty$. Let
$\bm{U}\in\mathcal{L}_{a}(\mathcal{B}(\mathcal{S}))$ and
$V\in\mathcal{L}_{b}(\mathcal{B}(\mathcal{S^{\prime}}))$ where
$\mathcal{B}(\mathcal{S})$ and $\mathcal{B}(\mathcal{S^{\prime}})$ denote the
$\sigma$\- fields. Then for any $1\leq a,b,c<\infty$ such that
$a^{-1}+b^{-1}+c^{-1}=1$,

$\displaystyle|E[\bm{U}V]-E[\bm{U}]E[V]|\leq
C_{\alpha}\|\bm{U}\|_{a}\|V\|_{b}[\alpha(\mathcal{S},\mathcal{S^{\prime}})]^{1/{c}},$
(S.1)

where $\|\bm{U}\|_{a}:=\|(\bm{U}^{T}\bm{U})^{1/2}\|_{a}$ and
$\alpha(\mathcal{S},\mathcal{S^{\prime}})=\sup\\{|P(AB)-P(A)P(B)|:A\in\mathcal{B}(S),B\in\mathcal{B}(S^{\prime})\\}$
and $C_{\alpha}$ is some constant.

###### Proof.

It is taken from Theorem 17.2.3 in Ibragimov (1975) and Lemma 1 in Deo (1973)
shared identical results. ∎

######  Lemma S.1.2.

Suppose that $|X|\leq c$ for a random variable $X$, then for any $\lambda>0$
we have

$\displaystyle
E(e^{\lambda(X-E(X))})\leq\exp\left\\{\lambda^{2}\sigma^{2}\left(\frac{e^{\lambda
c}-1-\lambda c}{(\lambda c)^{2}}\right)\right\\}$

where $\sigma^{2}=var(X)$.

###### Proof.

This is the Bennett’s inequality; please see Page 51 of Wainwright (2019). ∎

###  S.1.1 Results for $\widehat{\mu}(t)$

We now provide the uniform convergence result for the mean $\widehat{\mu}(t)$,
$t\in[0,1]$.

######  Lemma S.1.3.

Under conditions (I)[(i),(ii), (iii), (iv)(c)], (III)(i), (IV)(ii) and (V)(i),

$\displaystyle\underset{t\in[0,1]}{\sup}|\widehat{\mu}(t)-\mu(t)|$
$\displaystyle=O_{p}\left(h_{\mu}^{2}+\log(n)\sqrt{\frac{\sum_{i=1}^{n}N_{i}w_{i}^{2}}{h_{\mu}}+\sum_{i=1}^{n}N_{i}(N_{i}-1)w_{i}^{2}}\right).$

###### Proof.

The proof requires the following two results. The arguments analogous to
Theorem 5.1 in Zhang and Wang (2016) together with Lemma S.1.5 and Proposition
S.1.4 yield the desired result. ∎

######  Proposition S.1.4.

Suppose conditions (I)[(i),(ii), (iii), (iv)(c)], (III)(i), (IV)(ii) and
(V)(i) hold. Define

$\displaystyle S_{n}(t)$
$\displaystyle=\sum_{i=1}^{n}w_{i}\sum_{j=1}^{N_{i}}K_{h_{\mu}}(T_{ij}-t)\mathcal{U}_{ij},$

where $\mathcal{U}_{ij}=V_{ij}+U_{ij}$. Then it follows that

$\displaystyle\sup_{t}|S_{n}(t)|$
$\displaystyle=O_{P}\left(\log(n)\sqrt{\frac{\sum_{i=1}^{n}N_{i}w_{i}^{2}}{h_{\mu}}+\sum_{i=1}^{n}N_{i}(N_{i}-1)w_{i}^{2}}\right).$
(S.2)

Moreover, it holds that

$\displaystyle\sup_{t}\left|\sum_{i=1}^{n}w_{i}\sum_{j=1}^{N_{i}}K_{h_{\mu}}(T_{ij}-t)e_{ij}\right|=O_{P}\left(\log(n)\sqrt{\frac{\sum_{i=1}^{n}N_{i}w_{i}^{2}}{h_{\mu}}+\sum_{i=1}^{n}N_{i}(N_{i}-1)w_{i}^{2}}\right).$
(S.3)

###### Proof.

Statement (S.3) follows from Lemma 5 in Zhang and Wang (2016). It suffices to
prove Statement (S.2), and the proof is divided into three steps. For
simplicity in notation, we denote $h_{\mu}=h$ in the below proof.

Step 1: Consider the following decomposition:
$S_{n}(t)=S_{n}^{(1)}(t)+S_{n}^{(2)}(t)$, where

$\displaystyle S_{n}^{(1)}(t)$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}w_{i}\sum_{j=1}^{N_{i}}\\{K_{h}(T_{ij}-t)\mathcal{U}_{ij}I(|\mathcal{U}_{ij}|\leq
B_{n})-E[K_{h}(T_{ij}-t)\mathcal{U}_{ij}I(|\mathcal{U}_{ij}|\leq B_{n})]\\},$
$\displaystyle S_{n}^{(2)}(t)$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}w_{i}\sum_{j=1}^{N_{i}}\\{K_{h}(T_{ij}-t)\mathcal{U}_{ij}I(|\mathcal{U}_{ij}|>B_{n})-E[K_{h}(T_{ij}-t)\mathcal{U}_{ij}I(|\mathcal{U}_{ij}|>B_{n})]\\}.$

Let $D_{i}(\bm{s}_{i})=\sup_{t}|V(\bm{s}_{i},t)|+\sup_{t}|U_{i}(t)|$ and
$\sup_{\bm{s}}E[|D_{i}(\bm{s})|^{2+\varrho}]<\infty$ by Condition
(I)[(iv)(c)]. Hence, we see that

$\displaystyle\sup_{t}|S_{n}^{(2)}(t)|\leq
2h^{-1}\sum_{i=1}^{n}w_{i}\sum_{j=1}^{N_{i}}B_{n}^{-(1+\varrho)}|\mathcal{U}_{ij}|^{2+\varrho}\leq
2h^{-1}\sum_{i=1}^{n}w_{i}\sum_{j=1}^{N_{i}}B_{n}^{-(1+\varrho)}|D_{i}(\bm{s}_{i})|^{2+\varrho}.$

Taking expectation, it follows that

$\displaystyle E[\sup_{t}|S_{n}^{(2)}(t)|]\leq 2Ch^{-1}B_{n}^{-(1+\varrho)}.$
(S.4)

Step 2: Let $t_{1},t_{2},\ldots,t_{r}$ be an equidistance partition of
$[0,1]$, and it follows that

$\displaystyle\sup_{t\in[0,1]}|S_{n}^{(1)}(t)|\leq\sup_{1\leq i\leq
r}|S_{n}^{(1)}(t_{i})|+\sup_{|t-\widetilde{t}|\leq
1/r}|S_{n}^{(1)}(t)-S_{n}^{(1)}(\widetilde{t})|.$

The Lipschitz conditions of kernel implies that

$\displaystyle\sup_{|t-\widetilde{t}|\leq
1/r}|S_{n}^{(1)}(t)-S_{n}^{(1)}(\widetilde{t})|\leq\sup_{|t-\widetilde{t}|\leq
1/r}\sum_{i=1}^{n}w_{i}\sum_{j=1}^{N_{i}}\frac{C}{h^{2}}|t-\widetilde{t}||\mathcal{U}_{ij}|\leq
Ch^{-2}r^{-1}\sum_{i=1}^{n}w_{i}\sum_{j=1}^{N_{i}}|\mathcal{U}_{ij}|.$

Taking expectation, we show that

$\displaystyle\sup_{|t-\widetilde{t}|\leq
1/r}|S^{(1)}_{n}(t)-S_{n}^{(1)}(\widetilde{t})|=O_{P}(h^{-2}r^{-1}\sum_{i=1}^{n}w_{i}N_{i})=O_{P}(h^{-2}r^{-1}).$
(S.5)

Using Lemma S.1.5, it follows that

$\displaystyle\mathbb{P}(\sup_{1\leq i\leq
r}|S_{n}^{(1)}(t_{i})|>8\log(n)Q_{n})$ $\displaystyle\leq$ $\displaystyle
2r\exp\left\\{-\frac{\log^{2}(n)Q_{n}^{2}}{2Q_{n}^{2}+\frac{2}{3}CB_{n}n^{-1}h^{-1}p^{2}\delta_{n}^{-2}\log(n)Q_{n}}\right\\}+\frac{Cr\Delta_{n}^{2}}{p^{2}}\psi(Cp^{2}\delta_{n}^{-2},n)\phi(p))\sqrt{\frac{B_{n}n^{-1}h^{-1}\Delta_{n}^{2}\delta_{n}^{-2}}{\log(n)Q_{n}}}.$

Combining (S.4)-(LABEL:eq:proposition:uniform:rate:step1:eq3), we show that

$\displaystyle\sup_{t}|S_{n}(t)|=O_{P}\left(\log(n)Q_{n}+h^{-2}r^{-1}+h^{-1}B_{n}^{-(1+\varrho)}\right),$
(S.7)

provided the RHS of (LABEL:eq:proposition:uniform:rate:step1:eq3) goes to
zero.

Step 3: The convergence of RHS of
(LABEL:eq:proposition:uniform:rate:step1:eq3) is implied by the following
sufficient conditions:

$\displaystyle\log(r)-\log^{2}(n)$ $\displaystyle\to$ $\displaystyle-\infty,$
$\displaystyle\log(r)-\frac{\log(n)nQ_{n}h\delta_{n}^{2}}{B_{n}p^{2}}$
$\displaystyle\to$ $\displaystyle-\infty,$
$\displaystyle\frac{B_{n}\Delta_{n}^{6}\delta_{n}^{-2}r^{2}\psi^{2}(Cp^{2}\delta_{n}^{-2},n)\phi^{2}(p)}{\log(n)nQ_{n}hp^{4}}$
$\displaystyle\to$ $\displaystyle 0.$

Taking $r=h^{-2}Q_{n}^{-1}$, $B_{n}=(Q_{n}h)^{-\frac{1}{1+\varrho}}$, (S.7)
becomes

$\displaystyle\sup_{t}|S_{n}(t)|=O_{P}\left(\log(n)Q_{n}\right),$

and the sufficient conditions become

$\displaystyle\log(r)-\log^{2}(n)$ $\displaystyle\to$ $\displaystyle-\infty,$
$\displaystyle\log(r)-\frac{\log(n)nQ_{n}^{\frac{2+\varrho}{1+\varrho}}h^{\frac{2+\varrho}{1+\varrho}}\delta_{n}^{2}}{p^{2}}$
$\displaystyle\to$ $\displaystyle-\infty,$
$\displaystyle\frac{\Delta_{n}^{6}\delta_{n}^{-2}\psi^{2}(Cp^{2}\delta_{n}^{-2},n)\phi^{2}(p)}{\log(n)nQ_{n}^{\frac{4+3\varrho}{1+\varrho}}h^{\frac{6+5\varrho}{1+\varrho}}p^{4}}$
$\displaystyle\to$ $\displaystyle 0.$

Finally, using the fact that $\psi(i,j)\leq\min(i,j)$ and $\alpha(p)\leq
p^{-\rho}$, the above conditions can be implied by

$\displaystyle\frac{nQ_{n}^{\frac{2+\varrho}{1+\varrho}}h^{\frac{2+\varrho}{1+\varrho}}\delta_{n}^{2}}{p^{2}}\to\infty,\frac{\Delta_{n}^{6}\delta_{n}^{-2}}{\log(n)nQ_{n}^{\frac{4+3\varrho}{1+\varrho}}h^{\frac{6+5\varrho}{1+\varrho}}p^{2\rho}}\to
0.$

To guarantee existence such $p$, we need

$\displaystyle
nQ_{n}^{\frac{2+\varrho}{1+\varrho}}h^{\frac{2+\varrho}{1+\varrho}}\delta_{n}^{2}\to\infty,[nQ_{n}^{\frac{2+\varrho}{1+\varrho}}h^{\frac{2+\varrho}{1+\varrho}}]^{\rho}>>\frac{\Delta_{n}^{6}\delta_{n}^{-2}}{\log(n)nQ_{n}^{\frac{4+3\varrho}{1+\varrho}}h^{\frac{6+5\varrho}{1+\varrho}}},$

which are identical to

$\displaystyle
nQ_{n}^{\frac{2+\varrho}{1+\varrho}}h^{\frac{2+\varrho}{1+\varrho}}\delta_{n}^{2}\to\infty,\frac{n^{\rho+1}Q_{n}^{\frac{(2+\varrho)\rho+4+3\varrho}{1+\varrho}}h^{\frac{(2+\varrho)\rho+6+5\varrho}{1+\varrho}}\delta_{n}^{2}\log(n)}{\Delta_{n}^{6}}\to\infty.$

For example, when $\varrho\to\infty$, the conditions become

$\displaystyle
nQ_{n}h\delta_{n}^{2}\to\infty,\frac{n^{\rho+1}Q_{n}^{\rho+3}h^{\rho+5}\delta_{n}^{2}\log(n)}{\Delta_{n}^{6}}\to\infty.$

∎

######  Lemma S.1.5.

Suppose conditions (I)[(i),(ii), (iii), (iv)(c)], (III)(i), (IV)(ii) and
(V)[(i),(ii)] hold. Define
$S_{n}(t)=\sum_{i=1}^{n}w_{i}\sum_{j=1}^{N_{i}}K_{h}(T_{ij}-t)\mathcal{U}_{ij}$,
where $\mathcal{U}_{ij}:=\mathcal{U}(\bm{s}_{i},T_{ij})=V_{ij}+U_{ij}$. Then
the following holds for any $\epsilon,p>0$ and some constant $C>0$:

$\displaystyle\mathbb{P}(|S_{n}(t)|>8\epsilon)\leq
2\exp\left\\{-\frac{\epsilon^{2}}{2Q_{n}^{2}+\frac{2}{3}CB_{n}h^{-1}p^{2}\delta_{n}^{-2}\epsilon}\right\\}+\frac{C\Delta_{n}^{2}}{p^{2}}\psi(Cp^{2}\delta_{n}^{-2},n)\phi(p)\sqrt{\frac{B_{n}n^{-1}h^{-1}\Delta_{n}^{2}\delta_{n}^{2}}{\epsilon}}.$

where
$Q_{n}^{2}=\sum_{i=1}^{n}w_{i}^{2}N_{i}h^{-1}+\sum_{i=1}^{n}w_{i}^{2}N_{i}(N_{i}-1)$.

###### Proof.

The proof is divided into three steps.

Step 1: By definition of $\Delta_{n}$ and without loss of generality, we can
assume $\bm{s}_{1},\ldots,\bm{s}_{n}\in[0,C\Delta_{n}]^{2}:=\Omega_{n}$ for
some $C>0$. If not, we can shift the points. For each pair $i,k$, we define
rectangles

$\displaystyle A_{ku}^{(1)}$ $\displaystyle=$
$\displaystyle\\{\bm{s}=(s_{1},s_{2}):2(k-1)p\leq s_{1}<(2k-1)p,2(u-1)p\leq
s_{2}<(2u-1)p\\},$ $\displaystyle A_{ku}^{(2)}$ $\displaystyle=$
$\displaystyle\\{\bm{s}=(s_{1},s_{2}):2(k-1)p\leq s_{1}<(2k-1)p,(2u-1)p\leq
s_{2}<2up\\},$ $\displaystyle A_{ku}^{(3)}$ $\displaystyle=$
$\displaystyle\\{\bm{s}=(s_{1},s_{2}):(2k-1)p\leq s_{1}<2kp,2(u-1)p\leq
s_{2}<(2u-1)p\\},$ $\displaystyle A_{ku}^{(4)}$ $\displaystyle=$
$\displaystyle\\{\bm{s}=(s_{1},s_{2}):(2k-1)p\leq s_{1}<2kp,(2u-1)p\leq
s_{2}<2up\\}.$

Correspondingly, let us define

$\displaystyle V_{ku}^{(1)}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}I(\bm{s}_{i}\in
A_{ku}^{(1)})w_{i}\sum_{j=1}^{N_{i}}K_{h}(T_{ij}-t)\mathcal{U}_{ij},$
$\displaystyle V_{ku}^{(2)}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}I(\bm{s}_{i}\in
A_{ku}^{(2)})w_{i}\sum_{j=1}^{N_{i}}K_{h}(T_{ij}-t)\mathcal{U}_{ij},$
$\displaystyle V_{ku}^{(3)}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}I(\bm{s}_{i}\in
A_{ku}^{(3)})w_{i}\sum_{j=1}^{N_{i}}K_{h}(T_{ij}-t)\mathcal{U}_{ij},$
$\displaystyle V_{ku}^{(4)}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}I(\bm{s}_{i}\in
A_{ku}^{(4)})w_{i}\sum_{j=1}^{N_{i}}K_{h}(T_{ij}-t)\mathcal{U}_{ij}.$

It is not difficult to verify that

$\displaystyle
S_{n}(t)=\sum_{k=1}^{M}\sum_{u=1}^{M}(V_{ku}^{(1)}+V_{ku}^{(2)}+V_{ku}^{(3)}+V_{ku}^{(4)}),$
(S.8)

where $M$ is an integer such that $M\leq C\Delta_{n}/p$. Noting that

$|{V}_{ku}^{(1)}|\leq 2B_{n}h^{-1}\sum_{i=1}^{n}I(\bm{s}_{i}\in
A_{ku}^{(1)})w_{i}N_{i}\leq
B_{n}h^{-1}Cp^{2}\delta_{n}^{-2}\max_{i}w_{i}N_{i}\leq
CB_{n}n^{-1}h^{-1}p^{2}\delta_{n}^{-2},$

Lemma 1.2 of Bosq (2012) implies that there is a sequence of independent
variables $\\{W_{ku}\\}$ such that ${V}_{ku}^{(1)}$ and $W_{ku}$ share the
same distribution, and further

$\displaystyle\mathbb{P}(|{V}_{ku}^{(1)}-W_{ku}|>x)$ $\displaystyle\leq$
$\displaystyle
C\sqrt{\frac{B_{n}n^{-1}h^{-1}p^{2}\delta_{n}^{-2}}{x}}\alpha(A_{ku}^{(1)},\Omega_{n}/A_{ku}^{(1)})$
(S.9) $\displaystyle\leq$ $\displaystyle
C\sqrt{\frac{B_{n}n^{-1}h^{-1}p^{2}\delta_{n}^{-2}}{x}}\psi(Cp^{2}\delta_{n}^{-2},n)\phi(p).$

Step 2: By direct examination, we show that

$\displaystyle E(|W_{ku}|^{2})=E(|{V}_{ku}^{(1)}|^{2})$ $\displaystyle\leq$
$\displaystyle E\left\\{\left(\sum_{i=1}^{n}I(\bm{s}_{i}\in
A_{ku}^{(1)})w_{i}\sum_{j=1}^{N_{i}}K_{h}(T_{ij}-t)\mathcal{U}_{ij}\right)^{2}\right\\}$
$\displaystyle\leq$ $\displaystyle\sum_{i=1}^{n}I(\bm{s}_{i}\in
A_{ku}^{(1)})w_{i}^{2}Var\left(\sum_{j=1}^{N_{i}}K_{h}(T_{ij}-t)\mathcal{U}_{ij}\right)$
$\displaystyle+\sum_{0<d(\bm{s}_{i},\bm{s}_{i^{\prime}})\leq
c_{n}}Cov\left(I(\bm{s}_{i}\in
A_{ku}^{(1)})w_{i}\sum_{j=1}^{N_{i}}K_{h}(T_{ij}-t)\mathcal{U}_{ij},I(\bm{s}_{i^{\prime}}\in
A_{ku}^{(1)})w_{i^{\prime}}\sum_{j=1}^{N_{i}}K_{h}(T_{{i^{\prime}}j}-t)\mathcal{U}_{{i^{\prime}}j}\right)$
$\displaystyle+\sum_{d(\bm{s}_{i},\bm{s}_{i^{\prime}})>c_{n}}Cov\left(I(\bm{s}_{i}\in
A_{ku}^{(1)})w_{i}\sum_{j=1}^{N_{i}}K_{h}(T_{ij}-t)\mathcal{U}_{ij},I(\bm{s}_{i^{\prime}}\in
A_{ku}^{(1)})w_{i^{\prime}}\sum_{j=1}^{N_{i}}K_{h}(T_{{i^{\prime}}j}-t)\mathcal{U}_{{i^{\prime}}j}\right)$
$\displaystyle\leq$ $\displaystyle\sum_{i=1}^{n}I(\bm{s}_{i}\in
A_{ku}^{(1)})w_{i}^{2}\sum_{j=1}^{N_{i}}\sum_{j^{\prime}=1}^{N_{i}^{\prime}}E\left\\{|K_{h}(T_{ij}-t)\mathcal{U}_{ij}K_{h}(T_{{i}j^{\prime}}-t)\mathcal{U}_{{i}j^{\prime}}|\right\\}$
$\displaystyle+\sum_{0<d(\bm{s}_{i},\bm{s}_{i^{\prime}})\leq
c_{n}}I(\bm{s}_{i}\in A_{ku}^{(1)})I(\bm{s}_{i^{\prime}}\in
A_{ku}^{(1)})w_{i}w_{i^{\prime}}\sum_{j=1}^{N_{i}}\sum_{j^{\prime}=1}^{N_{i}^{\prime}}E\left\\{|K_{h}(T_{ij}-t)\mathcal{U}_{ij}K_{h}(T_{{i^{\prime}}j^{\prime}}-t)\mathcal{U}_{{i^{\prime}}j^{\prime}}|\right\\}$
$\displaystyle+\sum_{d(\bm{s}_{i},\bm{s}_{i^{\prime}})>c_{n}}I(\bm{s}_{i}\in
A_{ku}^{(1)})I(\bm{s}_{i^{\prime}}\in
A_{ku}^{(1)})w_{i}w_{i^{\prime}}\sum_{j=1}^{N_{i}}\sum_{j^{\prime}=1}^{N_{i}^{\prime}}E\left\\{|K_{h}(T_{ij}-t)\mathcal{U}_{ij}K_{h}(T_{{i^{\prime}}j^{\prime}}-t)\mathcal{U}_{{i^{\prime}}j^{\prime}}|\right\\}.$

Taking summation, we show that

$\displaystyle\sum_{k=1}^{M}\sum_{u=1}^{M}E(|W_{ku}|^{2})$ $\displaystyle\leq$
$\displaystyle\sum_{i=1}^{n}w_{i}^{2}\sum_{j=1}^{N_{i}}\sum_{j^{\prime}=1}^{N_{i}}E\left\\{|K_{h}(T_{ij}-t)\mathcal{U}_{ij}K_{h}(T_{{i}j^{\prime}}-t)\mathcal{U}_{{i}j^{\prime}}|\right\\}$
(S.10) $\displaystyle+\sum_{0<d(\bm{s}_{i},\bm{s}_{i^{\prime}})\leq
c_{n}}w_{i}w_{i^{\prime}}\sum_{j=1}^{N_{i}}\sum_{j^{\prime}=1}^{N_{i}^{\prime}}E\left\\{|K_{h}(T_{ij}-t)\mathcal{U}_{ij}K_{h}(T_{{i^{\prime}}j^{\prime}}-t)\mathcal{U}_{{i^{\prime}}j^{\prime}}|\right\\}$
$\displaystyle+\sum_{d(\bm{s}_{i},\bm{s}_{i^{\prime}})>c_{n}}w_{i}w_{i^{\prime}}\sum_{j=1}^{N_{i}}\sum_{j^{\prime}=1}^{N_{i}^{\prime}}E\left\\{|K_{h}(T_{ij}-t)\mathcal{U}_{ij}K_{h}(T_{{i^{\prime}}j^{\prime}}-t)\mathcal{U}_{{i^{\prime}}j^{\prime}}|\right\\}$
$\displaystyle:=$ $\displaystyle V_{1}+V_{2}+V_{3}.$

Notice that if $j\neq j^{\prime}$, it follows that

$\displaystyle
E\left\\{|K_{h}(T_{ij}-t)\mathcal{U}_{ij}K_{h}(T_{{i}j^{\prime}}-t)\mathcal{U}_{{i}j^{\prime}}|\right\\}$
$\displaystyle=$ $\displaystyle
E\left\\{\int_{0}^{1}\int_{0}^{1}\frac{1}{h^{2}}K\left(\frac{x_{1}-t}{h}\right)K\left(\frac{x_{2}-t}{h}\right)|\mathcal{U}(\bm{s}_{i},x_{1})\mathcal{U}(\bm{s}_{i},x2)|f_{T}(x_{1})f_{T}(x_{2})dx_{1}dx_{2}\right\\}$
$\displaystyle=$ $\displaystyle
E\left\\{\int_{-C}^{C}\int_{-C}^{C}K\left(y_{1}\right)K\left(y_{2}\right)|\mathcal{U}(\bm{s}_{i},t+hy_{1})\mathcal{U}(\bm{s}_{i},t+hy_{2})|f_{T}(t+hy_{1})f_{T}(t+hy_{2})dy_{1}dy_{2}\right\\}$
$\displaystyle\leq$ $\displaystyle
CE\left\\{\int_{-C}^{C}\int_{-C}^{C}|\mathcal{U}(\bm{s}_{i},t+hy_{1})\mathcal{U}(\bm{s}_{i},t+hy_{2})|dy_{1}dy_{2}\right\\}$
$\displaystyle\leq$ $\displaystyle
C\left\\{\int_{-C}^{C}\int_{-C}^{C}E(|\mathcal{U}(\bm{s}_{i},t+hy_{1})\mathcal{U}(\bm{s}_{i},t+hy_{2})|)dy_{1}dy_{2}\right\\}$
$\displaystyle\leq$ $\displaystyle
C\int_{-C}^{C}\sup_{\bm{s},t}E\left\\{|\mathcal{U}(\bm{s},t)|^{2}\right\\}dy\leq
C.$

Similarly, if $j=j^{\prime}$, we have

$\displaystyle
E\left\\{|K_{h}(T_{ij}-t)\mathcal{U}_{ij}K_{h}(T_{{i}j^{\prime}}-t)\mathcal{U}_{{i}j^{\prime}}|\right\\}$
$\displaystyle=$ $\displaystyle
E\left\\{K_{h}^{2}(T_{ij}-t)\mathcal{U}_{ij}^{2}\right\\}$ $\displaystyle=$
$\displaystyle
E\left\\{\int_{0}^{1}\frac{1}{h^{2}}K^{2}\left(\frac{x-t}{h}\right)|\mathcal{U}^{2}(\bm{s}_{i},x)|f_{T}(x)dx\right\\}$
$\displaystyle=$ $\displaystyle
E\left\\{\int_{-C}^{C}\frac{1}{h}K^{2}\left(y\right)\mathcal{U}^{2}(\bm{s}_{i},t+hy)f_{T}(x)dx\right\\}$
$\displaystyle\leq$ $\displaystyle
Ch^{-1}E\left\\{\int_{-C}^{C}\mathcal{U}^{2}(\bm{s}_{i},t+hy)dx\right\\}\leq
Ch^{-1}.$

By the same calculation, we show that

$\displaystyle
E\left\\{|K_{h}(T_{ij}-t)\mathcal{U}_{ij}K_{h}(T_{{i^{\prime}}j^{\prime}}-t)\mathcal{U}_{{i^{\prime}}j^{\prime}}|\right\\}$
$\displaystyle\leq$ $\displaystyle C,\textrm{if }i\neq i^{\prime},$
$\displaystyle
E\left\\{|K_{h}^{2+\varrho}(T_{ij}-t)|\mathcal{U}_{ij}|^{2+\varrho}\right\\}$
$\displaystyle\leq$ $\displaystyle
Ch^{-(1+\varrho)}E\left\\{\int_{-C}^{C}|\mathcal{U}(\bm{s}_{i},t+hy)|^{2+\varrho}dy\right\\}$
$\displaystyle\leq$ $\displaystyle Ch^{-(1+\varrho)}.$

Hence, if $i\neq i^{\prime}$, we have

$\displaystyle
E\left\\{|K_{h}(T_{ij}-t)\mathcal{U}_{ij}K_{h}(T_{{i^{\prime}}j^{\prime}}-t)\mathcal{U}_{{i^{\prime}}j^{\prime}}|\right\\}$
$\displaystyle\leq$ $\displaystyle
E^{\frac{1}{2+\varrho}}\left\\{|K_{h}^{2+\varrho}(T_{ij}-t)|\mathcal{U}_{ij}|^{2+\varrho}\right\\}E^{\frac{1}{2+\varrho}}\left\\{|K_{h}^{2+\varrho}(T_{i^{\prime}j^{\prime}}-t)|\mathcal{U}_{i^{\prime}j^{\prime}}|^{2+\varrho}\right\\}\phi^{\frac{\varrho}{2+\varrho}}(d(\bm{s}_{i},\bm{s}_{i^{\prime}}))$
$\displaystyle\leq$ $\displaystyle
Ch^{-\frac{2(1+\varrho)}{2+\varrho}}\phi^{\frac{\varrho}{2+\varrho}}(d(\bm{s}_{i},\bm{s}_{i^{\prime}})).$

Using the above bounds and the definitions of $V_{1},V_{2},V_{3}$ in (S.10),
we show that

$\displaystyle V_{1}$ $\displaystyle\leq$ $\displaystyle
C\sum_{i=1}^{n}w_{i}^{2}N_{i}h^{-1}+C\sum_{i=1}^{n}w_{i}^{2}N_{i}(N_{i}-1),$
$\displaystyle V_{2}$ $\displaystyle\leq$ $\displaystyle
C\sum_{0<d(\bm{s}_{i},\bm{s}_{i^{\prime}})\leq
c_{n}}w_{i}w_{i^{\prime}}N_{i}N_{i^{\prime}},$ $\displaystyle V_{3}$
$\displaystyle\leq$ $\displaystyle\sum_{d(\bm{s}_{i},\bm{s}_{i^{\prime}})\geq
c_{n}}w_{i}w_{i^{\prime}}N_{i}N_{i^{\prime}}h^{-\frac{2(1+\varrho)}{2+\varrho}}\phi^{\frac{\varrho}{2+\varrho}}(d(\bm{s}_{i},\bm{s}_{i^{\prime}}))$
$\displaystyle\leq$ $\displaystyle
C\sum_{m=c_{n}}^{\infty}n^{-2}h^{-\frac{2(1+\varrho)}{2+\varrho}}\phi^{\frac{\varrho}{2+\varrho}}(m)nm\delta_{n}^{-2}=Ch^{-\frac{2(1+\varrho)}{2+\varrho}}n^{-1}\delta_{n}^{-2}\sum_{m=c_{n}}^{\infty}m\phi^{\frac{\varrho}{2+\varrho}}(m).$

Hence, (S.10) and the rate conditions imply that

$\displaystyle\sum_{k=1}^{M}\sum_{u=1}^{M}E(|W_{ku}|^{2})\leq
V_{1}+V_{2}+V_{3}=Q_{n}^{2}$ (S.11)

and

$\displaystyle|W_{ku}|\leq 2B_{n}h^{-1}\sum_{i=1}^{n}I(\bm{s}_{i}\in
A_{ku}^{(1)})w_{i}N_{i}\leq
B_{n}h^{-1}Cp^{2}\delta_{n}^{-2}\max_{i}w_{i}N_{i}\leq
CB_{n}n^{-1}h^{-1}p^{2}\delta_{n}^{-2}.$ (S.12)

Step 3: For any $\epsilon>0$, simple calculation leads to

$\displaystyle\mathbb{P}(|\sum_{k=1}^{M}\sum_{u=1}^{M}V_{ku}^{(1)}|>2\epsilon)$
$\displaystyle\leq$
$\displaystyle\mathbb{P}(|\sum_{k=1}^{M}\sum_{u=1}^{M}{V}_{ku}^{(1)}|>2\epsilon,|{V}_{ku}^{(1)}-W_{ku}|\leq\epsilon/M^{2},\textrm{
for all
}k,u)+\sum_{k=1}^{M}\sum_{u=1}^{M}\mathbb{P}(|{V}_{ku}^{(1)}-W_{ku}|>\epsilon/M^{2})$
$\displaystyle\leq$
$\displaystyle\mathbb{P}(|\sum_{k=1}^{M}\sum_{u=1}^{M}W_{ku}|>\epsilon)+\sum_{k=1}^{M}\sum_{u=1}^{M}\mathbb{P}(|{V}_{ku}^{(1)}-W_{ku}|>\epsilon/M^{2})$
$\displaystyle:=$ $\displaystyle R_{1}+R_{2}.$

Using (S.11) and (S.12), Bernstein inequality implies that

$\displaystyle R_{1}\leq
2\exp\left\\{-\frac{\epsilon^{2}}{2Q_{n}^{2}+\frac{2}{3}CB_{n}h^{-1}p^{2}\delta_{n}^{-2}\epsilon}\right\\}.$

The definition of $M$ and (S.9) together show that

$\displaystyle R_{2}\leq
CM^{2}\sqrt{\frac{B_{n}n^{-1}h^{-1}M^{2}p^{2}\delta_{n}^{-2}}{\epsilon}}\psi(Cp^{2}\delta_{n}^{-2},n)\phi(p)\leq\frac{C\Delta_{n}^{2}}{p^{2}}\psi(Cp^{2}\delta_{n}^{-2},n)\phi(p))\sqrt{\frac{B_{n}n^{-1}h^{-1}\Delta_{n}^{2}\delta_{n}^{-2}}{\epsilon}}.$

Combining the above three inequalities, we conclude that

$\displaystyle\mathbb{P}(|\sum_{k=1}^{M}\sum_{u=1}^{M}V_{ku}^{(1)}|>2\epsilon)$
$\displaystyle\leq$ $\displaystyle
2\exp\left\\{-\frac{\epsilon^{2}}{2Q_{n}^{2}+\frac{2}{3}CB_{n}h^{-1}p^{2}\delta_{n}^{-2}\epsilon}\right\\}$
(S.13)
$\displaystyle+\frac{C\Delta_{n}^{2}}{p^{2}}\psi(Cp^{2}\delta_{n}^{-2},n)\phi(p)\sqrt{\frac{B_{n}n^{-1}h^{-1}\Delta_{n}^{2}\delta_{n}^{-2}}{\epsilon}}.$

Similar inequalities as (S.13) can be proved for
$V_{ku}^{(2)},V_{ku}^{(3)},V_{ku}^{(4)}$ by the same argument. Hence, we
complete the proof using union bound and (S.8).

∎

###  S.1.2 Results for $\widehat{R}(s;t_{1},t_{2})$

In this section, we provide point-wise and $L_{2}$ convergence rates for the
spatial covariance estimator $\widehat{R}(s;t_{1},t_{2})$. The following
theorem provides the point-wise bias and variance.

###### Proof of Theorem 4.1.

For convenience in notation, we denote
$K_{ij}(t_{1})=K_{h_{c}}(T_{ij}-t_{1})$,
$L_{ii^{\prime}}=L_{b}(\|\bm{s}_{i}-\bm{s}_{i^{\prime}}\|-\|\bm{s}_{0}\|)$,
$\mathcal{A}_{ij,i^{\prime}j^{\prime}}=(T_{ij},T_{i^{\prime}j^{\prime}},\|\bm{s}_{i}-\bm{s}_{i^{\prime}}\|)^{T}$,
$\bm{a}=(t_{1},t_{2},\|\bm{s}_{0}\|)^{T}$, and
$\bm{\alpha}=(\alpha_{1},\alpha_{2},\alpha_{3})^{T}$. Let
$C_{ij,i^{\prime}j^{\prime}}:=[Z_{ij}-\widehat{\mu}(T_{ij})][Z_{i^{\prime}j^{\prime}}-\widehat{\mu}(T_{i^{\prime}j^{\prime}})]$,
$1\leq i,i^{\prime}\leq n$, $1\leq j\leq N_{i}$, $1\leq j^{\prime}\leq
N_{i^{\prime}}$. Because of Lemma S.1.3, for simplicity, in the following
proof we directly work with $\mu(t)$ instead of $\widehat{\mu}(t)$.

Note that $\widehat{R}(\|\bm{s}_{0}\|;t_{1},t_{2})=\widehat{\alpha}_{0}$,
where

$\displaystyle(\widehat{\alpha}_{0},\widehat{\bm{\alpha}}^{T})^{T}$
$\displaystyle=$
$\displaystyle\underset{(\alpha_{0},\bm{\alpha})\in\mathbb{R}^{4}}{\text{arg
min}}\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v_{i,i^{\prime}}\sum_{j=1}^{N_{i}}\sum_{j^{\prime}=1}^{N_{i^{\prime}}}\bigg{\\{}C_{ij,i^{\prime}j^{\prime}}-\alpha_{0}-\bm{\alpha}^{T}(\mathcal{A}_{ij,i^{\prime}j^{\prime}}-\bm{a})\bigg{\\}}^{2}$
(S.14) $\displaystyle\qquad\qquad\times
K_{ij}(t_{1})K_{i^{\prime}j^{\prime}}(t_{2})L_{ii^{\prime}}.$

Define

$\displaystyle\mathcal{P}=\begin{bmatrix}P_{00}&\bm{P}_{10}^{T}\\\
\bm{P}_{10}&\bm{P}_{11}\end{bmatrix}$ (S.15)

where

$\displaystyle P_{00}$ $\displaystyle=\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v_{i,i^{\prime}}\sum_{j=1}^{N_{i}}\sum_{j^{\prime}=1}^{N_{i^{\prime}}}K_{ij}(t_{1})K_{i^{\prime}j^{\prime}}(t_{2})L_{ii^{\prime}},$
(S.16)

and

$\displaystyle\bm{P}_{lk}$ $\displaystyle=\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v_{i,i^{\prime}}\sum_{j=1}^{N_{i}}\sum_{j^{\prime}=1}^{N_{i^{\prime}}}(\mathcal{A}_{ij,i^{\prime}j^{\prime}}-\bm{a})^{l}((\mathcal{A}_{ij,i^{\prime}j^{\prime}}-\bm{a})^{T})^{k}K_{ij}(t_{1})K_{i^{\prime}j^{\prime}}(t_{2})L_{ii^{\prime}},$
(S.17)

for $l,k\in\\{(1,0),(1,1)\\}$. Then, we can write

$\displaystyle\widehat{\bm{\theta}}:=\begin{pmatrix}\widehat{\alpha}_{0}-\alpha_{0}\\\
\widehat{\bm{\alpha}}-\bm{\alpha}\end{pmatrix}=\mathcal{P}^{-1}\mathcal{Q},\qquad\mathcal{Q}=\begin{pmatrix}q_{0}\\\
\bm{q}\end{pmatrix}$

where

$\displaystyle q_{0}$ $\displaystyle=\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v_{i,i^{\prime}}\sum_{j=1}^{N_{i}}\sum_{j^{\prime}=1}^{N_{i^{\prime}}}\eta_{ij,i^{\prime}j^{\prime}}K_{ij}(t_{1})K_{i^{\prime}j^{\prime}}(t_{2})L_{ii^{\prime}},$

and

$\displaystyle\bm{q}$ $\displaystyle=\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v_{i,i^{\prime}}\sum_{j=1}^{N_{i}}\sum_{j^{\prime}=1}^{N_{i^{\prime}}}\eta_{ij,i^{\prime}j^{\prime}}(\mathcal{A}_{ij,i^{\prime}j^{\prime}}-\bm{a})K_{ij}(t_{1})K_{i^{\prime}j^{\prime}}(t_{2})L_{ii^{\prime}},$

for
$\eta_{ij,i^{\prime}j^{\prime}}=C_{ij,i^{\prime}j^{\prime}}-\alpha_{0}(\bm{a})-(\mathcal{A}_{ij,i^{\prime}j^{\prime}}-\bm{a})^{T}\bm{\alpha}$.
First, we consider the elements of (S.15). It follows from Conditions (I)[(i),
(ii)] and (III) that

$\displaystyle E\\{K_{ij}(t_{1})K_{i^{\prime}j^{\prime}}(t_{2})\\}$
$\displaystyle=\int K(u)K(v)f_{T}(t_{1}+uh_{c})f_{T}(t_{2}+vh_{c})dudv$
$\displaystyle=f_{T}(t_{1})f_{T}(t_{2})(1+o(1)).$

Therefore, using Conditions (II)[(ii), (iii), (iv)], and (III), and proceeding
similar to Lu and Tjøstheim (2014) we obtain that

$\displaystyle E\\{P_{00}\\}$
$\displaystyle=\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v_{i,i^{\prime}}\sum_{j=1}^{N_{i}}\sum_{j^{\prime}=1}^{N_{i^{\prime}}}L_{ii^{\prime}}E\\{K_{ij}(t_{1})K_{i^{\prime}j^{\prime}}(t_{2})\\}$
$\displaystyle=\int\int
f_{T}(t_{1})f_{T}(t_{2})L((\|\bm{s}_{1}-\bm{s}_{2}\|-\|\bm{s}_{0}\|)/b)b^{-1}g_{\bm{S}}(\bm{s}_{1})g_{\bm{S}}(\bm{s}_{2})d\bm{s}_{1}d\bm{s}_{2}(1+o(1)).$

Now consider a transformation
$(\|\bm{s}_{1}-\bm{s}_{2}\|-\|\bm{s}_{0}\|)/b=s$, with
$\bm{s}_{j}=(u_{j},v_{j})$ for $j=1,2$, where
$u_{1}=u_{2}+(\|\bm{s}_{0}+bs\|)\cos(\vartheta)$ and
$v_{1}=v_{2}+(\|\bm{s}_{0}+bs\|)\sin(\vartheta)$ for $0<\vartheta\leq 2\pi$.
Let $u_{2}^{*}=u_{2}+\|\bm{s}_{0}\|\cos(\vartheta)$,
$v_{2}^{*}=v_{2}+\|\bm{s}_{0}\|\sin(\vartheta)$. Note that
$d\bm{s}_{1}=du_{1}dv_{1}=b(\|\bm{s}_{0}\|+bs)dsd\vartheta$. Consequently,

$\displaystyle f_{T}(t_{1})f_{T}(t_{2})\int b^{-1}$ $\displaystyle
L(s)g_{\bm{S}}(u_{2}^{*}+bs\cos(\vartheta),v_{2}^{*}+bs\sin(\vartheta))g_{\bm{S}}(u_{2},v_{2})$
$\displaystyle\qquad\qquad\qquad\times b(\|\bm{s}_{0}\|+bs)dsd\vartheta
du_{2}dv_{2}(1+o_{p}(1))$
$\displaystyle=f_{T}(t_{1})f_{T}(t_{2})A_{0}^{*}(\|\bm{s}_{0}\|)(1+o_{p}(1)).$

Next we show that $\text{var}\\{P_{00}\\}\rightarrow 0$. It follows from
Conditions (I)(ii), (II), (III), and (IV)(iv) that

$\displaystyle\text{var}\\{P_{00}\\}$
$\displaystyle=\text{var}\left\\{\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v_{i,i^{\prime}}\sum_{j=1}^{N_{i}}\sum_{j^{\prime}=1}^{N_{i^{\prime}}}L_{ii^{\prime}}K_{ij}(t_{1})K_{i^{\prime}j^{\prime}}(t_{2})\right\\}$
$\displaystyle=O_{p}\left(h_{c}^{-2}b^{-1}\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v^{2}_{i,i^{\prime}}N_{i}N_{i^{\prime}}\right),$

which goes to zero based on the Condition (IV)(iv). Consequently,

$\displaystyle P_{00}$
$\displaystyle=f_{T}(t_{1})f_{T}(t_{2})A_{0}^{*}(\|\bm{s}_{0}\|)(1+o_{p}(1)).$

For convenience in notation, in the following expression, we denote
$f_{j}=f_{T}(t_{j})$, $f^{\prime}_{j}=\partial f_{T}(t_{j})/\partial t_{j}$,
for $j=1,2$, $A_{0}^{*}=A_{0}^{*}(\|\bm{s}_{0}\|)$,
$A_{1}^{*}=A_{1}^{*}(\|\bm{s}_{0}\|)$, and
$A_{10}^{*}=A_{1}^{*}+\frac{A_{0}^{*}}{\|\bm{s}_{0}\|}$. The calculations
analogous to above yield the following, ignoring higher order terms,

$\displaystyle\mathcal{P}$
$\displaystyle=\begin{pmatrix}f_{1}f_{2}A_{0}^{*}&h_{c}^{2}f^{\prime}_{1}f_{2}A_{0}^{*}&h_{c}^{2}f_{1}f^{\prime}_{2}A_{0}^{*}&\mu_{L,2}b^{2}A_{10}^{*}\\\
h_{c}^{2}f^{\prime}_{1}f_{2}A_{0}^{*}&h_{c}^{2}f_{1}f_{2}A_{0}^{*}&h_{c}^{4}f^{\prime}_{1}f^{\prime}_{2}A_{0}^{*}&h_{c}^{2}f^{\prime}_{1}f_{2}b^{2}\mu_{L,2}A_{10}^{*}\\\
h_{c}^{2}f_{1}f^{\prime}_{2}A_{0}^{*}&h_{c}^{4}f^{\prime}_{1}f^{\prime}_{2}A_{0}^{*}&h_{c}^{2}f_{1}f_{2}A_{0}^{*}&h_{c}^{2}f_{1}f^{\prime}_{2}b^{2}\mu_{L,2}A_{10}^{*}\\\
b^{2}\mu_{L,2}A_{10}^{*}&h_{c}^{2}f^{\prime}_{1}f_{2}b^{2}\mu_{L,2}A_{10}^{*}&h_{c}^{2}f_{1}f^{\prime}_{2}b^{2}\mu_{L,2}A_{10}^{*}&b^{2}\mu_{L,2}f_{1}f_{2}A_{0}^{*}\end{pmatrix}.$
(S.18)

Bias Calculation: Observe that, by Taylor expansion, the leading term in
$E\left\\{\eta_{ij,i^{\prime}j^{\prime}}|T_{ij},T_{i^{\prime}j^{\prime}}\right\\}$
is $1/2[(T_{ij}-t_{1})^{2}\partial^{2}R(\|\bm{s}_{0}\|;t_{1},t_{2})/\partial
t_{1}^{2}+(T_{i^{\prime}j^{\prime}}-t_{2})^{2}\partial^{2}R(\|\bm{s}_{0}\|;t_{1},t_{2})/\partial
t_{2}^{2}$ (∥si-si′∥-∥s0∥)2∂2R(∥s0∥;t1,t2)/∂∥s0∥]. We write

$\displaystyle E[q_{0}]$ $\displaystyle=\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v_{i,i^{\prime}}\sum_{j=1}^{N_{i}}\sum_{j^{\prime}=1}^{N_{i^{\prime}}}E[\eta_{ij,i^{\prime}j^{\prime}}K_{ij}(t_{1})K_{i^{\prime}j^{\prime}}(t_{2})L_{ii^{\prime}}]$
$\displaystyle=\frac{1}{2}\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}L_{ii^{\prime}}v_{i,i^{\prime}}N_{i}N_{i^{\prime}}\int\big{\\{}(T_{ij}-t_{1})^{2}\partial^{2}R/\partial
t_{1}^{2}+(T_{i^{\prime}j^{\prime}}-t_{2})^{2}\partial^{2}R/\partial
t_{2}^{2}$
$\displaystyle+(\|\bm{s}_{i}-\bm{s}_{i^{\prime}}\|-\|\bm{s}_{0}\|)^{2}\partial^{2}R/\partial\|\bm{s}_{0}\|^{2}\big{\\}}K_{ij}(t_{1})K_{i^{\prime}j^{\prime}}(t_{2})f_{T}(T_{ij})f_{T}(T_{i^{\prime}j^{\prime}})dT_{ij}dT_{ij^{\prime}}+o(1)$
$\displaystyle=\frac{1}{2}\bigg{\\{}h_{c}^{2}\big{\\{}\partial^{2}R/\partial
t_{1}^{2}+\partial^{2}R/\partial
t_{2}^{2}\big{\\}}+b^{2}\mu_{L,2}\partial^{2}R/\partial\|\bm{s}_{0}\|^{2}\bigg{\\}}f_{T}(t_{1})f_{T}(t_{2})A_{0}^{*}+o_{p}(h_{c}^{2}+b^{2}),$

where $\mu_{L,2}=\int u^{2}L(u)du$. Consequently,

$\displaystyle
E\\{\widehat{R}(\|\bm{s}_{0}\|;t_{1},t_{2})-R(\|\bm{s}_{0}\|;t_{1},t_{2})~{}|~{}\mathcal{X}_{T}\\}$
$\displaystyle=\frac{1}{2}h_{c}^{2}\big{\\{}\partial^{2}R/\partial
t_{1}^{2}+\partial^{2}R/\partial t_{2}^{2}\big{\\}}$
$\displaystyle\qquad\qquad+\frac{1}{2}b^{2}\mu_{L,2}\partial^{2}R/\partial\|\bm{s}_{0}\|^{2}+o_{p}(h_{c}^{2}+b^{2}),$

where $\mathcal{X}_{T}=\\{T_{ij},i=1,\ldots,n,j=1,\ldots,N_{i}\\}$.

Variance Calculation: We now compute the conditional variance for the
estimator $\widehat{R}(\|\bm{s}_{0}\|;t_{1},t_{2})$. By definition

$\displaystyle\text{var}(\widehat{R}(\|\bm{s}_{0}\|;t_{1},t_{2})|\mathcal{X}_{T})$
$\displaystyle=$
$\displaystyle\text{var}(\bm{e}_{1}^{T}\bm{\theta}|\mathcal{X}_{T})=\text{var}\left\\{P^{00}q_{0}+\bm{P}^{{10}^{T}}\bm{q}_{1}|\mathcal{X}_{T}\right\\},$
(S.19)

where $P^{00}$ and $\bm{P}^{10}$ are the elements of the first row in matrix
$\mathcal{P}^{-1}$ which is the inverse of (S.18). The leading term in the
variance of (S.19) is

$\displaystyle\text{var}(P^{00}q_{0}|\mathcal{X}_{T})$
$\displaystyle=(P^{00})^{2}\bigg{[}\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v_{i,i^{\prime}}^{2}\text{var}\bigg{(}\sum_{j=1}^{N_{i}}\sum_{j^{\prime}=1}^{N_{i^{\prime}}}\eta_{ij,i^{\prime}j^{\prime}}K_{ij}(t_{1})K_{i^{\prime}j^{\prime}}(t_{2})L_{ii^{\prime}}~{}|~{}\mathcal{X}_{T}\bigg{)}$
$\displaystyle+\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}\underset{(l,l^{\prime})\neq(i,i^{\prime})}{\sum_{l}\sum_{l^{\prime}\neq
l}}v_{i,i^{\prime}}v_{l,l^{\prime}}$
$\displaystyle\times\text{cov}\bigg{(}\sum_{j}^{N_{i}}\sum_{j^{\prime}}^{N_{i^{\prime}}}\eta_{ij,i^{\prime}j^{\prime}}K_{ij}(t_{1})K_{i^{\prime}j^{\prime}}(t_{2})L_{ii^{\prime}},\sum_{k}^{N_{l}}\sum_{k^{\prime}}^{N_{l^{\prime}}}\eta_{lk,l^{\prime}k^{\prime}}K_{lk}(t_{1})K_{l^{\prime}k^{\prime}}(t_{2})L_{ll^{\prime}}~{}|~{}\mathcal{X}_{T}\bigg{)}\bigg{]}$
$\displaystyle:=(P^{00})^{2}\\{(I)+(II)\\}.$ (S.20)

We show that the term (II) goes to zero as $n\rightarrow\infty$. Therefore the
leading term for the variance is obtained from (I).

Term (I): Observe that

var
$\displaystyle\bigg{(}\sum_{j=1}^{N_{i}}\sum_{j^{\prime}=1}^{N_{i^{\prime}}}\eta_{ij,i^{\prime}j^{\prime}}K_{ij}(t_{1})K_{i^{\prime}j^{\prime}}(t_{2})L_{ii^{\prime}}~{}|~{}\mathcal{X}_{T}\bigg{)}$
$\displaystyle=\sum_{j=1}^{N_{i}}\sum_{j^{\prime}=1}^{N_{i^{\prime}}}\text{var}\bigg{(}\eta_{ij,i^{\prime}j^{\prime}}K_{ij}(t_{1})K_{i^{\prime}j^{\prime}}(t_{2})L_{ii^{\prime}}~{}|~{}\mathcal{X}_{T}\bigg{)}$
$\displaystyle\qquad+\sum_{j=1}^{N_{i}}\sum_{j^{\prime}=1}^{N_{i^{\prime}}}\underset{(k,k^{\prime})\neq(j,j^{\prime})}{\sum_{k}\sum_{k^{\prime}}}\text{cov}\big{\\{}\eta_{ij,i^{\prime}j^{\prime}}K_{ij}(t_{1})K_{i^{\prime}j^{\prime}}(t_{2})L_{ii^{\prime}},\eta_{ik,i^{\prime}k^{\prime}}K_{ik}(t_{1})K_{i^{\prime}k^{\prime}}(t_{2})L_{ii^{\prime}}~{}|~{}\mathcal{X}_{T}\big{\\}}.$

Let $R(K_{1})=\int K_{1}^{2}(u)du$ and $R(L)=\int L^{2}(u)u$. Denote

$\displaystyle\Xi_{R}(t_{1},t_{2},t_{1},t_{2};\|\bm{s}_{i}-\bm{s}_{i^{\prime}}\|)=\text{cov}\big{\\{}\eta_{ij,i^{\prime}j^{\prime}},\eta_{ik,i^{\prime}k^{\prime}}|T_{ij}=t_{1},T_{i^{\prime}j^{\prime}}=t_{2},T_{ik}=t_{1},T_{i^{\prime}k^{\prime}}=t_{2}\\}.$

Under Conditions (I)[(i), (ii), (vi)], (II) and (III), and using the standard
arguments, we write the expected value of part (I) in (S.20) as

$\displaystyle E[(I)]=\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v^{2}_{i,i^{\prime}}N_{i}N_{i^{\prime}}h_{c}^{-2}b^{-1}R(L)R^{2}(K_{1})f_{T}(t_{1})f_{T}(t_{2})~{}\Xi_{R}(t_{1},t_{2};\|\bm{s}_{0}\|)A_{0}^{*}$
$\displaystyle\qquad+\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v^{2}_{i,i^{\prime}}N_{i}N_{i^{\prime}}(N_{i}-1)(N_{i^{\prime}}-1)b^{-1}R(L)f_{T}^{2}(t_{1})f_{T}^{2}(t_{2})~{}\Xi_{R}(t_{1},t_{2},t_{1},t_{2};\|\bm{s}_{0}\|)A_{0}^{*}$
$\displaystyle\qquad+\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v^{2}_{i,i^{\prime}}N_{i}N_{i^{\prime}}(N_{i^{\prime}}-1)h_{c}^{-1}b^{-1}R(L)R(K_{1})f_{T}(t_{1})f_{T}^{2}(t_{2})~{}\Xi_{R}(t_{1},t_{2},t_{2};\|\bm{s}_{0}\|)A_{0}^{*}$
$\displaystyle\qquad+\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v^{2}_{i,i^{\prime}}N_{i}N_{i^{\prime}}(N_{i}-1)h_{c}^{-1}b^{-1}R(L)R(K_{1})f_{T}^{2}(t_{1})f_{T}(t_{2})~{}\Xi_{R}(t_{1},t_{2},t_{1};\|\bm{s}_{0}\|)A_{0}^{*}+o_{p}(1).$

Consequently, the leading term in the conditional variance is as follows

$\displaystyle(P^{00})^{2}E[(I)]$
$\displaystyle=\frac{\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v^{2}_{i,i^{\prime}}N_{i}N_{i^{\prime}}}{h_{c}^{2}bA_{0}^{*}f_{T}(t_{1})f_{T}(t_{2})}R(L)R^{2}(K_{1})~{}\Xi_{R}(t_{1},t_{2};\|\bm{s}_{0}\|)$
$\displaystyle\qquad+\frac{\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v^{2}_{i,i^{\prime}}N_{i}N_{i^{\prime}}(N_{i^{\prime}}-1)}{h_{c}bA_{0}^{*}f_{T}(t_{1})}R(L)R(K_{1})\Xi_{R}(t_{1},t_{2},t_{2};\|\bm{s}_{0}\|)$
$\displaystyle\qquad+\frac{\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v^{2}_{i,i^{\prime}}N_{i}N_{i^{\prime}}(N_{i}-1)}{h_{c}bA_{0}^{*}f_{T}(t_{2})}R(L)R(K_{1})\Xi_{R}(t_{1},t_{2},t_{1};\|\bm{s}_{0}\|)$
$\displaystyle\qquad+\frac{\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}v^{2}_{i,i^{\prime}}N_{i}(N_{i}-1)N_{i^{\prime}}(N_{i^{\prime}}-1)}{bA_{0}^{*}}R(L)\Xi_{R}(t_{1},t_{2},t_{1},t_{2};\|\bm{s}_{0}\|)\bigg{\\}}$
$\displaystyle\qquad\qquad\qquad\times(1+o_{p}(1)).$

Term (II): The idea is to relate the conditional covariance term with
unconditional covariance terms using the following fact

$\displaystyle E[\text{cov}(A,B|C)]$
$\displaystyle=\text{cov}(A,B)-\text{cov}(E[A|C],E[B|C]),$ (S.21)

for random variables $A$, $B$, and $C$. For simplicity in notation, let us
denote

$\sum_{(l,l^{\prime})\neq(i,i^{\prime})}=\sum_{i=1}^{n}\sum_{i^{\prime}\neq
i}\underset{(l,l^{\prime})\neq(i,i^{\prime})}{\sum_{l}\sum_{l^{\prime}\neq
l}}.$

We write the term (II) in (S.20) as

$\displaystyle E\bigg{\\{}$
$\displaystyle\sum_{(l,l^{\prime})\neq(i,i^{\prime})}v_{i,i^{\prime}}v_{l,l^{\prime}}\text{cov}\bigg{(}\sum_{j}^{N_{i}}\sum_{j^{\prime}}^{N_{i^{\prime}}}\eta_{ij,i^{\prime}j^{\prime}}K_{ij}(t_{1})K_{i^{\prime}j^{\prime}}(t_{2})L_{ii^{\prime}},$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\sum_{k}^{N_{l}}\sum_{k^{\prime}}^{N_{l^{\prime}}}\eta_{lk,l^{\prime}k^{\prime}}K_{lk}(t_{1})K_{l^{\prime}k^{\prime}}(t_{2})L_{ll^{\prime}}~{}|~{}\mathcal{X}_{T}\bigg{)}\bigg{\\}}$
$\displaystyle=\sum_{(l,l^{\prime})\neq(i,i^{\prime})}v_{i,i^{\prime}}v_{l,l^{\prime}}\text{cov}\bigg{(}\sum_{j}^{N_{i}}\sum_{j^{\prime}}^{N_{i^{\prime}}}\eta_{ij,i^{\prime}j^{\prime}}K_{ij}(t_{1})K_{i^{\prime}j^{\prime}}(t_{2})L_{ii^{\prime}},\sum_{k}^{N_{l}}\sum_{k^{\prime}}^{N_{l^{\prime}}}\eta_{lk,l^{\prime}k^{\prime}}K_{lk}(t_{1})K_{l^{\prime}k^{\prime}}(t_{2})L_{ll^{\prime}}\bigg{)}$
$\displaystyle-\sum_{(l,l^{\prime})\neq(i,i^{\prime})}v_{i,i^{\prime}}v_{l,l^{\prime}}\text{cov}\bigg{(}\sum_{j}^{N_{i}}\sum_{j^{\prime}}^{N_{i^{\prime}}}E[\eta_{ij,i^{\prime}j^{\prime}}K_{ij}(t_{1})K_{i^{\prime}j^{\prime}}(t_{2})L_{ii^{\prime}}~{}|~{}\mathcal{X}_{T}],$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\sum_{k}^{N_{l}}\sum_{k^{\prime}}^{N_{l^{\prime}}}E[\eta_{lk,l^{\prime}k^{\prime}}K_{lk}(t_{1})K_{l^{\prime}k^{\prime}}(t_{2})L_{ll^{\prime}}~{}|~{}\mathcal{X}_{T}]\bigg{)}.$
(S.22)

From Lemma S.1.6 it directly follows that the first term in the right hand
side of the equality in (S.22) goes to zero. Recall that

$\displaystyle
E\left\\{\eta_{ij,i^{\prime}j^{\prime}}|\mathcal{X}_{T}\right\\}=1/2[(T_{ij}-t_{1})^{2}\partial^{2}R/\partial
t_{1}^{2}+(T_{i^{\prime}j^{\prime}}-t_{2})^{2}\partial^{2}R/\partial
t_{2}^{2}+(\|\bm{s}_{i}-\bm{s}_{i^{\prime}}\|-\|\bm{s}_{0}\|)^{2}\partial^{2}R/\partial\|\bm{s}_{0}\|].$

Since $T_{ij}$’s are independent, direct calculations yield that the second
term in the right side of equality in (S.22) goes to zero under condition
(IV)(iv). Hence, the result is proved.

∎

######  Lemma S.1.6.

Under the conditions of Theorem 4.1, it holds that

$\displaystyle\sum_{(l,l^{\prime})\neq(i,i^{\prime})}v_{i,i^{\prime}}v_{l,l^{\prime}}\text{cov}\bigg{\\{}\sum_{j=1}^{N_{i}}\sum_{j^{\prime}=1}^{N_{i^{\prime}}}\eta_{ij,i^{\prime}j^{\prime}}L_{ii^{\prime}}K_{ij}(t_{1})K_{i^{\prime}j^{\prime}}(t_{2}),\sum_{k=1}^{N_{l}}\sum_{k^{\prime}=1}^{N_{l^{\prime}}}\eta_{lk,l^{\prime}k^{\prime}}L_{ll^{\prime}}K_{lk}(t_{1})K_{l^{\prime}k^{\prime}}(t_{2})\bigg{\\}}$
(S.23)

goes to zero as $n\rightarrow\infty$.

###### Proof.

Since the random variable $\eta_{ij,ij^{\prime}}$ is not bounded, we use the
truncation argument in the following proof. Let $L_{n}=h^{-4/(4+\varrho)}$ for
some $\varrho>0$ defined in Condition (II)(i). Define
$\eta_{ij,i^{\prime}j^{\prime}}=\eta_{ij,i^{\prime}j^{\prime},1}+\eta_{ij,i^{\prime}j^{\prime},2}$
where
$\eta_{ij,i^{\prime}j^{\prime},1}:=\eta_{ij,i^{\prime}j^{\prime}}I_{\\{|\eta_{ij,i^{\prime}j^{\prime}}|\leq
L_{n}\\}}$ and
$\eta_{ij,i^{\prime}j^{\prime},2}:=\eta_{ij,i^{\prime}j^{\prime}}I_{\\{|\eta_{ij,i^{\prime}j^{\prime}}|>L_{n}\\}}$.
Let

$\displaystyle\chi_{ij,i^{\prime}j^{\prime},l}:=\eta_{ij,i^{\prime}j^{\prime},l}K_{ij}(t_{1})K_{i^{\prime}j^{\prime}}(t_{2})\text{
and
}\varpi_{ij,i^{\prime}j^{\prime},l}:=\chi_{ij,i^{\prime}j^{\prime},l}-E\chi_{ij,i^{\prime}j^{\prime},l},\qquad
l=1,2.$

Observe that
$\varpi_{ij,i^{\prime}j^{\prime}}=\varpi_{ij,i^{\prime}j^{\prime},1}+\varpi_{ij,i^{\prime}j^{\prime},2}$.
Therefore,

$\displaystyle
E\varpi_{ij,i^{\prime}j^{\prime}}\varpi_{lk,l^{\prime}k^{\prime}}$
$\displaystyle=E\varpi_{ij,i^{\prime}j^{\prime},1}\varpi_{lk,l^{\prime}k^{\prime},1}+E\varpi_{ij,i^{\prime}j^{\prime},1}\varpi_{lk,l^{\prime}k^{\prime},2}$
$\displaystyle\qquad\qquad+E\varpi_{ij,i^{\prime}j^{\prime},2}\varpi_{lk,l^{\prime}k^{\prime},1}+E\varpi_{ij,i^{\prime}j^{\prime},2}\varpi_{lk,l^{\prime}k^{\prime},2}.$
(S.24)

We now provide bounds for each term in (S.24). By Cauchy-Schwartz inequality,
note that,

$\displaystyle|E\varpi_{ij,i^{\prime}j^{\prime},1}$
$\displaystyle\varpi_{lk,l^{\prime}k^{\prime},2}|$
$\displaystyle\leq\\{E\chi_{ij,i^{\prime}j^{\prime},1}^{2}\\}^{1/2}\\{E\chi_{lk,l^{\prime}k^{\prime},2}^{2}\\}^{1/2}$
$\displaystyle\leq
M_{71}h_{c}^{-1}\\{E(E(|\eta_{lk,l^{\prime}k^{\prime}}|^{2}I_{\\{|\eta_{lk,l^{\prime}k^{\prime}}|>L_{n}\\}}|T_{lk},T_{l^{\prime}k^{\prime}})K^{2}_{lk}(t_{1})K^{2}_{l^{\prime}k^{\prime}}(t_{2})\\}^{1/2}$
$\displaystyle\leq
M_{71}h_{c}^{-1}\\{h_{c}^{-2}L_{n}^{-\varrho}E(|\eta_{lk,l^{\prime}k^{\prime}}|^{2+\varrho}|T_{lk},T_{l^{\prime}k^{\prime}})\\}^{1/2}$
$\displaystyle\leq
M_{72}h_{c}^{-2}L_{n}^{-\varrho/2}=M_{72}h_{c}^{-8/(4+\varrho)},$

where the last step follows from Condition (I)(v)(b) and $M_{71}$ and $M_{72}$
are some constants. Analogous arguments yield

$\displaystyle
E\varpi_{ij,i^{\prime}j^{\prime},2}\varpi_{lk,l^{\prime}k^{\prime},1}$
$\displaystyle\leq
M_{72}h_{c}^{-2}L_{n}^{-\varrho/2}=M_{72}h_{c}^{-8/(4+\varrho)}\qquad\text{
and }$ $\displaystyle
E\varpi_{ij,i^{\prime}j^{\prime},2}\varpi_{lk,l^{\prime}k^{\prime},2}$
$\displaystyle\leq
M_{72}h_{c}^{-2}L_{n}^{-\varrho}=M_{72}h_{c}^{2(\varrho-4)/(4+\varrho)}.$

Consider

$\displaystyle
E\varpi_{ij,i^{\prime}j^{\prime},1}\varpi_{lk,l^{\prime}k^{\prime},1}$
$\displaystyle=\text{cov}\\{\eta_{ij,i^{\prime}j^{\prime},1}K_{ij}(t_{1})K_{i^{\prime}j^{\prime}}(t_{2}),\eta_{lk,l^{\prime}k^{\prime},1}K_{lk}(t_{1})K_{l^{\prime}k^{\prime}}(t_{2})\\}$
$\displaystyle=E\\{\eta_{ij,i^{\prime}j^{\prime},1}K_{ij}(t_{1})K_{i^{\prime}j^{\prime}}(t_{2})\eta_{lk,l^{\prime}k^{\prime},1}K_{lk}(t_{1})K_{l^{\prime}k^{\prime}}(t_{2})\\}$
$\displaystyle\qquad\qquad-E\\{\eta_{ij,i^{\prime}j^{\prime},1}K_{ij}(t_{1})K_{i^{\prime}j^{\prime}}(t_{2}))E\\{\eta_{lk,l^{\prime}k^{\prime},1}K_{lk}(t_{1})K_{l^{\prime}k^{\prime}}(t_{2})\\}$
$\displaystyle=\int
K_{1}(u)K_{1}(u^{\prime})K_{1}(v)K_{1}(v^{\prime})\\{\Xi_{1,i,i^{\prime},l,l^{\prime}}(u,v,u^{\prime},v^{\prime})$
$\displaystyle\qquad-\Xi_{1,i,i^{\prime}}(u,v)\Xi_{1,l,l^{\prime}}(u^{\prime},v^{\prime})\\}f_{T}(u)f_{T}(v)f_{T}(u^{\prime})f_{T}(v^{\prime})dudvdu^{\prime}dv^{\prime},$

where

$\displaystyle\Xi_{1,i,i^{\prime},l,l^{\prime}}(u,v,u^{\prime},v^{\prime})$